text:
export AbstractSover, SolverBigM, SolverSelfDual
export pointE_2_linearidx, linearidx_2_pointE, linearize, get_symmetricprimalpt
export solve_bigM

abstract type AbstractSover end

mutable struct SolverBigM{T, matT} <: AbstractSover
    model::SDCOContext{T, matT}
    model_bigM::SDCOContext{T, matT}

    ## Big M method parameters
    M1::T
    M2::T
end

function SolverBigM(pb::SDCOContext{T}) where T
    e = ones(pb.c)

    ## Declare internal bigM model
    c̄ = PointE(pb.c.mats, [pb.c.vec..., T(0), T(0)])

    Ā = Vector{eltype(pb.A)}(undef, length(pb.A)+1)
    for (i, ai) in enumerate(pb.A)
        Ā[i] = PointE(ai.mats, [ai.vec..., 0, -dot(ai, e)])
    end

    Ā[end] = PointE(-e.mats, [-e.vec..., -T(1), T(pb.nc)])

    b̄ = [pb.b..., 0]

    return SolverBigM(pb, SDCOContext(c̄, Ā, b̄, options=pb.options), NaN, NaN)
end

"""
    update_parameters!(Mctx, M1, M2)

update the bigM model with the new M1, M2 parameters.
"""
function update_parameters!(Mctx::SolverBigM, M1, M2)
    Mctx.M1 = M1
    Mctx.M2 = M2

    e = ones(Mctx.model.c)

    Mctx.model_bigM.b[end] = -M1
    Mctx.model_bigM.c.vec[end] = M2 - dot(Mctx.model.c, e)
end

"""
    z̄, M1, M2 = derivestartingpoint(Mctx::SolverBigM, x0::PointE)

Derive the starting point, and adequate M parameters for building an extended
big M problem, starting from a point satisfying the primal linear constraint A(x)=b.
"""
function derivestartingpoint(Mctx::SolverBigM{T, matspT}, x0::PointE{T, matT}) where {T, matT, matspT}
    pb = Mctx.model

    e = ones(x0)

    ξ2 = max(0, -min(pb, x0)) + 1
    M1 = max(0, dot(e, x0)) + 1
    ξ1 = M1 - dot(e, x0)
    add!(x0, ξ2 * e)

    y0 = zeros(T, pb.m)
    η = max(0, -min(pb, pb.c - evaluate(pb.A, y0))) + 1

    s0 = pb.c - evaluate(pb.A, y0) + η * e
    M2 = max(0, dot(s0, e)) + 1
    σ1 = η
    σ2 = M2 - dot(s0, e)

    ## complete starting points for extended problem
    x̄0 = PointE(x0.mats, [x0.vec..., ξ1, ξ2])
    ȳ0 = [y0..., η]
    s̄0 = PointE(s0.mats, [s0.vec..., σ1, σ2])
    z̄ = PointPrimalDual(x̄0, ȳ0, s̄0)

    return z̄, M1, M2
end

"""
"""
function extractinitpbsol!(z::PointPrimalDual, z̄::PointPrimalDual)
    extractinitpbsol!(z.x, z̄.x)
    z.y = z̄.y[1:end-1]
    extractinitpbsol!(z.s, z̄.s)
    nothing
end

function extractinitpbsol!(x::PointE{T, matT}, x̄::PointE{T, matT}) where {T, matT}
    for (i, mat) in enumerate(x̄.mats)
        x.mats[i] = mat
    end
    x.vec = x̄.vec[1:end-2]
    nothing
end

"""
    get_symmetricprimalpt(A)

Compute a symmetric point satisfying the primal constraint A(x) - b = 0.
"""
function get_symmetricprimalpt(A::Vector{PointE{T, matT}}, b) where {T, matT}

    ## Write down linear system matrix
    M = linearize(A)

    ## Solve
    vecsol = M \ b

    ## Build feasible point
    x0 = PointE(first(A).dims, length(first(A).vec), T, Dense{T})
    for (linind, val) in  enumerate(vecsol)
        s, i, j = linearidx_2_pointE(x0, linind)

        if s == -1 && j == -1
            x0.vec[i] = val
        else
            x0.mats[s][i, j] = val
        end
    end

    ## Symmetrize point
    for (i, mat) in enumerate(x0.mats)
        x0.mats[i] = (mat + transpose(mat)) / 2
    end

    return x0
end

function pointE_2_linearidx(ai, s, i, j)
    dims, veclen = ai.dims, length(ai.vec)
    return sum(dims[1:(s-1)].^2) + (i-1)*dims[s] + j
end

function pointE_2_linearidx(ai, i)
    dims, veclen = ai.dims, length(ai.vec)
    return sum(dims.^2) + i
end

function linearidx_2_pointE(ai, i)
    dims, veclen = ai.dims, length(ai.vec)
    if i > sum(dims.^2)
        return -1, i - sum(dims.^2), -1
    else
        s = 1
        while i > sum(dims[1:s].^2)
            s+=1
        end
        matlinind = i - sum(dims[1:s-1].^2) -1 #-1 for 0 based indexing
        i = matlinind ÷ dims[s]
        j = matlinind % dims[s]
        return s, i+1, j+1
    end
end

function linearize(A)
    n = sum(A[1].dims.^2) + length(A[1].vec)

    Is = []
    Js = []
    Vs = Float64[]
    for (ctrind, ai) in enumerate(A)
        ## Matrices
        for (s, mat) in enumerate(ai.mats)
            for i=1:size(mat, 1), j=1:size(mat, 2)
                if mat[i, j] != 0
                    push!(Is, ctrind)
                    push!(Js, pointE_2_linearidx(ai, s, i, j))
                    push!(Vs, mat[i, j])
                end
            end
        end

        ## Linear part
        for (vecind, coeff) in enumerate(ai.vec)
            push!(Is, ctrind)
            push!(Js, pointE_2_linearidx(ai, vecind))
            push!(Vs, coeff)
        end
    end

    return sparse(Is, Js, Vs, length(A), n)
end


text:
@with_kw struct Network{N_spec, N_reac, N_integrate, N_neq, T}

    species_noneq::Vector{String}
    grRec::Bool
    SupTh::Bool
    ssC  ::Bool

    iCion     :: Int32
    iH2formGr :: Int32
    iH2diss   :: Int32
    iCOdiss   :: Int32
    iHpRecGr  :: Int32
    iHepRecGr :: Int32
    iCpRecGr  :: Int32
    iSipRecGr :: Int32

    fac_H ::SVector{N_spec,Int32}
    fac_He::SVector{N_spec,Int32}
    fac_C ::SVector{N_spec,Int32}
    fac_O ::SVector{N_spec,Int32}
    fac_Si::SVector{N_spec,Int32}
    fac_S ::SVector{N_spec,Int32}
    fac_N ::SVector{N_spec,Int32}
    fac_Mg::SVector{N_spec,Int32}
    fac_Fe::SVector{N_spec,Int32}
    charge::SVector{N_spec,Int32}

    typ  :: SVector{N_reac,Int32}
    ir1  :: SVector{N_reac,Int32}
    ir2  :: SVector{N_reac,Int32}
    ip1  :: SVector{N_reac,Int32}
    ip2  :: SVector{N_reac,Int32}
    ip3  :: SVector{N_reac,Int32}
    ip4  :: SVector{N_reac,Int32}
    pro3 :: Bool
    pro4 :: Bool

    alpha :: SVector{N_reac,T}
    beta  :: SVector{N_reac,T}
    gamma :: SVector{N_reac,T}

    #always present
    iH    :: Int32
    ielec :: Int32
    iH2   :: Int32
    iHp   :: Int32

    #optional elements
    iHe   :: Int32
    iC    :: Int32
    iO    :: Int32
    iSi   :: Int32
    iS    :: Int32
    iN    :: Int32
    iMg   :: Int32
    iFe   :: Int32
    iCO   :: Int32
    iHep  :: Int32
    iCp   :: Int32
    iSip  :: Int32

    ineq          :: SVector{N_neq,Int32}
    idx_integrate :: SVector{N_integrate,Int32}
end

function read_umist_line(umist_line)
    num = umist_line[1]
    typ = umist_line[2]
    R::Vector{String} = [umist_line[3], umist_line[4]]
    P::Vector{String} = [umist_line[5], umist_line[6], umist_line[7], umist_line[8]]
    #alpha, beta, gamma, Tl, Tu, ST, ACC =
    #    umist_line[10], umist_line[11], umist_line[12],
    #    umist_line[13], umist_line[14], umist_line[15], umist_line[16]
    alpha, beta, gamma =
        umist_line[10], umist_line[11], umist_line[12]
    return num, typ, R, P, alpha, beta, gamma
end

function print_reaction(num, typ, R, P, alpha, beta, gamma)
    @printf "(%4s)" num
    #@printf "[%3s]" typ

    @printf " [α,β,γ = "
    @printf "%10.2e," alpha
    @printf "%10.2e," beta
    @printf "%10.2e]" gamma

    @printf "%8s" R[1]
    @printf "     + %8s" R[2]
    if typ == "GR"
        @printf "(gr)-->%8s" P[1]
    else
        @printf "    -->%8s" P[1]
    end
    if(P[2]!="")
        @printf "     + %8s" P[2]
    end
    if(P[3]!="")
        @printf "     + %8s" P[3]
    end
    if(P[4]!="")
        @printf "     + %8s" P[4]
    end
    print("\n")
end

function print_reaction(num, typ, R, P)
    @printf "%7s" R[1]
    @printf "     + %7s" R[2]
    if typ == "GR"
        @printf " (gr)-->%7s" P[1]
    else
        @printf "     -->%7s" P[1]
    end
    if(P[2]!="")
        @printf "     + %7s" P[2]
    end
    if(P[3]!="")
        @printf "     + %7s" P[3]
    end
    if(P[4]!="")
        @printf "     + %7s" P[4]
    end
    print("\n")
end

function search_species_reactants(network, s1::String)
    reac=[]
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        if any(s1 .== R)
            print_reaction(n, t, R, P, a, b, c)
            push!(reac, n)
        end
    end
    reac
end

function search_species_reactants(network, s1::String, s2::String)
    reac=[]
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        if any(s1 .== R) && any(s2 .== R)
            @printf "(i=%.4d) " i
            print_reaction(n, t, R, P, a, b, c)
            push!(reac, n)
        end
    end
    reac
end

function search_species_products(network, s1::String)
    reac=[]
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        if any(s1 .== P)
            print_reaction(n, t, R, P, a, b, c)
            push!(reac, n)
        end
    end
    reac
end

function search_species_products(network, s1::String, s2::String)
    reac=[]
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        if any(s1 .== P) && any(s2 .== P)
            print_reaction(n, t, R, P, a, b, c)
            push!(reac, n)
        end
    end
    reac
end

function print_all_reactions(network)
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        #@printf "(%3d)  " i
        print_reaction(n, t, R, P, a, b, c)
    end
end

function get_all_species_in_network(network)
    species::Vector{String} = []
    for i in 1:size(network, 1)
        n, t, R, P, a, b, c = read_umist_line(network[i,:])
        #@printf "(%3d)  " i
        #print_reaction(n, t, R, P, a, b, c)
        [push!(species, R[j]) for j in eachindex(R)]
        [push!(species, P[k]) for k in eachindex(P)]
    end
    species = unique(species)
    filter!(s -> (s != "")&&(s != "PHOTON")&&(s != "CRPHOT")&&(s != "CRP"), species)
end

function breakdown_species(s::String)
    charge::Int32 = 0
    elements = String[]
    counts   = Int32[]
    mult::Int32 = 1
    skipflag::Bool = false

    if s == "e-"
        return -1, elements, counts
    end

    for i in length(s):-1:1
        if skipflag
            #@show skipflag
            skipflag = false
            continue
        end
        char = s[i]
        #@show char
        isuppercase(char)
        if char == '+'
            charge += 1
        elseif char == '-'
            charge -= 1
        elseif isletter(char)
            elem = ""
            if islowercase(char)
                elem = s[i-1:i]
                skipflag = true
            else
                elem = string(char)
            end
            push!(elements, elem)
            push!(counts, mult)
            mult = 1
        else
            mult = parse(Int32, char)
        end
    end
    return charge, elements, counts
end

#s="HSi3Mg2+"
#breakdown_species(s)

function get_unique_elements(elements::Vector{Vector{String}})
    unique_elements = String[]
    for i in 1:length(elements)
        for j in eachindex(elements[i])
            push!(unique_elements, elements[i][j])
        end
    end
    return unique!(unique_elements)
end

function get_idx_integrate(elements::Vector{Vector{String}}, species_noneq::Vector{String}, dict::Dict)
    N_spec = length(elements)
    unique_elements = get_unique_elements(elements)
    #const N_uni_elem = length(unique_elements_a)
    #const uni_elem = SVector{N_uni_elem}(unique_elements_a)
    #N_conserve = length(unique_elements) + 1  #+1 for electron

    #species_derived = ["e-", "H", "He", "C", "O", "Si", "S"]
    species_derived = unique_elements
    push!(species_derived, "e-")
    #@show species_derived

    N_neq = length(species_noneq)
    if N_neq > 0
        for i in 1:N_neq
            push!(species_derived, species_noneq[i])
        end
    end
    ineq = SVector{N_neq}([dict[species_noneq[i]] for i in 1:N_neq])
    #@show species_derived
    #@show ineq
    #@show all_species
    idx_derived = sort([dict[species_derived[i]] for i in eachindex(species_derived)])

    idx_integrate_a = filter!(s -> all(s .!= idx_derived), collect(1:N_spec))
    return idx_integrate_a, species_derived
