text:
function ∇ERI_2e4c(BS::BasisSet, iA, T::DataType = Float64)
    # Pre allocate output
    out = zeros(T, BS.nbas, BS.nbas, BS.nbas, BS.nbas, 3)
    return ∇ERI_2e4c!(out, BS, iA, T)
end

function ∇ERI_2e4c!(out, BS::BasisSet, iA, T::DataType = Float64)

    if size(out) != (BS.nbas, BS.nbas, BS.nbas, BS.nbas, 3)
        throw(DimensionMismatch("Size of the output array needs to be (N, N, N, N, 3)."))
    end

    A = BS.atoms[iA]

    # Get shell index (s0) where basis of the desired atom start
    s0 = 0
    for a in eachindex(BS.atoms)
        if BS.atoms[a] != A
            s0 += length(BS.basis[a])
        else
            break
        end
    end

    # Shell indexes for basis in the atom A
    Ashells = [(s0 + i -1) for i = 1:length(BS.basis[iA])]
    notAshells = Int[]
    for i = 1:BS.nshells
        if !((i-1) in Ashells)
            push!(notAshells, i-1)
        end
    end

    bas_per_shell = [Libcint.CINTcgtos_spheric(i-1, BS.lc_bas) for i = 1:BS.nshells]
    ao_offset = [sum(bas_per_shell[1:(i-1)]) for i = 1:BS.nshells]
    NBmax = maximum(bas_per_shell)
    buf = zeros(Cdouble, 3*NBmax^4)

    # Find unique (i,j,k,l) combinations given permutational symmetry
    unique_idx = NTuple{4,Int16}[]
    N = Int16(BS.nshells - 1)
    ZERO = zero(Int16)
    for i = ZERO:N
        for j = i:N
            for k = ZERO:N
                for l = k:N
                    if index2(i,j) < index2(k,l)
                        continue
                    end
                    push!(unique_idx, (i,j,k,l))
                end
            end
        end
    end

    for (i,j,k,l) in unique_idx

        x_in_A = [x in Ashells for x = (i,j,k,l)]

        # If no basis is centered on A, skip
        # If all basis are centered on A, skip
        if !any(x_in_A) || all(x_in_A)
            continue
        end

        Li = bas_per_shell[i+1]
        Lj = bas_per_shell[j+1]
        Lk = bas_per_shell[k+1]
        Ll = bas_per_shell[l+1]
        Lijkl = Li*Lj*Lk*Ll

        ioff = ao_offset[i+1]
        joff = ao_offset[j+1]
        koff = ao_offset[k+1]
        loff = ao_offset[l+1]

        I = (ioff+1):(ioff+Li)
        J = (joff+1):(joff+Lj)
        K = (koff+1):(koff+Lk)
        L = (loff+1):(loff+Ll)

        # [i'j|kl]
        if x_in_A[1]
            cint2e_ip1_sph!(buf, Cint.([i,j,k,l]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
            for q in 1:3
                r = (1+Lijkl*(q-1)):(q*Lijkl)
                ∇q = reshape(buf[r], Int(Li), Int(Lj), Int(Lk), Int(Ll))
                out[I,J,K,L,q] += -∇q
            end
        end

        # [ij'|kl]
        if x_in_A[2]
            cint2e_ip1_sph!(buf, Cint.([j,i,k,l]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
            for q in 1:3
                r = (1+Lijkl*(q-1)):(q*Lijkl)
                ∇q = reshape(buf[r], Int(Lj), Int(Li), Int(Lk), Int(Ll))
                out[I,J,K,L,q] += -permutedims(∇q, (2,1,3,4))
            end
        end

        # [ij|k'l]
        if x_in_A[3]
            cint2e_ip1_sph!(buf, Cint.([k,l,i,j]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
            for q in 1:3
                r = (1+Lijkl*(q-1)):(q*Lijkl)
                ∇q = reshape(buf[r], Int(Lk), Int(Ll), Int(Li), Int(Lj))
                out[I,J,K,L,q] += -permutedims(∇q, (3,4,1,2))
            end
        end

        # [ij|kl']
        if x_in_A[4]
            cint2e_ip1_sph!(buf, Cint.([l,k,i,j]), BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
            for q in 1:3
                r = (1+Lijkl*(q-1)):(q*Lijkl)
                ∇q = reshape(buf[r], Int(Ll), Int(Lk), Int(Li), Int(Lj))
                out[I,J,K,L,q] += -permutedims(∇q, (3,4,2,1))
            end
        end

        for q in 1:3
            if i != j && k != l && index2(i,j) != index2(k,l)
                # i,j permutation
                out[J, I, K, L, q] .= permutedims(out[I, J, K, L, q], (2,1,3,4))
                # k,l permutation
                out[I, J, L, K, q] .= permutedims(out[I, J, K, L, q], (1,2,4,3))

                # i,j + k,l permutatiom
                out[J, I, L, K, q] .= permutedims(out[I, J, K, L, q], (2,1,4,3))

                # ij, kl permutation
                out[K, L, I, J, q] .= permutedims(out[I, J, K, L, q], (3,4,1,2))
                # ij, kl + k,l permutation
                out[L, K, I, J, q] .= permutedims(out[I, J, K, L, q], (4,3,1,2))
                # ij, kl + i,j permutation
                out[K, L, J, I, q] .= permutedims(out[I, J, K, L, q], (3,4,2,1))
                # ij, kl + i,j + k,l permutation
                out[L, K, J, I, q] .= permutedims(out[I, J, K, L, q], (4,3,2,1))

            elseif k != l && index2(i,j) != index2(k,l)
                # k,l permutation
                out[I, J, L, K, q] .= permutedims(out[I, J, K, L, q], (1,2,4,3))
                # ij, kl permutation
                out[K, L, I, J, q] .= permutedims(out[I, J, K, L, q], (3,4,1,2))
                # ij, kl + k,l permutation
                out[L, K, I, J, q] .= permutedims(out[I, J, K, L, q], (4,3,1,2))

            elseif i != j && index2(i,j) != index2(k,l)
                # i,j permutation
                out[J, I, K, L, q] .= permutedims(out[I, J, K, L, q], (2,1,3,4))

                # ij, kl permutation
                out[K, L, I, J, q] .= permutedims(out[I, J, K, L, q], (3,4,1,2))
                # ij, kl + i,j permutation
                out[K, L, J, I, q] .= permutedims(out[I, J, K, L, q], (3,4,2,1))

            elseif i != j && k != l
                # i,j permutation
                out[J, I, K, L, q] .= permutedims(out[I, J, K, L, q], (2,1,3,4))
                # k,l permutation
                out[I, J, L, K, q] .= permutedims(out[I, J, K, L, q], (1,2,4,3))

                # i,j + k,l permutatiom
                out[J, I, L, K, q] .= permutedims(out[I, J, K, L, q], (2,1,4,3))
            elseif index2(i,j) != index2(k,l)
                # ij, kl permutation
                out[K, L, I, J, q] .= permutedims(out[I, J, K, L, q], (3,4,1,2))
            end
        end
    end

    return out
end

function ∇sparseERI_2e4c(BS::BasisSet, iA, T::DataType = Float64)

    A = BS.atoms[iA]

    # Number of unique integral elements
    N = Int((BS.nbas^2 - BS.nbas)/2) + BS.nbas
    N = Int((N^2 - N)/2) + N

    # Pre allocate output
    ∇x = zeros(T, N)
    ∇y = zeros(T, N)
    ∇z = zeros(T, N)
    indexes = Array{NTuple{4,Int16}}(undef, N)
    mask = repeat([false], N)

    # Get shell index (s0) where basis of the desired atom start
    s0 = 0
    for a in eachindex(BS.atoms)
        if BS.atoms[a] != A
            s0 += length(BS.basis[a])
        else
            break
        end
    end

    # Shell indexes for basis in the atom A
    # C indexing
    Ashells = [(s0 + i -1) for i = 1:length(BS.basis[iA])]
    notAshells = Int[]
    for i = 1:BS.nshells
        if !((i-1) in Ashells)
            push!(notAshells, i-1)
        end
    end

    # Pre compute a list of number of basis for each shell (2l +1)
    bas_per_shell = [Libcint.CINTcgtos_spheric(i-1, BS.lc_bas) for i = 1:BS.nshells]
    NBmax = maximum(bas_per_shell)

    # Offset list for each shell, used to map shell index to AO index
    ao_offset = [sum(bas_per_shell[1:(i-1)]) - 1 for i = 1:BS.nshells]

    # Unique shell pairs with i < j
    num_ij = Int((BS.nshells^2 - BS.nshells)/2) + BS.nshells

    # Save ij pairs
    ij_vals = Array{NTuple{2,Int32}}(undef, num_ij)
    lim = Int32(BS.nshells - 1)
    for i = UnitRange{Int32}(zero(Int32),lim)
        for j = UnitRange{Int32}(i, lim)
            idx = index2(i,j) + 1
            ij_vals[idx] = (i,j)
        end
    end

    buf_arrays = [zeros(Cdouble, 3*NBmax^4) for _ = 1:Threads.nthreads()]

    # i,j,k,l => Shell indexes starting at zero
    # I, J, K, L => AO indexes starting at one
    for ij in eachindex(ij_vals)
    #@sync for ij in eachindex(ij_vals)
        #Threads.@spawn begin
        #@inbounds begin
            buf = buf_arrays[Threads.threadid()]
            i,j = ij_vals[ij]
            Li, Lj = bas_per_shell[i+1], bas_per_shell[j+1]
            Lij = Li*Lj
            ioff = ao_offset[i+1]
            joff = ao_offset[j+1]
            for kl in ij:num_ij
                k,l = ij_vals[kl]

                # If all basis are centered at A, or none is, the derivative is zero
                x_in_A = [x in Ashells for x = (i,j,k,l)]
                if !any(x_in_A) || all(x_in_A)
                    continue
                end

                Lk, Ll = bas_per_shell[k+1], bas_per_shell[l+1]
                Lijk = Lij*Lk
                Lijkl = Lijk*Ll
                koff = ao_offset[k+1]
                loff = ao_offset[l+1]

                # NOTE: Using loops instead of array operations could make this more efficient
                # Compute ERI
                bufx = zeros(Cdouble, Int(Li), Int(Lj), Int(Lk), Int(Ll))
                bufy = zeros(Cdouble, Int(Li), Int(Lj), Int(Lk), Int(Ll))
                bufz = zeros(Cdouble, Int(Li), Int(Lj), Int(Lk), Int(Ll))
                if x_in_A[1]
                    cint2e_ip1_sph!(buf, [i,j,k,l], BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                    bufx += reshape(buf[1:Lijkl], Int(Li), Int(Lj), Int(Lk), Int(Ll))
                    bufy += reshape(buf[Lijkl+1:2*Lijkl], Int(Li), Int(Lj), Int(Lk), Int(Ll))
                    bufz += reshape(buf[2*Lijkl+1:3*Lijkl], Int(Li), Int(Lj), Int(Lk), Int(Ll))
                end

                if x_in_A[2]
                    cint2e_ip1_sph!(buf, [j,i,k,l], BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                    bufx += permutedims(reshape(buf[1:Lijkl],           Int(Lj), Int(Li), Int(Lk), Int(Ll)), (2,1,3,4))
                    bufy += permutedims(reshape(buf[Lijkl+1:2*Lijkl],   Int(Lj), Int(Li), Int(Lk), Int(Ll)), (2,1,3,4))
                    bufz += permutedims(reshape(buf[2*Lijkl+1:3*Lijkl], Int(Lj), Int(Li), Int(Lk), Int(Ll)), (2,1,3,4))
                end

                if x_in_A[3]
                    cint2e_ip1_sph!(buf, [k,l,i,j], BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                    bufx += permutedims(reshape(buf[1:Lijkl],           Int(Lk), Int(Ll), Int(Li), Int(Lj)), (3,4,1,2))
                    bufy += permutedims(reshape(buf[Lijkl+1:2*Lijkl],   Int(Lk), Int(Ll), Int(Li), Int(Lj)), (3,4,1,2))
                    bufz += permutedims(reshape(buf[2*Lijkl+1:3*Lijkl], Int(Lk), Int(Ll), Int(Li), Int(Lj)), (3,4,1,2))
                end

                if x_in_A[4]
                    cint2e_ip1_sph!(buf, [l,k,i,j], BS.lc_atoms, BS.natoms, BS.lc_bas, BS.nbas, BS.lc_env)
                    bufx += permutedims(reshape(buf[1:Lijkl],           Int(Ll), Int(Lk), Int(Li), Int(Lj)), (3,4,2,1))
                    bufy += permutedims(reshape(buf[Lijkl+1:2*Lijkl],   Int(Ll), Int(Lk), Int(Li), Int(Lj)), (3,4,2,1))
                    bufz += permutedims(reshape(buf[2*Lijkl+1:3*Lijkl], Int(Ll), Int(Lk), Int(Li), Int(Lj)), (3,4,2,1))
                end

                ### This block aims to retrieve unique elements within buf and map them to AO indexes
                # is, js, ks, ls are indexes within the shell e.g. for a p shell is = (1, 2, 3)
                # bl, bkl, bjkl are used to map the (i,j,k,l) index into a one-dimensional index for buf
                # That is, get the correct integrals for the AO quartet.
                for ls = 1:Ll
                    L = loff + ls
                    bl = Lijk*(ls-1)
                    for ks = 1:Lk
                        K = koff + ks
                        L < K ? break : nothing

                        # L ≥ K
                        KL = (L * (L + 1)) >> 1 + K
                        bkl = Lij*(ks-1) + bl
                        for js = 1:Lj
                            J = joff + js
                            bjkl = Li*(js-1) + bkl
                            for is = 1:Li
                                I = ioff + is
                                J < I ? break : nothing

                                IJ = (J * (J + 1)) >> 1 + I

                                idx = index2(IJ,KL) + 1
                                ∇x[idx] = -bufx[is + bjkl]
                                ∇y[idx] = -bufy[is + bjkl]
                                ∇z[idx] = -bufz[is + bjkl]
                                if KL ≥ IJ
                                    indexes[idx] = (I, J, K, L)
                                else
                                    indexes[idx] = (K, L, I, J)
                                end
                                mask[idx] = true
                            end
                        end
                    end
                end
            end
        #end #inbounds
        #end #spawn
    end #sync
    return indexes[mask], ∇x[mask], ∇y[mask], ∇z[mask]
end


text:
using Documenter, BitcoinPrimitives

makedocs(sitename="BitcoinPrimitives",
         doctest=true)

ext: jl


text:
##########################################
#   1-dimensional latent and 1PL pars    #
##########################################
"""
```julia
__probability(latent_val::Float64, parameters::Parameters1PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 1PL model at `latent_val` point.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent_val::Float64`** : Required. The point in the latent space in which compute the probability.
- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function __probability(latent_val::Float64, parameters::Parameters1PL)
   return 1 / (1 + _exp_c(parameters.b - latent_val))::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters1PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 1PL model at `Latent1D` point.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _probability(latent::Latent1D, parameters::Parameters1PL)
   return __probability(latent.val, parameters)::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters1PL, g_item::Vector{Float64}, g_latent::Vector{Float64})
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 1PL model at `Latent1D` point.
It updates the gradient vectors if they are not empty.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters1PL`** : Required. A 1-parameter logistic parameters object.
- **`g_item::Vector{Float64}`** : Optional. Vector of gradients with respect to item parameters.
- **`g_latent::Vector{Float64}`** : Optional. Vector of gradients with respect to the 1-dimensional latent.

# Output

A `Float64` scalar. Updates gradient vectors.
"""
function _probability(
    latent::Latent1D,
    parameters::Parameters1PL,
    g_item::Vector{Float64},
    g_latent::Vector{Float64},
)
    p = _probability(latent, parameters)

    if size(g_item, 1) > 0
        g_item .= p * (1 - p)
    end

    if size(g_latent, 1) > 0
        g_latent .= -p * (1 - p)
    end

    return p::Float64
end

##########################################
#   1-dimensional latent and 2PL pars    #
##########################################
"""
```julia
__probability(latent_val::Float64, parameters::Parameters2PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 2PL model at `latent_val` point.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent_val::Float64`** : Required. The point in the latent space in which compute the probability.
- **`parameters::Parameters2PL`** : Required. A 2-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function __probability(latent_val::Float64, parameters::Parameters2PL)
    return _sig_c(parameters.a * (latent_val - parameters.b))::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters2PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 2PL model at `Latent1D` point.
It follows the parametrization ``a(θ - b)``.

# Arguments
- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters2PL`** : Required. A 2-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _probability(latent::Latent1D, parameters::Parameters2PL)
    return __probability(latent.val, parameters)::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters2PL, g_item::Vector{Float64}, g_latent::Vector{Float64})
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 2PL model at `Latent1D` point.
It updates the gradient vectors if they are not empty.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters2PL`** : Required. A 2-parameter logistic parameters object.
- **`g_item::Vector{Float64}`** : Optional. Vector of gradients with respect to item parameters.
- **`g_latent::Vector{Float64}`** : Optional. Vector of gradients with respect to the 1-dimensional latent.

# Output

A `Float64` scalar. Updates the gradient vectors.

# Example

_Probability of a correct response is computed for each examinee and each item._

"""
function _probability(
    latent::Latent1D,
    parameters::Parameters2PL,
    g_item::Vector{Float64},
    g_latent::Vector{Float64},
)
    p = _probability(latent, parameters)

    if size(g_item, 1) > 0 || size(g_latent, 1) > 0
        p1p = p * (1 - p)
        if size(g_item, 1) > 0
            g_item .= [(latent.val - parameters.b) * p1p, -parameters.a * p1p]
        end

        if size(g_latent, 1) > 0
            g_latent .= parameters.a * p1p
        end
    end

    return p::Float64
end

##########################################
#   1-dimensional latent and 3PL pars    #
##########################################
"""
```julia
__probability(latent_val::Float64, parameters::Parameters3PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 3PL model at `latent_val` point.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent_val::Float64`** : Required. The point in the latent space in which compute the probability.
- **`parameters::Parameters3PL`** : Required. A 3-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function __probability(latent_val::Float64, parameters::Parameters3PL)
    return parameters.c +
    (1 - parameters.c) * (1 / (1 + _exp_c(-parameters.a * (latent_val - parameters.b))))::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters3PL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 3PL model at `Latent1D` point.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters3PL`** : Required. A 3-parameter logistic parameters object.

# Output

A `Float64` scalar.
"""
function _probability(latent::Latent1D, parameters::Parameters3PL)
    return __probability(latent.val, parameters)::Float64
end

"""
```julia
_probability(latent::Latent1D, parameters::Parameters3PL,  g_item::Vector{Float64}, g_latent::Vector{Float64})
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 3PL model at `Latent1D` point.
It updates the gradient vectors if they are not empty.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::Latent1D`** : Required. A 1-dimensional `Latent1D` latent variable.
- **`parameters::Parameters3PL`** : Required. A 3-parameter logistic parameters object.
- **`g_item::Vector{Float64}`** : Optional. Vector of gradients with respect to item parameters.
- **`g_latent::Vector{Float64}`** : Optional. Vector of gradients with respect to the 1-dimensional latent.

# Output

A `Float64` scalar. Updates the gradient vectors.
"""
function _probability(
    latent::Latent1D,
    parameters::Parameters3PL,
    g_item::Vector{Float64},
    g_latent::Vector{Float64},
)
    p = _probability(latent, parameters)

    if size(g_item, 1) > 0 || size(g_latent, 1) > 0
        q1c = (1 - p) / (1 - parameters.c)

        if size(g_item, 1) > 0
            g_item .= [
                (latent.val - parameters.b) * q1c * (p - parameters.c),
                -parameters.a * q1c * (p - parameters.c),
                q1c,
            ]
        end

        #by Kim's book
        if size(g_latent, 1) > 0
            g_latent .= parameters.a * (p - parameters.c) * q1c
        end

    end

    return p::Float64
end

##########################################
#   N-dimensional latent and NPL pars    #
##########################################
"""
```julia
__probability(latent_vals::Vector{Float64}, parameters::ParametersNPL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 3PL model at `latent_vals` points.
N-dimensional latent.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent_vals::Vector{Float64}`** : Required. The points in the latent space in which compute the probability.
- **`parameters::ParametersNPL`** : Required. A 2-parameter logistic parameters object suitbale for N-dimensional latent.

# Output

A `Float64` scalar.
"""
function __probability(latent_vals::Vector{Float64}, parameters::ParametersNPL)
    return _sig_c(latent_vals' * parameters.a - parameters.b)::Float64
end

"""
```julia
_probability(latent::LatentND, parameters::ParametersNPL)
```

# Description

It computes the probability (ICF) of a correct response for item `parameters` under the 3PL model at `Latent1D` point.
It updates the gradient vectors if they are not empty.
It follows the parametrization ``a(θ - b)``.

# Arguments

- **`latent::LatentND`** : Required. A N-dimensional latent variable.
- **`parameters::ParametersNPL`** : Required. A 2-parameter logistic parameters object suitable for a N-dimensional latent.

# Output

A `Float64` scalar.
"""
function _probability(latent::LatentND, parameters::ParametersNPL)
    return _sig_c(latent.val' * parameters.a - parameters.b)::Float64
end

##########################################
#   Externals  for examinees and items   #
##########################################

"""
```julia
probability(examinee::AbstractExaminee, item::AbstractItem)
```

# Description

It computes the probability (ICF) that an `examinee` answers correctly at `item`.

# Arguments

- **`examinee::AbstractExaminee`** : Required. An Examinee.
- **`item::AbstractItem`** : Required. An Item.

# Output

A `Float64` scalar.
"""
function probability(examinee::AbstractExaminee, item::AbstractItem)
    return _probability(examinee.latent, item.parameters)::Float64
end

"""
```julia
probability(examinees::Vector{<:AbstractExaminee}, items::Vector{<:AbstractItem})
```

# Description

It computes the probability (ICF) that a vector of `examinees` answers correctly at a vector of `items`.

# Arguments

- **`examinees::Vector{<:AbstractExaminee}`** : Required. A vector of `Examinee`s.
- **`items::Vector{<:AbstractItem}`** : Required. A vector of `Item`s.

# Output

A I x N `Float64` matrix.
"""
function probability(examinees::Vector{<:AbstractExaminee}, items::Vector{<:AbstractItem})
    return mapreduce(e -> map(i -> probability(e, i), items), hcat, examinees)::Matrix{Float64}
end

"""
```julia
probability(items::Vector{<:AbstractItem}, examinees::Vector{<:AbstractExaminee})
```

# Description

It computes the probability (ICF) that a vector of `examinees` answers correctly at a vector of `items`.

# Arguments

- **`examinees::Vector{<:AbstractExaminee}`** : Required. A vector of `Examinee`s.
- **`items::Vector{<:AbstractItem}`** : Required. A vector of `Item`s.

# Output

A I x N `Float64` matrix.
"""
function probability(items::Vector{<:AbstractItem}, examinees::Vector{<:AbstractExaminee})
    return probability(examinees, items)::Matrix{Float64}
end

##########################################
#         Externals for matrix data      #
##########################################
"""
```julia
probability(parameters_matrix::Matrix{Float64}, latent_matrix::Matrix{Float64})
```

# Description

It computes the probability (ICF) of a correct response for item parameters and latents values provided in matrix form.
Not suitable for 3PL models, for such a kind of model use `probability_3PL()`.
It follows the parametrization ``a \theta - b``.

# Arguments

- **`parameters_matrix::Matrix{Float64}`** : Required. A ``I x (\text{n_latents}+ 1)`` matrix with item parameters. intercept (b) must be in first column, latents coefficients (``a_j``) in next columns ``(2, ..., \text{n_latents} + 1)``.
- **`latent_matrix::Matrix{Float64}`** : Required. A ``\text{n_latents} \times N`` or ``(\text{n_latents} + 1) x N`` matrix with latents values.

# Outputs

A ``I x N`` `Float64` matrix.
"""
function probability(parameters_matrix::Matrix{Float64}, latents_matrix::Matrix{Float64})
    if size(latents_matrix, 1) == (size(parameters_matrix, 2)-1)
        latents_matrix = vcat(.-ones(Float64, size(latents_matrix, 2))', latents_matrix)
    elseif size(latents_matrix, 1) != size(parameters_matrix, 2)
        error("The latents matrix must have n_pars-1 columns (the intercept is automatically added) or n_pars columns.")
    end
    return _sig_c.(_gemmblasAB(parameters_matrix, latents_matrix))::Matrix{Float64}
end

"""
```julia
probability_3PL(parameters_matrix::Matrix{Float64}, latent_matrix::Matrix{Float64})
```

# Description

Only for models which has guessing parameter (c) in last row of parameters_matrix. It computes the probability (ICF) of a correct response for item parameters and latents values provided in matrix form.
It follows the parametrization ``a \theta - b``.

# Arguments

- **`parameters_matrix::Matrix{Float64}`** : Required. A ``I x (\text{n_latents}+ 1)`` matrix with item parameters. intercept (b) must be in first column, latents coefficients (``a_j``) in next columns ``(2, ..., \text{n_latents} + 1)``.
- **`latent_matrix::Matrix{Float64}`** : Required. A ``\text{n_latents} \times N`` matrix with latents values.

# Output

A ``I x N`` `Float64` matrix.
"""
function probability_3PL(
    parameters_matrix::Matrix{Float64},
    latents_matrix::Matrix{Float64},
)
    parameters_core = parameters_matrix[:, 1:(end-1)]
    if size(latents_matrix, 1) == (size(parameters_core, 2)-1)
        latents_matrix = vcat(ones(Float64, size(latents_matrix, 2))', latents_matrix)
    elseif size(latents_matrix, 1) != size(parameters_core, 2)
        error("The latents matrix must have n_pars-2 columns (the intercept is automatically added) or n_pars-1 columns.")
    end
    ret = _sig_c.(_gemmblasAB(parameters_core, latents_matrix))
    return _matrix_rows_vec(ret, parameters_matrix[:, end], (x, y) -> y + (1 - y) * x)::Matrix{Float64}
end


text:
@testset "Anonymous Functions" begin

    add_curry(n) = (x) -> x + n

    add1 = mocked(add_curry(1))

    # make sure that normal call works
    @test add1(1) == 2

    # able to apply patch?
    apply(add1 => (x) -> x + 10) do
        @test add1(1) == 11
    end

    # test conditional patch
    apply(add1 => (x) -> x > 5 ? Fallback() : x) do
        @test add1(1) == 1  # no fallback
        @test add1(6) == 7  # fallback
    end

    # anonymous functions have no multiple dispatch so we can patch
    # with anything! This should not be a normal use case, however.
    apply(add1 => (x,y,z) -> x + y + z) do
        @test_throws MethodError add1(1)
        @test add1(1, 2, 3) == 6
    end

end


text:
module brml

push!(LOAD_PATH, joinpath(pwd(), "src"))

# push all subdirectories from src
map(d -> push!(LOAD_PATH, joinpath(pwd(), "src", d)),
    filter(d -> isdir(joinpath("src", d)), readdir("src")))

push!(LOAD_PATH, joinpath(pwd(), "Demos"))

println("\nAvailable Demos:\n")
demos=readdir("Demos")
for i=1:length(demos)
#    try
        if demos[i][1:4]=="demo" && demos[i][end-2:end]==".jl"
            println(demos[i])
            reload(demos[i])
        end
#    end
end

using Reexport

@reexport using BrmlFunctions, HMM, PotentialInference, GraphPlot

end


text:
# define our own method to avoid type piracy with Base.showable
_showable(a::AbstractVector{<:MIME}, x) = any(m -> showable(m, x), a)
_showable(m, x) = showable(m, x)

"""
A vector of MIME types (or vectors of MIME types) that IJulia will try to
render. IJulia will try to render every MIME type specified in the first level
of the vector. If a vector of MIME types is specified, IJulia will include only
the first MIME type that is renderable (this allows for the expression of
priority and exclusion of redundant data).

For example, since "text/plain" is specified as a first-child of the array,
IJulia will always try to include a "text/plain" representation of anything that
is displayed. Since markdown and html are specified within a sub-vector, IJulia
will always try to render "text/markdown", and will only try to render
"text/html" if markdown isn't possible.
"""
const ijulia_mime_types = Vector{Union{MIME, AbstractVector{MIME}}}([
    MIME("text/plain"),
    MIME("image/svg+xml"),
    [MIME("image/png"),MIME("image/jpeg")],
    [
        MIME("text/markdown"),
        MIME("text/html"),
    ],
    MIME("text/latex"),
])

"""
MIME types that when rendered (via stringmime) return JSON data. See
`ijulia_mime_types` for a description of how MIME types are selected.

This is necessary to embed the JSON as is in the displaydata bundle (rather than
as stringify'd JSON).
"""
const ijulia_jsonmime_types = Vector{Union{MIME, Vector{MIME}}}([
    [[MIME("application/vnd.vegalite.v$n+json") for n in 4:-1:2]...,
    [MIME("application/vnd.vega.v$n+json") for n in 5:-1:3]...],
    MIME("application/vnd.dataresource+json"), MIME("application/vnd.plotly.v1+json")
])

register_mime(x::Union{MIME, Vector{MIME}})= push!(ijulia_mime_types, x)
register_mime(x::AbstractVector{<:MIME}) = push!(ijulia_mime_types, Vector{Mime}(x))
register_jsonmime(x::Union{MIME, Vector{MIME}}) = push!(ijulia_jsonmime_types, x)
register_jsonmime(x::AbstractVector{<:MIME}) = push!(ijulia_jsonmime_types, Vector{Mime}(x))

# return a String=>Any dictionary to attach as metadata
# in Jupyter display_data and pyout messages
metadata(x) = Dict()

"""
Generate the preferred MIME representation of x.

Returns a tuple with the selected MIME type and the representation of the data
using that MIME type.
"""
function display_mimestring(mime_array::Vector{MIME}, x)
    for m in mime_array
        if _showable(m, x)
            return display_mimestring(m, x)
        end
    end
    error("No displayable MIME types in mime array.")
end

display_mimestring(m::MIME, x) = (m, limitstringmime(m, x))

# text/plain output must have valid Unicode data to display in Jupyter
function display_mimestring(m::MIME"text/plain", x)
    s = limitstringmime(m, x)
    return m, (isvalid(s) ? s : "(binary data)")
end

"""
Generate the preferred json-MIME representation of x.

Returns a tuple with the selected MIME type and the representation of the data
using that MIME type (as a `JSONText`).
"""
function display_mimejson(mime_array::Vector{MIME}, x)
    for m in mime_array
        if _showable(m, x)
            return display_mimejson(m, x)
        end
    end
    error("No displayable MIME types in mime array.")
end

display_mimejson(m::MIME, x) = (m, JSON.JSONText(limitstringmime(m, x)))

"""
Generate a dictionary of `mime_type => data` pairs for all registered MIME
types. This is the format that Jupyter expects in display_data and
execute_result messages.
"""
function display_dict(x)
    data = Dict{String, Union{String, JSONText}}()
    for m in ijulia_mime_types
        try
            if _showable(m, x)
                mime, mime_repr = display_mimestring(m, x)
                data[string(mime)] = mime_repr
            end
        catch
            if m == MIME("text/plain")
                rethrow() # text/plain is required
            end
        end
    end

    for m in ijulia_jsonmime_types
        try
            if _showable(m, x)
                mime, mime_repr = display_mimejson(m, x)
                data[string(mime)] = mime_repr
            end
        catch
        end
    end

    return data

end

# queue of objects to display at end of cell execution
const displayqueue = Any[]

# remove x from the display queue
function undisplay(x)
    i = findfirst(isequal(x), displayqueue)
    i !== nothing && splice!(displayqueue, i)
    return x
end

import Base: ip_matches_func

function show_bt(io::IO, top_func::Symbol, t, set)
    # follow PR #17570 code in removing top_func from backtrace
    eval_ind = findlast(addr->ip_matches_func(addr, top_func), t)
    eval_ind !== nothing && (t = t[1:eval_ind-1])
    Base.show_backtrace(io, t)
end

# wrapper for showerror(..., backtrace=false) since invokelatest
# doesn't support keyword arguments.
showerror_nobt(io, e, bt) = showerror(io, e, bt, backtrace=false)

# return the content of a pyerr message for exception e
function error_content(e, bt=catch_backtrace();
                       backtrace_top::Symbol=SOFTSCOPE[] ? :softscope_include_string : :include_string,
                       msg::AbstractString="")
    tb = map(String, split(sprint(show_bt,
                                        backtrace_top,
                                        bt, 1:typemax(Int)),
                                  "\n", keepempty=true))

    ename = string(typeof(e))
    evalue = try
        # Peel away one LoadError layer that comes from running include_string on the cell
        isa(e, LoadError) && (e = e.error)
        sprint((io, e, bt) -> invokelatest(showerror_nobt, io, e, bt), e, bt; context=InlineIOContext(stderr))
    catch
        "SYSTEM: show(lasterr) caused an error"
    end
    pushfirst!(tb, evalue) # fperez says this needs to be in traceback too
    if !isempty(msg)
        pushfirst!(tb, msg)
    end
    Dict("ename" => ename, "evalue" => evalue,
                 "traceback" => tb)
end

#######################################################################


text:
@enum Output begin
    json = 1
    yaml = 2
    text = 3
    table = 4
end

# TODO: Should we make region an enum as well?
mutable struct AWSConfig
    creds::AWSCredentials
    region::String
    output::Output
end

function AWSConfig(;
        profile=nothing,
        creds=AWSCredentials(profile=profile),
        region=get(ENV, "AWS_DEFAULT_REGION", "us-east-1"),
        output=json
    )
    return AWSConfig(creds, region, output)
end

function Base.show(io::IO, config::AWSConfig)
    println(io, "
    AWSConfig(
        creds=AWSCredentials(
            access_key=$(config.creds.access_key_id)
            secret_key=$(config.creds.secret_key[1:4])...
        ),
        region=$(config.region),
        output=$(config.output)
    )")
end

max_stars_repo_name: mattBrzezinski/Julia-AWS-Status

max_issues_repo_name: mattBrzezinski/Julia-AWS-Status

max_forks_repo_name: mattBrzezinski/Julia-AWS-Status


text:
function Uniform(prog::Program, name::String)
	u = ccall( (:glGetUniformLocation, lib), Uniform,
		(Program, Ptr{GLchar}), prog, bytestring(name))
	if u.location < 0
		GetError()
		error("uniform ", name, " not found")
	else
		return u
	end
end

import Base.write
importall ImmutableArrays

for (T,t) in ((Float32, "f"), (Int32, "i"), (Uint32, "ui"))
	@eval begin
		# Write single vector
		write(u::Uniform, val::$T) =
			ccall(($(string("glUniform1", t)), lib), Void,
				(Uniform, $T), u, val)
		write(u::Uniform, val::Vector2{$T}) =
			ccall(($(string("glUniform2", t, "v")), lib), Void,
				(Uniform, GLsizei, Ptr{Vector2{$T}}), u, 1, &val)
		write(u::Uniform, val::Vector3{$T}) =
			ccall(($(string("glUniform3", t, "v")), lib), Void,
				(Uniform, GLsizei, Ptr{Vector3{$T}}), u, 1, &val)
		write(u::Uniform, val::Vector4{$T}) =
			ccall(($(string("glUniform4", t, "v")), lib), Void,
				(Uniform, GLsizei, Ptr{Vector4{$T}}), u, 1, &val)

		# Write matrix
		write(u::Uniform, mat::Matrix2x2{$T}) =
			ccall(($(string("glUniformMatrix2", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix2x2{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix2x3{$T}) =
			ccall(($(string("glUniformMatrix2x3", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix2x3{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix2x4{$T}) =
			ccall(($(string("glUniformMatrix2x4", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix2x4{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix3x2{$T}) =
			ccall(($(string("glUniformMatrix3x2", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix3x2{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix3x3{$T}) =
			ccall(($(string("glUniformMatrix3", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix3x3{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix3x4{$T}) =
			ccall(($(string("glUniformMatrix3x4", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix3x4{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix4x2{$T}) =
			ccall(($(string("glUniformMatrix4x2", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix4x2{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix4x3{$T}) =
			ccall(($(string("glUniformMatrix4x3", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix4x3{$T}}),
				u, 1, false, &mat)
		write(u::Uniform, mat::Matrix4x4{$T}) =
			ccall(($(string("glUniformMatrix4", t, "v")), lib), Void,
				(Uniform, GLsizei, GLboolean, Ptr{Matrix4x4{$T}}),
				u, 1, false, &mat)
	end
end

write(u::Uniform, val::Bool) = write(u, convert(Uint32, val))


text:
module DesertIslandDisk

export  @tostring
export  @remindme
export  printschema
export  df_transpose
export  read_partitioned_dataframe_csv
export  browserview_html
export  stop_script_exection_and_return_to_repl_or_shell

greet() = print("Hello World!")

include("tostring.jl")
include("remindme.jl")
include("printschema.jl")
include("df_transpose.jl")
include("repl_script_end_gracefully.jl")
include("read_partitioned_dataframe_csv.jl")
include("browserview_html.jl")

end # module

max_stars_repo_name: mahiki/DesertIslandDisk

max_issues_repo_name: mahiki/DesertIslandDisk

max_forks_repo_name: mahiki/DesertIslandDisk


text:
using Thyme, Compose, Gadfly

module Scallop

using Thyme, Compose, Gadfly

export @scalloptest

include("backtest/backtest.jl")
include("book/book.jl")
include("chart/chart.jl")
include("metrics/metrics.jl")
include("underlying/underlying.jl")
include("walkforward/walkforward.jl")

include("scalloptest.jl")

end


text:
struct EvaluatedMatrices
    A1::Array{Float64,2}
    A2::Array{Float64,2}
    A3::Array{Float64,2}
    A4::Array{Float64,2}
    B1::Array{Float64,2}
    B2::Array{Float64,2}
    B3::Array{Float64,2}
    B4::Array{Float64,2}
    B5::Array{Float64,2}
    C1_states::Array{Float64,2}
    C2_states::Array{Float64,2}
    C_shocks::Array{Float64,2}
    C_constant::Array{Float64,2}
end

struct StaticSub
    YA1::Array{Float64,2}
    YA2::Array{Float64,2}
    YA3::Array{Float64,2}
    YA4::Array{Float64,2}
    YB1::Array{Float64,2}
    YB2::Array{Float64,2}
    YB3::Array{Float64,2}
    YB4::Array{Float64,2}
    YC1::Array{Float64,2}
    YC2::Array{Float64,2}
end

struct ExoSolution
    W::Array{Float64,2}
    Z::Array{Float64,2}
end

struct DynamicSolution
    X::Array{Float64,2}
    J::Array{Float64,2}
end

struct StaticSolution
    Y::Array{Float64,2}
end

struct ShocksSolution
    orderedME::Array{Float64,2}
end

struct ConstantSolution
    orderedMC::Array{Float64,2}
end

struct RSDSGESolution
    MX::Array{Float64,3}
    ME::Array{Float64,3}
    MC::Array{Float64,3}
end

function EvaluatedMatrices(model,r)
    A1 = createstacked(model,model.perturbationsystem.A1,r)
    A2 = createstacked(model,model.perturbationsystem.A2,r)
    A3 = createstacked(model,model.perturbationsystem.A3,r)
    A4 = createstacked(model,model.perturbationsystem.A4,r)
    B1 = createsummed(model,model.perturbationsystem.B1,r)
    B2 = createsummed(model,model.perturbationsystem.B2,r)
    B3 = createsummed(model,model.perturbationsystem.B3,r)
    B4 = createsummed(model,model.perturbationsystem.B4,r)
    B5 = createsummed(model,model.perturbationsystem.B5,r)
    C1 = createsummed(model,model.perturbationsystem.C1_states,r)
    C2 = createsummed(model,model.perturbationsystem.C2_states,r)
    C_shocks = createsummed(model,model.perturbationsystem.C_shocks,r)
    C_constant = createstacked(model,model.perturbationsystem.C_constant,r)
    C_constant_p = createstacked(model,model.perturbationsystem.C_constant_p,r)
    multiplybyhats!(C_constant,model,r,true)
    multiplybyhats!(C_constant_p,model,r,false)
    # Collapse the constant matrices
    C_constant += C_constant_p
    C_constant = sum(C_constant,2)
    return EvaluatedMatrices(A1,A2,A3,A4,B1,B2,B3,B4,B5,C1,C2,C_shocks,C_constant)
end

function evaluatematrices(model)
    out = Array{EvaluatedMatrices}(0)
    for r in 1:model.meta.numregimes
        push!(out,EvaluatedMatrices(model,r))
    end
    return out
end

function createarglistvalues(model,r,rp)
    arg = Array{Float64}(1)
    arg = [model.transmatrix.values[r,rp]] # Transition prob first
    append!(arg,model.steadystate.values) # Add in the steady state values
    # Non-switching vars first
    append!(arg,model.parameters.values_withbar[r, .!model.parameters.isswitching])
    # For parameters that affect SS, just need the bar version
    append!(arg,model.parameters.values_withbar[r,model.parameters.affectsSS])
    # Finally, for parameters that switch but _don't_ affect SS, need both R and and RP version
    append!(arg,model.parameters.values_withbar[r,model.parameters.isswitching .& .!model.parameters.affectsSS])
    append!(arg,model.parameters.values_withbar[rp,model.parameters.isswitching .& .!model.parameters.affectsSS])
    return arg
end

function evaluate(model,fmatrix,r,rp)
    args = createarglistvalues(model,r,rp)
    out = zeros(Float64,size(fmatrix)...)
    if nnz(fmatrix) > 0
        R,C,F = findnz(fmatrix)
        for (r,c,f) in zip(R,C,F)
            out[r,c] = f(args...)
        end
    end
    return out
end

function createsummed(model,fmatrix,r)
    out = zeros(Float64,size(fmatrix)...)
    if !isempty(fmatrix)
        for rp in 1:model.meta.numregimes
            out += evaluate(model,fmatrix,r,rp)
        end
    end
    return out
end

function createstacked(model,fmatrix,r)
    if !isempty(fmatrix)
        out = zeros(size(fmatrix,1),0)
        for rp in 1:model.meta.numregimes
            out = hcat(out,evaluate(model,fmatrix,r,rp))
        end
    else
        out = zeros(size(fmatrix,1),size(fmatrix,2)*model.meta.numregimes)
    end
    return out
end

function multiplybyhats!(C_constant,model,r,today)
    # This function goes through the C_constant array and multiplies each element by it's relevant hat
    # since that is the derivative of the parameter wrt perturbationconstant (and the element already in the matrix
    # is the derivative of the system wrt the parameter
    for rp in 1:model.meta.numregimes
        i = 0
        for p in 1:model.meta.numparameters
            if model.parameters.affectsSS[p]
                i+=1
                x = (rp-1)*countnz(model.parameters.affectsSS) + i
                for eq in 1:size(C_constant,1)
                    if today
                        C_constant[eq,x] = C_constant[eq,x]*model.parameters.decomposition[p].hat[r]
                    else
                        C_constant[eq,x] = C_constant[eq,x]*model.parameters.decomposition[p].hat[rp]
                    end
                end
            end
        end
    end
    return
end

function substatics(model,m,r)
    A1 = m.A1[model.equations.isstatic,:]
    A2 = m.A2[model.equations.isstatic,:]
    A3 = m.A3[model.equations.isstatic,:]
    A4 = m.A4[model.equations.isstatic,:]
    B1 = m.B1[model.equations.isstatic,:]
    B2 = m.B2[model.equations.isstatic,:]
    B3 = m.B3[model.equations.isstatic,:]
    B4 = m.B4[model.equations.isstatic,:]
    B5 = m.B5[model.equations.isstatic,:]
    C1 = m.C1_states[model.equations.isstatic,:]
    C2 = m.C2_states[model.equations.isstatic,:]
    return StaticSub(-B5\A1,-B5\A2,-B5\A3,-B5\A4,-B5\B1,-B5\B2,-B5\B3,-B5\B4,-B5\C1,-B5\C2)
end

function substatics(model,m)
    solutions = Array{StaticSub}(model.meta.numregimes)
    for r in 1:model.meta.numregimes
        solutions[r] = substatics(model,m[r],r)
    end
    return solutions
end

function solveexos(model,m,staticsol,r)
    # These shouldn't be needed
    A1 = m.A1[model.equations.isexo,:]
    A2 = m.A2[model.equations.isexo,:]
    A3 = m.A3[model.equations.isexo,:]
    A4 = m.A4[model.equations.isexo,:]
    B1 = m.B1[model.equations.isexo,:]
    B2 = m.B2[model.equations.isexo,:]

    B3 = m.B3[model.equations.isexo,:]
    B4 = m.B4[model.equations.isexo,:]
    B5 = m.B5[model.equations.isexo,:]
    C1 = m.C1_states[model.equations.isexo,:]
    C2 = m.C2_states[model.equations.isexo,:]

    A1 += B5*staticsol.YA1
    A2 += B5*staticsol.YA2
    A3 += B5*staticsol.YA3
    A4 += B5*staticsol.YA4
    B1 += B5*staticsol.YB1
    B2 += B5*staticsol.YB2
    B3 += B5*staticsol.YB3
    B4 += B5*staticsol.YB4
    C1 += B5*staticsol.YC1
    C2 += B5*staticsol.YC2

    # These should all be non zero. This is just double checking that. I'll remove these ops for production, since
    # they carry performance cost.
    if countnz(A1)>0
        error("Non-zero A1 entries in exo block")
    end
    if countnz(A2)>0
        error("Non-zero A2 entries in exo block")
    end
    if countnz(A3)>0
        error("Non-zero A3 entries in exo block")
    end
    if countnz(A4)>0
        error("Non-zero A4 entries in exo block")
    end
    if countnz(B1)>0
        error("Non-zero B1 entries in exo block")
    end
    if countnz(B3)>0
        println(full(B3))
        error("Non-zero B3 entries in exo block")
    end

    sol = -hcat(B2,B4)\hcat(C1,C2)
    W = sol[1:countnz(model.perturbationsystem.isW),:]
    Z = sol[(1+countnz(model.perturbationsystem.isW)):end,:]
    return ExoSolution(W,Z)
end

function solveexos(model,m,staticsubs)
    solutions = Array{ExoSolution}(model.meta.numregimes)
    for r in 1:model.meta.numregimes
        solutions[r] = solveexos(model,m[r],staticsubs[r],r)
    end
    return solutions
end

function arrangeexomatrices(model,exosols)
    nz = countnz(model.perturbationsystem.isZ)
    nw = countnz(model.perturbationsystem.isW)
    nx = countnz(model.perturbationsystem.isX)
    Wz = zeros(Float64,0,nz)
    Zz = zeros(Float64,0,nz)
    for r in 1:model.meta.numregimes
        Wz = vcat(Wz,exosols[r].W[:,(1+nx):end])
        Zz = vcat(Zz,exosols[r].Z[:,(1+nx):end])
    end
    return Wz,Zz
end

function createsystem(model,m,staticsol,exosol,Wz,Zz,dynamics,r)
    A1 = m.A1[dynamics,:]
    A2 = m.A2[dynamics,:]
    A3 = m.A3[dynamics,:]
    A4 = m.A4[dynamics,:]
    B1 = m.B1[dynamics,:]
    B2 = m.B2[dynamics,:]
    B3 = m.B3[dynamics,:]
    B4 = m.B4[dynamics,:]
    B5 = m.B5[dynamics,:]
    C1 = m.C1_states[dynamics,:]
    C2 = m.C2_states[dynamics,:]

    A1 += B5*staticsol.YA1
    A2 += B5*staticsol.YA2
    A3 += B5*staticsol.YA3
    A4 += B5*staticsol.YA4
    B1 += B5*staticsol.YB1
    B2 += B5*staticsol.YB2
    B3 += B5*staticsol.YB3
    B4 += B5*staticsol.YB4
    C1 += B5*staticsol.YC1
    C2 += B5*staticsol.YC2

    return  A1*(model.perturbationsystem.Jx*model.perturbationsystem.Xs[:,:,r] + model.perturbationsystem.Jz*exosol.Z) + A2*(Wz*exosol.Z) + A3*(model.perturbationsystem.Xx*model.perturbationsystem.Xs[:,:,r] + model.perturbationsystem.Xz*exosol.Z) + A4*(Zz*exosol.Z) + B1*model.perturbationsystem.Js[:,:,r] + B2*exosol.W + B3*model.perturbationsystem.Xs[:,:,r] + B4*exosol.Z + hcat(C1, C2)
end

function createsystem(model,m,statics,exos,Wx,Wz)
    dynamics = .!model.equations.isstatic .& .!model.equations.isexo
    system = Array{SymPy.Sym}(countnz(dynamics)*model.meta.numregimes,model.meta.numstates)
    for r in 1:model.meta.numregimes
        system[(1+(r-1)*countnz(dynamics)):r*countnz(dynamics),:] = createsystem(model,m[r],statics[r],exos[r],Wx,Wz,dynamics,r)
    end
    return reshape(system,length(system)) # The reshape is needed because sympy can only take 1d arrays. It doesn't change anything important.
end

function markequations(model,block)
    dynamics = .!model.equations.isstatic .& .!model.equations.isexo
    forblock = falses(0,model.meta.numstates)
    removefromsystem = falses(0,model.meta.numstates)
    keepeq = falses(countnz(dynamics),model.meta.numstates)
    keepeq[block.equations[dynamics],block.states] = true
    removeeq = falses(countnz(dynamics),model.meta.numstates)
    removeeq[block.equations[dynamics],:] = true
    for r in 1:model.meta.numregimes
        forblock = vcat(forblock,keepeq)
        removefromsystem = vcat(removefromsystem,removeeq)
    end
    return forblock,removefromsystem
end

function replaceblocks(model,block,system)
    # It can still depend on any or all of the exos, so we don't sub any of those out.
    tosubout = model.perturbationsystem.Xs[block.vars[model.perturbationsystem.isX],.!block.states,:]
    tosubout = reshape(tosubout,length(tosubout))
    tosubout2 = model.perturbationsystem.Js[block.vars[model.perturbationsystem.isJ],.!block.states,:]
    tosubout = vcat(tosubout,reshape(tosubout2,length(tosubout2)))

    sympyzero = Sym(0.0)
    zeroarray = Array{SymPy.Sym}(length(tosubout))
    fill!(zeroarray,sympyzero)
    replacearg = createpairedtuplelist(tosubout,zeroarray)
    if !isempty(replacearg)>0
        return squeeze(SymPy.subs(system,replacearg...),2),Dict{SymPy.Sym,SymPy.Sym}(zip(tosubout,zeroarray)) # I squeeze because subs returns a 2d with singleton trailing for some reason
    else
        return system,Dict{SymPy.Sym,SymPy.Sym}()
    end
end

function separateblocks(model,system)
    restricteddecisionrules = Dict{SymPy.Sym,SymPy.Sym}() # To hold the decision rule solutions that are restricted to zero because of block exogeneity
    blocks = Array{Array{SymPy.Sym,1}}(0)
    usedinblock = falses(length(system))
    for block in model.perturbationsystem.blocks
        blockwhich,removewhich = markequations(model,block)
        blockwhich = reshape(blockwhich,length(system))
        removewhich = reshape(removewhich,length(system))
        usedinblock = usedinblock .| removewhich
        # Now, replace any decision rules for states that don't exist in the block by zero
        blocksystem,blockdecisionrules = replaceblocks(model,block,system[blockwhich])
        push!(blocks,blocksystem)
        merge!(restricteddecisionrules,blockdecisionrules)
    end
    return blocks,system[.!usedinblock],restricteddecisionrules
end

function subsdict(target,dictionary)
    subsargs = ()
    if !isempty(dictionary) > 0
        for (key,val) in dictionary
            subsargs = (subsargs...,(key,val))
        end
        return SymPy.subs(target,subsargs...)
    else
        return target
    end
end

function solveblocks(model,blocks,system,restricteddecisionrules,verbose::Bool,method)
    blockdict = Dict{SymPy.Sym,SymPy.Sym}() # To hold the solutions to all the blocks
    if length(blocks)>0
        printlnif(verbose,"Solving the $(length(blocks)) model sub-blocks:")
        for (i,block) in enumerate(blocks)
            printif(verbose,"Block $(i) of $(length(blocks)) ($(length(block)) equations)...")
            if method == 1
                sol = solve_withfallback(block)
            elseif method == 2
                sol = numericsolve(reshape(block,length(block)))
            else
                error("Option $(method) is not a valid solution method.")
            end
            if !isa(sol,Dict)
                error("Handling multiple solutions for blocks not yet implemented.")
            end
            merge!(blockdict,sol)
            printlnif(verbose,"done.")
        end
        printif(verbose,"Subbing block solutions into the full system...")
        system = subsdict(system,blockdict)
        system = subsdict(system,restricteddecisionrules)
        if ndims(system) > 1 && !isempty(system)
             # Squeeze because subs returns trailing singleton
            system = squeeze(system,2)
        end
        printlnif(verbose,"done.")
    end
    return system,blockdict
end

function convert2rationals(system)
    for (i,expression) in enumerate(system)
        system[i]=SymPy.nsimplify(expression)
    end
    return system
end

function systemissolvable(system)
    if length(system) == length(SymPy.free_symbols(system))
        return true
    else
        println(system)
        println(SymPy.free_symbols(system))
        error("Poorly specified system. You have $(length(system)) equations for $(length(SymPy.free_symbols(system))) unknowns.")
    end
end

function solve_withfallback(system)
    # This function wraps solve to first try solving using floats. If that fails because sympy is sometimes inaccurate with floats
    # it then converts to rationals and tries again.
    if !isempty(system)
        if systemissolvable(system)
            try
                solution = SymPy.solve(system,method="f5b")
                return solution
            catch
                warn("Solve failed. Converting to rationals and trying again.")
                system = convert2rationals(system)
                solution = SymPy.solve(system,method="f5b")
                return solution
            end
        end
    else
        # Sub-blocks account for the whole system. This is a no-op
        return Dict{SymPy.Sym,SymPy.Sym}()
    end
end

function lambdifysystem(system,args)
    out = Array{Function}(length(system))
    for (i,expr) in enumerate(system)
        out[i] = SymPy.lambdify(expr,args)
    end
    return out
end

function evaluatearray(x,farray)
    out = similar(x)
    for (i,f) in enumerate(farray)
        out[i] = f(x...)
    end
    return out
end

function numericsolve(system)
    if !isempty(system)
        symbols = SymPy.free_symbols(system)
        tosolve = lambdifysystem(system,symbols)
        fwrap = x->evaluatearray(x,tosolve)
        nlout = NLsolve.nlsolve(NLsolve.not_in_place(fwrap), zeros(length(symbols)))
        if !NLsolve.converged(nlout)
            error("Numerical solver failed to converge.")
        end
        dynamicsolutions = Dict{SymPy.Sym,Float64}(zip(symbols,nlout.zero))
        return dynamicsolutions
    else
        # Sub-blocks account for the whole system. This is a no-op
        return Dict{SymPy.Sym,SymPy.Sym}()
    end
end

function hasimaginary(sol)
    for val in values(sol)
        if imag(val) > 1e-12
            return true
        end
    end
    return false
end

function dict2array(dict,indices)
    out = Array{valtype(dict)}(size(indices))
    for i in eachindex(indices)
        out[i] = dict[indices[i]]
    end
    return out
end

function formatdynamicsolution(model,solution::Dict)
    out = Array{DynamicSolution}(model.meta.numregimes)
    for r in 1:model.meta.numregimes
        Xsols = dict2array(solution,model.perturbationsystem.Xs[:,:,r])
        Jsols = dict2array(solution,model.perturbationsystem.Js[:,:,r])
        out[r] = DynamicSolution(Xsols,Jsols)
    end
    return out
end

function arrangesolutionmatrices(model,exosolution,dynamicsolution)
    Xx = Array{Float64}(0,countnz(model.perturbationsystem.isX[model.vars.isstate]))
    Xz = Array{Float64}(0,countnz(model.perturbationsystem.isZ[model.vars.isstate]))
    Jx = Array{Float64}(0,countnz(model.perturbationsystem.isX[model.vars.isstate]))
    Jz = Array{Float64}(0,countnz(model.perturbationsystem.isZ[model.vars.isstate]))
    for r in 1:model.meta.numregimes
        Xx = vcat(Xx,dynamicsolution[r].X[:,model.perturbationsystem.isX[model.vars.isstate]])
        Jx = vcat(Jx,dynamicsolution[r].J[:,model.perturbationsystem.isX[model.vars.isstate]])
        Xz = vcat(Xz,dynamicsolution[r].X[:,model.perturbationsystem.isZ[model.vars.isstate]])
        Jz = vcat(Jz,dynamicsolution[r].J[:,model.perturbationsystem.isZ[model.vars.isstate]])
    end
    return Xx,Xz,Jx,Jz
end

function solvestatics(model,staticsols,exosols,dynamicsols,Xx,Xz,Jx,Jz,Wz,Zz)
    # This function uses the exo and dynamic solutions to return the actual solutions for static variables, rather
    # than 'solutions' that just express them as a function of other variables
    out = Array{StaticSolution}(model.meta.numregimes)
    for (r,static,exo,dynamic) in zip(1:model.meta.numregimes,staticsols,exosols,dynamicsols)
        out[r] = StaticSolution(static.YA1*(Jx*dynamic.X + Jz*exo.Z) + static.YA2*(Wz*exo.Z) + static.YA3*(Xx*dynamic.X + Xz*exo.Z) + static.YA4*(Zz*exo.Z) + static.YB1*dynamic.J + static.YB2*exo.W + static.YB3*dynamic.X + static.YB4*exo.Z + hcat(static.YC1, static.YC2))
    end
    return out
end

function convert2decision(model,dynamics,exos,statics)
    decision = zeros(Float64,model.meta.numvars,model.meta.numvars,model.meta.numregimes)
    Xindices = model.vars.states_indices[model.perturbationsystem.isX[model.vars.isstate]]
    Zindices = model.vars.states_indices[model.perturbationsystem.isZ[model.vars.isstate]]
    indices = vcat(Xindices,Zindices)

    for r in 1:model.meta.numregimes
        decision[model.perturbationsystem.isY,indices,r] = statics[r].Y
        decision[model.perturbationsystem.isJ,indices,r] = dynamics[r].J
        decision[model.perturbationsystem.isX,indices,r] = dynamics[r].X
        decision[model.perturbationsystem.isZ,indices,r] = exos[r].Z
        decision[model.perturbationsystem.isW,indices,r] = exos[r].W
    end

    return decision
end

function solveshocks(model,m,Xx,Xz,Jx,Jz,Wz,Zz,r)
    return ShocksSolution(-hcat(m.B1, m.B2, m.A1*Jx+m.A3*Xx+m.B3, m.A1*Jz+m.A2*Wz+m.A3*Xz+m.A4*Zz+m.B4, m.B5)\m.C_shocks)
end

function rearrangeshocks(model,ordered)
    # J,W,X,Z,Y
    shockssol = zeros(Float64,model.meta.numvars,model.meta.numshocks,model.meta.numregimes)
    for r in 1:model.meta.numregimes
        offset = 1
        shockssol[model.perturbationsystem.isJ,:,r] = ordered[r].orderedME[offset:(offset-1+countnz(model.perturbationsystem.isJ)),:]
        offset += countnz(model.perturbationsystem.isJ)
        shockssol[model.perturbationsystem.isW,:,r] = ordered[r].orderedME[offset:(offset-1+countnz(model.perturbationsystem.isW)),:]
        offset += countnz(model.perturbationsystem.isW)
        shockssol[model.perturbationsystem.isX,:,r] = ordered[r].orderedME[offset:(offset-1+countnz(model.perturbationsystem.isX)),:]
        offset += countnz(model.perturbationsystem.isX)
        shockssol[model.perturbationsystem.isZ,:,r] = ordered[r].orderedME[offset:(offset-1+countnz(model.perturbationsystem.isZ)),:]
        offset += countnz(model.perturbationsystem.isZ)
        shockssol[model.perturbationsystem.isY,:,r] = ordered[r].orderedME[offset:(offset-1+countnz(model.perturbationsystem.isY)),:]
    end
    return shockssol
end

function solveshocks(model,m,Xx,Xz,Jx,Jz,Wz,Zz)
    out = Array{ShocksSolution}(0)
    for r in 1:model.meta.numregimes
        push!(out,solveshocks(model,m[r],Xx,Xz,Jx,Jz,Wz,Zz,r))
    end
    # Now rearrange it into the same order as the variables are listed in the model
    shockssol = rearrangeshocks(model,out)
    return shockssol
end

function solveconstant(model,m,Xx,Xz,Jx,Jz,Wz,Zz,r)
    return ConstantSolution(-hcat(m.B1, m.B2, m.A1*Jx+m.A3*Xx+m.B3, m.A1*Jz+m.A2*Wz+m.A3*Xz+m.A4*Zz+m.B4, m.B5)\m.C_constant)
end

function rearrangeconstant(model,ordered)
    # J,W,X,Z,Y
    constantsol = zeros(Float64,model.meta.numvars,1,model.meta.numregimes)
    for r in 1:model.meta.numregimes
        offset = 1
        constantsol[model.perturbationsystem.isJ,1,r] = ordered[r].orderedMC[offset:(offset-1+countnz(model.perturbationsystem.isJ)),1]
        offset += countnz(model.perturbationsystem.isJ)
        constantsol[model.perturbationsystem.isW,1,r] = ordered[r].orderedMC[offset:(offset-1+countnz(model.perturbationsystem.isW)),1]
        offset += countnz(model.perturbationsystem.isW)
        constantsol[model.perturbationsystem.isX,1,r] = ordered[r].orderedMC[offset:(offset-1+countnz(model.perturbationsystem.isX)),1]
        offset += countnz(model.perturbationsystem.isX)
        constantsol[model.perturbationsystem.isZ,1,r] = ordered[r].orderedMC[offset:(offset-1+countnz(model.perturbationsystem.isZ)),1]
        offset += countnz(model.perturbationsystem.isZ)
        constantsol[model.perturbationsystem.isY,1,r] = ordered[r].orderedMC[offset:(offset-1+countnz(model.perturbationsystem.isY)),1]
    end
    return constantsol
end

function addssvalues!(constantsol,model)
    for r in 1:model.meta.numregimes
        constantsol[:,1,r] += model.steadystate.values
    end
    return
end

function solveconstant(model,m,Xx,Xz,Jx,Jz,Wz,Zz)
    out = Array{ConstantSolution}(0)
    for r in 1:model.meta.numregimes
        push!(out,solveconstant(model,m[r],Xx,Xz,Jx,Jz,Wz,Zz,r))
    end
    # Now rearrange it into the same order as the variables are listed in the model
    constantsol = rearrangeconstant(model,out)
    # Add steady state values
    addssvalues!(constantsol,model)
    return constantsol
end

function printif(verbose,string)
    if verbose
        print(string)
    end
    return
end

function printlnif(verbose,string)
    if verbose
        println(string)
    end
    return
end

function solve(model::RSDSGEModel,verbose::Bool = true; method::Int=1)
    printif(verbose,"Evaluating coefficient matrices...")
    m = evaluatematrices(model)
    printlnif(verbose,"done.")

    printif(verbose,"Subbing out static variables...")
    statics = substatics(model,m)
    printlnif(verbose,"done.")

    printif(verbose,"Solving exogenous equations...")
    exosolutions = solveexos(model,m,statics)
    Wz,Zz = arrangeexomatrices(model,exosolutions)
    printlnif(verbose,"done.")

    printif(verbose,"Constructing system of quadratic equations...")
    system = createsystem(model,m,statics,exosolutions,Wz,Zz)
    printlnif(verbose,"done")

    blocks,system,restricteddecisionrules = separateblocks(model,system)
    system,blockdict = solveblocks(model,blocks,system,restricteddecisionrules,verbose,method)
    merge!(blockdict,restricteddecisionrules)

    printif(verbose,"Solving complete system of quadratic equations ($(length(system)) equations)...")
    if method == 1
        dynamicsolutions = solve_withfallback(system)
    elseif method == 2
        dynamicsolutions = numericsolve(system)
    else
        error("Option $(method) is not a valid solution method.")
    end

    printlnif(verbose,"done")

    if isa(dynamicsolutions,Dict)
        dynamicsolutions = [dynamicsolutions]
    end
    println("Found $(length(dynamicsolutions)) solution(s).")
    solutions = Array{RSDSGESolution}(0)
    for (i,sol) in enumerate(dynamicsolutions)
        # Ignore any imaginary solutions
        if !hasimaginary(sol)
            # Merge in the block dict solutions
            sol = merge(sol,blockdict)
            # Here's where I probably want to test for stability using FRWZ MSS method?

            # Format into nicer form (i.e., not a dictionary)
            dynamics = formatdynamicsolution(model,sol)
            # Arrange the states solution in Xx, Xz etc
            Xx,Xz,Jx,Jz = arrangesolutionmatrices(model,exosolutions,dynamics)
            # Back out the statics
            staticsol = solvestatics(model,statics,exosolutions,dynamics,Xx,Xz,Jx,Jz,Wz,Zz)
            # Format into decision rule matrix
            decisions = convert2decision(model,dynamics,exosolutions,staticsol)
            # Solve shocks
            shocks = solveshocks(model,m,Xx,Xz,Jx,Jz,Wz,Zz)
            # Solve constant
            constant = solveconstant(model,m,Xx,Xz,Jx,Jz,Wz,Zz)
            # Wrap into solution type and push into array
            push!(solutions,RSDSGESolution(decisions,shocks,constant))
        else
            println("Solution number $(i) has imaginary components, ignoring...")
        end
    end
    return solutions
end


text:
using .Circuits
using QXZoo
using QXTns

export convert_to_tnc

"""
_convert_to_tnc(circ::QXZoo.Circuit.Circ; kwargs...)

Function to convert a QXZoo circuit to a QXTools tensor network circuit
"""
function _convert_to_tnc(circ::QXZoo.Circuit.Circ; kwargs...)
    tnc = TensorNetworkCircuit(circ.num_qubits)
    for gate in circ.circ_ops
        mat_elements = gate_matrix(gate)
        qubits = gate_qubits(gate)
        push!(tnc, qubits, mat_elements; kwargs...)
    end
    tnc
end

"""
convert_to_tnc(circ::QXZoo.Circuit.Circ;
               input::Union{String, Nothing}=nothing,
               output::Union{String, Nothing}=nothing,
               no_input::Bool=false,
               no_output::Bool=false,
               kwargs...)

Function to convert a QXZoo circuit to a QXTools tensor network circuit
"""
function convert_to_tnc(circ::QXZoo.Circuit.Circ;
                        input::Union{String, Nothing}=nothing,
                        output::Union{String, Nothing}=nothing,
                        no_input::Bool=false,
                        no_output::Bool=false,
                        kwargs...)
    tnc = _convert_to_tnc(circ; kwargs...)
    if !no_input add_input!(tnc, input) end
    if !no_output add_output!(tnc, output) end
    tnc
end