text:
"""
    Fermi.CoupledCluster.cc_update_energy(T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4}, moints::IntegralHelper{T,O},
                          alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute CC energy from T1 and T2 amplitudes constructed with a set of restricted orbitals.
NOTE: It does not include off-diagonal Fock contributions: See `od_cc_update_energy`.
"""
function cc_update_energy(T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4}, moints::IntegralHelper{T,E,O},
                          alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    Vovov = moints["OVOV"]
    CC_energy = zero(T)
    TWO = T(2)
    @tensoropt (k=>x, l=>x, c=>100x, d=>100x)  begin
        B[l,c,k,d] := TWO*T2[k,l,c,d]
        B[l,c,k,d] -= T1[l,c]*T1[k,d]
        B[l,c,k,d] -= T2[l,k,c,d]
        CC_energy += B[l,c,k,d]*Vovov[k,c,l,d]
        CC_energy += TWO*T1[l,c]*T1[k,d]*Vovov[l,c,k,d]
    end

    return CC_energy
end

"""
    Fermi.CoupledCluster.od_cc_update_energy(T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4}, moints::IntegralHelper{T,O},
                          alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute contributions from off-diagonal Fock terms to the CC energy from T1 and T2 amplitudes constructed with a set of restricted orbitals.
See also: `cc_update_energy`
"""
function od_cc_update_energy(T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4}, moints::IntegralHelper{T,Chonky,O},
                          alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    f = moints["Fia"]
    TWO = T(2)
    CC_energy = cc_update_energy(T1, T2, moints, alg)
    @tensor CC_energy += TWO*f[k,c]*T1[k,c]

    return CC_energy
end
"""
    Fermi.CoupledCluster.cc_update_T1!(newT1::AbstractArray{T,2}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,O},
                      alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute DᵢₐTᵢₐ from old T1 and T2 amplitudes. Final updated T1 amplitudes can be obtained by applying the denominator 1/Dᵢₐ.
NOTE: It does not include off-diagonal Fock contributions: See `od_cc_update_T1`.
"""
function cc_update_T1!(newT1::AbstractArray{T,2}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,E,O},
                      alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = moints["OOOO"], moints["OOOV"], moints["OOVV"], moints["OVOV"], moints["OVVV"], moints["VVVV"]
    TWO = T(2)
    FOUR = T(4)
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT1[i,a] -= T1[k,c]*Voovv[i,k,c,a]
        newT1[i,a] += TWO*T1[k,c]*Vovov[k,c,i,a]
        newT1[i,a] -= T2[k,i,c,d]*Vovvv[k,d,a,c]
        newT1[i,a] += TWO*T2[i,k,c,d]*Vovvv[k,d,a,c]
        newT1[i,a] += -TWO*T2[k,l,a,c]*Vooov[k,i,l,c]
        newT1[i,a] += T2[l,k,a,c]*Vooov[k,i,l,c]
        newT1[i,a] += -TWO*T1[k,c]*T1[l,a]*Vooov[l,i,k,c]
        newT1[i,a] -= T1[k,c]*T1[i,d]*Vovvv[k,d,a,c]
        newT1[i,a] += TWO*T1[k,c]*T1[i,d]*Vovvv[k,c,a,d]
        newT1[i,a] += T1[k,c]*T1[l,a]*Vooov[k,i,l,c]
        newT1[i,a] += -TWO*T1[k,c]*T2[i,l,a,d]*Vovov[l,c,k,d]
        newT1[i,a] += -TWO*T1[k,c]*T2[l,i,a,d]*Vovov[k,c,l,d]
        newT1[i,a] += T1[k,c]*T2[l,i,a,d]*Vovov[l,c,k,d]
        newT1[i,a] += -TWO*T1[i,c]*T2[l,k,a,d]*Vovov[l,c,k,d]
        newT1[i,a] += T1[i,c]*T2[l,k,a,d]*Vovov[k,c,l,d]
        newT1[i,a] += -TWO*T1[l,a]*T2[i,k,d,c]*Vovov[k,c,l,d]
        newT1[i,a] += T1[l,a]*T2[i,k,c,d]*Vovov[k,c,l,d]
        newT1[i,a] += T1[k,c]*T1[i,d]*T1[l,a]*Vovov[l,c,k,d]
        newT1[i,a] += -TWO*T1[k,c]*T1[i,d]*T1[l,a]*Vovov[k,c,l,d]
        newT1[i,a] += FOUR*T1[k,c]*T2[i,l,a,d]*Vovov[k,c,l,d]
    end
end

"""
    Fermi.CoupledCluster.od_cc_update_T1!(newT1::AbstractArray{T,2}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,O},
                      alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute non-diagonal Fock contributions to DᵢₐTᵢₐ from old T1 and T2 amplitudes. Final updated T1 amplitudes
can be obtained by applying the denominator 1/Dᵢₐ.
See also: `cc_update_T1`.
"""
function od_cc_update_T1!(newT1::AbstractArray{T,2}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,E,O},
                      alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    # Include non-RHF terms
    fov = moints["Fia"]
    foo = moints["Fij"]
    fvv = moints["Fab"]
    TWO = T(2)
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT1[i,a] += fov[i,a]
        newT1[i,a] -= foo[i,k]*T1[k,a]
        newT1[i,a] += fvv[c,a]*T1[i,c]
        newT1[i,a] -= fov[k,c]*T1[i,c]*T1[k,a]
        newT1[i,a] += TWO*fov[k,c]*T2[i,k,a,c]
        newT1[i,a] -= fov[k,c]*T2[k,i,a,c]
    end

    # Include RHF terms
    cc_update_T1!(newT1, T1, T2, moints, alg)
end