end
function construct_index_arrays(reduced_network::Matrix{Any}, dict::Dict)

    N_reac = size(reduced_network, 1)
    num_a = zeros(Int32, N_reac)
    typ_a = zeros(Int32, N_reac)

    ir1_a = zeros(Int32, N_reac)
    ir2_a = zeros(Int32, N_reac)
    ip1_a = zeros(Int32, N_reac)
    ip2_a = zeros(Int32, N_reac)
    ip3_a = zeros(Int32, N_reac)
    ip4_a = zeros(Int32, N_reac)

    alpha_a = zeros(N_reac)
    beta_a  = zeros(N_reac)
    gamma_a = zeros(N_reac)

    iCion     = 0
    iH2formGr = 0
    iH2diss   = 0
    iCOdiss   = 0
    iHpRecGr  = 0
    iHepRecGr = 0
    iCpRecGr  = 0
    iSipRecGr = 0

    for i in 1:N_reac
        num_a[i], t, R, P, alpha_a[i], beta_a[i], gamma_a[i] = read_umist_line(reduced_network[i,:])
        if any(R .== "C+") && any(R .== "e-")
            #alpha_a[i] = 4.67e-12
            #gamma_a[i] = 3.76
            #@show alpha_a[i], beta_a[i], gamma_a[i]
        end
        if any(R .== "H-") && any(R .== "PHOTON")
            #alpha_a[i] = 55.8e-10
            #gamma_a[i] = 3.76
            #@show alpha_a[i], beta_a[i], gamma_a[i]
        end
        if any(R .== "PHOTON") && any(R .== "C") #C photoionization
            iCion = i
            #@show alpha_a[i], gamma_a[i]
        end
        if t=="GR" && any(R .== "H") && any(P .== "H2") #H2 photodissociation
            iH2formGr = i
        end
        if any(R .== "H2") && any(R .== "PHOTON") && any(P .== "H") #H2 photodissociation
            iH2diss = i
        end
        if any(R .== "CO") && any(R .== "PHOTON") && any(P .== "C") && any(P .== "O") #CO photodissociation
            iCOdiss = i
        end
        if t=="GR" && any(R .== "H+") && any(R .== "e-") && any(P .== "H") && any(P .== "PHOTON") #H+ recomb. on grains
            iHpRecGr = i
        end
        if t=="GR" && any(R .== "He+") && any(R .== "e-") && any(P .== "He") && any(P .== "PHOTON") #He+ recomb. on grains
            iHepRecGr = i
        end
        if t=="GR" && any(R .== "C+") && any(R .== "e-") && any(P .== "C") && any(P .== "PHOTON") #C+ recomb. on grains
            iCpRecGr = i
        end
        if t=="GR" && any(R .== "Si+") && any(R .== "e-") && any(P .== "Si") && any(P .== "PHOTON") #H+ recomb. on grains
            iSipRecGr = i
        end

        if any(R .== "PHOTON")
            typ_a[i] = 1
        elseif any(R .== "CRP")
            typ_a[i] = 2
        elseif any(R .== "CRPHOT")
            typ_a[i] = 3
        else
            typ_a[i] = 0
        end

        #overwrite if it's a grain-process
        if t == "GR"
            typ_a[i] = 4
        end

        #overwrite if it's an ion-neutral process
        #if SupTh == true && t == "IN"
        #    if any(R .== "C+") && any(R .== "H2") && any(P .== "CH+") && any(P .== "H")
        #        typ_a[i] = 5
        #        @show i,R,P
        #    end
        #end

        ir1_a[i], ir2_a[i] = dict[R[1]], dict[R[2]]
        ip1_a[i], ip2_a[i], ip3_a[i], ip4_a[i] = dict[P[1]], dict[P[2]], dict[P[3]], dict[P[4]]
    end
    return typ_a, ir1_a, ir2_a, ip1_a, ip2_a, ip3_a, ip4_a, alpha_a, beta_a, gamma_a,
    iCion, iH2formGr, iH2diss, iCOdiss, iHpRecGr, iHepRecGr, iCpRecGr, iSipRecGr
end
function setup_conservation_arrays(elem::String, elements::Vector{Vector{String}}, counts::Vector{Vector{Int32}})
    N_spec = length(elements)
    idx_spec = Int32[]
    count_spec = Int32[]
    fac = zeros(Int32, N_spec)

    for i in eachindex(elements)
        for j in eachindex(elements[i])
            #@show i, j, elements[i][j]
            if elements[i][j] == elem
                push!(idx_spec, i)
                push!(count_spec, counts[i][j])
                fac[i] = counts[i][j]
                #@show counts[i][j], all_species[i]
            end
        end
    end
    fac
end

function initialize_chemistry_network(all_species::Vector{String};
    species_noneq::Vector{String} = String[],
    grRec::Bool = true,
    SupTh::Bool = false,
    ssC  ::Bool = false,
    T::DataType = Float64
    )
    N_spec = length(all_species)

    all_species_2 = vcat(all_species , ["", "PHOTON", "CRP", "CRPHOT"]);

    #read UMIST table
    umist = readdlm((@__DIR__)*"/RATE12.dist.txt",':')

    reac_num=[]
    for i in 1:size(umist,1)
        #println(i)
        n, t, R, P, a, b, c = read_umist_line(umist[i,:])
        if any(R[1].==all_species_2) && any(R[2].==all_species_2) && any(P[1].==all_species_2) && any(P[2].==all_species_2) && any(P[3].==all_species_2) && any(P[4].==all_species_2)
            #print_reaction(n, t, R, P, a, b, c)
            push!(reac_num,n)
        end
    end
    reduced_network = umist[reac_num,:]

    #add 7 more reactions (H2 & CO dissociation & grain-processes)
    more_reactions = readdlm((@__DIR__)*"/user-rates.txt",':')
    reac_num=[]
    for i in 1:size(more_reactions,1)
        #println(i)
        n, t, R, P, a, b, c = read_umist_line(more_reactions[i,:])
        if any(R[1].==all_species_2) && any(R[2].==all_species_2) && any(P[1].==all_species_2) && any(P[2].==all_species_2) && any(P[3].==all_species_2) && any(P[4].==all_species_2)
            #print_reaction(n, t, R, P, a, b, c)
            push!(reac_num,n)
        end
    end
    #@show reac_num
    more_reactions = more_reactions[reac_num,:]
    reduced_network = vcat(reduced_network, more_reactions);
    print_all_reactions(reduced_network)

    N_reac = size(reduced_network, 1)
    #@show N_reac, N_spec

    dict = Dict(all_species .=> collect(1:N_spec) );
    dict[""] = dict["PHOTON"] = dict["CRP"] = dict["CRPHOT"] = 0

    charge = @. Int32(getindex(breakdown_species(all_species) , 1))
    elements = @. getindex(breakdown_species(all_species) , 2)
    counts = @. getindex(breakdown_species(all_species) , 3)

    idx_integrate, species_derived = get_idx_integrate(elements, species_noneq, dict)
    N_integrate = length(idx_integrate)

    N_neq = length(species_noneq)
    ineq = SVector{N_neq,Int32}([dict[species_noneq[i]] for i in 1:N_neq])

    typ, ir1, ir2, ip1, ip2, ip3, ip4, alpha, beta, gamma,
    iCion, iH2formGr, iH2diss, iCOdiss, iHpRecGr, iHepRecGr, iCpRecGr, iSipRecGr = construct_index_arrays(reduced_network, dict)

    net = Network{N_spec, N_reac, N_integrate, N_neq, Float64}(
    grRec = grRec,
    SupTh = SupTh,
    ssC   = ssC,
    iCion     = iCion    ,
    iH2formGr = iH2formGr,
    iH2diss   = iH2diss  ,
    iCOdiss   = iCOdiss  ,
    iHpRecGr  = iHpRecGr ,
    iHepRecGr = iHepRecGr,
    iCpRecGr  = iCpRecGr ,
    iSipRecGr = iSipRecGr,
    species_noneq = species_noneq,
    fac_H  = setup_conservation_arrays("H" , elements, counts),
    fac_He = setup_conservation_arrays("He", elements, counts),
    fac_C  = setup_conservation_arrays("C" , elements, counts),
    fac_O  = setup_conservation_arrays("O" , elements, counts),
    fac_Si = setup_conservation_arrays("Si", elements, counts),
    fac_S  = setup_conservation_arrays("S" , elements, counts),
    fac_N  = setup_conservation_arrays("N" , elements, counts),
    fac_Mg = setup_conservation_arrays("Mg", elements, counts),
    fac_Fe = setup_conservation_arrays("Fe", elements, counts),
    charge = charge,
    typ    = typ  ,
    ir1    = ir1  ,
    ir2    = ir2  ,
    ip1    = ip1  ,
    ip2    = ip2  ,
    ip3    = ip3  ,
    ip4    = ip4  ,
    pro3   = any(ip3 .!= 0),
    pro4   = any(ip4 .!= 0),
    alpha  = alpha,
    beta   = beta ,
    gamma  = gamma,
    #always present
    iH = dict["H"],
    ielec = dict["e-"],
    iH2 = dict["H2"],
    iHp = dict["H+"],
    iHe = any(all_species.=="He") ? dict["He"] : 0,
    iC  = any(all_species.=="C" ) ? dict["C" ] : 0,
    iO  = any(all_species.=="O" ) ? dict["O" ] : 0,
    iSi = any(all_species.=="Si") ? dict["Si"] : 0,
    iS  = any(all_species.=="S" ) ? dict["S" ] : 0,
    iN  = any(all_species.=="N" ) ? dict["N" ] : 0,
    iMg = any(all_species.=="Mg" ) ? dict["Mg" ] : 0,
    iFe = any(all_species.=="Fe" ) ? dict["Fe" ] : 0,
    iCO  = any(all_species.=="CO" ) ? dict["CO" ] : 0,
    iHep = any(all_species.=="He+") ? dict["He+"] : 0,
    iCp = any(all_species.=="C+") ? dict["C+"] : 0,
    iSip = any(all_species.=="Si+") ? dict["Si+"] : 0,
    #N_neq = N_neq,
    #N_integrate = N_integrate,
    idx_integrate = idx_integrate,
    ineq = ineq
    )
    return net, dict
end

@with_kw struct AbundTotal{T}
    abC :: T = 1.4e-4
    abO :: T = 3.2e-4
    # abSi = 0.0
    # abSi = 1.5e-5
    abSi :: T = 1.7e-6 #depleted abundance
    # abS = 1e-6
    # abN = 7.6e-5
    # abFe = 1.6e-6
    # abMg = 1.4e-5
    abS :: T = 0.0
    abN :: T = 0.0
    abFe :: T = 0.0
    abMg :: T = 0.0
end

import Base.+
+(a::AbundTotal{T}, b::AbundTotal{T}) where {T<:Real} =
AbundTotal{T}(
a.abC  + b.abC ,
a.abO  + b.abO ,
a.abSi + b.abSi,
a.abS  + b.abS ,
a.abN  + b.abN ,
a.abFe + b.abFe,
a.abMg + b.abMg)

import Base.-
-(a::AbundTotal{T}, b::AbundTotal{T}) where {T<:Real} =
AbundTotal{T}(
a.abC  - b.abC ,
a.abO  - b.abO ,
a.abSi - b.abSi,
a.abS  - b.abS ,
a.abN  - b.abN ,
a.abFe - b.abFe,
a.abMg - b.abMg)

import Base.*
*(a::AbundTotal, c::T) where {T<:Real} = AbundTotal(c*a.abC, c*a.abO, c*a.abSi, c*a.abS, c*a.abN, c*a.abFe, c*a.abMg)
*(c::T, a::AbundTotal) where {T<:Real} = *(a::AbundTotal, c::T)
import Base./
/(a::AbundTotal, c::T) where {T<:Real} = AbundTotal(a.abC/c, a.abO/c, a.abSi/c, a.abS/c, a.abN/c, a.abFe/c, a.abMg/c)

#=
function test()
    all_species =
    String[
    "H", "H-", "H2", "H+", "H2+", "H3+", "e-", "He", "He+", "HeH+",
    "C", "C+", "CO", "HCO+",
    "O", "O+", "OH", "OH+", "H2O+", "H3O+", "H2O",
    "O2", "CO+", "O2+",
    "CH2", "CH2+", "CH", "CH+", "CH3+",
    "Si+", "Si"
    #"SiO", "SiO+", "SiOH+",
    #"S", "S+",
    ]

    species_noneq=String[]
    #net::Network{31,286,25,0,Float64}, dict = initialize_chemistry_network(all_species, species_noneq)
    net, dict = initialize_chemistry_network(all_species, species_noneq)

    #@code_warntype initialize_chemistry_network(all_species, species_noneq)
    @unpack iH2, iCO, iC, fac_H, fac_C, fac_O, charge = net
    abtot = AbundTotal()
    @unpack abC_s, abO_s, abSi_s = abtot

    net
end
@code_warntype test()
#n = test()
=#

###### specific reactions
#=
reactions = Int32[]

rea_H = [6155, 75, 4921, 6093, 489, 731, 730, 2604, 1278, 1277, 2895, 733, 406, 408, 434, 6162, 492]
rea_He = [734, 3431, 491, 458, 2677, 2610, 2678, 3065, 6094,6166, 3517, 6173, 967]
rea_O = [2686, 2951, 2689, 2662, 1296, 1295, 1293, 1294, 1266, 1264, 1265, 1427, 5643, 5199,
1646, 4029, 408, 1613, 2789, 2767, 2904, 378, 5196, 1643, 1642, 5705, 948, 405, 1421, 5412]
rea_CH = [2852, 2647, 6095, 2648, 1175, 1176, 1177, 1158, 1161, 1160, 1159, 3055, 5370, 5374]
rea_CO = [5199, 1646, 488, 2654, 1349, 5310, 5231, 5232, 64, 1349, 823, 6158, 707, 194, 876]
#rea_photo = [5997, 5915, 5989, 5887, 5888, 5902]
rea_suprathermal = [2618, 2189]

reactions = vcat(rea_H, rea_He, rea_O, rea_CH, rea_CO, rea_suprathermal)
unique!(reactions);

reduced_network = umist[reactions,:]
const all_species = get_all_species_in_network(reduced_network);

#add photoprocesses automatically
charge_a = @. Int32(getindex(breakdown_species(all_species) , 1));
photo_reactions = []
for i in 1:size(all_species,1)
    #if charge_a[i] <= 0
    if true
        r = search_species_reactants(umist, "PHOTON", all_species[i])
        if all_species[i] != "CO" #will do CO shielding manually
            global photo_reactions = vcat(r, photo_reactions)
        end
    end
end

#add photoprocesses to the network
all_reactions = vcat(reactions, photo_reactions)
unique!(all_reactions);
reduced_network = umist[all_reactions,:];

const N_reac = size(reduced_network, 1)
const all_species = get_all_species_in_network(reduced_network);
const N_spec = size(all_species,1)
#println(all_species)
#print_all_reactions(reduced_network)
#@show N_reac, N_spec
=#

######

ext: jl


text:
include("src/main.jl")

BUSES, LINES, GENERATORS, rPTDF =  load_case_data(datafile="basecase")
# basecas = 15 Bus example from A. Papavasiliou, 2016, "Analysis of Distribution Locational Marginal Prices"

ld_objective, ld_bus_results, ld_line_results = Solve_LinDist(BUSES, LINES, GENERATORS)
info("Results of standard lindist. Objective: $ld_objective")
println(ld_bus_results)
println(ld_line_results)


text:
# Autogenerated wrapper script for Elfutils_jll for x86_64-linux-musl
export libasm, libdw, libelf

using Zlib_jll
using Bzip2_jll
using XZ_jll
using argp_standalone_jll
using fts_jll
using obstack_jll
JLLWrappers.@generate_wrapper_header("Elfutils")
JLLWrappers.@declare_library_product(libasm, "libasm.so.1")
JLLWrappers.@declare_library_product(libdw, "libdw.so.1")
JLLWrappers.@declare_library_product(libelf, "libelf.so.1")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, Bzip2_jll, XZ_jll, argp_standalone_jll, fts_jll, obstack_jll)
    JLLWrappers.@init_library_product(
        libasm,
        "lib/libasm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdw,
        "lib/libdw.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libelf,
        "lib/libelf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type AbstractIterativeSampler end

export gradient

function _sample_state!(is::AbstractIterativeSampler)
    samples        = batches(is.samples)

    # Monte-Carlo sampling
    σ_old = first(samples)

    init_sampler!(is.sampler, is.bnet, σ_old, is.sampler_cache)
    for σ_next=samples
        !samplenext!(σ_next, σ_old,
                     is.sampler, is.bnet, is.sampler_cache) && break
        σ_old = σ_next
    end
end

function _center_gradient!(is::AbstractIterativeSampler)
    # Center the gradient so that it has zero-average
    # do it for every block of gradients / cogradients
    for (∇vec_avg, ∇vals_vec)=zip(is.∇vec_avg, is.∇vals_vec)
        workers_mean!(∇vec_avg, ∇vals_vec, is.parallel_cache) # MPI
        ∇vals_vec .-= ∇vec_avg
    end
end

function gradient(is::AbstractIterativeSampler, iter; kwargs...)
    loss, grad_before_precond = sample!(is; kwargs...)
    grad = precondition!(grad_before_precond, is.precond_alg, iter)
    return loss, grad
end


text:
"""
Store a "hit record" for a ray intersecting a geometric object.
"""
struct HitRecord
    # intersection point
    point::Vector
    # surface normal at intersection point
    normal::Vector
    # reference to material
    material

    function HitRecord(point, normal, material)
        @assert(abs(norm(normal) - 1) < 1e-11, "hit record normal must be normalized")
        new(point, normal, material)
    end
end


text:
import Statistics, GLM, DataFrames

"""
    datagrid(df::DataFrames.DataFrame; cols=:all, n::Int=10, kwargs...)

Create a reference grid of data.

# Arguments
- `df`: A `DataFrame` (can be a `StatsModel`).
- `cols`: The target columns. The rest will be maintained "fixed".
- `n`: For numeric targets columns, what desired length (controls the spacing).
- `fix_num`: How to fix the numeric variables. Can be a function or a number.
- `fix_fac`: How to fix the factors. Should be a String indicating an existing level.

# Examples
```julia
df = simulate_data_correlation([[0.2, 0.5], [0.4, 0.2]])

grid = datagrid(df, n=3)
```
"""
function datagrid(df::DataFrames.DataFrame; cols=:all, n::Int=10, kwargs...)
    if cols == :all
        cols = names(df)
    elseif isa(cols, Vector{String})
        cols = Symbol.(cols)
    elseif isa(cols, Vector{Int})
        cols = names(df)[cols]
    elseif !isa(cols, Vector{Symbol})
        throw(ArgumentError("`cols` must be Vector{String}, Vector{Int}, Vector{Symbol} or :all."))
    end

    # Split data
    df_target = df[cols]
    df_to_keep_fixed = df[setdiff(names(df), cols)]

    # Create grid
    grid = colwise(x -> datagrid(x, n=n), df_target)
    grid = vec(collect(Base.product(grid...)))
    grid = DataFrame(collect.([ getindex.(grid,t) for t in 1:length(grid[1]) ]), names(df_target))

    # Add fixed
    df_to_keep_fixed = keep_fixed(df_to_keep_fixed, n=nrow(grid); kwargs...)
    grid = hcat(grid, df_to_keep_fixed)

    return grid
end

function datagrid(model::StatsModels.DataFrameRegressionModel; cols=:all, n::Int=10, kwargs...)
    data = model_data(model)
    # Select only Predictors
    data = data[values(model_description(model))["Predictors"]]
    grid = datagrid(data, cols=cols, n=n; kwargs...)
    return grid
end

# datagrid: Vectors ------------------------------------------------------------------------

function datagrid(X::AbstractVector{<:Union{Real, Missing}}; n::Int=10)
    X = skipmissing(X)
    X = collect(range(minimum(X), stop=maximum(X), length=n))
end

function datagrid(X::DataFrames.CategoricalVector; n=nothing)
    X = levels(X)
end

function datagrid(X::AbstractVector{<:Union{String, Missing}}; n=nothing)
    X = unique(skipmissing(X))
end

# keep_fixed -------------------------------------------------------------------------------

function keep_fixed(X::AbstractVector{<:Union{Real, Missing}}; fix_num=Statistics.mean, kwargs...)
    if isa(fix_num, Number)
        return fix_num
    else
        X = skipmissing(X)
        fix_num(X)
    end
end

function keep_fixed(X::DataFrames.CategoricalVector; fix_fac=1, kwargs...)
    if isa(fix_fac, Int)
        return levels(X)[fix_fac]
    elseif isa(fix_fac, String)
        if fix_fac in levels(X)
            return fix_fac
        else
            throw(ArgumentError("`fix_fac` must be an existing level."))
        end
    else
        throw(ArgumentError("`fix_fac` must be an Int or a String."))
    end
end

function keep_fixed(X::AbstractVector{<:Union{String, Missing}}; fix_fac=1, kwargs...)
    if isa(fix_fac, Int)
        X = skipmissing(X)
        return unique(X)[fix_fac]
    elseif isa(fix_fac, String)
        return fix_fac
    else
        throw(ArgumentError("`fix_fac` must be an Int or a String."))
    end
end

function keep_fixed(df::DataFrames.DataFrame; fix_num=Statistics.mean, fix_fac=1, n::Int=1, kwargs...)
    fixed = colwise(x -> keep_fixed(x, fix_num=fix_num, fix_fac=fix_fac; kwargs...), df)
    fixed = DataFrame([repeat([x], n) for x in fixed], names(df))
    return fixed
end


text:
# download input files from: https://gist.github.com/isaacsas/7648fe065f1884307a2906a0a48ea12b
# then add ReactionNetworkImporter package (unregistered)
#]
# add https://github.com/isaacsas/ReactionNetworkImporter.jl.git

using DiffEqBase, DiffEqBiological,Plots,OrdinaryDiffEq, Sundials
using ReactionNetworkImporter
using TimerOutputs

# parameters
networkname = "tester"
tf = 10.

# input files, specify path
#speciesf = "SOME-PATH-TO/BCR_pop.txt"
#rxsf = "SOME-PATH-TO/BCR_rxn.txt"

const to = TimerOutput()
reset_timer!(to)

# get the reaction network
@timeit to "netgen" rn,initialpop = get_rxnetwork_simple(RSSAFile(), networkname, speciesf, rxsf; printrxs = false)
@timeit to "addodes" addodes!(rn; build_jac=false, build_symfuncs=false)
@timeit to "ODEProb" oprob = ODEProblem(rn,convert.(Float64,initialpop),(0.,tf))
show(to)

# note solvers run _much_ faster the second time
reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false,calck=false); end; show(to)
reset_timer!(to); @timeit to "Rodas4" begin sol = solve(oprob,Rodas4(autodiff=false),dense=false); end;

reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=false); end; show(to)
reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(linear_solver=:GMRES),dense=false); end; show(to)

# note, restart REPL to see timing with first pass overhead for the next solver
reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(),dense=false); end; show(to)
reset_timer!(to); @timeit to "CVODE_BDF" begin sol = solve(oprob,CVODE_BDF(),dense=false); end; show(to)

max_issues_repo_name: isaacsas/ParseRxNetwork


text:
"""
    function workflow()
---
This records the framework of the package.
It can also serve as a template pipeline.

## Parameter example
- The keys must be of these names

```julia
    Parameters = Dict(:nSire=> 100,
                      :nDam => 200, # male:female = 1:2
                      :nSib => 2nsib,
                      :nC7e => nsib, # number of sib to be challenged
                      :nG8n => 10, # number of generations
                      :nQTL => [1000, 1000],
                      :h²   => [.5, .5],
                      :p8e  => .5, # percentage of affected in the binary trait
                      :e19e => .8, # edit_successsful_rate
                      :w4t  => weight # weight on the binary trait EBV
                      )
    par = (; Parameters...)     # named tuple.  contents as above
```
"""
function workflow(; debug = true)
    #copyright()
    if_normalization()
end

"""
    function data_collection()
This function collects some real data for downstream simulation.  It
1. remove duplicates
2. phase with `beagle`.
3. serialize the data to disk
4. also save a `jld2` file
"""
function data_collection()
    @time gt600()           # -> c.vcf.gz
    @time base = vcf2dic(joinpath(dat_dir, "run/c.vcf.gz"))
    @time serialize(joinpath(dat_dir, "run/ns.ser"), base)
    @save joinpath(dat_dir, "run/base.jld2") base
end

"""
    function a_simple_snp_blup_test()
This function using a real dataset from Nofima, which has 600 ID and 50,793 loci.
It consists 4 scenarios:
1. SNP-BLUP with 012 genotypes.
2. SNP-BLUP with normalized genotypes.
"""
function if_normalization()
    @load "dat/run/base.jld2" base
    nb = begin                  # a smaller new base
        loci = sort(randperm(size(base.hap)[1])[1:400])
        (hap = base.hap[loci, :], pos = base.pos[loci])
    end
    p = mean(nb.hap, dims=2)
    qtl = sim_QTL(nb, 350)[1]       # one QTL
    gt = hap2gt(nb.hap)             # → nlc × nid
    Q = qtl_gt(nb.hap, qtl.pos)     # QTL genotypes
    tbv = Q'qtl.effect
    z = zeros(size(gt))
    copyto!(z, gt)

    # snp-blup with no normalization
    lhs, rhs = begin
        x = [ones(600) z]
        x'x, x'tbv
    end
    esnp = lhs \ rhs
    ebv = gt * esnp[2:end]
    cor(ebv, tbv), norm(ebv - tbv)
    # returns usually close to (1, 0), meaning not necessary to normalize genotypes.
end

function tst_gwas()
    # Simulation with the base population
    nqtl = 30
    @load "dat/run/base.jld2" base
    base
    h² = .75
    qtl = sim_QTL(base, 30)[1]
    _, pt = phenotype(base.hap, qtl, h²)
    gt = hap2gt(base.hap)
    gt, pt, qtl
end


text:
using Combinatorics
using Serialization

function init_mp(PG::Pair_ER)
    zeros(2*ne(PG.N1.G),2*ne(PG.N2.G))
end

function sum_combi(neig_i1::Array{Int64,1},neig_i2::Array{Int64,1},k::Int64,
    mp_previous::Array{Float64,2},gen_i1::Vector{Vector{Int64}},gen_i2::Vector{Vector{Int64}},perm_k::Vector{Vector{Int64}})
    if  k == 0
        return Float64(0)
    end
    sum_inj = -Inf#Float64(0)
    for s1 in gen_i1
        for s2 in gen_i2
            for ps2 in perm_k
                sum_inj = log_sum_exp(sum_inj, sum(mp_previous[neig_i1[s1[j]],neig_i2[s2[ps2[j]]]] for j in 1:k))
                #sum_inj += prod(exp(mp_previous[neig_i1[s1[j]],neig_i2[s2[ps2[j]]]]) for j in 1:k)
            end
        end
    end
    sum_inj
    #log(sum_inj)