"""
    Fermi.CoupledCluster.cc_update_T2!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,O},
                    alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute Dᵢⱼₐᵦ⋅Tᵢⱼₐᵦ from old T1 and T2 amplitudes. Final updated T2 amplitudes can be obtained by applying the denominator 1/Dᵢⱼₐᵦ.
NOTE: It does not include off-diagonal Fock contributions: See `od_cc_update_T2`.
"""
function cc_update_T2!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,E,O},
                    alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    Voooo, Vooov, Voovv, Vovov, Vovvv, Vvvvv = moints["OOOO"], moints["OOOV"], moints["OOVV"], moints["OVOV"], moints["OVVV"], moints["VVVV"]

    # Include (VV,VV) term. This is the only part that changes with Density Fitting
    cc_update_T2_v4_term!(newT2, T1, T2, moints, alg)

    TWO = T(2)
    FOUR = T(4)
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        newT2[i,j,a,b] += Vovov[i,a,j,b]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*Voooo[i,k,j,l]
        newT2[i,j,a,b] += T2[k,l,a,b]*Voooo[i,k,j,l]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,a]*Vovvv[k,c,b,d]
        newT2[i,j,a,b] -= T1[i,c]*T1[j,d]*T1[k,b]*Vovvv[k,d,a,c]
        newT2[i,j,a,b] += T1[i,c]*T1[k,a]*T1[l,b]*Vooov[l,j,k,c]
        newT2[i,j,a,b] += T1[j,c]*T1[k,a]*T1[l,b]*Vooov[k,i,l,c]
        newT2[i,j,a,b] += T2[k,l,a,c]*T2[i,j,d,b]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += -TWO*T2[i,k,a,c]*T2[l,j,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += -TWO*T2[l,k,a,c]*T2[i,j,d,b]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,d,b]*Vovov[l,c,k,d]
        newT2[i,j,a,b] += T2[i,k,a,c]*T2[l,j,b,d]*Vovov[l,c,k,d]
        newT2[i,j,a,b] += -TWO*T2[i,k,a,c]*T2[j,l,b,d]*Vovov[l,c,k,d]
        newT2[i,j,a,b] += T2[k,i,a,c]*T2[l,j,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += -TWO*T2[k,i,a,c]*T2[j,l,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T2[i,j,a,c]*T2[l,k,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += -TWO*T2[i,j,a,c]*T2[k,l,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T2[k,j,a,c]*T2[i,l,d,b]*Vovov[l,c,k,d]
        newT2[i,j,a,b] += FOUR*T2[i,k,a,c]*T2[j,l,b,d]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T2[i,j,d,c]*T2[l,k,a,b]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T1[k,a]*T1[l,b]*Vovov[k,c,l,d]
        newT2[i,j,a,b] += T1[i,c]*T1[j,d]*T2[l,k,a,b]*Vovov[l,c,k,d]
        newT2[i,j,a,b] += T1[k,a]*T1[l,b]*T2[i,j,d,c]*Vovov[l,c,k,d]
        P_OoVv[i,j,a,b] := T1[j,c]*Vovvv[i,a,c,b]
        P_OoVv[i,j,a,b] -= T1[k,b]*Vooov[j,k,i,a]
        P_OoVv[i,j,a,b] -= T2[k,i,a,c]*Vovov[k,c,j,b]
        P_OoVv[i,j,a,b] -= T1[i,c]*T1[k,a]*Vovov[k,c,j,b]
        P_OoVv[i,j,a,b] -= T1[i,c]*T1[k,b]*Voovv[j,k,c,a]
        P_OoVv[i,j,a,b] += TWO*T2[i,k,a,c]*Vovov[k,c,j,b]
        P_OoVv[i,j,a,b] -= T2[i,k,a,c]*Voovv[j,k,c,b]
        P_OoVv[i,j,a,b] -= T2[k,j,a,c]*Voovv[i,k,c,b]
        P_OoVv[i,j,a,b] += -TWO*T1[l,b]*T2[i,k,a,c]*Vooov[l,j,k,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[k,i,a,c]*Vooov[l,j,k,c]
        P_OoVv[i,j,a,b] -= T1[j,c]*T2[i,k,d,b]*Vovvv[k,c,a,d]
        P_OoVv[i,j,a,b] -= T1[j,c]*T2[k,i,a,d]*Vovvv[k,d,b,c]
        P_OoVv[i,j,a,b] -= T1[j,c]*T2[i,k,a,d]*Vovvv[k,c,b,d]
        P_OoVv[i,j,a,b] += T1[j,c]*T2[l,k,a,b]*Vooov[l,i,k,c]
        P_OoVv[i,j,a,b] += T1[l,b]*T2[i,k,a,c]*Vooov[k,j,l,c]
        P_OoVv[i,j,a,b] -= T1[k,a]*T2[i,j,d,c]*Vovvv[k,d,b,c]
        P_OoVv[i,j,a,b] += T1[k,a]*T2[i,l,c,b]*Vooov[l,j,k,c]
        P_OoVv[i,j,a,b] += TWO*T1[j,c]*T2[i,k,a,d]*Vovvv[k,d,b,c]
        P_OoVv[i,j,a,b] -= T1[k,c]*T2[i,j,a,d]*Vovvv[k,d,b,c]
        P_OoVv[i,j,a,b] += TWO*T1[k,c]*T2[i,j,a,d]*Vovvv[k,c,b,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T2[i,l,a,b]*Vooov[k,j,l,c]
        P_OoVv[i,j,a,b] += -TWO*T1[k,c]*T2[i,l,a,b]*Vooov[l,j,k,c]
        P_OoVv[i,j,a,b] += T2[j,k,c,d]*T2[i,l,a,b]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += -TWO*T1[k,c]*T1[j,d]*T2[i,l,a,b]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[j,d]*T2[i,l,a,b]*Vovov[l,c,k,d]
        P_OoVv[i,j,a,b] += -TWO*T1[k,c]*T1[l,a]*T2[i,j,d,b]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += T1[k,c]*T1[l,a]*T2[i,j,d,b]*Vovov[l,c,k,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,b,d]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += -TWO*T1[i,c]*T1[k,a]*T2[j,l,b,d]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[k,a]*T2[l,j,d,b]*Vovov[l,c,k,d]
        P_OoVv[i,j,a,b] += T1[i,c]*T1[l,b]*T2[k,j,a,d]*Vovov[k,c,l,d]
        P_OoVv[i,j,a,b] += -TWO*T2[i,k,d,c]*T2[l,j,a,b]*Vovov[k,c,l,d]

        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
end

function cc_update_T2_v4_term!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,Chonky,O},
                    alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

    Vvvvv = moints["VVVV"]
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        τ[i,j,a,b] := T2[i,j,a,b] + T1[i,a]*T1[j,b]
        newT2[i,j,a,b] += τ[i,j,c,d]*Vvvvv[c,a,d,b]
    end
end

function cc_update_T2_v4_term!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,E,O},
                    alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractDFERI, O<:AbstractRestrictedOrbitals}

    Bvv = moints["BVV"]
    @tensor τ[i,j,a,b] := T2[i,j,a,b] + T1[i,a]*T1[j,b]

    o_size = size(T1, 1)
    v_size = size(T1, 2)
    cdb = FermiMDzeros(T, v_size, v_size, v_size)
    newT2a = FermiMDzeros(T, o_size, o_size, v_size)
    for a = 1:v_size
        Ba = @views Bvv[:,:,a]
        @tensor cdb[c,d,b] = Ba[Q,c]*Bvv[Q,d,b]
        @tensor newT2a[i,j,b] = τ[i,j,c,d]*cdb[c,d,b]
        newT2[:,:,a,:] += newT2a
    end
end

"""
    Fermi.CoupledCluster.od_cc_update_T2!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,O},
                    alg::RCCSDa) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Compute non-diagonal Fock contribution to Dᵢⱼₐᵦ⋅Tᵢⱼₐᵦ from old T1 and T2 amplitudes. Final updated T2 amplitudes
can be obtained by applying the denominator 1/Dᵢⱼₐᵦ. See also: `cc_update_T2`.
"""
function od_cc_update_T2!(newT2::AbstractArray{T,4}, T1::AbstractArray{T,2}, T2::AbstractArray{T,4}, moints::IntegralHelper{T,E,O},
                    alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    # Include non-RHF terms
    fov = moints["Fia"]
    foo = moints["Fij"]
    fvv = moints["Fab"]
    ONE = one(T)
    @tensoropt (i=>x, j=>x, k=>x, l=>x, a=>10x, b=>10x, c=>10x, d=>10x) begin
        P_OoVv[i,j,a,b] := -ONE*foo[i,k]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] += fvv[c,a]*T2[i,j,c,b]
        P_OoVv[i,j,a,b] -= fov[k,c]*T1[i,c]*T2[k,j,a,b]
        P_OoVv[i,j,a,b] -= fov[k,c]*T1[k,a]*T2[i,j,c,b]
        newT2[i,j,a,b] += P_OoVv[i,j,a,b] + P_OoVv[j,i,b,a]
    end
    #Include RHF terms
    cc_update_T2!(newT2, T1, T2, moints, alg)
end

"""
    Fermi.CoupledCluster.update_amp!(newT1::AbstractArray{T,2}, newT2::Array{T,4}, T1::Array{T, 2}, T2::Array{T, 4}, moints::IntegralHelper{T,O},
                     alg::A) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Computes new T1 and T2 amplitudes from old ones. It assumes arbitrary restricted orbitals.
"""
function update_amp!(newT1::AbstractArray{T,2}, newT2::AbstractArray{T,4}, T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4},
                    moints::IntegralHelper{T,E,O}, alg::RCCSDa) where {T<:AbstractFloat, E<:AbstractERI, O<:AbstractRestrictedOrbitals}

    # Clean the arrays
    fill!(newT1, 0.0)
    fill!(newT2, 0.0)

    # Get new amplitudes
    od_cc_update_T1!(newT1, T1, T2, moints, alg)
    od_cc_update_T2!(newT2, T1, T2, moints, alg)

    # Orbital energies line
    if haskey(moints.cache, "D1")
        d = moints["D1"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        d = FermiMDArray([ϵo[i]-ϵv[a] for i=eachindex(ϵo), a=eachindex(ϵv)])
        moints["D1"] = d
    end

    if haskey(moints.cache, "D2")
        D = moints["D2"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        D = FermiMDArray([ϵo[i]+ϵo[j]-ϵv[a]-ϵv[b] for i=eachindex(ϵo), j=eachindex(ϵo), a=eachindex(ϵv), b=eachindex(ϵv)])
        moints["D2"] = D
    end

    newT1 ./= d
    newT2 ./= D
end

"""
    Fermi.CoupledCluster.update_amp!(newT1::AbstractArray{T,2}, newT2::Array{T,4}, T1::Array{T, 2}, T2::Array{T, 4}, moints::IntegralHelper{T,O},
                     alg::A) where {T<:AbstractFloat, O<:AbstractRestrictedOrbitals}

Computes new T1 and T2 amplitudes from old ones. It assumes Restricted Hartree-Fock reference.
"""
function update_amp!(newT1::AbstractArray{T,2}, newT2::AbstractArray{T,4}, T1::AbstractArray{T, 2}, T2::AbstractArray{T, 4}, moints::IntegralHelper{T,E,RHFOrbitals},
                     alg::RCCSDa) where {T<:AbstractFloat,E<:AbstractERI}

    # Clean the arrays
    fill!(newT1, 0.0)
    fill!(newT2, 0.0)

    # Get new amplitudes
    cc_update_T1!(newT1, T1, T2, moints, alg)
    cc_update_T2!(newT2, T1, T2, moints, alg)

    # Orbital energies line
    if haskey(moints.cache, "D1")
        d = moints["D1"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        d = FermiMDArray([ϵo[i]-ϵv[a] for i=eachindex(ϵo), a=eachindex(ϵv)])
        moints["D1"] = d
    end

    if haskey(moints.cache, "D2")
        D = moints["D2"]
    else
        Fd = moints["Fd"]
        ndocc = moints.molecule.Nα
        frozen = Options.get("drop_occ")
        inac = Options.get("drop_vir")
        ϵo = Fd[(1+frozen:ndocc)]
        ϵv = Fd[(1+ndocc):end-inac]

        D = FermiMDArray([ϵo[i]+ϵo[j]-ϵv[a]-ϵv[b] for i=eachindex(ϵo), j=eachindex(ϵo), a=eachindex(ϵv), b=eachindex(ϵv)])
        moints["D2"] = D
    end

    # Orbital energies line
    d, D = moints["D1"], moints["D2"]

    newT1 ./= d
    newT2 ./= D
end


text:
using ModelingToolkit
using DataFrames
using CSV
using Tables

# State variable and parameters
@variables t
@variables Rtg13A_c(t) Rtg13I_c(t) Rtg13A_n(t) Rtg13I_n(t) Rtg3I_c(t) Rtg3A_c(t) Rtg3A_n(t) Rtg3I_n(t) Rtg1_c(t) Rtg1_n(t)
@variables s(t) Rtg2I_c(t) Rtg2A_c(t)
@variables Mks(t) Rtg2Mks_c(t) Bmh(t) BmhMks(t)
@parameters ΣRtg1 ΣRtg2 ΣRtg3 ΣMks ΣBmh
@parameters ksV = 1.0 ksD = 1.0 n_S = 1.0 k2I = 1.0 mul_S = 1.0
@parameters k2M = 1.0 kn2M = 1.0 kBM = 1.0 knBM = 1.0
@parameters k13I = 1.0 k13IV = 1.0 k13ID = 1.0 k3A_c = 1.0 k3I_c = 1.0 k3I_n = 1.0
@parameters k13_c = 1.0 kn13_c = 1.0 k13_n = 1.0 kn13_n = 1.0
@parameters k1in = 1.0 k1out = 1.0 k3inA = 1.0 k3outA = 1.0 k3inI = 1.0 k3outI = 1.0

# Distributions of RTG1 and RTG3 proteins
rtg13_nucleus(u) = u[Rtg13I_n] + u[Rtg13A_n]
rtg13_cytosol(u) = u[Rtg13I_c] + u[Rtg13A_c]
rtg3_nucleus(u) = rtg13_nucleus(u) + u[Rtg3A_n] + u[Rtg3I_n]
rtg3_cytosol(u) = rtg13_cytosol(u) + u[Rtg3A_c] + u[Rtg3I_c]
rtg1_nucleus(u) = rtg13_nucleus(u) + u[Rtg1_n]
rtg1_cytosol(u) = rtg13_cytosol(u) + u[Rtg1_c]

# Mapping from symbols to MTK parameters
const symbol2mtk = Dict(Symbol(k) => k for k in (ΣRtg1, ΣRtg2, ΣRtg3, ΣMks, ΣBmh, ksV, ksD, n_S, k2I, mul_S, k2M, kn2M, kBM, knBM, k13I, k13IV, k13ID, k3A_c, k3I_c, k3I_n, k13_c, kn13_c, k13_n, kn13_n, k1in, k1out, k3inA, k3outA, k3inI, k3outI))

# Utility functions

"""Loads a CSV file and returns a RowTable (An tuple of NamedTuple for each row)"""
load_data(filename) = Tables.rowtable(dropmissing(CSV.read(filename, DataFrame), disallowmissing=true))

"""Load boolean conditions for model validation."""
load_conditions(filename="boolean_table_RTG13.csv") = load_data(joinpath(@__DIR__, "data", filename))

"""Load model parameters"""
function load_parameters(filename="solution_rtgMTK_optim.csv")
    paramsets = load_data(joinpath(@__DIR__, "data", filename))
    map(paramsets) do p
        Dict(symbol2mtk[k] => i for (k, i) in pairs(p))
    end
end

ext: jl


text:
#! /nfs/zfs-student-3/users/2013/mngomane/.brew/bin/julia

include("tests/shell.jl")

max_stars_repo_name: mngomane/projects

max_issues_repo_name: mngomane/projects

max_forks_repo_name: mngomane/projects


text:
using RR; using Cxx
using RR: current_session
cxx"""
#include <cxxabi.h>
#include "llvm/MC/SubtargetFeature.h"
#include "llvm/MC/MCObjectFileInfo.h"
#include "llvm/MC/MCContext.h"
#include "llvm/MC/MCDisassembler.h"
#include "llvm/MC/MCInst.h"
#include "llvm/MC/MCInstPrinter.h"
#include "llvm/MC/MCAsmInfo.h"
#include "llvm/Support/TargetSelect.h"
#include "llvm/Support/TargetRegistry.h"
#include "llvm/MC/MCInstrInfo.h"
using namespace llvm;
"""

include(Pkg.dir("RR","scratch","disassembler.jl"))

using DWARF.CallFrameInfo
using Gallium.Unwinder: find_fde
using ObjFileBase: handle
using Gallium: Location
using AbstractTrees

import Gallium.GlibcDyldModules: compute_entry_ptr

current_vm() = current_vm(timeline)

#=
stack_remap = compute_remap()

current_vm(timeline) = Gallium.TransparentRemap(current_task(current_session(timeline)), stack_remap::Vector{Gallium.Remap})
current_vm() = current_vm(timeline)
=#

function get_insts(session, stack)
    stack = isa(stack, Gallium.NativeStack) ? stack.stack[end] : stack
    base, mod = Gallium.find_module(session, modules, UInt(stack.ip))
    modrel = UInt(UInt(stack.ip)-base)
    if isnull(mod.xpdata)
        loc, fde = find_fde(mod, modrel)
        seekloc = loc
        cie = realize_cie(fde)
        nbytes = UInt(CallFrameInfo.fde_range(fde, cie))
    else
        entry = Gallium.Unwinder.find_seh_entry(mod, modrel)
        loc = entry.start
        # Need to translate from virtual to file addresses. Hardcode 0xa00 for
        # now.
        seekloc = loc - 0xa00
        nbytes = entry.stop - entry.start
    end
    if ObjFileBase.isrelocatable(handle(mod))
        # For JIT frames, base is the start of .text, so we need to add that
        # offset back
        text = first(filter(x->sectionname(x)==
            ObjFileBase.mangle_sname(handle(mod),"text"),Sections(handle(mod))))
        seekloc += sectionoffset(text)
    end
    seek(handle(mod), seekloc)
    insts = read(handle(mod), UInt8, nbytes)
    base, loc, insts
end

function Gallium.continue!(session::RR.ReplayTimeline)
    RR.single_step!(timeline)
    while true
        res = RR.continue!(timeline)
        stop_task = icxx"$res.break_status.task;"
        (focus_tid == 0 || (stop_task != 0 && focus_tid == icxx"$stop_task->rec_tid;")) && break
    end
end

global focus_tid = 0

RunDebugREPL() = ASTInterpreter.RunDebugREPL(compute_stack(modules))


text:
include("main.jl")

REWRITE_DATA = true

BUSES, LINES, GENERATORS, rPTDF = load_case_data()
ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity()

# Test dependency on variance estimation confidence
casename = "sample_number_variation"
case_path = "results/$(casename)"

if REWRITE_DATA && isdir(case_path)
    rm(case_path, recursive=true)
end

if !(isdir(case_path)) || REWRITE_DATA
    # Only run if there are no saved results
    mkdir(case_path)

    testvalues = [20, 50, 100, 500, 1000, 5000]
    # var_conf = [0.1, 0.05, 0.01, 0.005, 0.001]
    var_conf = [0.05]

    sample_size = 750
    # Create global set of true observations to make tests comparable
    OBSERVATIONS = rand(ERROR_DIST_TRUE, sample_size)

    cctest_dict = DataFrame(xi= Any[], eta=Any[], v_violation=Any[], f_violation=Any[], median_costs=Any[], mean_costs=Any[], min_costs=Any[], max_costs=Any[], quantile10_cost=Any[], quantile90_cost=Any[], median_delta=Any[], mean_delta=Any[], v_violation_sample=Any[])
    ddtest_dict = DataFrame(xi = Any[], eta=Any[], v_violation=Any[], f_violation=Any[], median_costs=Any[], mean_costs=Any[], min_costs=Any[], max_costs=Any[], quantile10_cost=Any[], quantile90_cost=Any[], median_delta=Any[], mean_delta=Any[], v_violation_sample=Any[])

    solvetimes_dr = []

    for v in testvalues
    for vc in var_conf
        info("Runnig dr for N = $v and ξ=$vc")

        ξ = vc
        ERROR_VARIANCES_TRUE, ERROR_DIST_TRUE, ERROR_VARIANCES_WC, ERROR_VARIANCES_SAMPLE = init_stochasticity(sample_N=v)

        objective, bus_results, line_results, solvetimecc = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_SAMPLE)
        cctest = test_results(bus_results, line_results, objective; samples = sample_size, message="Standard CC", print_on=false)

        ddobjective, ddbus_results, ddline_results, solvetimedr = Solve_DD_DR_CCLinDist(BUSES, LINES, GENERATORS, ERROR_VARIANCES_WC)
        ddtest = test_results(ddbus_results, ddline_results, ddobjective; samples = sample_size, message="DD CC", print_on=false)

        push!(cctest_dict, vcat([vc, v], cctest)) # Redundant Data Saving for easier plotting
        push!(ddtest_dict, vcat([vc, v], ddtest))
    end
    end

    writetable("$(case_path)/cctest_dict.csv", cctest_dict)
    writetable("$(case_path)/ddtest_dict.csv", ddtest_dict)

else
    # If there is saved data load it
    info("Data already exists at $case_path")
    cctest_dict = readtable("$(case_path)/cctest_dict.csv")
    ddtest_dict = readtable("$(case_path)/ddtest_dict.csv")
end


text:
#using Revise
using JuMP, EAGO

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 1000,
                                                    "iteration_limit" => 100000,

                                                    "cp_depth" => -1,
                                                    "cp_repetitions" => -1,
                                                    "cp_forward_reverse_limit" => 2,

                                                    "cut_min_iterations" => 2,
                                                    "cut_max_iterations" => 2,
                                                    "objective_cut_on" => true,
                                                    "subgrad_tighten" => true,

                                                    "obbt_depth" => 4,
                                                    "obbt_repetitions" => 4,
                                                    "obbt_aggressive_on" => true,

                                                    "upper_bounding_depth" => 4,

                                                    "fbbt_lp_depth" => 100000,
                                                    "fbbt_lp_repetitions" => 3))

# ----- Variables ----- #
x_Idx = Any[1, 2, 3, 4, 5, 6]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[5], 0.0)
JuMP.set_lower_bound(x[1], 0.0)
JuMP.set_lower_bound(x[4], 0.0)
JuMP.set_lower_bound(x[2], 0.0)
JuMP.set_lower_bound(x[6], 0.0)
JuMP.set_lower_bound(x[3], 0.0)
JuMP.set_upper_bound(x[1], 1.0)
JuMP.set_upper_bound(x[2], 1.0)
JuMP.set_upper_bound(x[3], 1.0)
JuMP.set_upper_bound(x[4], 1.0)
JuMP.set_upper_bound(x[5], 1.0)
JuMP.set_upper_bound(x[6], 1.0)

# ----- Objective ----- #
@NLobjective(m, Min, -exp(-(10* ((-0.1312)+x[1])^2+0.05* ((-0.1696)+x[2])^2+17* ((-0.5569)+x[3])^2+3.5* ((-0.0124)+x[4])^2+1.7* ((-0.8283)+x[5])^2+8* ((-0.5886)+x[6])^2))+1.2*exp(-(0.05* ((-0.2329)+x[1])^2+10* ((-0.4135)+x[2])^2+17* ((-0.8307)+x[3])^2+0.1* ((-0.3736)+x[4])^2+8* ((-0.1004)+x[5])^2+14* ((-0.9991)+x[6])^2))+3*exp(-(3* ((-0.2348)+x[1])^2+3.5* ((-0.1451)+x[2])^2+1.7* ((-0.3522)+x[3])^2+10* ((-0.2883)+x[4])^2+17* ((-0.3047)+x[5])^2+8* ((-0.665)+x[6])^2))+3.2*exp(-(17* ((-0.4047)+x[1])^2+8* ((-0.8828)+x[2])^2+0.05* ((-0.8732)+x[3])^2+10* ((-0.5743)+x[4])^2+0.1* ((-0.1091)+x[5])^2+14* ((-0.0381)+x[6])^2)))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
println("run time: $run_time")


text:
import AverageShiftedHistograms: ash, xy, pdf, Kernels

## ---------- ENERGY MAPS ----------
"""
    EnergyMap

Energy map for Local Discriminant Basis. Current available types are:
- [`TimeFrequency`](@ref)
- [`ProbabilityDensity`](@ref)
"""
abstract type EnergyMap end

@doc raw"""
    TimeFrequency <: EnergyMap

An energy map based on time frequencies, a measure based on the differences of derived
quantities from projection ``Z_i``, such as mean class energies or cumulants.

**See also:** [`EnergyMap`](@ref), [`ProbabilityDensity`](@ref), [`Signatures`](@ref)
"""
struct TimeFrequency <: EnergyMap end

@doc raw"""
    ProbabilityDensity <: EnergyMap

An energy map based on probability density, a measure based on the differences among the
pdfs of ``Z_i``. Since we do not know the true density functions of the coefficients, the
PDFs are estimated using the Average Shifted Histogram (ASH).

**See also:** [`EnergyMap`](@ref), [`TimeFrequency`](@ref), [`Signatures`](@ref)
"""
struct ProbabilityDensity <: EnergyMap end

@doc raw"""
    Signatures <: EnergyMap

An energy map based on signatures, a measure that uses the Earth Mover's Distance (EMD) to
compute the discriminating measure of a coordinate. Signatures provide us with a fully
data-driven representation, which can be efficiently used with EMD. This representation is
more efficient than a histogram and is able to represent complex data structure with fewer
samples.

Here, a signature for the coefficients in the ``j``-th level, ``k``-th node, ``l``-th index
of class ``c`` is defined as

``s_{j,k,l}^{(c)} = \{(\alpha_{i;j,k,l}^{(c)}, w_{i;j,k,l}^{(c)})\}_{i=1}^{N_c}``

where ``\alpha_{i;j,k,l}^{(c)}`` and ``w_{i;j,k,l}^{(c)}`` are the expansion coefficients
and weights at location ``(j,k,l)`` for signal ``i`` of class ``c`` respectively. Currently,
the two valid types of weights are `:equal` and `:pdf`.

# Parameters
- `weight::Symbol`: (Default: `:equal`) Type of weight to be used to compute ``w_{i;j,k,l}^{(c)}``. Available
    methods are `:equal` and `:pdf`.

**See also:** [`EnergyMap`](@ref), [`TimeFrequency`](@ref), [`ProbabilityDensity`](@ref)
"""
struct Signatures <: EnergyMap
    weight::Symbol
    # Constructor
    Signatures(weight=:equal) = weight ∈ [:equal, :pdf] ? new(weight) :
        throw(ValueError("Invalid weight type. Valid weight types are :equal and :pdf."))
end

@doc raw"""
    energy_map(Xw, y, method)

Computes the energy map based on the decomposed input signals `Xw` and corresponding labels
`y`.

# Arguments
- `Xw::AbstractArray{S} where S<:AbstractFloat`: Decomposed 1D or 2D signals.
- `y::AbstractVector{T} where T`: Corresponding labels to the signals in `Xw`.
- `method::EnergyMap`: Type of energy map to compute. Supported types are:
    - `TimeFrequency()`
    - `ProbabilityDensity()`
    - `Signatures()`

# Returns
- `Γ`: Computed energy map. Depending on the input `method`, the data structure of `Γ` is:
    - `TimeFrequency()` ``\Longrightarrow`` `Array{S,3}` for 1D signals or `Array{S,4}` for
      2D signals.
    - `ProbabilityDensity()` ``\Longrightarrow`` `Array{S,4}` for 1D signals or `Array{S,5}`
      for 2D signals.
    - `Signatures(weight = :equal)`: ``\Longrightarrow`` Vector{NamedTuple{(:coef, :weight),
      Tuple{Array{S}, S}}}
    - `Signatures(weight = :pdf`: ``\Longrightarrow`` Vector{NamedTuple{(:coef, :weight),
      Tuple{Array{S}, Array{S}}}}

!!! tip
    When setting `method = Signatures()` of any type of weight, the output is a vector of
    named tuples, where each named tuple contains the coefficients and weights of all the
    signals of the same label.

# Examples
```julia
using Wavelets, WaveletsExt

X, y = generateclassdata(ClassData(:tri, 5, 5, 5))
Xw = wpdall(X, wavelet(WT.haar))

energy_map(Xw, y, TimeFrequency())
energy_map(Xw, y, ProbabilityDensity())
energy_map(Xw, y, Signatures())
```

**See also:** [`EnergyMap`](@ref). [`TimeFrequency`](@ref), [`ProbabilityDensity`](@ref)
"""
function energy_map(Xw::AbstractArray{S}, y::AbstractVector{T}, method::TimeFrequency) where
                   {S<:AbstractFloat, T}
    # --- Sanity check ---
    N = ndims(Xw)
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    sz = size(Xw)[1:(N-2)]
    L = size(Xw, N-1)
    Nx = size(Xw, N)
    # parameter checking
    @assert 3 ≤ N ≤ 4
    @assert Nx == Ny
    @assert nc > 1
    @assert 1 ≤ L-1 ≤ maxtransformlevels(min(sz...))

    # --- Construct normalized energy map for each class ---
    Γ = Array{S, N}(undef, (sz...,L,nc))
    map_size = prod([sz...,L])                      # number of elements per class of energy map
    slice_dim = N==3 ? [1] : [1,2]                  # dimension for norm computation of each slice of signal
    for (i,cᵢ) in enumerate(c)
        idx = findall(yᵢ -> yᵢ==cᵢ, y)
        @inbounds x = N==3 ? Xw[:,1,idx] : Xw[:,:,1,idx]      # Set of original signals of class cᵢ
        @inbounds xw = N==3 ? Xw[:,:,idx] : Xw[:,:,:,idx]     # Set of decomposed signals of class cᵢ
        norm_sum = sum(mapslices(xᵢ -> norm(xᵢ,2)^2, x, dims = slice_dim))
        en = sum(xw.^2, dims=N)
        start_index = (i-1)*map_size                # start_index for current class in energy map
        for j in eachindex(en)
            @inbounds Γ[start_index+j] = en[j] / norm_sum
        end
    end
    return Γ
end

function energy_map(Xw::AbstractArray{S}, y::AbstractVector{T},
                    method::ProbabilityDensity) where {S<:AbstractFloat, T}
    # --- Sanity check ---
    N = ndims(Xw)
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    sz = size(Xw)[1:(N-2)]
    L = size(Xw, N-1)
    Nx = size(Xw, N)
    # parameter checking
    @assert 3 ≤ N ≤ 4
    @assert Nx == Ny
    @assert nc > 1
    @assert 1 ≤ L-1 ≤ maxtransformlevels(min(sz...))

    # --- Construct empirical probability density for each coefficent of each class ---
    nbins = ceil(Int, (30*Nx)^(1/5))    # number of bins/histogram
    mbins = ceil(Int, 100/nbins)        # number of histograms M/nbins, M=100 is arbitrary
    pdf_len = (nbins+1)*mbins           # vector length of empirical pdf
    slice_size = prod([sz...,L])        # number of elements per class of energy map slice
    map_size = prod([sz...,L,pdf_len])  # dimension for norm computation of each slice of signal
    Γ = Array{Float64,N+1}(undef, (sz..., L, pdf_len, nc))
    for (i,cᵢ) in enumerate(c)      # iterate over each class
        idx = findall(yᵢ -> yᵢ==cᵢ, y)
        @inbounds xw = N==3 ? Xw[:,:,idx] : Xw[:,:,:,idx]  # wavelet packet for class cᵢ
        for j in 1:slice_size
            @inbounds z = @view Xw[j:slice_size:end]    # Coefficients at index j for all signals
            @inbounds zᵢ = @view xw[j:slice_size:end]   # Coefficients at index j for signals in class cᵢ
            # ASH parameter setup
            σ = std(z)
            s = 0.5
            δ = (maximum(z)-minimum(z)+σ)/((nbins+1)*mbins-1)
            rng = range(minimum(z)-s*σ, step=δ, length=(nbins+1)*mbins)
            # Empirical PDF
            epdf = ash(zᵢ, rng=rng, m=mbins, kernel=Kernels.triangular)
            start_index = (i-1)*map_size    # start_index for current class in energy map
            end_index = i*map_size          # end index for current class in energy map
            @inbounds _, Γ[(start_index+j):slice_size:end_index] = xy(epdf) # copy distribution into Γ
        end
    end
    return Γ
end

function energy_map(Xw::AbstractArray{S}, y::AbstractVector{T}, method::Signatures) where
                   {S<:AbstractFloat, T}
    # --- Sanity check ---
    N = ndims(Xw)
    c = unique(y)       # unique classes
    nc = length(c)      # number of classes
    Ny = length(y)
    sz = size(Xw)[1:(N-2)]
    L = size(Xw, N-1)
    Nx = size(Xw, N)
    # parameter checking
    @assert 3 ≤ N ≤ 4
    @assert Nx == Ny
    @assert nc > 1
    @assert 1 ≤ L-1 ≤ maxtransformlevels(min(sz...))
    @assert method.weight ∈ [:equal, :pdf]

    # --- Form signatures in a structure of a named tuple ---
    Γ = method.weight==:equal ?
        Vector{NamedTuple{(:coef, :weight), Tuple{Array{S}, S}}}(undef, nc) :      # equal weights
        Vector{NamedTuple{(:coef, :weight), Tuple{Array{S}, Array{S}}}}(undef, nc) # pdf-based weights
    slice_size = prod([sz...,L])        # number of elements per class of energy map slice
    for (i, cᵢ) in enumerate(c)
        idx = findall(yᵢ -> yᵢ==cᵢ, y)
        xw = N==3 ? Xw[:,:,idx] : Xw[:,:,:,idx] # wavelet packet for class cᵢ
        Nc = length(idx)                        # number of data in class cᵢ
        if method.weight == :equal
            w = 1/Nc
        else
            nbins = ceil(Int, (30*Nx)^(1/5)) # number of bins/histogram
            mbins = ceil(Int, 100/nbins)     # number of histograms M/nbins, M=100 is arbitrary

            # compute weights
            w = Array{S,N}(undef, (sz...,L,Nc))
            for j in 1:slice_size
                z = @view xw[j:slice_size:end]  # coefficients at j for cᵢ signals
                # ASH parameter setup
                σ = std(z); s = 0.5
                δ = (maximum(z)-minimum(z)+σ)/((nbins+1)*mbins-1)
                rng = range(minimum(z)-s*σ, step=δ, length=(nbins+1)*mbins)
                # Empirical PDF
                epdf = ash(z, rng=rng, m=mbins, kernel=Kernels.triangular)
                for k in 1:Nc
                    start_index = (k-1)*slice_size
                    w[start_index+j] = pdf(epdf, z[k])
                end
            end
        end
        Γ[i] = (coef = xw, weight = w)
    end
    return Γ
end


text:
""" Wiring diagrams as a symmetric monoidal category.

This module provides a high-level categorical interface to wiring diagrams,
building on the low-level imperative interface and the operadic interface. It
also defines data types and functions to represent diagonals, codiagonals,
duals, caps, cups, daggers, and other gadgets in wiring diagrams.
"""
module MonoidalDirectedWiringDiagrams
export Ports, Junction, PortOp, BoxOp, functor, permute,
  implicit_mcopy, implicit_mmerge, junctioned_mcopy, junctioned_mmerge,
  junction_diagram, add_junctions, add_junctions!, rem_junctions, merge_junctions,
  junction_caps, junction_cups, junctioned_dunit, junctioned_dcounit

using AutoHashEquals
using LightGraphs

using ...GAT, ...Theories
import ...Theories: dom, codom, id, compose, ⋅, ∘,
  otimes, ⊗, munit, braid, σ, oplus, ⊕, mzero, swap,
  mcopy, delete, Δ, ◊, mmerge, create, ∇, □, dual, dunit, dcounit, mate, dagger,
  plus, zero, coplus, cozero, meet, join, top, bottom, trace
import ...Syntax: functor, head
using ..DirectedWiringDiagrams
import ..DirectedWiringDiagrams: Box, WiringDiagram, input_ports, output_ports
import ..UndirectedWiringDiagrams: add_junctions!, junction_diagram

# Categorical interface
#######################

# Ports as objects
#-----------------

""" A list of ports.

The objects in categories of wiring diagrams.
"""
@auto_hash_equals struct Ports{Theory,Value}
  ports::Vector{Value}
  Ports{T}(ports::Vector{V}) where {T,V} = new{T,V}(ports)
end
Ports(ports::Vector) = Ports{Any}(ports)

# Iterator interface.
Base.iterate(A::Ports, args...) = iterate(A.ports, args...)
Base.keys(A::Ports) = keys(A.ports)
Base.length(A::Ports) = length(A.ports)
Base.eltype(A::Ports{T,V}) where {T,V} = V

# Indexing interface.
Base.getindex(A::Ports, i::Int) = A.ports[i]
Base.getindex(A::Ports{T}, i::UnitRange) where T = Ports{T}(A.ports[i])
Base.firstindex(A::Ports) = 1
Base.lastindex(A::Ports) = length(A)

Base.cat(A::Ports{T}, B::Ports{T}) where T = Ports{T}([A.ports; B.ports])
Base.reverse(A::Ports{T}) where T = Ports{T}(reverse(A.ports))

Box(value, inputs::Ports, outputs::Ports) =
  Box(value, collect(inputs), collect(outputs))
Box(inputs::Ports, outputs::Ports) = Box(collect(inputs), collect(outputs))

WiringDiagram(value, inputs::Ports{T}, outputs::Ports{T}) where T =
  WiringDiagram{T}(value, collect(inputs), collect(outputs))

WiringDiagram(inputs::Ports{T}, outputs::Ports{T}) where T =
  WiringDiagram{T}(collect(inputs), collect(outputs))

dom(d::WiringDiagram{T}) where T = Ports{T}(input_ports(d))
codom(d::WiringDiagram{T}) where T = Ports{T}(output_ports(d))

# Symmetric monoidal category
#----------------------------

""" Wiring diagrams as a symmetric monoidal category.

Extra structure, such as copying or merging, can be added to wiring diagrams in
different ways, but wiring diagrams always form a symmetric monoidal category in
the same way.
"""
@instance SymmetricMonoidalCategory{Ports, WiringDiagram} begin
  @import dom, codom

  function id(A::Ports)
    f = WiringDiagram(A, A)
    add_wires!(f, ((input_id(f),i) => (output_id(f),i) for i in eachindex(A)))
    return f
  end

  function compose(f::WiringDiagram, g::WiringDiagram; unsubstituted::Bool=false)
    if length(codom(f)) != length(dom(g))
      # Check only that f and g have the same number of ports.
      # The port types will be checked when the wires are added.
      error("Incompatible domains $(codom(f)) and $(dom(g))")
    end
    h = WiringDiagram(dom(f), codom(g))
    fv = add_box!(h, f)
    gv = add_box!(h, g)
    add_wires!(h, ((input_id(h),i) => (fv,i) for i in eachindex(dom(f))))
    add_wires!(h, ((fv,i) => (gv,i) for i in eachindex(codom(f))))
    add_wires!(h, ((gv,i) => (output_id(h),i) for i in eachindex(codom(g))))
    unsubstituted ? h : substitute(h, [fv,gv])
  end

  otimes(A::Ports, B::Ports) = cat(A, B)
  munit(::Type{Ports}) = Ports([])

  function otimes(f::WiringDiagram, g::WiringDiagram; unsubstituted::Bool=false)
    h = WiringDiagram(otimes(dom(f),dom(g)), otimes(codom(f),codom(g)))
    m, n = length(dom(f)), length(codom(f))
    fv = add_box!(h, f)
    gv = add_box!(h, g)
    add_wires!(h, (input_id(h),i) => (fv,i) for i in eachindex(dom(f)))
    add_wires!(h, (input_id(h),i+m) => (gv,i) for i in eachindex(dom(g)))
    add_wires!(h, (fv,i) => (output_id(h),i) for i in eachindex(codom(f)))
    add_wires!(h, (gv,i) => (output_id(h),i+n) for i in eachindex(codom(g)))
    unsubstituted ? h : substitute(h, [fv,gv])
  end

  function braid(A::Ports, B::Ports)
    h = WiringDiagram(otimes(A,B), otimes(B,A))
    m, n = length(A), length(B)
    add_wires!(h, ((input_id(h),i) => (output_id(h),i+n) for i in 1:m))
    add_wires!(h, ((input_id(h),i+m) => (output_id(h),i) for i in 1:n))
    h
  end
end

munit(::Type{Ports{T}}) where T = Ports{T}([])
munit(::Type{Ports{T,V}}) where {T,V} = Ports{T}(V[])

# Unbiased version of braiding (permutation).

function permute(A::Ports, σ::Vector{Int}; inverse::Bool=false)
  @assert length(A) == length(σ)
  B = Ports([ A[σ[i]] for i in eachindex(σ) ])
  if inverse
    f = WiringDiagram(B, A)
    add_wires!(f, ((input_id(f),σ[i]) => (output_id(f),i) for i in eachindex(σ)))
  else
    f = WiringDiagram(A, B)
    add_wires!(f, ((input_id(f),i) => (output_id(f),σ[i]) for i in eachindex(σ)))
  end
  return f
end

# Functors
#---------

""" Apply functor in a category of wiring diagrams.

Defined by compatible mappings of ports and boxes.
"""
function functor(d::WiringDiagram, f_ports, f_box; contravariant::Bool=false, kw...)
  functor_impl = contravariant ? contravariant_functor : covariant_functor
  functor_impl(d, f_ports, f_box; kw...)
end

function covariant_functor(d::WiringDiagram, f_ports, f_box)
  result = WiringDiagram(f_ports(dom(d)), f_ports(codom(d)))
  add_boxes!(result, (f_box(box(d, v)) for v in box_ids(d)))
  add_wires!(result, wires(d))
  result
end

function contravariant_functor(d::WiringDiagram, f_ports, f_box;
                               monoidal_contravariant::Bool=false)
  result = WiringDiagram(f_ports(codom(d)), f_ports(dom(d)))
  add_boxes!(result, (f_box(box(d, v)) for v in box_ids(d)))

  nports = (port::Port) -> length(
    (port.kind == InputPort ? input_ports : output_ports)(d, port.box))
  map_port = (port::Port) -> Port(
    if (port.box == input_id(d)) output_id(d)
    elseif (port.box == output_id(d)) input_id(d)
    else port.box end,
    port.kind == InputPort ? OutputPort : InputPort,
    monoidal_contravariant ? nports(port) - port.port + 1 : port.port
  )
  add_wires!(result, map(wires(d)) do wire
    Wire(wire.value, map_port(wire.target), map_port(wire.source))
  end)
  result
end

# Diagonals and codiagonals
#--------------------------

""" Implicit copy in wiring diagram.

Copies are represented by multiple outgoing wires from a single port and
deletions by no outgoing wires.
"""
function implicit_mcopy(A::Ports, n::Int)
  f = WiringDiagram(A, otimes(repeat([A], n)))
  m = length(A)
  add_wires!(f, ((input_id(f),i) => (output_id(f),i+m*(j-1))
                  for i in 1:m for j in 1:n))
  return f
end

""" Implicit merge in wiring diagram.

Merges are represented by multiple incoming wires into a single port and
creations by no incoming wires.
"""
function implicit_mmerge(A::Ports, n::Int)
  f = WiringDiagram(otimes(repeat([A],n)), A)
  m = length(A)
  add_wires!(f, ((input_id(f),i+m*(j-1)) => (output_id(f),i)
                  for i in 1:m for j in 1:n))
  return f
end

""" Explicit copy in wiring diagram.

Copies and deletions are represented by junctions (boxes of type `Junction`).
"""
junctioned_mcopy(A::Ports, n::Int; kw...) = junction_diagram(A, 1, n; kw...)

""" Explicit merge in wiring diagram.

Merges and creations are represented by junctions (boxes of type `Junction`).
"""
junctioned_mmerge(A::Ports, n::Int; kw...) = junction_diagram(A, n, 1; kw...)

# Implicit diagonals and codiagonals are the default in untyped wiring diagrams.
mcopy(A::Ports{Any}, n::Int) = implicit_mcopy(A, n)
mmerge(A::Ports{Any}, n::Int) = implicit_mmerge(A, n)

# Default implementation for biased (co)diagonals fall backs to unbiased
# (co)diagonals, if they are defined.
mcopy(A::Ports) = mcopy(A, 2)
delete(A::Ports) = mcopy(A, 0)
mmerge(A::Ports) = mmerge(A, 2)
create(A::Ports) = mmerge(A, 0)
plus(A::Ports) = plus(A, 2)
zero(A::Ports) = plus(A, 0)
coplus(A::Ports) = coplus(A, 2)
cozero(A::Ports) = coplus(A, 0)

# Cartesian category
#-------------------

mcopy(A::Ports{MonoidalCategoryWithDiagonals}, n::Int) = implicit_mcopy(A, n)

mcopy(A::Ports{CartesianCategory}, n::Int) = implicit_mcopy(A, n)

# Cocartesian category
#---------------------

mmerge(A::Ports{MonoidalCategoryWithCodiagonals}, n::Int) = implicit_mmerge(A, n)

mmerge(A::Ports{CocartesianCategory}, n::Int) = implicit_mmerge(A, n)

# Biproduct category
#-------------------

# The coherence laws relating diagonal to codiagonal do not hold for general
# bidiagonals, so an explicit representation is needed.
mcopy(A::Ports{MonoidalCategoryWithBidiagonals}, n::Int) = junctioned_mcopy(A, n)
mmerge(A::Ports{MonoidalCategoryWithBidiagonals}, n::Int) = junctioned_mmerge(A, n)

mcopy(A::Ports{BiproductCategory}, n::Int) = implicit_mcopy(A, n)
mmerge(A::Ports{BiproductCategory}, n::Int) = implicit_mmerge(A, n)

# Dagger category
#----------------

dagger(f::WiringDiagram{DaggerSymmetricMonoidalCategory}) =
  functor(f, identity, dagger, contravariant=true)

# Compact closed category
#------------------------

junctioned_dunit(A::Ports) = junction_caps(A, otimes(dual(A),A))
junctioned_dcounit(A::Ports) = junction_cups(A, otimes(A,dual(A)))

dual(A::Ports{CompactClosedCategory}) = dual_ports(A)
dunit(A::Ports{CompactClosedCategory}) = junctioned_dunit(A)
dcounit(A::Ports{CompactClosedCategory}) = junctioned_dcounit(A)
mate(f::WiringDiagram{CompactClosedCategory}) =
  functor(f, dual, mate, contravariant=true, monoidal_contravariant=true)

dual(A::Ports{DaggerCompactCategory}) = dual_ports(A)
dunit(A::Ports{DaggerCompactCategory}) = junctioned_dunit(A)
dcounit(A::Ports{DaggerCompactCategory}) = junctioned_dcounit(A)
dagger(f::WiringDiagram{DaggerCompactCategory}) =
  functor(f, identity, dagger, contravariant=true)
mate(f::WiringDiagram{DaggerCompactCategory}) =
  functor(f, dual, mate, contravariant=true, monoidal_contravariant=true)

# Traced monoidal category
#-------------------------

""" Trace (feedback loop) in a wiring diagram.
"""
function trace(X::Ports, f::WiringDiagram; unsubstituted::Bool=false)
  n = length(X)
  @assert dom(f)[1:n] == X && codom(f)[1:n] == X
  A, B = dom(f)[n+1:end], codom(f)[n+1:end]
  h = WiringDiagram(A, B)
  fv = add_box!(h, f)
  add_wires!(h, (fv,i) => (fv,i) for i in 1:n)
  add_wires!(h, (input_id(h),i) => (fv,i+n) for i in eachindex(A))
  add_wires!(h, (fv,i+n) => (output_id(h),i) for i in eachindex(B))
  unsubstituted ? h : substitute(h)
end

const TracedMon = TracedMonoidalCategory

trace(X::Ports{TracedMon}, A::Ports{TracedMon},
      B::Ports{TracedMon}, f::WiringDiagram{TracedMon}) = trace(X, f)

# Bicategory of relations
#------------------------

const BiRel = BicategoryRelations

mcopy(A::Ports{BiRel}, n::Int) = junctioned_mcopy(A, n)
mmerge(A::Ports{BiRel}, n::Int) = junctioned_mmerge(A, n)

dunit(A::Ports{BiRel}) = junction_caps(A)
dcounit(A::Ports{BiRel}) = junction_cups(A)

dagger(f::WiringDiagram{BiRel}) =
  functor(f, identity, dagger, contravariant=true)

meet(f::WiringDiagram{BiRel}, g::WiringDiagram{BiRel}) =
  compose(mcopy(dom(f)), otimes(f,g), mmerge(codom(f)))
top(A::Ports{BiRel}, B::Ports{BiRel}) = compose(delete(A), create(B))

# Abelian bicategory of relations
#--------------------------------

const AbBiRel = AbelianBicategoryRelations

# Additive notation. FIXME: Use @instance.
oplus(f::WiringDiagram{AbBiRel}, g::WiringDiagram{AbBiRel}) = otimes(f,g)
⊕(f::WiringDiagram{AbBiRel}, g::WiringDiagram{AbBiRel}) = oplus(f,g)
mzero(::Type{T}) where T <: Ports{AbBiRel} = munit(T)
swap(A::Ports{AbBiRel}, B::Ports{AbBiRel}) = braid(A, B)

mcopy(A::Ports{AbBiRel}, n::Int) = junctioned_mcopy(A, n; op=:times)
mmerge(A::Ports{AbBiRel}, n::Int) = junctioned_mmerge(A, n; op=:times)

dunit(A::Ports{AbBiRel}) = junction_caps(A; op=:times)
dcounit(A::Ports{AbBiRel}) = junction_cups(A; op=:times)

plus(A::Ports{AbBiRel}, n::Int) = junctioned_mmerge(A, n, op=:plus)
coplus(A::Ports{AbBiRel}, n::Int) = junctioned_mcopy(A, n, op=:plus)

dagger(f::WiringDiagram{AbBiRel}) =
  functor(f, identity, dagger, contravariant=true)

meet(f::WiringDiagram{AbBiRel}, g::WiringDiagram{AbBiRel}) =
  compose(mcopy(dom(f)), oplus(f,g), mmerge(codom(f)))
join(f::WiringDiagram{AbBiRel}, g::WiringDiagram{AbBiRel}) =
  compose(coplus(dom(f)), oplus(f,g), plus(codom(f)))
top(A::Ports{AbBiRel}, B::Ports{AbBiRel}) = compose(delete(A), create(B))
bottom(A::Ports{AbBiRel}, B::Ports{AbBiRel}) = compose(cozero(A), zero(B))

# Junctions
###########

""" Junction node in a wiring diagram.

Junction nodes are used to explicitly represent copies, merges, deletions,
creations, caps, and cups.
"""
struct Junction{Op,Value} <: AbstractBox
  value::Value
  input_ports::Vector
  output_ports::Vector
  Junction{Op}(value::Value, inputs::Vector, outputs::Vector) where {Op,Value} =
    new{Op,Value}(value, inputs, outputs)
end

Junction(args...) = Junction{nothing}(args...)
Junction{Op}(value, ninputs::Int, noutputs::Int) where Op =
  Junction{Op}(value, repeat([value], ninputs), repeat([value], noutputs))

head(junction::Junction{Op}) where Op = Op

Base.:(==)(j1::Junction, j2::Junction) =
  head(j1) == head(j2) && j1.value == j2.value &&
  input_ports(j1) == input_ports(j2) && output_ports(j1) == output_ports(j2)

""" Wiring diagram with a junction node for each of the given ports.
"""
junction_diagram(A::Ports, nin::Int, nout::Int; op=nothing) =
  junction_diagram(Junction{op}, A, nin, nout)

function junction_diagram(make_junction, A::Ports, nin::Int, nout::Int)
  f = WiringDiagram(otimes(repeat([A], nin)), otimes(repeat([A], nout)))
  m = length(A)
  for (i, value) in enumerate(A)
    v = add_box!(f, make_junction(
      value, repeat([value], nin), repeat([value], nout)))
    add_wires!(f, ((input_id(f),i+m*(j-1)) => (v,j) for j in 1:nin))
    add_wires!(f, ((v,j) => (output_id(f),i+m*(j-1)) for j in 1:nout))
  end
  return f
end

""" Wiring diagram of nested cups made out of junction nodes.
"""
junction_cups(A::Ports; kw...) = junction_cups(A, cat(A,reverse(A)); kw...)
junction_cups(A::Ports, inputs::Ports; op=nothing) =
  junction_cups(Junction{op}, A, inputs)

function junction_cups(make_junction, A::Ports, inputs::Ports)
  @assert length(inputs) == 2*length(A)
  f = WiringDiagram(inputs, munit(typeof(inputs)))
  m, ports = length(A), collect(inputs)
  for (i, value) in enumerate(A)
    j1, j2 = i, 2m-i+1 # Outer cups to inner cups
    v = add_box!(f, make_junction(value, ports[[j1,j2]], empty(ports)))
    add_wires!(f, [(input_id(f),j1) => (v,1), (input_id(f),j2) => (v,2)])
  end
  return f
end

""" Wiring diagram of nested caps made out of junction nodes.
"""
junction_caps(A::Ports; kw...) = junction_caps(A, cat(reverse(A),A); kw...)
junction_caps(A::Ports, outputs::Ports; op=nothing) =
  junction_caps(Junction{op}, A, outputs)

function junction_caps(make_junction, A::Ports, outputs::Ports)
  @assert length(outputs) == 2*length(A)
  f = WiringDiagram(munit(typeof(outputs)), outputs)
  m, ports = length(A), collect(outputs)
  for (i, value) in enumerate(A)
    j1, j2 = m-i+1, m+i # Inner caps to outer caps
    v = add_box!(f, make_junction(value, empty(ports), ports[[j1,j2]]))
    add_wires!(f, [(v,1) => (output_id(f),j1), (v,2) => (output_id(f),j2)])
  end
  return f
end

""" Add junction nodes to wiring diagram.

Transforms from the implicit to the explicit representation of diagonals and
codiagonals. This operation is inverse to `rem_junctions`.
"""
function add_junctions(d::WiringDiagram)
  add_junctions!(copy(d))
end
function add_junctions!(d::WiringDiagram)
  add_output_junctions!(d, input_id(d))
  add_input_junctions!(d, output_id(d))
  for v in box_ids(d)
    add_input_junctions!(d, v)
    add_output_junctions!(d, v)
  end
  return d
end

function add_input_junctions!(d::WiringDiagram, v::Int)
  for (port, port_value) in enumerate(input_ports(d, v))
    wires = in_wires(d, v, port)
    nwires = length(wires)
    if nwires != 1
      rem_wires!(d, wires)
      jv = add_box!(d, Junction(port_value, nwires, 1))
      sources = sort!([ wire.source for wire in wires ])
      add_wire!(d, Port(jv, OutputPort, 1) => Port(v, InputPort, port))
      add_wires!(d, [ source => Port(jv, InputPort, i)
                      for (i, source) in enumerate(sources) ])
    end
  end
end

function add_output_junctions!(d::WiringDiagram, v::Int)
  for (port, port_value) in enumerate(output_ports(d, v))
    # TODO: Sort out wires.
    wires = out_wires(d, v, port)
    nwires = length(wires)
    if nwires != 1
      rem_wires!(d, wires)
      jv = add_box!(d, Junction(port_value, 1, nwires))
      targets = sort!([ wire.target for wire in wires ])
      add_wire!(d, Port(v, OutputPort, port) => Port(jv, InputPort, 1))
      add_wires!(d, [ Port(jv, OutputPort, i) => target
                      for (i, target) in enumerate(targets) ])
    end
  end
end

""" Remove junction nodes from wiring diagram.

Transforms from the explicit to the implicit representation of diagonals and
codiagonals. This operation is inverse to `add_junctions`.
"""
function rem_junctions(d::WiringDiagram; op=nothing)
  junction_ids = filter(v -> box(d,v) isa Junction{op}, box_ids(d))
  junction_diagrams = map(junction_ids) do v
    junction = box(d,v)::Junction
    inputs, outputs = input_ports(junction), output_ports(junction)
    wiring = WiringDiagram(inputs, outputs)
    add_wires!(wiring, (input_id(wiring), src) => (output_id(wiring), tgt)
                        for src in eachindex(inputs), tgt in eachindex(outputs))
    wiring
  end
  substitute(d, junction_ids, junction_diagrams)
end

""" Merge adjacent junction nodes into single junctions.
"""
function merge_junctions(d::WiringDiagram; op=nothing)
  junction_ids = filter(v -> box(d,v) isa Junction{op}, box_ids(d))
  junction_graph, vmap = induced_subgraph(DiGraph(graph(d)), junction_ids)
  for edge in edges(junction_graph)
    # Only merge junctions with equal values.
    if box(d, vmap[src(edge)]).value != box(d, vmap[dst(edge)]).value
      rem_edge!(junction_graph, edge)
    end
  end
  components = [ [vmap[v] for v in component]
    for component in weakly_connected_components(junction_graph)
    if length(component) > 1 ]
  values = [ box(d, first(component)).value for component in components ]
  encapsulate(d, components;
    discard_boxes=true, values=values, make_box=Junction{op})
end

# Operations on ports and boxes
###############################

""" Port value wrapping another value.

Represents unary operations on ports in wiring diagrams.
"""
struct PortOp{Op}
  value::Any
end

head(::PortOp{Op}) where Op = Op

Base.:(==)(op1::PortOp, op2::PortOp) =
  head(op1) == head(op2) && op1.value == op2.value

""" Box wrapping another box.

Represents unary operations on boxes in wiring diagrams.
"""
struct BoxOp{op} <: AbstractBox
  box::AbstractBox
end

head(::BoxOp{Op}) where Op = Op
input_ports(op::BoxOp) = input_ports(op.box)
output_ports(op::BoxOp) = output_ports(op.box)

Base.:(==)(op1::BoxOp, op2::BoxOp) =
  head(op1) == head(op2) && op1.box == op2.box

# Duals
#------

const DualPort = PortOp{:dual}

dual_port(x) = DualPort(x)
dual_port(dual::DualPort) = dual.value

dual_ports(ports::Vector) = [ dual_port(x) for x in Iterators.reverse(ports) ]
dual_ports(ports::Ports{T}) where T = Ports{T}(dual_ports(collect(ports)))

# Adjoints
#---------

const DaggerBox = BoxOp{:dagger}

input_ports(dagger::DaggerBox) = output_ports(dagger.box)
output_ports(dagger::DaggerBox) = input_ports(dagger.box)

dagger(box::Box) = DaggerBox(box)
dagger(dagger::DaggerBox) = dagger.box
dagger(junction::Junction{Op}) where Op = Junction{Op}(
  junction.value, output_ports(junction), input_ports(junction))

const MateBox = BoxOp{:mate}

input_ports(mate::MateBox) = dual_ports(output_ports(mate.box))
output_ports(mate::MateBox) = dual_ports(input_ports(mate.box))

mate(box::Box) = MateBox(box)
mate(mate::MateBox) = mate.box

# Assume that mates and daggers commute, as in a dagger compact category.
# Normalize to apply mates before daggers.
dagger(mate::MateBox) = DaggerBox(mate)
mate(dagger::DaggerBox) = dagger(mate(dagger.box))

end


text:
######################################################################
# Errors used throughout P2PNet
# -----
# Licensed under Apache License 2.0

export CancellationError
struct CancellationError <: Exception
    msg::String
end
CancellationError() = CancellationError("")

export ConnectionError
struct ConnectionError <: Exception
    msg::String
    cause::Vector{Exception}
end
ConnectionError(msg::AbstractString = "") = ConnectionError(msg, [])
function Base.showerror(io::IO, err::ConnectionError)
    print(io, "Connection error")
    if !isempty(err.msg)
        print(io, ": $(err.msg)")
    end

    if !isempty(err.cause)
        println(io)
        println(io, "$(length(err.cause)) failures: ")

        for sub ∈ err.cause
            showerror(io, sub)
            # TODO: Print stacktrace please
        end
    end
end

export NotFoundError
struct NotFoundError <: Exception
    msg::String
end
NotFoundError() = NotFoundError("")

export NotImplementedError
struct NotImplementedError <: Exception
    msg::String
end
NotImplementedError() = NotImplementedError("")
function Base.show(io::IO, err::NotImplementedError)
    print(io, "Not implemented")
    if !isempty(err.msg)
        print(io, ": ")
        print(io, err.msg)
    end
end

export ParityError
struct ParityError <: Exception
    msg::String
end
ParityError() = ParityError("")

export StateError
struct StateError <: Exception
    msg::String
end
StateError() = StateError("")

export @notimplemented
macro notimplemented()
    :(throw(NotImplementedError()))
end


text:
"""global menu configuration parameters"""
CONFIG = Dict()

"""

    config( <see arguments> )

Keyword-only function to configure global menu parameters

# Arguments
 - `charset::Symbol=:na`: ui characters to use (`:ascii` or `:unicode`); overridden by other arguments
 - `cursor::Char='>'|'→'`: character to use for cursor
 - `up_arrow::Char='^'|'↑'`: character to use for up arrow
 - `down_arrow::Char='v'|'↓'`: character to use for down arrow
 - `checked::String="[X]"|"✓"`: string to use for checked
 - `unchecked::String="[ ]"|"⬚")`: string to use for unchecked
 - `scroll::Symbol=:na`: If `:wrap` then wrap the cursor around top and bottom, if :`nowrap` do not wrap cursor
 - `supress_output::Bool=false`: For testing. If true, menu will not be printed to console.
"""
function config(;charset::Symbol = :na,
                scroll::Symbol = :na,
                cursor::Char = '\0',
                up_arrow::Char = '\0',
                down_arrow::Char = '\0',
                checked::String = "",
                unchecked::String = "",
                supress_output::Union{Void, Bool}=nothing)

    if !(charset in [:na, :ascii, :unicode])
        error("charset should be :ascii or :unicode, recieved $charset")
    end

    if !(scroll in [:na, :wrap, :nowrap])
        error("scroll must be :wrap or :nowrap, recieved $scroll")
    end

    if scroll == :wrap
        CONFIG[:scroll_wrap] = true
    elseif scroll == :nowrap
        CONFIG[:scroll_wrap] = false
    end

    if charset == :ascii
        cursor = '>'
        up_arrow = '^'
        down_arrow = 'v'
        checked = "[X]"
        unchecked = "[ ]"
    elseif charset == :unicode
        cursor = '→'
        up_arrow = '↑'
        down_arrow = '↓'
        checked = "✓"
        unchecked = "⬚"
    end

    if isa(supress_output, Bool)
        CONFIG[:supress_output] = supress_output
    end

    cursor != '\0'      && (CONFIG[:cursor] = cursor)
    up_arrow != '\0'    && (CONFIG[:up_arrow] = up_arrow)
    down_arrow != '\0'  && (CONFIG[:down_arrow] = down_arrow)
    checked != ""       && (CONFIG[:checked] = checked)
    unchecked != ""     && (CONFIG[:unchecked] = unchecked)

    # Don't return anything
    nothing
end

# Set up defaults
config(charset=:ascii, scroll=:nowrap, supress_output=false)


text:
## Copyright (c) 2013 Miguel Bazdresch
##
## This file is distributed under the 2-clause BSD License.

# Gaston test framework
#
# There are two kinds of tests:
#   - Tests that should error -- things that Gaston should not allow, such as
#     building a figure with an invalid plotstyle.
#   - Tests that should succeed -- things that should work without producing
#     an error.

using Gaston
if VERSION >= v"0.7-"
	using Test
else
	using Base.Test
end

@testset "Pass expected" begin
	closeall()
    set(terminal="null")
	# figures
	@test figure() == 1
	@test figure() == 2
	@test figure(4) == 4
	@test closefigure(4) == 4
	@test closefigure() == 2
	@test closeall() == 1
	@test closeall() == 0
	# plots
	@test plot(1:10) == 1
	@test plot(1:10,handle=2) == 2
	@test plot(1:10,handle=4) == 4
	@test closeall() == 3
	@test plot(sin.(-3:0.01:3), legend = "sine", plotstyle = "lines",
		color = "blue", marker = "ecircle", linewidth = 2, pointsize = 1.1,
		title = "test plot 1", xlabel = "x", ylabel = "y",
		box = "inside horizontal left top", axis ="loglog") == 1
	@test plot(1:10,xrange = "[2:3]") == 1
	@test plot(1:10,xrange = "[-1.1:3.4]") == 1
	@test plot(1:10,xrange = "[:3.4]") == 1
	@test plot(1:10,xrange = "[3.4:]") == 1
	@test plot(1:10,xrange = "[3.4:*]") == 1
	@test plot(1:10,xrange = "[*:3.4]") == 1
	@test plot(3,4,plotstyle="points",pointsize=3,xrange="[2.95:3.05]",yrange="[3.95:4.045]") == 1
	@test plot(rand(10).+im.*rand(10)) == 1
	@test plot(3+4im,plotstyle="points",pointsize=3,xrange="[2.95:3.05]",yrange="[3.95:4.045]") == 1
	@test begin
		err = Gaston.ErrorCoords(rand(40))
		plot(1:40,err=err,plotstyle="errorbars");
	end == 1
	@test begin
		err = Gaston.ErrorCoords(rand(40))
		plot(1:40,err=err,plotstyle="errorlines");
	end == 1
	@test begin
		err = Gaston.ErrorCoords(rand(40),rand(40))
		plot(1:40,err=err,plotstyle="errorbars");
	end == 1
	@test begin
		err = Gaston.ErrorCoords(rand(40),rand(40))
		plot(1:40,err=err,plotstyle="errorlines");
	end == 1
	@test begin
		fin = Gaston.FinancialCoords(0.1*rand(10),0.1*rand(10),0.1*rand(10),0.1*rand(10))
		plot(1:10,financial=fin,plotstyle="financebars");
	end == 1
	# histograms
	@test histogram(rand(1000)) == 1
	@test histogram(randn(1000), bins=100, norm=1, color="blue", linewidth=2,
        title="test histogram", xlabel="x", ylabel="y",
        box="inside horizontal left top") == 1
    # imagesc
	z = rand(5,6)
	@test imagesc(z,title="test imagesc 1",xlabel="xx",ylabel="yy") == 1
	@test imagesc(1:6,1:5,z,title="test imagesc 3",xlabel="xx",ylabel="yy") == 1
	# surf
	@test surf(rand(10,10)) == 1
	@test surf(rand(10,10)) == 1
	@test surf(0:9,2:11,rand(10,10)) == 1
    @test surf(0:9,2:11,(x,y)->x*y) == 1
	@test surf(0:9,2:11,(x,y)->x*y,title="test",plotstyle="pm3d") == 1
	# printfigure
	set(outputfile="$(tempdir())/gastontest")
	@test begin
		plot(1:10)
		printfigure()
	end == 1
	@test begin
		plot(1:10,handle=2)
		printfigure(handle=2,term="png")
		closefigure()
	end == 2
	@test begin
		plot(1:10)
		printfigure(term="eps")
	end == 1
	@test printfigure(term="pdf") == 1
	@test begin
		set(print_size="640,480")
		printfigure(term="svg")
	end == 1
    if !is_windows()
        @test printfigure(term="gif") == 1 # This test does not pass in Windows
    end
	# build a multiple-plot figure manually
	@test begin
		ac = Gaston.AxesConf(title="T")
		x1, exp_pdf = Gaston.hist(randn(10000),25)
		exp_pdf .= exp_pdf./(step(x1)*sum(exp_pdf))
		exp_cconf = Gaston.CurveConf(plotstyle="boxes",
									 color="blue",
									 legend="E")
		exp_curve = Gaston.Curve(x1,exp_pdf,exp_cconf)
		x2 = -5:0.05:5
		theo_pdf = @. 1/sqrt(2π)*exp((-x2^2)/2)
		theo_cconf = Gaston.CurveConf(color="black",legend="T")
		theo_curve = Gaston.Curve(x2,theo_pdf,theo_cconf)
		figure(1)
		Gaston.push_figure!(1,ac,exp_curve,theo_curve)
		Gaston.llplot()
	end == nothing
	# set
	@test set(legend="A") == nothing
	@test set(plotstyle="linespoints") == nothing
	@test set(color="red") == nothing
	@test set(marker="ecircle") == nothing
	@test set(linewidth=3) == nothing
	@test set(pointsize=3) == nothing
	@test set(palette="gray") == nothing
	@test set(title="A") == nothing
	@test set(xlabel="A") == nothing
	@test set(ylabel="A") == nothing
	@test set(zlabel="A") == nothing
	@test set(fill="solid") == nothing
	@test set(grid="on") == nothing
    if !is_windows()
        @test set(terminal="x11") == nothing # This test does not pass in Windows
    end
	@test set(outputfile="A") == nothing
	@test set(print_color="red") == nothing
	@test set(print_fontface="A") == nothing
	@test set(print_fontscale=1) == nothing
	@test set(print_linewidth=3) == nothing
	@test set(print_size="10,10") == nothing
	closeall()
end

@testset "Failure expected" begin
	closeall()
    set(terminal="null")
	# figure-related
	@test_throws ErrorException figure("invalid")
	@test_throws ErrorException figure(1.0)
    @test_throws ErrorException figure(1:2)
	@test_throws ErrorException closefigure(-1)
	@test_throws ErrorException closefigure("invalid")
	@test_throws ErrorException closefigure(1.0)
	@test_throws ErrorException closefigure(1:2)
	# plot
	@test_throws AssertionError plot(0:10,0:11)
	@test_throws AssertionError plot(0:10,legend=0)
	@test_throws AssertionError plot(0:10,plotstyle="invalid")
	@test_throws AssertionError plot(0:10,marker="invalid")
	@test_throws AssertionError plot(0:10,marker=0)
	@test_throws AssertionError plot(0:10,linewidth="b")
	@test_throws AssertionError plot(0:10,linewidth=im)
	@test_throws AssertionError plot(0:10,pointsize="b")
	@test_throws AssertionError plot(0:10,pointsize=im)
	@test_throws AssertionError plot(0:10,title=0)
	@test_throws AssertionError plot(0:10,xlabel=0)
	@test_throws AssertionError plot(0:10,ylabel=0)
	@test_throws AssertionError plot(0:10,axis="invalid")
	@test_throws AssertionError plot(1:10,xrange = "2:3")
	@test_throws AssertionError plot(1:10,yrange = "ab")
	f = Gaston.FinancialCoords([1,2],[1,2],[1,2],[1,2])
	@test_throws AssertionError plot(1:10,financial=f)
	er = Gaston.ErrorCoords([0.1,0.1])
	@test_throws AssertionError plot(1:10,err=er)
	# plot!
	plot(1:10)
	@test_throws AssertionError plot!(0:10,legend=0)
	@test_throws MethodError plot!(0:10,axis="loglog")
	# imagesc
	z = rand(5,6)
	@test_throws AssertionError imagesc(1:5,1:7,z)
	# histogram
	@test_throws MethodError histogram(0:10+im*0:10)
	# set
	@test_throws AssertionError set(legend=3)
	@test_throws AssertionError set(plotstyle="A")
	@test_throws AssertionError set(color=3)
	@test_throws AssertionError set(marker="xyz")
	@test_throws AssertionError set(linewidth="A")
	@test_throws AssertionError set(pointsize="A")
	@test_throws AssertionError set(title=3)
	@test_throws AssertionError set(xlabel=3)
	@test_throws AssertionError set(ylabel=3)
	@test_throws AssertionError set(zlabel=3)
	@test_throws AssertionError set(fill="red")
	@test_throws AssertionError set(grid="xyz")
	@test_throws AssertionError set(terminal="x12")
	@test_throws AssertionError set(outputfile=3)
	@test_throws AssertionError set(print_color=3)
	@test_throws AssertionError set(print_fontface=3)
	@test_throws AssertionError set(print_fontscale="1")
	@test_throws AssertionError set(print_linewidth="3")
	@test_throws AssertionError set(print_size=10)
	closeall()
end

@testset "linestyle tests" begin
    closeall()
    set(terminal="null")
    @test plot(1:10) == 1 # solid
    @test plot(1:10, linestyle="") == 1 # solid
    @test plot(1:10, linestyle="-") == 1 # dashed
    @test plot(1:10, linestyle=".") == 1 # dotted
    @test plot(1:10, linestyle="_") == 1 # em-dashed
    @test plot(1:10, linestyle="- ._. ") == 1 # complex pattern
    @test_throws AssertionError plot(1:10, linestyle=" ") # only spaces not allowed
    @test_throws AssertionError plot(1:10, linestyle="-=") # = is not allowed in pattern
    closeall()
end


text:
#!/usr/bin/env julia
#=
This is to test Beagle usage.
When imputing, Beagle can either take two separate files, one as reference,
and the other as genotype file to be imputed. Or, these two files can be
merged as one as the genotype file.

Here I used the 4,796 ID who were genotyped with 17k chips.  Randomly sample
100 ID to be masked, using the rest as reference.
=#

using DelimitedFiles, Plots, Statistics
dat = readdlm("num.txt")
dat = reshape(dat, 4, :)'
n = length(dat[:, 1])
mtwo = mean(dat[:, 1])
mone=mean(dat[:,3])

plot(dat[:, 1], label="Two", dpi=300)
plot!([1, n], [mtwo, mtwo], label="Mean two", dpi=300)
plot!(dat[:, 3], label="One", dpi=300)
plot!([1, n], [mone, mone], label="Mean one", dpi=300)

savefig("one-vs-two.png")

max_stars_repo_name: xijiang/Refactered-NSG

max_issues_repo_name: xijiang/Refactered-NSG

max_forks_repo_name: xijiang/Refactered-NSG


text:
"""
    Ziegler-Biersack-Littmark (ZBL) Potential

See https://docs.lammps.org/pair_zbl.html
"""
mutable struct ZBL <: Potential
    ε0::Float64      # Electrical permittivity of vacuum.
    e::Float64       # Electron charge.
    zi::Float64      # Atomic number of the atom i.
    zj::Float64      # Atomic number of the atom j.
    a::Float64
    rcutfac::Float64 # Rcut
end

"""
    ZBL(params::Dict)

Creates a ZBL potential.
"""
function ZBL(params::Dict)
    # Creates the ZBL model
    ε0 = params["ε0"]
    e = params["e"]
    zi = params["zi"]
    zj = params["zj"]
    a = 0.46850 / (zi^0.23 + zj^0.23)
    rcutfac = params["rcutfac"]
    return ZBL(ε0, e, zi, zj, a, rcutfac)
end

"""
    potential_energy(p::ZBL, r::Position, args...)

Calculates the potential energy of ZBL with a switching function.
"""
function potential_energy(p::ZBL, r::Position, args...)
    return  zbl(p, norm(r)) + S(p, norm(r))
end

"""
    zbl(p::ZBL, r::Float64)

ZBL function without the switching function S.
"""
function zbl(p::ZBL, r::Float64)
    return 1 / (4.0 * π * p.ε0) * p.zi * p.zj * p.e^2 / r * ϕ(p, r / p.a)
end

"""
    S(r)

Switching function that ramps the energy, force, and curvature smoothly to zero
between an inner and outer cutoff. Here, the inner and outer cutoff are the same
for all pairs of atom types.
"""
function S(p::ZBL, r::Float64)
    return r <= p.rcutfac ? 0.5 * (cos(pi * r / p.rcutfac) + 1.0) : 0
end

# TODO: check this function: "exp" or "e"?
"""
    ϕ(p::ZBL, x::Float64)

Auxiliary function, necessary to calculate ZBL.
"""
function ϕ(p::ZBL, x::Float64)
    return   0.18175 * exp(-3.1998*x)  + 0.50986 * exp(-0.94229*x) +
             0.28022 * exp(-0.40290*x) + 0.2817  * exp(-0.20162*x)
end


text:
using Printf

function output_result(stepnum, Qbase, cellxmax, cellymax, specific_heat_ratio, out_file_front, out_ext, out_dir, Rd, nval)

    stepnum = string(stepnum)
    while length(stepnum) < 6
        stepnum = "0"*stepnum
    end

    fff = out_dir*"/"*out_file_front*stepnum*out_ext
    open(fff,"w") do f
        write(f,"result:rho[kg/m^3], u[m/s], v[m/s], p[Pa], T[K]\n")
        for i in 2:cellxmax-1
            for j in 2:cellymax-1
                for l in 1:nval
                    a = @sprintf("%8.8e", Qbase[i,j,l])
                    write(f, a*" ")
                end
                T = Qbase[i,j,4]/(Qbase[i,j,1]*Rd)
                a = @sprintf("%8.8e", T)
                write(f, a*"\n")
            end
        end
    end
    println("\nwrite "*fff)
end

function reset_write(fwrite)
    open( fwrite, "w" ) do f
    end
end

function output_physicaltime(fwrite, t, dt)
    a = string(t)
    b = @sprintf("%8.8f", t*dt)
    open(fwrite, "a") do f
        write(f, "\n")
        write(f, "step = "*a)
        write(f, "     ")
        write(f, "time[s] = "*b)
        write(f, "\n")
        write(f, "\n")
    end
end

function output_innertime(fwrite, tau, norm2, nval)
    a = string(tau)
    open(fwrite, "a") do f
        write(f, a)
        write(f, " ")
        for l in 1:nval
            a = @sprintf("%8.8e", norm2[l])
            write(f, a)
            write(f, " ")
        end
        write(f, "\n")
    end
end

function output_fin(fwrite, start_t, end_t, nt, dt, in_nt, cellxmax, cellymax)

    etime = end_t - start_t     # 経過時間
    outtime = ["temp"]
    outtime[1] = string(etime)        # 文字列への変換

    et = replace(outtime[1], " milliseconds"=>"")
    et = parse(Float64, et)           # 整数への変換(ミリ秒)

    et_h = floor(et / 10^3 / 60 /60)  # n時間の計算
    et = et - et_h * (10^3 * 60 * 60) # et - n hour
    et_m = floor(et / 10^3 / 60 )     # n分の計算
    et = et - et_m * (10^3 * 60)      # et - n min
    et_s = floor(et / 10^3)           # n秒の計算
    et = et - et_s * (10^3)           # et -n sec

    et_h = @sprintf( "%02d", et_h)    # 文字列への変換
    et_m = @sprintf( "%02d", et_m)    # 文字列への変換
    et_s = @sprintf( "%02d", et_s)    # 文字列への変換
    et = et_h*":"*et_m*":"*et_s       # 結合

    numcell = string(cellxmax * cellymax)
    dt      = string(dt)
    nt      = string(nt)
    in_nt   = string(in_nt)

    open(fwrite, "a") do f
        write(f, "\n")
        write(f, "\n")
        write(f, "---------------------------- \n")
        write(f, "  fin program  ")
        write(f, "\n")
        write(f, "  number of cell            = " * numcell)
        write(f, "\n")
        write(f, "  dt                        = " * dt)
        write(f, "\n")
        write(f, "  number of time iteration  = " * nt)
        write(f, "\n")
        write(f, "  number of inner iteration = " * in_nt)
        write(f, "\n")
        write(f, "  elapse time (hh:mm:ss)    = " * et)
        write(f, "\n")
        write(f, "---------------------------- \n")
    end
end

max_stars_repo_name: hide-dog/test

max_issues_repo_name: hide-dog/test

max_forks_repo_name: hide-dog/test


text:
using SegmentIntersections, Test

include("point_tests.jl")
include("segment_tests.jl")
include("brute_tests.jl")
include("bentley_ottmann/event_tests.jl")
include("bentley_ottmann/queue_tests.jl")
include("bentley_ottmann/status_tests.jl")
include("bentley_ottmann/results_tests.jl")


text:
using ViewerGL
GL = ViewerGL

function allexamples()
    include("./Convex.jl")
    include("./Convex2d.jl")
    include("./Grid.jl")
    include("./Mesh.jl")
    include("./numbering.jl")
    include("./Points.jl")
    include("./Polygon.jl")
    include("./Polyhedron.jl")
    include("./Text.jl")
end

allexamples()


text:
module Raman
import Luna.PhysData: c, ε_0, ħ, k_B, roomtemp, amg
import Luna.PhysData: raman_parameters
import Luna.Maths: planck_taper

abstract type AbstractRamanResponse end

# make Raman responses broadcast like a scalar
Broadcast.broadcastable(R::AbstractRamanResponse) = Ref(R)

"""
    hrpre(R::AbstractRamanResponse, t)

Get the pre (without damping) response function at time `t`.

"""
function hrpre end

"""
    hrdamp(R::AbstractRamanResponse, ρ)

Get the damping (dephasing) constant `τ2` for density `̢ρ`.

"""
function hrdamp end

"""
    (R::AbstractRamanResponse)(t, ρ)

Get the full response function at time `t` and density `̢ρ`.

"""
function (R::AbstractRamanResponse)(t, ρ)
    hrpre(R, t) * exp(-t/hrdamp(R, ρ))
end

struct RamanRespSingleDampedOscillator{Tτ2} <: AbstractRamanResponse
    K::Float64 # overall scale factor
    Ω::Float64 # frequency
    τ2ρ::Tτ2 # coherence dephasing time function density to τ2: ρ -> τ2
end

"""
    RamanRespNormedSingleDampedOscillator(K, Ω, τ2)

Construct a simple normalised single damped oscillator model with scale factor `K`,
angular frequency `Ω` and density independent coherence time `τ2`.

The scale factor `K` is applied after normalising the integral of the response function
to unity.
"""
function RamanRespNormedSingleDampedOscillator(K, Ω, τ2)
    K *= (Ω^2 + 1/τ2^2)/Ω # normalise SDO model to unity integral, scaled by prefactor K.
    τ2ρ = let τ2=τ2
        ρ -> τ2
    end
    RamanRespSingleDampedOscillator(K, Ω, τ2ρ)
end

function hrpre(R::RamanRespSingleDampedOscillator, t)
    t > 0.0 ? R.K*sin(R.Ω*t) : 0.0
end

hrdamp(R::RamanRespSingleDampedOscillator, ρ) = R.τ2ρ(ρ)

"""
    RamanRespVibrational(Ωv, dαdQ, μ; τ2=nothing, Bρ=nothing, Aρ=nothing)

Construct a molecular vibrational Raman model (single damped oscillator).

# Arguments
- `Ωv::Real`: vibrational frequency [rad/s]
- `dαdQ::Real`: isotropic averaged polarizability derivative [m^2]
- `μ::Real`: reduced molecular mass [kg]
- `τ2::Real=nothing`: coherence time [s]
- `Bρ::Real=nothing` : density dependent broadening coefficient [Hz/amagat]
- `Aρ::Real=nothing` : self diffusion coefficient [Hz amagat]

Only one of `τ2` or `Bρ` should be specified.
If `Bρ` is specified then `Aρ` must be too.

# References
- Full model description: To be published, Yingying paper.
- We followed closely: Phys. Rev. A, vol. 92, no. 6, p. 063828, Dec. 2015,
  But note that that paper uses weird units, and we converted it to SI for
  the above reference.
"""
function RamanRespVibrational(Ωv, dαdQ, μ; τ2=nothing, Bρ=nothing, Aρ=nothing)
    K = (4π*ε_0)^2*dαdQ^2/(4μ*Ωv)
    τ2ρ = if isnothing(Bρ)
        isnothing(τ2) && error("one of `τ2` or `Bρ` must be specified")
        let τ2=τ2
            ρ -> τ2
        end
    else
        !isnothing(τ2) && error("only one of `τ2` or `Bρ` must be specified")
        isnothing(Aρ) && error("if `Bρ` is specified you must also specify `Aρ`")
        let Bρ=Bρ, Aρ=Aρ
            ρ -> 1/(pi*(Aρ/(ρ/amg) + Bρ*ρ/amg))
        end
    end
    RamanRespSingleDampedOscillator(K, Ωv, τ2ρ)
end

# TODO: we assume here that all rotational levels have same linewidth
struct RamanRespRotationalNonRigid{TR, Tτ2} <: AbstractRamanResponse
    Rs::TR # List of Raman responses
    τ2ρ::Tτ2 # coherence dephasing time function density to τ2: ρ -> τ2
end

"""
    RamanRespRotationalNonRigid(B, Δα, τ2, qJodd, qJeven;
                                D=0.0, minJ=0, maxJ=50, temp=roomtemp,
                                τ2=nothing, Bρ=nothing, Aρ=nothing)

Construct a rotational nonrigid rotor Raman model.

# Arguments
- `B::Real`: the rotational constant [1/m]
- `Δα::Real`: molecular polarizability anisotropy [m^3]
- `qJodd::Integer`: nuclear spin parameter for odd `J`
- `qJeven::Integer`: nuclear spin parameter for even `J`
- `D::Real=0.0`: centrifugal constant [1/m]
- `minJ::Integer=0`: J value to start at
- `maxJ::Integer=50`: J value to sum until
- `temp::Real=roomtemp`: temperature
- `τ2::Real=nothing`: coherence time [s]
- `Bρ::Real=nothing` : density dependent broadening coefficient [Hz/amagat]
- `Aρ::Real=nothing` : self diffusion coefficient [Hz amagat]

Only one of `τ2` or `Bρ` should be specified.
If `Bρ` is specified then `Aρ` must be too.

# References
- Full model description: To be published, Yingying paper.
- We followed closely: Phys. Rev. A, vol. 92, no. 6, p. 063828, Dec. 2015,
  But note that that paper uses weird units, and we converted it to SI for
  the above reference.
"""
function RamanRespRotationalNonRigid(B, Δα, qJodd::Int, qJeven::Int;
    D=0.0, minJ=0, maxJ=50, temp=roomtemp, τ2=nothing, Bρ=nothing, Aρ=nothing)
    J = minJ:maxJ # range of J values to start with
    EJ = @. 2π*ħ*c*(B*J*(J + 1) - D*(J*(J + 1))^2) # energy of each J level
    # limit J range to those which have monotonic increasing energy
    mJ = findfirst(x -> x < 0.0, diff(EJ))
    if !isnothing(mJ)
        if length(minJ:mJ) <= 2 # need at least 1 pair of J levels
            error("Raman rigid rotation model cannot sum over levels")
        end
        J = J[1:mJ]
        EJ = EJ[1:mJ]
    end
    # get nuclear degeneracy
    DJ = map(x -> isodd(x) ? qJodd : qJeven, J)
    # get population of each level
    ρ = @. DJ*(2*J + 1)*exp(-EJ/(k_B*temp))
    ρ ./= sum(ρ)
    # we need pairs of J, J+2 levels, so limit J to just the starting levels
    J = J[1:end-2]
    # angular frequency of each J, J+2 pair
    Ω = (EJ[3:end] .- EJ[1:end-2])./ħ
    # absolute Raman prefactor
    K = -(4π*ε_0)^2*2/15*Δα^2/ħ
    τ2ρ = if isnothing(Bρ)
              isnothing(τ2) && error("one of `τ2` or `Bρ` must be specified")
              let τ2=τ2
                  ρ -> τ2
              end
          else
              !isnothing(τ2) && error("only one of `τ2` or `Bρ` must be specified")
              isnothing(Aρ) && error("if `Bρ` is specified you must also specify `Aρ`")
              let Bρ=Bρ, Aρ=Aρ
                  ρ -> 1/(pi*(Aρ/(ρ/amg) + Bρ*ρ/amg))
              end
          end
    Rs = [RamanRespSingleDampedOscillator((K*(J[i] + 1)*(J[i] + 2)/(2*J[i] + 3)
                                            *(ρ[i+2]/(2*J[i] + 5) - ρ[i]/(2*J[i] + 1))),
                                          Ω[i], τ2ρ) for i=1:length(J)]
    RamanRespRotationalNonRigid(Rs, τ2ρ)
end

hrpre(R::RamanRespRotationalNonRigid, t) = sum(hrpre(Ri, t) for Ri in R.Rs)

hrdamp(R::RamanRespRotationalNonRigid, ρ) = R.τ2ρ(ρ)

struct CombinedRamanResponse
    Rs::Vector{Any} # list of Raman responses
    t::Vector{Float64} # time grid
    w::Vector{Float64} # filter window
    hpres::Vector{Vector{Float64}} # pre Raman responses for each R in Rs
end

function CombinedRamanResponse(t, Rs)
    tt = collect(0:(length(t) - 1)) .* (t[2] - t[1])
    hpres = [hrpre.(R, tt) for R in Rs]
    w = planck_taper(tt, -tt[end], -tt[end]*0.7, tt[end]*0.7, tt[end])
    CombinedRamanResponse(Rs, tt, w, hpres)
end

function (R::CombinedRamanResponse)(ht::AbstractVector, ρ)
    fill!(ht, 0.0)
    for i=1:length(R.Rs)
        ht .+= R.hpres[i] .* exp.(-R.t ./ hrdamp.(R.Rs[i], ρ))
    end
    ht .*= R.w
end

(R::CombinedRamanResponse)(t::Number, ρ) = sum(Ri(t, ρ) for Ri in R.Rs)

"""
    molecular_raman_response(rp; kwargs...)

Get the Raman response function for the Raman parameters in named tuple `rp`.

# Keyword Arguments
- `rotation::Bool = true`: whether to include the rotational contribution
- `vibration::Bool = true`: whether to include the vibrational contribution
- `minJ::Integer = 0`: the minimum rotational quantum number to include
- `maxJ::Integer = 50`: the maximum rotational quantum number to include
- `temp::Real = roomtemp`: the temperature
"""
function molecular_raman_response(t, rp; rotation=true, vibration=true, minJ=0, maxJ=50, temp=roomtemp)
    Rs = []
    if rotation && (rp.rotation != :none)
        if rp.rotation != :nonrigid
            throw(DomainError(rp.rotation, "Unknown Rotational Raman model $(rp.rotation)"))
        end
        if haskey(rp, :Bρr)
            hr = RamanRespRotationalNonRigid(rp.B, rp.Δα, rp.qJodd, rp.qJeven, Bρ=rp.Bρr, Aρ=rp.Aρr,
                                            D=rp.D, minJ=minJ, maxJ=maxJ, temp=temp)
        else
            hr = RamanRespRotationalNonRigid(rp.B, rp.Δα, rp.qJodd, rp.qJeven, τ2=rp.τ2r,
                                            D=rp.D, minJ=minJ, maxJ=maxJ, temp=temp)
        end
        push!(Rs, hr)
    end
    if vibration  && (rp.vibration != :none)
        if rp.vibration != :sdo
            throw(DomainError(rp.rotation, "Unknown Vibrational Raman model $(rp.vibration)"))
        end
        if haskey(rp, :Bρv)
            hv = RamanRespVibrational(rp.Ωv, rp.dαdQ, rp.μ, Bρ=rp.Bρv, Aρ=rp.Aρv)
        else
            hv = RamanRespVibrational(rp.Ωv, rp.dαdQ, rp.μ, τ2=rp.τ2v)
        end
        push!(Rs, hv)
    end
    CombinedRamanResponse(t, Rs)
end

"""
    raman_response(material; kwargs...)

Get the Raman response function for `material`.

For details on the keyword arguments see [`molecular_raman_response`](@ref).
"""
function raman_response(t, material; kwargs...)
    rp = raman_parameters(material)
    if rp.kind == :molecular
        return molecular_raman_response(t, rp; kwargs...)
    elseif rp.kind == :normedsdo
        return CombinedRamanResponse(t, [RamanRespNormedSingleDampedOscillator(rp.K, rp.Ω, rp.τ2)])
    else
        throw(DomainError(rp.kind, "Unknown Raman model $(rp.kind)"))
    end
end

end


text:
using Test
using SolidStateDetectors
using StaticArrays

import SolidStateDetectors.ConstructiveSolidGeometry as CSG
import SolidStateDetectors.ConstructiveSolidGeometry: Geometry, Dictionary

T = Float64

default_units = SolidStateDetectors.default_unit_tuple()
no_translations = (rotation = one(SMatrix{3, 3, T, 9}), translation = zero(CartesianVector{T}))

@testset "Test primitive parsing" begin
    @testset "Cone" begin
        for bot in (2.0, Dict("from" => 1.0, "to" => 2.0)),
            top in (2.0, 1.0, Dict("from" => 1.0, "to" => 2.0), Dict("from" => 2.0, "to" => 4.0)),
            φ in ((30,180),(0,360))

            dict = Dict("difference" => [
                Dict("cone"   => Dict(
                    "r"       => Dict("bottom" => bot, "top" => top),
                    "phi"     => Dict("from" => "$(φ[1])°", "to" => "$(φ[2])°"),
                    "h"       => 1.0))
                for i in 1:2
            ])

            c = Geometry(T, dict, default_units, no_translations)

            # Conversion from Geometry -> Dict and Dict -> Geometry should result in the same geometry
            output = Dictionary(c)
            name = collect(keys(output["difference"][1]))[1] # "tube" or "cone"
            if haskey(output["difference"][1][name], "phi")
                @test CSG._parse_value(T, output["difference"][1][name]["phi"]["to"], default_units.angle) == CSG._parse_value(T, "$(φ[2])°", default_units.angle)
            end
            # Internally, the primitive stores the offset of phi in a rotation, but it should not be part of the config file
            @test !haskey(output["difference"][1], "rotation")
            @test c.a.rotation ≈ Geometry(T, output, default_units, no_translations).a.rotation

            # No warnings or errors when decomposing the Cones into surfaces
            @test_nowarn CSG.surfaces(c.a)
            @test_nowarn CSG.surfaces(c.b)

            # Check if all Cones are saved the right way
            @test c.a isa CSG.Cone
            @test c.b isa CSG.Cone
        end

        # Test different constructor versions and dictionary construction
        cone1 = @inferred CSG.Cone(CSG.ClosedPrimitive,r=1f0, φ = π, hZ = 1f0, origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
        cone2 = @inferred CSG.Cone{Float32}(r=1.0, φ = π)
        @test cone1 === cone2
        dict = Dict("tube"   => Dict(
                "r"       => 1.0,
                "h"       => 1.0))
        cone4 = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(cone4)
        @test dict == output

        ## Test in method for several cones
        # Non transformed geometries
        @test !in(CartesianPoint{Float32}(0,-1,0),cone1)
        @test in(CartesianPoint{Float64}(0,-1,0),cone4)
        cone_open = @inferred CSG.Cone(CSG.OpenPrimitive, r=1f0,hZ=1f0, origin = zero(CartesianPoint{Float32}),rotation = one(SMatrix{3, 3, Float32, 9}))
        @test !in(CartesianPoint{Float32}(1,0,0),cone_open)
        # Boundary behaviour of closed and open primitives for transformed
        # geometries (origin, rotation)
        cone_closed_trafo = @inferred CSG.Cone(CSG.ClosedPrimitive,r=1.0, hZ=3.0, origin=CartesianPoint{Float32}(0,0,1),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))#90 degree rotation arond y-axis
        cone_open_trafo = @inferred CSG.Cone(CSG.OpenPrimitive,r=1.0, hZ=3.0, origin=CartesianPoint{Float32}(0,0,1),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        @test in(CartesianPoint{Float64}(3,0,1),cone_closed_trafo)
        @test !in(CartesianPoint{Float64}(3,0,1),cone_open_trafo)
        tol=1e-8
        @test !in(CartesianPoint{Float64}(3+tol,0,1),cone_closed_trafo)
        @test in(CartesianPoint{Float64}(3-tol,0,1),cone_open_trafo)
        #Cones with r being any valid combination
        cone_empty = @inferred CSG.Cone{Float64}(r=((1,2),(3,5)), hZ=2, rotation = SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        @test in(CartesianPoint{Float64}(1,0,2.5+tol), cone_empty)
        @test !in(CartesianPoint{Float64}(1,0,2.5-tol), cone_empty)
        @inferred CSG.Cone{Float64}(r=(nothing,(3,5)))
        @inferred CSG.Cone{Float64}(r=((3,5),nothing))
        @inferred CSG.Cone{Float64}(r=(1,2))
        @inferred CSG.Cone{Float64}(r=(nothing,(3,5)))
        @inferred CSG.Cone{Float64}(r=((1,),(2,)))
        @inferred CSG.Cone(r=(nothing,(3,5)))
        @inferred CSG.Cone(r=((3,5),nothing))
        @inferred CSG.Cone(r=(1,2))
        @inferred CSG.Cone(r=(nothing,(3,5)))
        @inferred CSG.Cone(r=((1,),(2,)))
        #Test where phi is a Tuple
        tuple_cone = @inferred CSG.Cone{Float64}(φ=(π/4,3*π/4), rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        rot_cone = @inferred CSG.Cone{Float64}(φ=π/2, rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0) * SMatrix{3}(cos(π/4),sin(π/4),0,-sin(π/4),cos(π/4),0,0,0,1))
        @test tuple_cone ==rot_cone
    end
    @testset "Torus" begin
        for r_tube in (2.0, Dict("from" => 1.0, "to" => 2.0)),
            φ in ((30,180),(0,360)),
            θmin in 0:45:90,
            θmax in 135:45:270

            dict = Dict("difference" => [
                Dict( "torus" => Dict(
                    "r_torus" => 10.0,
                    "r_tube"  => r_tube,
                    "phi"     => Dict("from" => "$(φ[1])°", "to" => "$(φ[2])°"),
                    "theta"   => Dict("from" => "$(θmin)°", "to" => "$(θmax)°")))
                for i in 1:2
            ])
            t = Geometry(T, dict, default_units, no_translations)

            # Conversion from Geometry -> Dict and Dict -> Geometry should result in the same geometry
            output = Dictionary(t)
            if haskey(output["difference"][1]["torus"], "phi")
                @test CSG._parse_value(T, output["difference"][1]["torus"]["phi"]["to"], default_units.angle) == CSG._parse_value(T, "$(φ[2])°", default_units.angle)
            end
            # Internally, the primitive stores the offset of phi in a rotation, but it should not be part of the config file
            @test !haskey(output["difference"][1], "rotation")
            @test t.a.rotation ≈ Geometry(T, output, default_units, no_translations).a.rotation

            # No warnings or errors when decomposing the Torus into surfaces
            @test_nowarn CSG.surfaces(t.a)
            @test_nowarn CSG.surfaces(t.b)

            # Check if all Torus are saved the right way
            @test t.a isa CSG.Torus
            @test t.b isa CSG.Torus
        end

        # Test different constructor versions and dictionary construction
        torus1 = @inferred CSG.Torus(CSG.ClosedPrimitive, r_torus = 1f0, r_tube = 1f0)
        torus2 = @inferred CSG.Torus{Float32}(r_torus = 1.0, r_tube = 1.0)
        @test torus1 == torus2
        dict = Dict("torus"   => Dict(
                "r_torus"       => 1.0,
                "r_tube"        => 1.0))
                # "phi"           => π,
                # "theta"         => (π,2π)))
        torus4 = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(torus4)
        @test dict == output

        ## Test in method for several tori
        # Non transformed geometries
        @test !in(CartesianPoint{Float32}(0,0,0.5),torus1)
        @test in(CartesianPoint{Float32}(0,0,0),torus1)
        torus1_open = @inferred CSG.Torus(CSG.OpenPrimitive, r_torus = 1f0, r_tube = 1f0)
        @test in(CartesianPoint{Float32}(0,2,0),torus1)
        @test !in(CartesianPoint{Float32}(0,2,0),torus1_open)
        # Boundary behaviour of closed and open primitives for transformed
        # geometries (origin, rotation)
        torus_closed_trafo = @inferred CSG.Torus(CSG.ClosedPrimitive,r_torus=2.0, r_tube=3.0, φ = π, origin=CartesianPoint{Float32}(0,0,1),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        torus_open_trafo = @inferred CSG.Torus(CSG.OpenPrimitive,r_torus=2.0, r_tube=3.0, φ = π, origin=CartesianPoint{Float32}(0,0,1),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))

        @test in(CartesianPoint{Float64}(1,0,-1),torus_closed_trafo)
        @test !in(CartesianPoint{Float64}(1,0,-1),torus_open_trafo)
        @test in(CartesianPoint{Float64}(0,0,0),torus_closed_trafo)
        @test !in(CartesianPoint{Float64}(0,0,0),torus_open_trafo)
        tol = 1e-8
        @test !in(CartesianPoint{Float64}(0,-tol,0),torus_closed_trafo)
        @test in(CartesianPoint{Float64}(0,+tol,0),torus_open_trafo)
        #Test all valid keyword argument types
        @inferred CSG.Torus{Float64}(φ=nothing,θ=nothing)
        CSG.Torus{Float64}(φ=π,θ=(π,2π)) #inferred does not work here, due to evaluation of _get_conemantle_type(θ) during runtime
        #Test where phi is a Tuple
        tuple_torus = @inferred CSG.Torus{Float64}(φ=(π/4,3*π/4),θ=nothing,rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        rot_torus = @inferred CSG.Torus{Float64}(φ=π/2,θ=nothing,rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0) * SMatrix{3}(cos(π/4),sin(π/4),0,-sin(π/4),cos(π/4),0,0,0,1))
        @test tuple_torus ==rot_torus
    end
    @testset "Ellipsoid" begin
        ellip1 = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1f0,origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
        ellip2 = @inferred CSG.Ellipsoid{Float32}(r=1.0)
        @test ellip1 === ellip2

        dict = Dict("sphere"   => Dict(
                "r"       => 1.0))
        ellip3 = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(ellip3)
        @test dict == output

        ellip_open = @inferred CSG.Ellipsoid(CSG.OpenPrimitive,r=1f0,origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
        @test in(CartesianPoint{Float32}(1,0,0),ellip1)
        @test !in(CartesianPoint{Float32}(1,0,0),ellip_open)

        ellip_closed_trafo = @inferred CSG.Ellipsoid(CSG.ClosedPrimitive,r=1.0, origin=1/sqrt(3)*CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
        ellip_open_trafo = @inferred CSG.Ellipsoid(CSG.OpenPrimitive,r=1.0, origin=1/sqrt(3)*CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0.5,sqrt(3)/2,0,-sqrt(3)/2,0.5,0,0,0,1))
        @test !in(CartesianPoint{Float64}(-1e-8,0,0),ellip_closed_trafo)
        @test in(CartesianPoint{Float64}(1e-8,0,0),ellip_open_trafo)
    end
    @testset "Box" begin
        box1 = @inferred CSG.Box(CSG.ClosedPrimitive,hX=1f0, hY=2f0, hZ=1f0, origin = zero(CartesianPoint{Float16}),rotation = one(SMatrix{3, 3, Float16, 9}))
        box2 = @inferred CSG.Box{Float32}(hX=1.0, hY=2f0, hZ=1f0)
        @test box1 === box2

        dict = Dict("box"   => Dict(
                "hX"       => 1.0,
                "hY"       => 42.0,
                "hZ"       => 3.14))
        box3 = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(box3)
        @test dict == output

        box_closed_trafo = @inferred CSG.Box(CSG.ClosedPrimitive, hX=2., origin=CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0,0,1,0,1,0,-1,0,0))
        box_open_trafo = @inferred CSG.Box(CSG.OpenPrimitive, hX=2., origin=CartesianPoint{Float32}(1,1,1),rotation=SMatrix{3}(0,0,1,0,1,0,-1,0,0))
        tol=1e-8
        @test in(CartesianPoint{Float64}(1,1,3),box_closed_trafo)
        @test !in(CartesianPoint{Float64}(1,1,3),box_open_trafo)
        @test !in(CartesianPoint{Float64}(1,1,3+tol),box_closed_trafo)
        @test in(CartesianPoint{Float64}(1,1,3-tol),box_open_trafo)
    end
    @testset "RegularPrism" begin
        prism1 = @inferred CSG.RegularPrism{3}(CSG.ClosedPrimitive,r=1f0, hZ = 2f0)
        prism2 = @inferred CSG.RegularPrism{3,Float32}(r=1.0, hZ = 2f0)
        @test prism1 === prism2

        dict = Dict("TriangularPrism"   => Dict(
                "r"       => 1.0,
                "h"        => 1.0))
        prism = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(prism)
        @test dict == output

        dict = Dict("QuadranglePrism"   => Dict(
                "r"       => 1.0,
                "z"      =>1.0))
        prism = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(prism)
        output["QuadranglePrism"]["z"]=output["QuadranglePrism"]["h"]#Workaround for tests, since dictionary uses the key "h", since "z" is outdated
        delete!(output["QuadranglePrism"], "h")
        @test dict == output

        dict = Dict("PentagonalPrism"   => Dict(
                "r"       => 1.0,
                "z"      =>(1.0,2.0)))
        prism = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(prism)
        output["PentagonalPrism"]["z"]=output["PentagonalPrism"]["h"]
        delete!(output["PentagonalPrism"], "h")
        #@test dict == output

        dict = Dict("HexagonalPrism"   => Dict(
                "r"       => 1.0,
                "h"      =>1.0))
        prism = Geometry(T,dict,default_units,no_translations)
        output = Dictionary(prism)
        @test dict == output

        # Test in() method
        @test in(CartesianPoint{Float32}(0,0,0),prism1)
        @test !in(CartesianPoint{Float32}(1.5,0,0),prism1)
        tol = 1e-8
        prism_closed = @inferred CSG.RegularPrism{4}(CSG.ClosedPrimitive; r =1.,hZ =1., origin=CartesianPoint{Float32}(0,0,5),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        prism_open = @inferred CSG.RegularPrism{5}(CSG.OpenPrimitive; r =1.,hZ =1., origin=CartesianPoint{Float32}(0,0,5),rotation=SMatrix{3}(0,0,-1,0,1,0,1,0,0))
        @test in(CartesianPoint{Float64}(0,0,4),prism_closed)
        @test !in(CartesianPoint{Float64}(0,0,4),prism_open)
        @test !in(CartesianPoint{Float64}(0,0,4-tol),prism_closed)
        @test in(CartesianPoint{Float64}(0,0,4+tol),prism_open)

        #Test different prism
        @inferred CSG.RegularPrism{3,Float32}()
        @inferred CSG.RegularPrism{4,Float32}()
        @inferred CSG.RegularPrism{5,Float32}()
        @inferred CSG.RegularPrism{6,Float32}()
    end
end


text:
immutable Chi <: ContinuousUnivariateDistribution
    df::Float64
end