end

function m_passing(PG::Pair_ER,e1::Int64,e2::Int64,mp_previous::Array{Float64,2})
    # trees T_1(e1) and T_2(e2)
    dic_comb = PG.dic_comb
    dic_perm = PG.dic_perm
    d1 = PG.N1.deg[e1]
    d2 = PG.N2.deg[e2]
    deg_min = min(d1,d2)
    if deg_min == 1
        output = psi_log(PG,0,d1-1,d2-1)
    else
        neig_i1 = PG.N1.neig[e1]
        neig_i2 = PG.N2.neig[e2]
        #gen_i1,gen_i2,perm_k = dic_comb[(d1-1,0)], dic_comb[(d2-1,0)], dic_perm[0]
        output = -Inf#psi_log(PG,0,d1-1,d2-1) + sum_combi(neig_i1,neig_i2,0,mp_previous, gen_i1, gen_i2, perm_k)#Float64(0)
        for k in 0:(deg_min-1)
            gen_i1,gen_i2,perm_k = dic_comb[(d1-1,k)], dic_comb[(d2-1,k)], dic_perm[k]
            output = log_sum_exp(output, psi_log(PG,k,d1-1,d2-1) + sum_combi(neig_i1,neig_i2,k,mp_previous, gen_i1, gen_i2, perm_k))
        end
    end
    output
end

@inline function compute_next!(PG::Pair_ER,mp_previous::Array{Float64,2},mp_next::Array{Float64,2},ne1::Int64,ne2::Int64)
    @inbounds for e1 in 1:ne1
        @inbounds for e2 in 1:ne2
            mp_next[e1,e2] = m_passing(PG,e1,e2,mp_previous)
        end
    end
    nothing
end

function run_bp(PG::Pair_ER,n_iter::Int64;mp_previous=nothing,verbose=false)
    if mp_previous === nothing
        mp_previous = init_mp(PG)
    end
    mp_next = similar(mp_previous)
    mp_opt = similar(mp_previous)
    (ne1,ne2) = size(mp_previous)
    res = zeros(6,n_iter)
    perf = 0
    last_iter = 0
    ov_opt = 0
    for i in 1:n_iter
        compute_next!(PG,mp_previous,mp_next,ne1,ne2)
        (mp_previous, mp_next) = (mp_next, mp_previous)
        M_loglr = create_matrix_lr(PG,mp_previous)
        #M_loglr = log.(M_loglr)
        ov1, ov2, v1,v2 = eval_M(PG, M_loglr)
        m1,m2 = eval_edges(PG, M_loglr)
        new_perf = (m1+m2)/2
        if new_perf > perf
            perf = new_perf
            mp_opt = deepcopy(mp_previous)
            last_iter = i
            ov_opt = (ov1+ov2)/2
        end
        res[1,i] = ov1
        res[2,i] = ov2
        res[3,i] = v1
        res[4,i] = v2
        res[5,i] = m1
        res[6,i] = m2
        if verbose
            println(i, " | ", m1, " | ", m2, " | " ,ov1," | ",ov2, " | ", v1, " | ", v2)
        end
        if isnan(v1) || isnan(v2)
            break
        end
    end
    mp_opt#res, mp_opt, last_iter, perf, ov_opt/PG.n
end

function lr_root(PG::Pair_ER,mp,d1,d2,neig_i1,neig_i2)
    dic_comb = PG.dic_comb
    dic_perm = PG.dic_perm
    deg_min = min(d1,d2)
    gen_i1,gen_i2,perm_k = dic_comb[(d1,0)], dic_comb[(d2,0)], dic_perm[0]
    sum = psi_log(PG,0,d1,d2) + sum_combi(neig_i1,neig_i2,0,mp, gen_i1, gen_i2, perm_k)
    for k in 1:deg_min
        gen_i1,gen_i2,perm_k = dic_comb[(d1,k)], dic_comb[(d2,k)], dic_perm[k]
        sum = log_sum_exp(sum, psi_log(PG,k,d1,d2) + sum_combi(neig_i1,neig_i2,k,mp, gen_i1, gen_i2, perm_k))
    end
    sum
end

function create_matrix_lr(PG::Pair_ER,mp)
    v1 = vertices(PG.N1.G)
    v2 = vertices(PG.N2.G)
    matrix_lr = zeros(length(v1),length(v2))
    deg1, neig1 = neighbors_edges(PG.N1.G)
    deg2, neig2 = neighbors_edges(PG.N2.G)
    for i in v1
        for j in v2
            matrix_lr[i,j] = lr_root(PG,mp,deg1[i],deg2[j],neig1[i],neig2[j])
        end
    end
    matrix_lr
end

function init_mp_biased(PG::Pair_ER,bias=10)
    mp_init = ones(2*ne(PG.N1.G),2*ne(PG.N2.G))
    for (i,j) in eachrow(PG.edge_G1G2)
        mp_init[i,j] = bias
    end
    mp_init
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    neighbourhoodOffsets(radius::Float64)

Return every (unique) `CartesianIndex{2}` in a `radius`-sized circle around `index`.
"""
function neighbourhood(index::CartesianIndex{2}, radius::Float64,
                       settings::Settings = defaultSettings)
    offsets = neighbourhoodOffsets(radius)
    neighbours = (i -> i + index).(offsets)
    return unique(correctCoords(neighbours, settings))
end

"""
    neighbourhoodOffsets(radius::Float64)

Return every `CartesianIndex{2}` in a `radius`-sized circle around (0,0).
"""
function neighbourhoodOffsets(radius::Float64)
    r_int = ceil(Int, radius)
    offsets = (CartesianIndices((r_int,r_int) .* 2 .+ 1) .-
               CartesianIndex(r_int+1, r_int+1))[:]
    filter!(i -> sum(i.I.^2) <= radius^2, offsets)
    return offsets
end

"""
    directNeighbours(ind::CartesianIndex{2})

Return the eight direct neighbours of index `ind`.
"""
function directNeighbours(ind::CartesianIndex{2}, settings::Settings = defaultSettings)
    neighbours = (CartesianIndices((3,3)) .- (CartesianIndex(2,2) - ind))[:]
    # remove the index itself
    deleteat!(neighbours, 5)
    return correctCoords(neighbours, settings)
end

"""
    correctCoords(coords::AbstractVector{CartesianIndex{2}})

Remove (if the lattice is a plane) or wrap (if the lattice is a toroid)
coordinates outside of the lattice size.
"""
function correctCoords(coords::AbstractVector{CartesianIndex{2}},
                       settings::Settings = defaultSettings)
    if settings.toroid
        return wrapCoordsOnToroid(coords, settings)
    else
        return removeCoordsOutsideBounds(coords, settings)
    end
end

"""
    removeCoordsOutsideBounds(coords::AbstractVector{CartesianIndex{2}})

Remove coordinates outside of the lattice size.
"""
function removeCoordsOutsideBounds(coords::AbstractVector{CartesianIndex{2}},
                                   settings::Settings = defaultSettings)
    filter(i -> all((1,1) .<= i.I .<= settings.latticeSize), coords)
end

"""
    wrapCoordsOnToroid(coords::AbstractVector{CartesianIndex{2}},

Wrap coordinates outside of the lattice size around the torus.
"""
function wrapCoordsOnToroid(coords::AbstractVector{CartesianIndex{2}},
                            settings::Settings = defaultSettings)
    # like mod, but returns y instead of 0
    mod_replace_zero(x,y) = if (m = mod(x,y)) == 0 y else m end
    map(i -> CartesianIndex(mod_replace_zero.(i.I,settings.latticeSize)), coords)
end

"""
    latticeDistance(a::CartesianIndex{2}, b::CartesianIndex{2},

Euclidean distance on the lattice.
"""
function latticeDistance(a::CartesianIndex{2}, b::CartesianIndex{2},
                         settings::Settings = defaultSettings)
    diff = abs.((a - b).I)
    if settings.toroid
        diff = min.(diff, settings.latticeSize .- diff)
    end
    return sqrt(sum(diff))
end

@inline neighbourhoodKernel(kernel::Symbol, radius::Float64) =
    neighbourhoodKernel(Val(kernel), radius)

"""
    neighbourhoodKernel(::Val{:cone}, radius::Float64)

Returns a linear neighbourhood decay function.
"""
function neighbourhoodKernel(::Val{:cone}, radius::Float64)
    (dist::Float64) -> (radius - dist)/radius
end

"""
    neighbourhoodKernel(::Val{:gauss}, radius::Float64)

Returns a gaussian neighbourhood decay function.
"""
function neighbourhoodKernel(::Val{:gauss}, radius::Float64)
    (dist::Float64) -> exp(-(dist/radius)^2/2)
end

"""
    neighbourhoodKernel(::Val{:mexhat}, radius::Float64)

Returns a "Mexican hat" (Ricker wavelet) neighbourhood decay function.
"""
function neighbourhoodKernel(::Val{:mexhat}, radius::Float64)
    (dist::Float64) -> begin
        square = (dist/radius)^2
        (1 - square) * exp(-square/2)
    end
end


text:
"""
    OptSummary{T<:AbstractFloat}

Summary of an `NLopt` optimization. Idea and code taken from
[`MixedModels`](https://github.com/JuliaStats/MixedModels.jl).
`T` is the type of the function argument(s) and of the function value.
"""
mutable struct OptSummary{T<:AbstractFloat}
    "copy of initial param values in the optimization"
    initial::Vector{T}
    "lower bounds on the parameter values"
    lowerbd::Vector{T}
    "relative tolerance on the function value f(x), as in NLopt"
    ftol_rel::T
    "absolute tolerance on the function value f(x), as in NLopt"
    ftol_abs::T
    "relative tolerance on the argument value x, as in NLopt"
    xtol_rel::T
    "absolute tolerance on the argument value x"
    xtol_abs::Vector{T}
    "initial step sizes, as in NLopt"
    initial_step::Vector{T}
    "maximum number of function evaluations, as in NLopt"
    maxfeval::Int
    "copy of the final parameter values from the optimization"
    final::Vector{T}
    "final value of the objective function"
    fmin::T
    "number of function evaluations"
    feval::Int
    "name of the optimizer used, as a `Symbol`"
    algorithm::Symbol
    "return value, as a `Symbol`"
    returnvalue::Symbol
end
function OptSummary(initial::Vector{T}, lowerbd::Vector{T}, algorithm::Symbol;
        ftol_rel::T = zero(T),
        ftol_abs::T = zero(T),
        xtol_rel::T = zero(T),
        xtol_abs::Vector{T} = fill(zero, length(initial)),
        initial_step::Vector{T} = T[]) where {T<:AbstractFloat}
    OptSummary(initial, lowerbd,
        ftol_rel, ftol_abs, xtol_rel, xtol_abs,
        initial_step,
        -1, # maxeval: 0 or negative for no limit
        copy(initial), # final x
        T(Inf), # final value (will be finite after minimization)
        -1, # number of evals
        algorithm,
        :FAILURE)
end

function NLopt.Opt(optsum::OptSummary)
    lb = optsum.lowerbd
    np = length(lb) # number of parameters to optimize
    opt = NLopt.Opt(optsum.algorithm, np)
    NLopt.ftol_rel!(opt, optsum.ftol_rel) # relative criterion on objective
    NLopt.ftol_abs!(opt, optsum.ftol_abs) # absolute criterion on objective
    NLopt.xtol_rel!(opt, optsum.xtol_rel) # relative criterion on parameters
    NLopt.xtol_abs!(opt, optsum.xtol_abs) # absolute criterion on parameters
    NLopt.lower_bounds!(opt, lb)
    NLopt.maxeval!(opt, optsum.maxfeval)
    if isempty(optsum.initial_step)
        optsum.initial_step = NLopt.initial_step(opt, similar(lb))
    else
        NLopt.initial_step!(opt, optsum.initial_step)
    end
    return opt
end


text:
function turing_inference(
    prob::DiffEqBase.DEProblem,
    alg,
    t,
    data,
    priors;
    likelihood_dist_priors = [InverseGamma(2, 3)],
    likelihood = (u,p,t,σ) -> MvNormal(u, σ[1]*ones(length(u))),
    num_samples=1000, sampler = Turing.NUTS(0.65),
    syms = [Turing.@varname(theta[i]) for i in 1:length(priors)],
    kwargs...,
)
    N = length(priors)
    Turing.@model mf(x, ::Type{T} = Float64) where {T <: Real} = begin
        theta = Vector{T}(undef, length(priors))
        for i in 1:length(priors)
            theta[i] ~ NamedDist(priors[i], syms[i])
        end
        σ = Vector{T}(undef, length(likelihood_dist_priors))
        for i in 1:length(likelihood_dist_priors)
            σ[i] ~ likelihood_dist_priors[i]
        end
        p_tmp = remake(prob, u0 = convert.(T, (prob.u0)), p = theta)
        _saveat = t === nothing ? Float64[] : t
        sol_tmp = solve(p_tmp, alg; saveat = _saveat, kwargs...)

        if sol_tmp isa DiffEqBase.AbstractEnsembleSolution
            failure = any((s.retcode != :Success for s in sol_tmp)) && any((s.retcode != :Terminated for s in sol_tmp))
        else
            failure = sol_tmp.retcode != :Success && sol_tmp != :Terminated
        end

        if failure
            @logpdf() = -Inf
            return
        end
        if sol_tmp isa DiffEqBase.AbstractNoTimeSolution
            res = sol_tmp.u
            x ~ likelihood(res, theta, Inf, σ)
        else
            u = sol_tmp.u
            for i = 1:length(t)
                res = sol_tmp.u[i]
                _t = sol_tmp.t[i]
                x[:,i] ~ likelihood(res, theta, _t, σ)
            end
        end
        return
    end

    # Instantiate a Model object.
    model = mf(data)
    chn = sample(model, sampler, num_samples)
    return chn
end


text:
"""
Streaming Interface of Automa.jl.

NOTE: This module is still experimental. The behavior may change without
deprecations.
"""
module Stream

import Automa
import TranscodingStreams: TranscodingStream

"""
    @mark()

Mark at the current position.

Note: `mark(stream)` doesn't work as expected because the reading position is
not updated while scanning the stream.
"""
macro mark()
    esc(:(buffer.markpos = p))
end

"""
    @markpos()

Get the markerd position.
"""
macro markpos()
    esc(:(buffer.markpos))
end

"""
    @relpos(pos)

Get the relative position of the absolute position `pos`.
"""
macro relpos(pos)
    esc(:(@assert buffer.markpos > 0; $(pos) - buffer.markpos + 1))
end

"""
    @abspos(pos)

Get the absolute position of the relative position `pos`.
"""
macro abspos(pos)
    esc(:(@assert buffer.markpos > 0; $(pos) + buffer.markpos - 1))
end

"""
    generate_reader(funcname::Symbol, machine::Automa.Machine; kwargs...)

Generate a streaming reader function of the name `funcname` from `machine`.

The generated function consumes data from a stream passed as the first argument
and executes the machine with filling the data buffer.

This function returns an expression object of the generated function.  The user
need to evaluate it in a module in which the generated function is needed.

# Keyword Arguments
- `arguments`: Additional arguments `funcname` will take (default: `()`).
    The default signature of the generated function is `(stream::TranscodingStream,)`,
    but it is possible to supply more arguments to the signature with this keyword argument.
- `context`: Automa's codegenerator (default: `Automa.CodeGenContext()`).
- `actions`: A dictionary of action code (default: `Dict{Symbol,Expr}()`).
- `initcode`: Initialization code (default: `:()`).
- `loopcode`: Loop code (default: `:()`).
- `returncode`: Return code (default: `:(return cs)`).

The generated code looks like this:
```julia
function {funcname}(stream::TranscodingStream, {arguments}...)
    buffer = stream.state.buffer1
    data = buffer.data
    {declare variables used by the machine}
    {initcode}
    @label __exec__
    {set up the variables and the data buffer}
    {execute the machine}
    {loopcode}
    if cs ≤ 0 || p > p_eof ≥ 0
        @label __return__
        {returncode}
    end
    @goto __exec__
end
```
"""
function generate_reader(
        funcname::Symbol,
        machine::Automa.Machine;
        arguments::Tuple=(),
        context::Automa.CodeGenContext=Automa.CodeGenContext(),
        actions::Dict{Symbol,Expr}=Dict{Symbol,Expr}(),
        initcode::Expr=:(),
        loopcode::Expr=:(),
        returncode::Expr=:(return cs))
    if returncode.head != :return
        returncode = Expr(:return, returncode)
    end
    functioncode = :(function $(funcname)(stream::$(TranscodingStream)) end)
    for arg in arguments
        push!(functioncode.args[1].args, arg)
    end
    functioncode.args[2] = quote
        buffer = stream.state.buffer1
        data = buffer.data
        $(Automa.generate_init_code(context, machine))
        $(initcode)

        @label __exec__
        if p_eof ≥ 0 || eof(stream)
            p_eof = buffer.marginpos - 1
        end
        p = buffer.bufferpos
        p_end = buffer.marginpos - 1
        $(Automa.generate_exec_code(context, machine, actions))
        Base.skip(stream, p - buffer.bufferpos)

        $(loopcode)

        if cs ≤ 0 || p > p_eof ≥ 0
            @label __return__
            $(returncode)
        end
        @goto __exec__
    end
    return functioncode
end

end  # module


text:
## in evalSR.jl:
#
#   member-assignment validieren
#
## in fctSR.jl den Gewichtungs-Parameter
#
#   AIC Validierung?
#   (oder innere-/äussere Distanzen?)
#
#

max_stars_repo_name: PGSch/SuperResolution

max_issues_repo_name: PGSch/SuperResolution

max_forks_repo_name: PGSch/SuperResolution


text:
function testmul(A,x)
    tmul!(similar(x), A, x)
end

# y .= A * x
tmul!(y::AbstractVector, A::SparseMatrixCSCSymmHerm, x::AbstractVector) = tmul!(y, A, x, 1, 0)

# C .= α * A * B + β * C
function tmul!(C::StridedVecOrMat{T}, sA::SparseMatrixCSCSymmHerm, B::StridedVecOrMat,
              α::Number, β::Number) where T

    fuplo = sA.uplo == 'U' ? nzrangeup : nzrangelo
    _tmul!(fuplo, C, sA, B, T(α), T(β))
end

_mul!(::Val{'U'}, fadj, C, sA, B, α, β) = _mul!(nzrangeup, fadj, C, sA, B, α, β)
_mul!(::Val{'L'}, fadj, C, sA, B, α, β) = _mul!(nzrangelo, fadj, C, sA, B, α, β)

function _tmul!(nzrang::Function, C::StridedVecOrMat{T}, sA, B, α, β) where T
    A = sA.data
    n = size(A, 2)
    m = size(B, 2)
    n == size(B, 1) == size(C, 1) && m == size(C, 2) || throw(DimensionMismatch())
    rv = rowvals(A)
    nzv = nonzeros(A)
    let z = T(0), sumcol=z, αxj=z, aarc=z, α = α #, row
    # row = 0 # zero(eltype(rv))
    if β != 1
        β != 0 ? rmul!(C, β) : fill!(C, z)
    end
    @inbounds for k = 1:m
        for col1 = 1:n
            col = col1
            αxj = B[col,k] * α
            sumcol = z
            for j = nzrang(A, col)
                row = rv[j]
                aarc = nzv[j]
                if row == col
                    sumcol += (sA isa Hermitian ? real : identity)(aarc) * B[row,k]
                else #if (row < col) == ul
                    C[row,k] += aarc * αxj
                    sumcol += (sA isa Hermitian ? adjoint : transpose)(aarc) * B[row,k]
                end
            end
            C[col,k] += α * sumcol
        end
    end
    end
    C
end

#nzrangeup(A, i) = nzrangeup(A, rowvals(A), i)
#nzrangelo(A, i) = nzrangelo(A, rowvals(A), i)

nzrangeup(A, i) = nzrangeup(A, nzrange(A, i), i)
nzrangelo(A, i) = nzrangelo(A, nzrange(A, i), i)
function nzrangeup(A, r::AbstractUnitRange, i)
    r1 = r.start; r2 = r.stop
    rv = rowvals(A)
    r1:searchsortedlast(rowvals(A), i, r1, r2, Forward)
    # @inbounds r2 < r1 || rv[r2] <= i ? r : r1:searchsortedlast(rv, i, r1, r2, Forward)
end
function nzrangeup(A, r::AbstractVector{<:Integer}, i)
    view(r, 1:searchsortedlast(view(rowvals(A), r), i))
end
function nzrangelo(A, r::AbstractUnitRange, i)
    r1 = r.start; r2 = r.stop
    rv = rowvals(A)
    # searchsortedfirst(rv, i, r1, r2, Forward):r2
    @inbounds r2 < r1 || rv[r1] >= i ? r : searchsortedfirst(rv, i, r1, r2, Forward):r2
end
function nzrangelo(A, r::AbstractVector{<:Integer}, i)
    view(r, searchsortedfirst(view(rowvals(A), r), i):length(r))
end

function spmatmul_orig(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
                  sortindices::Symbol = :sortcols) where {Tv,Ti}
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)

    @inbounds begin
        ip = 1
        xb = zeros(Ti, mA)
        x  = zeros(Tv, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k] != i
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = i
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
                        for vp in colptrC[i]:(ip - 1)
                nzvalC[vp] = x[rowvalC[vp]]
            end
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = SparseArrays.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C
end

function spmatmul_mod(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
                  sortindices::Symbol = :sortcols) where {Tv,Ti}
    mA, nA = size(A)
    mB, nB = size(B)
    nA==mB || throw(DimensionMismatch())

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval
    colptrB = B.colptr; rowvalB = B.rowval; nzvalB = B.nzval
    # TODO: Need better estimation of result space
    # nnzC = min(mA*nB, length(nzvalA) + length(nzvalB))
    nnzC = estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)

    @inbounds begin
        ip = 1
        xb = fill(false, mA)
        x  = Vector{Tv}(undef, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                resize!(rowvalC, nnzC + max(nnzC,mA))
                resize!(nzvalC, nnzC + max(nnzC,mA))
                nnzC = length(nzvalC)
            end
            colptrC[i] = ip
            for jp in colptrB[i]:(colptrB[i+1] - 1)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in colptrA[j]:(colptrA[j+1] - 1)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if !xb[k]
                        rowvalC[ip] = k
                        ip += 1
                        xb[k] = true
                        x[k] = nzC
                    else
                        x[k] += nzC
                    end
                end
            end
            for vp in colptrC[i]:(ip - 1)
                nzvalC[vp] = x[rowvalC[vp]]
                xb[rowvalC[vp]] = false
            end
        end
        colptrC[nB+1] = ip
    end

    deleteat!(rowvalC, colptrC[end]:length(rowvalC))
    deleteat!(nzvalC, colptrC[end]:length(nzvalC))

    # The Gustavson algorithm does not guarantee the product to have sorted row indices.
    Cunsorted = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    C = SparseArrays.sortSparseMatrixCSC!(Cunsorted, sortindices=sortindices)
    return C
end

# Gustavsen's matrix multiplication algorithm revisited
function spmatmul_alt(A::SparseMatrixCSC{Tv,Ti}, B::SparseMatrixCSC{Tv,Ti};
                  sortindices::Symbol = :sortcols) where {Tv,Ti}

    mA, nA = size(A)
    nB = size(B, 2)
    nA == size(B, 1) || throw(DimensionMismatch())

    rowvalA = rowvals(A); nzvalA = nonzeros(A)
    rowvalB = rowvals(B); nzvalB = nonzeros(B)
    nnzC = estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)
    nzpercol = nnzC ÷ max(nB, 1)

    @inbounds begin
        ip = 1
        #x  = Vector{Tv}(undef, mA)
        xb = fill(false, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                nnzC += max(mA, nnzC>>2)
                resize!(rowvalC, nnzC)
                resize!(nzvalC, nnzC)
            end
            ip0 = ip
            colptrC[i] = ip
            for jp in nzrange(B, i)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in nzrange(A, j)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k]
                        nzvalC[k+ip0-1] += nzC
                    else
                        nzvalC[k+ip0-1] = nzC
                        xb[k] = true
                        ip += 1
                    end
                end
            end
            if ip > ip0
                vp = ip0
                for k = 1:mA
                    if xb[k]
                        rowvalC[vp] = k
                        nzvalC[vp] = nzvalC[k+ip0-1]
                        xb[k] = false
                        vp += 1
                    end
                end
            end
        end
        colptrC[nB+1] = ip
    end

    resize!(rowvalC, ip - 1)
    resize!(nzvalC, ip - 1)

    # This modification of Gustavson algorithm has sorted row indices
    SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
end
# Gustavsen's matrix multiplication algorithm revisited
function spmatmul(A::SparseMatrixCSCUnion{Tv,Ti}, B::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
    mA, nA = size(A)
    nB = size(B, 2)
    nA == size(B, 1) || throw(DimensionMismatch())

    rowvalA = rowvals(A); nzvalA = nonzeros(A)
    rowvalB = rowvals(B); nzvalB = nonzeros(B)
    nnzC = max(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10, mA)
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)
    nzpercol = nnzC ÷ max(nB, 1)

    @inbounds begin
        ip = 1
        xb = fill(false, mA)
        x = Vector{Tv}(undef, mA)
        for i in 1:nB
            if ip + mA - 1 > nnzC
                nnzC += max(mA, nnzC>>2)
                resize!(rowvalC, nnzC)
                resize!(nzvalC, nnzC)
            end
            colptrC[i] = ip0 = ip
            k0 = ip - 1
            for jp in nzrange(B, i)
                nzB = nzvalB[jp]
                j = rowvalB[jp]
                for kp in nzrange(A, j)
                    nzC = nzvalA[kp] * nzB
                    k = rowvalA[kp]
                    if xb[k]
                        #nzvalC[k+k0] += nzC
                        x[k] += nzC
                    else
                        #nzvalC[k+k0] = nzC
                        x[k] = nzC
                        xb[k] = true
                        rowvalC[ip] = k
                        ip += 1
                    end
                end
            end
            nnzi = ip - ip0
            if nnzi > 0
                if prefer_sort(ip-k0, mA)
                    alg = nnzi <= Base.Sort.SMALL_THRESHOLD+1 ? InsertionSort : QuickSort
                    #sort!(rowvalC, ip0, ip-1, QuickSort, Base.Order.Forward)
                    sort!(rowvalC, ip0, ip-1, alg, Base.Order.Forward)
                    for vp = ip0:ip-1
                        k = rowvalC[vp]
                        xb[k] = false
                        #nzvalC[vp] = nzvalC[k+k0]
                        nzvalC[vp] = x[k]
                    end
                elseif nnzi == mA+2
                    for k = 1:mA
                        xb[k] = false
                        rowvalC[k+k0] = k
                    end
                else
                    for k = 1:mA
                        if xb[k]
                            xb[k] = false
                            rowvalC[ip0] = k
                            #nzvalC[ip0] = nzvalC[k+k0]
                            nzvalC[ip0] = x[k]
                            ip0 += 1
                        end
                    end
                end
            end
        end
        colptrC[nB+1] = ip
    end

    resize!(rowvalC, ip - 1)
    resize!(nzvalC, ip - 1)

    # This modification of Gustavson algorithm has sorted row indices
    C = SparseMatrixCSC(mA, nB, colptrC, rowvalC, nzvalC)
    return C
end

# estimated number of non-zeros in matrix product
# it is assumed, that the non-zero indices are distributed independently and uniformly
# in both matrices. Over-estimation is possible if that is not the case.
function estimate_mulsize(m::Integer, nnzA::Integer, n::Integer, nnzB::Integer, k::Integer)
    p = (nnzA / (m * n)) * (nnzB / (n * k))
    p >= 1 ? m*k : p > 0 ? Int(ceil(-expm1(log1p(-p) * n)*m*k)) : 0 # (1-(1-p)^n)*m*k
end

# determine if sort! shall be used or the whole column be scanned
# based on empirical data on i7-3610QM CPU
# which measured runtimes of the scanning and the sorting loops of the algorithm.
# The parameters 6 and 3 might be modified.
prefer_sort(nz::Integer, m::Integer) = m > 6 && 3 * ilog2(nz) * nz < m

# minimal number of bits required to represent integer; ilog2(n) >= log2(n)
ilog2(n::Integer) = sizeof(n)<<3 - leading_zeros(n)

function spmatmul(aA::Adjoint{Tv,<:SparseMatrixCSCUnion{Tv,Ti}}, B::SparseMatrixCSCUnion{Tv,Ti}) where {Tv,Ti}
    A = parent(aA)
    mA, nA = size(A)
    nB = size(B, 2)
    nA == size(B, 1) || throw(DimensionMismatch())

    rowvalA = rowvals(A); nzvalA = nonzeros(A)
    rowvalB = rowvals(B); nzvalB = nonzeros(B)
    nnzC = max(estimate_mulsize(mA, nnz(A), nA, nnz(B), nB) * 11 ÷ 10, mA)
    colptrC = Vector{Ti}(undef, nB+1)
    rowvalC = Vector{Ti}(undef, nnzC)
    nzvalC = Vector{Tv}(undef, nnzC)

    @inbounds begin
        akku = Tv(0)
        ip = 1
        x = Vector{Tv}(undef, mA)
        xb = fill(false, mA)
        for i in 1:nB
            if ip + nA - 1 > nnzC
                nnzC += max(mA, nnzC>>2)
                resize!(rowvalC, nnzC)
                resize!(nzvalC, nnzC)
            end
            colptrC[i] = ip
            rangeB = nzrange(B, i)
            isempty(rangeB) && continue
            for j = rangeB
                rB = rowvalB[j]
                x[rB] = nzvalB[j]
                xb[rB] = true
            end
            for k in 1:nA
                valid = false
                for j = nzrange(A, k)
                    rB = rowvalA[j]
                    if xb[rB]
                        nzv = nzvalA[j] * x[rB]
                        if valid
                            akku += nzv
                        else
                            valid = true
                            akku = nzv
                        end
                    end
                end
                if valid
                    rowvalC[ip] = k
                    nzvalC[ip] = akku
                    ip += 1
                end
            end
            for j in rangeB
                rB = rowvalB[j]
                xb[rB] = false
            end
        end
        colptrC[nB+1] = ip
    end

    resize!(rowvalC, ip - 1)
    resize!(nzvalC, ip - 1)
    C = SparseMatrixCSC(nA, nB, colptrC, rowvalC, nzvalC)
    return C
end


text:
###########################################################################################################
###########################################################################################################
############################## OVERLAP LIKE G USING RECCURRENCE RELATIONS #################################
###########################################################################################################
###########################################################################################################
export AuxiliaryGr
###########################################################################################################
function AuxiliaryGr(n1::arb, q::Int, n2::arb, n3::arb, p1::arb, p2::arb, p3::arb, lim::Int)
    zero = parent(n1)(0)
    one = parent(n1)(1)
    two = parent(n1)(2)
    three = parent(n1)(3)
    four = parent(n1)(4)

    lowq = Zero(lim)
    upq1 = q
    upq2 = 2*q

    Grlrm = zeros(RF, upq2 - lowq + 2, upq2 - lowq + 2)
    Grlrp = zeros(RF, upq2 - lowq + 2, upq2 - lowq + 2)
    Grh = zeros(RF, upq2 - lowq + 2, upq2 - lowq + 2)

    GR = zeros(RF, upq2 - lowq + 1)

    Grlrm[1, 1] = real(AuxiliaryGrlm(n1, RF(n2 + 2*q), n3, p1, p3, p2))
    Grlrm[1, 2] = real(AuxiliaryGrlm(n1, n2 + 2*q, n3 + 10//10, p1, p3, p2))

    Grlrp[1, 1] = real(AuxiliaryGrlp(n1, n3, n2 + 2*q, p1, p3, p2))
    Grlrp[2, 1] = real(AuxiliaryGrlp(n1, n3 + 10//10, n2 + 2*q, p1, p3, p2))

    Grh[1, 1] = real(AuxiliaryGrh(n1, n2 + 2*q, n3, p1, p2, p3, lim))
    Grh[1, 2] = real(AuxiliaryGrh(n1, n2 + 2*q, n3 + 10//10, p1, p2, p3, lim))

    GR[1] = AuxiliaryG(n1, 0, n2 + 2*q, n3, p1, p2, p3, lim)

    #GR[2] = (
    #((p2 + p3) // (p2 - p3)) * ((n3 + 1) // (n2)) * GR[1] +
    #((p2) // (p2 - p3)) * ((1) // (n2)) * (
    #Grlrm[1, 2] - Grlrp[2, 1]
    #) +
    #((p3) // (p2 - p3)) * Grh[1, 2]
    #)

    for s1 in lowq : upq2 - 1
        s2 = s1 - lowq + 1

        Grlrm[s2 + 1, s2 + 1] = (
        ((p2) // (n2 + 2*q - s2 + 1)) * Grlrm[s2, s2 + 1]
        -
        ((n3 + s2) // (n2 + 2*q - s2 + 1)) * Grlrm[s2, s2]
        )

        Grlrm[s2 + 1, s2 + 2] = Grlrm[s2, s2 + 1] - two * Grlrm[s2 + 1, s2 + 1]
        #################
        #################
        Grlrp[s2 + 1, s2 + 1] = (
        ((p2) // (n2 + 2*q - s2 + 1)) * Grlrp[s2 + 1, s2]
        -
        ((n3 + s2) // (n2 + 2*q - s2 + 1)) * Grlrp[s2, s2]
        )

        Grlrp[s2 + 2, s2 + 1] = Grlrp[s2 + 1, s2] + two * Grlrp[s2 + 1, s2 + 1]
        #################
        #################
        Grh[s2 + 1, s2 + 1] = (
        ((p3) // (n2 + 2*q - s2 + 1)) * Grh[s2, s2 + 1]
        +
        ((n3 + s2) // (n2 + 2*q - s2 + 1)) * Grh[s2, s2]
        )

        Grh[s2 + 1, s2 + 2] = -Grh[s2, s2 + 1] + two * Grh[s2 + 1, s2 + 1]
        #################
        #################
        GR[s2 + 1] = (
        ((p2 + p3) // (p2 - p3)) * ((n3 + s2) // (n2 + 2*q - s2 + 1)) * GR[s2] +
        ((p2) // (p2 - p3)) * ((1) // (n2 + 2*q - s2 + 1)) * (
        Grlrm[s2, s2 + 1] - Grlrp[s2 + 1, s2]
        ) +
        ((p3) // (p2 - p3)) * ((1) // (n2 + 2*q - s2 + 1)) * Grh[s2, s2 + 1]
        )
    end

    GRR = zeros(RF, upq1 - lowq + 1)
    for s1 in lowq : upq1
        GRR[s1 + 1] = Power((-10//10),(s1)) * Binomial(q,s1) * GR[2*s1 + 1]
    end
    ((RF(10//10))//(RF(20//10)^(RF(2 * q)))) * sum(GRR)
end


text:
"""
    PopGen.read(infile::String; kwargs...)
Wraps the individual file importers to read a file in as a `PopData` object. File type is
inferred from the file extension (case insensitive): \n

| File Format         | Extensions             | Docstring     |
| :------------------ | :--------------------- | :------------ |
| delimited           | `.csv`, `.txt`, `.tsv` | `?delimited`  |
| genepop             | `.gen`, `.genepop`     | `?genepop`    |
| structure           | `.str`, `.structure`   | `?structure`  |
| plink               | `.bed`, `.ped`  | `?plink`  |
| variant call format (vcf) | `.vcf`, `.vcf.gz`| `?vcf`  |
| variant call format (bcf) | `.bcf`, `.bcf.gz`| `?bcf`  |

This function uses the same keyword arguments (and defaults) as the file importing
functions it wraps; please see their respective docstrings in the Julia help console.
for specific usage details (e.g. `?genepop`).

## Examples
```
PopGen.read("cavernous_assfish.gen", digits = 3)

PopGen.read("juglans_nigra.vcf")
```
"""
function read(infile::String; kwargs...)
    ext = splitext(infile)[2] |> lowercase

    if occursin(r".vcf$", infile) | occursin(r".vcf.gz$", infile)
        return vcf(infile; kwargs...)

    elseif occursin(r".bcf$", infile) | occursin(r".bcf.gz$", infile)
        return bcf(infile; kwargs...)

    elseif ext in [".gen", ".genepop"]
        return genepop(infile;kwargs...)

    elseif ext in [".csv", ".txt", ".tsv"]
        return delimited(infile; kwargs...)

    elseif ext in [".str", ".structure"]
        return structure(infile; kwargs...)

    elseif ext in [".bed", ".ped"]
       return plink(infile; kwargs...)

    else
        @error "File type not recognized by filename extension \n delimited: .csv | .tsv | .txt \n genepop: .gen | .genepop \n structure: .str | .structure \n variant call formant: .bcf | .vcf \n plink: .bed | .ped | .fam | .bim | .map"
    end
end

"""
    PopGen.write(data::PopData, filename::String, kwargs...)
    PopGen.write(data::PopData; filename::String, kwargs...)

Writes `PopData` to a specified file type inferred from the extension of `filename = ` (case insensitive). Additional keyword
arguments `kwargs...` are specific to the intended file type, and are listed in the docstrings of the specific
file writer with the format `?filetype`. For example, to find the appropriate keywords for a conversion
to Genepop format, call up the `?genepop` docstring.

| File Format | Extensions             | Docstring          |
| :---------- | :--------------------- | :----------------- |
| genepop     | `.gen`, `.genepop`     | ?genepop   |
| delimited   | `.csv`, `.txt`, `.tsv` | ?delimited |
| structure   | `.str`, `.structure`   | ?structure |

### Example
```
cats = @nancycats;
fewer_cats = omit(cats, name = samplenames(cats)[1:10]);
PopGen.write(fewer_cats, filename = "filtered_nancycats.gen", digits = 3, format = "h")
```
"""
function Base.write(data::PopData; filename::String, kwargs...)
    ext = split(filename, ".")[end] |> lowercase
    if ext in ["gen", "genepop"]
        genepop(data, filename = filename; kwargs...)
    elseif ext in ["str", "structure"]
        structure(data, filename = filename; kwargs...)
    elseif ext in ["csv", "txt", "tsv"]
        delimited(data, filename = filename; kwargs...)
    else
        @error "File type not recognized by filename extension. Please see the docstring"
    end
end

function Base.write(data::PopData, filename::String, kwargs...)
    ext = split(filename, ".")[end] |> lowercase
    if ext in ["gen", "genepop"]
        genepop(data, filename = filename; kwargs...)
    elseif ext in ["str", "structure"]
        structure(data, filename = filename; kwargs...)
    elseif ext in ["csv", "txt", "tsv"]
        delimited(data, filename = filename; kwargs...)
    else
        @error "File type not recognized by filename extension. Please see the docstring"
    end
end


text:
function broadcastmsg(roomname::String, msg)
    if !haskey(RoomDict, roomname)
        #Shouldn't Happen
        @warn("Called broadcastmsg with non-existent $(roomname) with $(msg)")
        return
    end
    return broadcastmsg(RoomDict[roomname], msg)
end

function broadcastmsg(listofws::Set{WebSocket}, msg;
                      cond = x -> true)
    for ws in listofws
        if cond(ws)
            writeguarded(ws, msg)
        end
    end
    nothing
end

function errormsg(ws, errortype)
    writeguarded(ws, "{ \"responsetype\" : \"$(errortype)-error\" }")
    nothing
end


text:
struct LMHMM
        hmm::HMM
        lm::LM # RNN based language model
end

import Base.+
+(hmm::HMM, lm::LM) = LMHMM(hmm, lm)

LMHMM() = Knet.load(KongYiji.load("lmhmm.jld2"), "lmhmm")

struct Ts
        val::Tv
        ic::Int
        ip::Int
        bp::Union{Ts, Nothing}
        hc
end

import Base.isless
isless(a::Ts, b::Ts) = a.val < b.val

function (m::LMHMM)(xs::Vector{String}; beams=100, o...)
        nc_max = mapreduce(ncodeunits, max, xs)
        np = length(m.hmm.tags)
        beams = max(beams, np)
        #Th = Heap{Ts{typeof(hiddens(m.lm))}}
        Th = Heap{Ts}

        dp = Vector{Th}(undef, nc_max + 1)
        for i = 1:nc_max + 1; dp[i] = Th(beams); end
        for i = 1:np; push!(dp[1], Ts(0, 1, i, nothing, (0,0))); end

        return [m(x, dp; o...) for x in xs]
end

function (m::LMHMM)(x::String; o...)
        return m([x]; o...)[1]
end

function (m::LMHMM)(x::String, dp::T; recover=true, withpos=false) where T
        @show x
        origin = x
        x = mask(x)
        chrs = codeunits(x)
        vtxs = collect(eachmatch(m.hmm.dict, chrs))
        sort!(vtxs)
        nv, nc, np = length(vtxs), length(chrs), size(dp, 1)
        pv = 1
        preic = 0
        for i = 2:nc + 1; empty!(dp[i]); end
        for ic = 1:nc + 1
                #@show ic, length(dp[ic])
                while pv <= nv && vtxs[pv].s < ic; pv += 1 end

                if !isempty(dp[ic])
                        preic = ic
                elseif (ic == nc + 1 || pv <= nv && vtxs[pv].s == ic)
                        jc = preic
                        for js in dp[jc]
                                is = next(m, js, 0, ic)
                                push!(dp[ic], is)
                        end
                        #empty!(dp[preic]) #TODO
                end

                if pv > nv || ic < vtxs[pv].s || ic == nc + 1; continue; end

                while pv <= nv && vtxs[pv].s == ic
                        jc = ic + length(vtxs[pv])
                        for is in dp[ic]
                                js = next(m, is, vtxs[pv].i, jc)
                                push!(dp[jc], js)
                        end
                        pv += 1
                end
        end

        t = findmax(dp[nc + 1].h.valtree) |> first
        cur = t
        words = String[]
        while !isnothing(cur.bp)
                pre = cur.bp
                push!(words, x[pre.ic:prevind(x, cur.ic)])
                cur = pre
        end
        reverse!(words)
        if recover; words = demask(words, origin); end
        ret = words

        if withpos
                cur = t
                postags = String[]
                while !isnothing(cur.bp)
                        pre = cur.bp
                        push!(postags, m.hmm.postags[cur.ip])
                        cur = pre
                end
                reverse!(postags)
                ret = [join(p, '/') for p in zip(words, postags)]
        end

        return ret
end

function next(m::LMHMM, s::Ts, iw::Int, ic::Int)
        val = logph2v(m.hmm, s.ip, iw)
        reset!(m.lm, s.hc)
        #@show typeof(s.hc)
        probp = m.lm(s.ip) |> vec
        #@show typeof(probp)
        ip = findmax(probp) |> last
        #@show typeof(ip)

        hc = hiddens(m.lm) |> deepcopy
        Ts(val, ic, ip, s, hc)
end

#### Utils

function ==(a::LMHMM, b::LMHMM)
        return all(fname -> getfield(a, fname) == getfield(b, fname), fieldnames(LMHMM))
end


text:
# LookupResult
dispose(x::LookupResult) = clang_LookupResult_dispose(x)

function clear(x::LookupResult, kind::CXLookupNameKind)
    @check_ptrs x
    return clang_LookupResult_clear(x, kind)
end

function dump(x::LookupResult)
    @check_ptrs x
    return clang_LookupResult_dump(x)
end

function empty(x::LookupResult)
    @check_ptrs x
    return clang_LookupResult_empty(x)
end

function getLookupName(x::LookupResult)
    @check_ptrs x
    return DeclarationName(clang_LookupResult_getLookupName(x))
end

function setLookupName(x::LookupResult, dn::DeclarationName)
    @check_ptrs x
    return clang_LookupResult_setLookupName(x, dn)
end

function getRepresentativeDecl(x::LookupResult)
    @check_ptrs x
    return NamedDecl(clang_LookupResult_getRepresentativeDecl(x))
end


text:
"""
    Xtal.lattice_sanity_check(M::AbstractMatrix)

Runs checks on a matrix intended to represent basis vectors of a crystal unit cell. Returns
nothing, but warns if the cell vectors form a left-handed coordinate system, and throws an
`AssertionError` if the cell vectors are not linearly independent.
"""
@inline function lattice_sanity_check(M::AbstractMatrix{<:Real})
    # Skip this check for lattices that are zero (meaning unspecified basis)
    if iszero(M)
        return nothing
    end
    @assert det(M) != 0 string(
        "cell vectors are not linearly independent.\n",
        "Matrix contents: ", M
    )
    det(M) < 0  && @warn "cell vectors form a left-handed coordinate system."
    return nothing
end

# TODO: can we use views on vectors of vectors to get a matrix?
function lattice_sanity_check(vs::AbstractVector{<:AbstractVector{<:Real}})
    return lattice_sanity_check(hcat(vs...))
end

"""
    BasisVectors{D}

Collection of `D` basis vectors spanning `D`-dimensional space.

Internally, a `BasisVectors{D}` is represented as an `SVector{D,SVector{D,Float64}}`. This avoids a
length declaration that would be required by an `SMatrix`.

`BasisVectors` can be indexed similarly to a matrix. However, there is one major difference:
indexing with a single value returns an `SVector{D,Float64}`, not a value as would happen with most
subtypes of `AbstractMatrix`.
"""
struct BasisVectors{D}
    vs::SVector{D,SVector{D,Float64}}
    # TODO: there's probably a nicer way to handle this, but whatever
    # Inner constructor for all matrices
    function BasisVectors{D}(M::AbstractMatrix{<:Real}) where D
        lattice_sanity_check(M)
        return new{D}([M[:,n] for n in 1:size(M,2)])
    end
    # And for all vectors of vectors
    function BasisVectors{D}(M::AbstractVector{AbstractVector{<:Real}}) where D
        lattice_sanity_check(M)
        return new{D}([M[:,n] for n in 1:size(M,2)])
    end
    # If SVectors are used directly, this one can be called
    function BasisVectors(vs::SVector{D,<:SVector{D,<:Real}}) where D
        lattice_sanity_check(vs)
        return new{D}(vs)
    end
    # Same for SMatrices
    function BasisVectors(M::SMatrix{D,D,<:Real}) where D
        lattice_sanity_check(M)
        return new{D}(SVector{D,SVector{D,Float64}}(M[:,n] for n in 1:D))
    end
end

BasisVectors(vs::Vararg{AbstractVector{<:Real}, D}) where D = BasisVectors{D}(hcat(vs...))
BasisVectors{D}(vs::AbstractVector{<:AbstractVector}) where D = BasisVectors{D}(hcat(vs...))

# TODO: might this be better off as an inner constructor?
function BasisVectors(M::AbstractMatrix{<:Real})
    return BasisVectors(SVector{D,SVector{D,Float64}}(M[:,n] for n in 1:D))
end

# This should get a vector
Base.getindex(b::BasisVectors, ind) = b.vs[ind]
# This should treat BasisVectors like matrices
Base.getindex(b::BasisVectors, i1, i2) = b[i2][i1]

vectors(b::BasisVectors) = b.vs
matrix(b::BasisVectors{D}) where D = SMatrix{D,D,Float64}(b[m,n] for m in 1:D, n in 1:D)

# This is needed for broadcasting
Base.length(b::BasisVectors) = length(b.vs)
Base.iterate(b::BasisVectors, state) = iterate(b.vs, state)
Base.iterate(b::BasisVectors) = iterate(b, 1)

# TODO: does this make sense?
Base.convert(::Type{SMatrix{D,D,Float64}}, b::BasisVectors{D}) where D = matrix(b)
# TODO: do we need other kinds of conversions?

Base.zero(::Type{BasisVectors{D}}) where D = BasisVectors(zeros(SVector{D,SVector{D,Float64}}))
Base.zero(::BasisVectors{D}) where D = BasisVectors(zeros(SVector{D,SVector{D,Float64}}))
Base.zeros(::Type{BasisVectors{D}}) where D = BasisVectors(zeros(SVector{D,SVector{D,Float64}}))

# Definitions for multiplication/division by a scalar
# TODO: there's probably a more efficient way to do this
Base.:*(s::Number, b::BasisVectors) = BasisVectors(matrix(b) * s)
Base.:*(b::BasisVectors, s::Number) = s * b
Base.:/(b::BasisVectors, s::Number) = BasisVectors(matrix(b) / s)

# And multiplication/division by vectors
Base.:*(b::BasisVectors, v::AbstractVecOrMat) = matrix(b) * v
Base.:*(v::AbstractVecOrMat, b::BasisVectors) = v * matrix(b)
Base.:\(b::BasisVectors, v::AbstractVecOrMat) = matrix(b) \ v

"""
    dual(b::BasisVectors)

Generates the dual lattice defined by a set of basis vectors.

"Dual" refers to the basis generated with the inverse transpose. This does not include factors of τ
that are used crystallographically (τ = 2π).
"""
dual(b::BasisVectors) = BasisVectors(inv(transpose(matrix(b))))

"""
    RealLattice{D}

Describes a real space crystal lattice with primitive and conventional basis vectors.
"""
struct RealLattice{D} <: AbstractLattice{D}
    prim::BasisVectors{D}
    conv::BasisVectors{D}
    # Inner constructor should take any AbstractMatrix{<:Real} as input
    function RealLattice{D}(
        prim::BasisVectors{D},
        conv::BasisVectors{D}
    ) where D
        # Perform checks on the lattice pairs
        @assert all(x -> x - round(Int, x) < TOL_DEF, matrix(prim)\matrix(conv)) string(
            "The larger set of basis vectors is not expressed in terms of integer multiples of",
            "the smaller set of basis vectors."
        )
        new(prim, conv)
    end
end

"""
    ReciprocalLattice{D}

Describes a reciprocal space crystal lattice with primitive and conventional basis vectors.
"""
struct ReciprocalLattice{D} <: AbstractLattice{D}
    prim::BasisVectors{D}
    conv::BasisVectors{D}
    # Inner constructor should take any AbstractMatrix{<:Real} as input
    function ReciprocalLattice{D}(
        prim::BasisVectors{D},
        conv::BasisVectors{D}
    ) where D
        # Perform checks on the lattice pairs
        @assert all(x -> x - round(Int, x) < TOL_DEF, matrix(conv)\matrix(prim)) string(
            "The smaller set of basis vectors is not expressed in terms of integer reciprocals of",
            "the larger set of basis vectors."
        )
        return new(prim, conv)
    end
end

# Get primitive and conventional lattices
# This is the preferred way of doing so
"""
    prim(l::AbstractLattice{D}) -> BasisVectors{D}

Returns the primitive lattice in a `RealLattice` or a `ReciprocalLattice`.
"""
prim(l::AbstractLattice) = l.prim

"""
    conv(l::AbstractLattice{D}) -> BasisVectors{D}

Returns the conventional lattice in a `RealLattice` or a `ReciprocalLattice`.
"""
conv(l::AbstractLattice) = l.conv

# Define the constructor for vectors of vectors
function (Type{AbstractLattice{D}})(
    prim::AbstractVector{<:AbstractVector{<:Real}},
    conv::AbstractVector{<:AbstractVector{<:Real}},
) where D
    return RealLattice{D}(BasisVectors{D}(prim), BasisVectors{D}(conv))
end

# Conversions between real and reciprocal lattices
# It's critical that the `RealLattice` and `ReciprocalLattice` constructors are idempotent
# so that conversion can be completely seamless
RealLattice(latt::RealLattice) = latt
ReciprocalLattice(latt::ReciprocalLattice) = latt

"""
    ReciprocalLattice(latt::RealLattice{D})

Converts a real lattice to its corresponding reciprocal lattice.
"""
function ReciprocalLattice(latt::RealLattice{D}) where D
    return ReciprocalLattice{D}(dual(prim(latt))*2π, dual(conv(latt)*2π))
end

"""
    RealLattice(latt::ReciprocalLattice{D})

Converts a reciprocal lattice to its corresponding real lattice.
"""
function RealLattice(latt::ReciprocalLattice{D}) where D
    return RealLattice{D}(dual(prim(latt))/2π, dual(conv(latt)/2π))
end

Base.convert(::Type{RealLattice}, latt::AbstractLattice) = RealLattice(latt)
Base.convert(::Type{ReciprocalLattice}, latt::AbstractLattice) = ReciprocalLattice(latt)
#=
function convert(::Type{RealLattice{D}}, latt::AbstractLattice{D}) where D
    return RealLattice(latt)
end

function convert(::Type{ReciprocalLattice{D}}, latt::AbstractLattice{D}) where D
    return ReciprocalLattice(latt)
end
=#

"""
    Xtal.lattice_pair_generator_3D(M::AbstractMatrix; prim=false, ctr=:P)

Generates a pair of 3D lattices that are related by common crystallographic transformations.

Returns an `NTuple{2, BasisVectors{3}` with the first matrix containing the primitive basis and the
second containing the conventional basis.
"""
function lattice_pair_generator_3D(M::AbstractMatrix; prim=false, ctr=:P)
    # Without special centering just return the pair of basis vectors
    ctr == :P && return BasisVectors{3}.((M,M))
    # Figure out which type of cell the input matrix defines
    # TODO: Try to ensure this operations favors a c axis pointing along the z direction
    if prim
        Mp = M
        Mc = inv(REDUCTION_MATRIX_3D[ctr]) * M
    else
        Mc = M
        Mp = REDUCTION_MATRIX_3D[ctr] * M
    end
    return BasisVectors{3}.((Mp,Mc))
end

# It appears these next two docstrings are broken!
#=
"""
    RealLattice{3}(M::AbstractMatrix{<:Real}; prim=false, ctr=:P)

Generates a 3-dimensional real space lattice given a set of conventional or primitive basis
vectors and centering information.

By default, inputs are assumed to describe a conventional cell.
"""
=#
function RealLattice{3}(M::AbstractMatrix{<:Real}; prim=false, ctr=:P)
    return RealLattice{3}(lattice_pair_generator_3D(M, prim=prim, ctr=ctr)...)
end

"""
    cell_lengths(M::AbstractMatrix) -> Vector{Float64}

Returns the lengths of the constituent vectors in a matrix representing cell vectors.
"""
cell_lengths(M::AbstractMatrix) = [norm(M[:,n]) for n = 1:size(M,2)]
# TODO: perhaps this is not necessary anymore
# But removing it might break the API
cell_lengths(b::BasisVectors) = cell_lengths(matrix(b))

"""
    lengths(b::BasisVectors) -> Vector{Float64}

Returns the lengths of the constituent vectors in a matrix representing cell vectors.
"""
lengths(b::BasisVectors) = cell_lengths(matrix(b))

"""
    cell_volume(M::AbstractMatrix) -> Float64

Returns the volume of a unit cell defined by a matrix. This volume does not carry the sign
(negative for cells that do not follow the right hand rule).
"""
cell_volume(M::AbstractMatrix) = abs(det(M))
# TODO: perhaps this is not necessary anymore
# But removing it might break the API
cell_volume(b::BasisVectors) = cell_volume(matrix(b))

"""
    volume(b::BasisVectors) -> Float64

Returns the volume of a unit cell defined by a matrix. This volume does not carry the sign
(negative for cells that do not follow the right hand rule).
"""
volume(b::BasisVectors) = cell_volume(matrix(b))

"""
    volume(l::AbstractLattice; primitive=true) -> Float64

Returns the volume of a lattice. By default, the primitive cell volume is used, but the
conventional cell volume may be calculated with `primitive=false`.
"""
volume(l::AbstractLattice; primitive::Bool=true) = volume(primitive ? prim(l) : conv(l))

"""
    Xtal.generate_pairs(D::Integer) -> Vector{NTuple{2,Int}}

Generate pairs of integers up to `D` in ascending order.
"""
function generate_pairs(D::Integer)
    out = Vector{NTuple{2,Int}}(undef, Int(D*(D-1)/2))
    c = 0
    for a = 1:D
        for b = (a+1):D
            c += 1
            out[c] = (a,b)
        end
    end
    return out
end

"""
    Xtal.generate_pairs(::Type{Val{D}}) -> SVector{D*(D-1)/2, NTuple{2,Int}}

Generate pairs of integers up to `D` in ascending order in an `SVector`.
"""
function generate_pairs(::Type{Val{D}}) where D
    N = Int(D*(D-1)/2)
    return SVector{N,NTuple{2,Int}}(generate_pairs(D))
end

"""
    cell_angle_cos(M::AbstractMatrix)

Generates the cosines of the unit cell angles.

The angles are generated in the correct order [α, β, γ] for 3-dimensional cells. This is achieved
by reversing the output of `generate_pairs()`. For crystals with more spatial dimensions, this
may lead to unexpected results.
"""
function cell_angle_cos(M::AbstractMatrix)
    dimpairs = reverse(generate_pairs(size(M,1)))
    return [dot(M[:,a], M[:,b])/(norm(M[:,a])*norm(M[:,b])) for (a,b) in dimpairs]
end

"""
    cell_angle_rad(M::AbstractMatrix)

Returns the angles (in radians) between each pair of basis vectors.
"""
cell_angle_rad(M::AbstractMatrix) = acos.(cell_angle_cos(M))
cell_angle_deg(M::AbstractMatrix) = acosd.(cell_angle_cos(M))

"""
    lengths(L::AbstractLattice{D}; prim=false) -> SVector{D,Float64}

Returns the cell vector lengths. By default, returns the lengths of the conventional cell vectors.
"""
function lengths(L::AbstractLattice{D}; prim=false) where D
    prim ? M = prim(L) : M = conv(L)
    return SVector{D,Float64}(cell_lengths(M))
end

"""
    angles_deg(L::AbstractLattice{D}; prim=false) -> SVector{D,Float64}

Returns the cell vector angles in degrees.
"""
function angles_deg(L::AbstractLattice{D}; prim=false) where D
    prim ? M = prim(L) : M = conv(L)
    return SVector{D,Float64}(cell_angles_deg(M))
end

"""
   angles_deg(L::AbstractLattice{D}; prim=false) -> SVector{D,Float64}

Returns the cell vector angles in radians.
"""
function angles_rad(L::AbstractLattice{D}; prim=false) where D
    prim ? M = prim(L) : M = conv(L)
    return SVector{D,Float64}(cell_angles_rad(M))
end

"""
    lattice2D(a::Real, b::Real, γ::Real) -> BasisVectors{2}

Constructs a set of basis vectors in an `SMatrix` that correspond to a unit cell in 2D with the
same length and angle parameters (in degrees).

By default, the b-vector is oriented along y. This selection corresponds to the default orientation
chosen by `lattice3D()`.
"""
function lattice2D(a::Real, b::Real, γ::Real)
    return BasisVectors(SMatrix{2,2,Float64}(a*sind(γ), a*cosd(γ), 0, b))
end

"""
    lattice3D(a::Real, b::Real, c::Real, α::Real, β::Real, γ::Real) -> BasisVectors{3}

Constructs a set of basis vectors in an `SMatrix` that correspond to a unit cell in 3D with the
same length and angle parameters (in degrees).

By default, the b-vector is oriented along y, and the a-vector is chosen to be perpendicular to
z, leaving the c-vector to freely vary. This selection allows for the most convenient orientation
of symmetry operations.
"""
function lattice3D(a::Real, b::Real, c::Real, α::Real, β::Real, γ::Real)
    c1 = c*(cosd(β) - cosd(γ)*cosd(α))/sind(γ)
    c2 = c*cosd(α)
    M = SMatrix{3,3,Float64}(a*sind(γ), a*cosd(γ), 0,
                                     0,        b,  0,
                                c1, c2, sqrt(c^2 - (c1^2 + c2^2)))
    return BasisVectors(M)
end

# TODO: Update and document this function a bit more.
# It works, but that really isn't enough for me.
# I think this was pulled directly from the WaveTrans source code
function maxHKLindex(M::AbstractMatrix{<:Real}, ecut::Real; c = CVASP)
    # I think the parts below convert a set of basis vectors into their reciprocals
    # But we already have ways to do that with `dual(::BasisVectors)`
    #--------------------------------------------------------------------------#
    cosines = cell_angle_cos(M)
    crosses = [cross(M[:,a], M[:,b]) for (a,b) in zip([2,3,1], [3,1,2])]
    triples = [dot(M[:,n], crosses[n])/(norm(crosses[n])*norm(M[:,n])) for n in 1:3]
    sines = hcat([[sqrt(1-c^2), sqrt(1-c^2) ,t] for (c,t) in zip(cosines, triples)]...)
    #--------------------------------------------------------------------------#
    nbmax = sqrt(c*ecut) ./ [norm(M[:,a])*sines[a,b] for a in 1:3, b in 1:3] .+ 1
    return floor.(Int, vec(maximum(nbmax, dims=1)))
end

# Assume that the basis vectors are defined in reciprocal space (??)
maxHKLindex(b::BasisVectors{3}, ecut::Real; c = CVASP) = maxHKLindex(matrix(b), ecut, c = c)

"""
    Xtal.maxHKLindex(L::AbstractLattice, ecut::Real; prim=true, c = CVASP)

Determines the maximum integer values of the reciprocal lattice vectors needed to store data out
to a specific energy cutoff for a 3D lattice.

By default, the energy cutoff is assumed to be in units of eV, the reciprocal lattice vector
lengths are assumed to be in angstroms, and the value of c (2m/ħ^2) is taken from VASP's default
value (which is incorrect!). Different values of c may be used for different units.

The functionality implemented here was taken from WaveTrans:
https://www.andrew.cmu.edu/user/feenstra/wavetrans/
"""
function maxHKLindex(L::AbstractLattice{3}, ecut::Real, prim=true, c = CVASP)
    M = ReciprocalLattice{3}(prim ? prim(L) : conv(L))
    return maxHKLindex(M, ecut, c=c)
end


text:
""" Upper bound computation (prefix is UB for Upper Bound) """

# All functions are public :

"""
Returns an object `parentToChild` : the already known points that won't need to be computed again.

Points we already know are the ones for which the variable `index` is already set to `value`.

Returns parentToChild
"""
function UB_computeParentToChild(prob::BiOKP, UB::UpperBound, index::Int, value::Int)

    CONFIG.debug && DEBUG_UB(prob, UB)

    # known solutions that go from the parent to the child node
	knownSols = nil(Sol)
	knownSolsLastElt = nil(Sol)

    # indicies of the known solutions inside the parent UB.sols
    previousIndicies = nil(Int)
    previousIndiciesLastElt = nil(Int)

    sols = UB.sols

    lenLB = length(sols)

    # attributes of the ParentToChild object : are lexicographic sols known by the child ?
    leftKnown = false
    rightKnown = false

    # right sol to memorize
    rightSol = Sol()

    iter = 1
    atLeastOneSolKnown = false # length(knownSols) >= 1
	while sols != nil(Sol)

		if !atLeastOneSolKnown && sols.head.x[index] == value

            # length(knownSols) >= 1 so atLeastOneSolKnown becomes true
			atLeastOneSolKnown = true

            # append the known solution to knownSols
			knownSols = cons(sols.head, nil(Sol))
			knownSolsLastElt = knownSols

            # append the index of the known solution to previousIndicies
            previousIndicies = cons(iter, nil(Int))
			previousIndiciesLastElt = previousIndicies

            rightSol = sols.head

            # update leftKnown and rightKnown booleans
            if iter == 1 # the left sol is a known solution for the child
                leftKnown = true
            elseif iter == lenLB
                rightKnown = true
            end

		elseif sols.head.x[index] == value

            # append the known solution to knownSols
			knownSolsLastElt.tail = cons(sols.head, nil(Sol))
			knownSolsLastElt = knownSolsLastElt.tail

            # append the index of the known solution to previousIndicies
            previousIndiciesLastElt.tail = cons(iter, nil(Int))
			previousIndiciesLastElt = previousIndiciesLastElt.tail

            rightSol = sols.head

            # update leftKnown and rightKnown booleans
            if iter == 1 # the left sol is a known solution for the child
                leftKnown = true
            elseif iter == lenLB
                rightKnown = true
            end

		end

		sols = sols.tail
        iter += 1
	end

    parentToChild = ParentToChild(knownSols,previousIndicies,rightSol,leftKnown,rightKnown)

    CONFIG.debug && DEBUG_parentToChild(prob, parentToChild, index=index, value=value)

	return parentToChild
end

"""
Computes the segments associated with the given `UB.sols` as a DualSet
"""
function UB_computeSegments!(prob::BiOKP, UB::UpperBound)

    CONFIG.debug && DEBUG_UB(prob, UB)

    sols = UB.sols

	nbConstraints = length(sols) + 1

	A = zeros(Float64, nbConstraints, 2)
	b = zeros(Float64, nbConstraints)

	max1 = sols.head.y[1]
	max2 = sols.head.y[2]

	iter = 1

	while sols.tail != nil(Sol)

        solL = sols.head
        solR = sols.tail.head

		A[iter, 1] = solL.y[2] - solR.y[2]
		A[iter, 2] = solR.y[1] - solL.y[1]

		b[iter] = (solL.y[2] - solR.y[2]) * solR.y[1] + (solR.y[1] - solL.y[1]) * solR.y[2]

		if solL.y[1] > max1
			max1 = solL.y[1]
		end
		if solR.y[1] > max1
			max1 = solR.y[1]
		end
		if solL.y[2] > max2
			max2 = solL.y[2]
		end
		if solR.y[2] > max2
			max2 = solR.y[2]
		end

		sols = sols.tail
		iter += 1
	end

    # add the constraint extreme horizontal and vertical to finish the polytope

    # vertical (on the right)
	A[nbConstraints-1, 1] = 1
	A[nbConstraints-1, 2] = 0

	b[nbConstraints-1] = max1

    # horizontal (on the left)
	A[nbConstraints, 1] = 0
	A[nbConstraints, 2] = 1

	b[nbConstraints] = max2

    UB_setSegments!(UB, DualSet(A, b)) # update segments

    CONFIG.debug && DEBUG_UB(prob, UB)

end


text:
"""
    BinomialLikelihood

Binomial likelihood with fixed number of trials equal to `n`.

```math
    p(y|f) = Binomial(y | n, f)
```
On calling, this would return a Binomial distribution with p=`f` and n=`n`.
"""
struct BinomialLikelihood
    n::Int
end

(l::BinomialLikelihood)(f::Real) = Binomial(l.n, logistic(f))

(l::BinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(Binomial.(l.n, logistic.(fs)))