text:
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                         SLATER-TYPE ORBİTALS ONE ELECTRON MOLECULAR INTEGRALS                           #
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
export TwoCenterOverlap, CTwoCenterOverlap, VTwoCenterOverlap, STwoCenterOverlap, TwoCenterOverlapRec,
TwoCenterNucAttractABA, CTwoCenterNucAttractABA, VTwoCenterNucAttractABA, STwoCenterNucAttractABA,
TwoCenterNucAttractAAB, CTwoCenterNucAttractAAB, VTwoCenterNucAttractAAB, VTwoCenterNucAttractAAB,
TwoCenterKinEnergy, CTwoCenterKinEnergy, VTwoCenterKinEnergy, STwoCenterKinEnergy
###########################################################################################################
###########################################################################################################
####################################### TWO CENTER OVERLAP INTEGRALS ######################################(43)
################################################### ALIGNED MOLECULAR COORDINATE SYSTEM
function TwoCenterOverlap(n1::arb, l1::Int, n2::arb, l2::Int, λ::Int, ρ::arb, τ::arb,lim::Int)
    Overlap = RF(0)
    for s1 in (EvenQ(l1-λ) ? -λ : -λ + 1)  : 2 : l1
        for s2 in (EvenQ(l2+λ) ? λ : λ + 1)  : 2 : l2
            for s3 in 0 : s1 + s2
                Overlap += SphPCG(s3, s1, s2, l1, λ, l2, λ, λ)*
                AuxiliaryG(RF(10//10), s3, RF(n1 - s1), RF(n2 - s2), RF(10//10), ρ, ρ * τ, lim)
            end
        end
    end
    return Power(-1,l2 + λ) * SlaterON(n1, n2, ρ, τ) * Overlap
end

function CTwoCenterOverlap(n1::arb, l1::Int, n2::arb, l2::Int, λ::Int, ρ::arb, τ::arb,lim::Float64)
    Overlap = RF(0)
    for s1 in (EvenQ(l1-λ) ? -λ : -λ + 1)  : 2 : l1
        for s2 in (EvenQ(l2+λ) ? λ : λ + 1)  : 2 : l2
            for s3 in 0 : s1 + s2
                Overlap += SphPCG(s3, s1, s2, l1, λ, l2, λ, λ)*
                CuhreAuxiliaryG(10//10, s3, n1 - s1, n2 - s2, 10//10, ρ, ρ * τ, lim)[1]
            end
        end
    end
    return Power(-1,l2 + λ) * SlaterON(n1, n2, ρ, τ) * Overlap
end

function VTwoCenterOverlap(n1::arb, l1::Int, n2::arb, l2::Int, λ::Int, ρ::arb, τ::arb,lim::Float64)
    Overlap = RF(0)
    for s1 in (EvenQ(l1-λ) ? -λ : -λ + 1)  : 2 : l1
        for s2 in (EvenQ(l2+λ) ? λ : λ + 1)  : 2 : l2
            for s3 in 0 : s1 + s2
                Overlap += SphPCG(s3, s1, s2, l1, λ, l2, λ, λ)*
                VegasAuxiliaryG(10//10, s3, n1 - s1, n2 - s2, 10//10, ρ, ρ * τ, lim)[1]
            end
        end
    end
    return Power(-1,l2 + λ) * SlaterON(n1, n2, ρ, τ) * Overlap
end

function STwoCenterOverlap(n1::arb, l1::Int, n2::arb, l2::Int, λ::Int, ρ::arb, τ::arb,lim::Float64)
    Overlap = RF(0)
    for s1 in (EvenQ(l1-λ) ? -λ : -λ + 1)  : 2 : l1
        for s2 in (EvenQ(l2+λ) ? λ : λ + 1)  : 2 : l2
            for s3 in 0 : s1 + s2
                Overlap += SphPCG(s3, s1, s2, l1, λ, l2, λ, λ)*
                SuaveAuxiliaryG(10//10, s3, n1 - s1, n2 - s2, 10//10, ρ, ρ * τ, lim)[1]
            end
        end
    end
    return Power(-1,l2 + λ) * SlaterON(n1, n2, ρ, τ) * Overlap
end

function TwoCenterOverlapRec(n1::arb, l1::Int, n2::arb, l2::Int, λ::Int, ρ::arb, τ::arb,lim::Int)
    OverlapRec = RF(0)
    for s1 in (EvenQ(l1-λ) ? -λ : -λ + 1)  : 2 : l1
        for s2 in (EvenQ(l2+λ) ? λ : λ + 1)  : 2 : l2
            for s3 in 0 : s1 + s2
                OverlapRec += SphPCG(s3, s1, s2, l1, λ, l2, λ, λ)*
                AuxiliaryGr(RF(10//10), s3, RF(n1 - s1), RF(n2 - s2), RF(10//10), ρ, ρ * τ, lim)
            end
        end
    end
    return SlaterON(n1, n2, ρ, τ)*OverlapRec
end
###########################################################################################################
function TwoCenterOverlap(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Int)

    TwoCenterOverlap(n1, l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), lim)
end

function CTwoCenterOverlap(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    CTwoCenterOverlap(n1, l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), lim)
end

function VTwoCenterOverlap(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    VTwoCenterOverlap(n1, l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), lim)
end

function STwoCenterOverlap(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    STwoCenterOverlap(n1, l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), lim)
end

function TwoCenterOverlapRec(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Int)

    TwoCenterOverlapRec(n1, l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), lim)
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
################################################### NON-ALIGNED MOLECULAR COORDINATE SYSTEM
function  TwoCenterOverlap(
n1::arb, l1::Int, m1::Int,
n2::arb, l2::Int, m2::Int,
ρ::arb, τ::arb, theta::arb, phi::arb,
x::Int, lim::Int)

    NOverlap = RF(0)

    for s in 0 : min(l1,l2)
        NOverlap += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        TwoCenterOverlap(n1, l1, n2, l2, s, ρ, τ,lim))
    end
    NOverlap
end

function  CTwoCenterOverlap(
    n1::arb, l1::Int, m1::Int,
    n2::arb, l2::Int, m2::Int,
    ρ::arb, τ::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    NOverlap = RF(0)

    for s in 0 : min(l1,l2)
        NOverlap += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        CTwoCenterOverlap(n1, l1, n2, l2, s, ρ, τ,lim))
    end
    NOverlap
end

function  VTwoCenterOverlap(
    n1::arb, l1::Int, m1::Int,
    n2::arb, l2::Int, m2::Int,
    ρ::arb, τ::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    NOverlap = RF(0)

    for s in 0 : min(l1,l2)
        NOverlap += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        VTwoCenterOverlap(n1, l1, n2, l2, s, ρ, τ,lim))
    end
     NOverlap
end

function  STwoCenterOverlap(
    n1::arb, l1::Int, m1::Int,
    n2::arb, l2::Int, m2::Int,
    ρ::arb, τ::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    NOverlap = RF(0)

    for s in 0 : min(l1,l2)
        NOverlap += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        STwoCenterOverlap(n1, l1, n2, l2, s, ρ, τ,lim))
    end
    NOverlap
end

function  TwoCenterOverlapRec(
    n1::arb, l1::Int, m1::Int,
    n2::arb, l2::Int, m2::Int,
    ρ::arb, τ::arb, theta::arb, phi::arb,
    x::Int, lim::Int)
        NOverlap = RF(0)

    for s in 0 : min(l1,l2)
        NOverlap += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        TwoCenterOverlapRec(n1, l1, n2, l2, s, ρ, τ,lim))
    end
        NOverlap
end

function TwoCenterOverlapRec(
    n1::arb, l1::Int, m1:: Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, x::Int, lim::Int)

    TwoCenterOverlapRec(n1, l1, m1, n2, l2, m2, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), theta, phi, x, lim)
end
###########################################################################################################
###########################################################################################################
################################# TWO CENTER NUCLEAR ATTRACTION INTEGRALS #################################(43)
################################################### ALIGNED MOLECULAR COORDINATE SYSTEM
function  TwoCenterNucAttractABA(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Int)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = TwoCenterOverlap(n1 - RF(1), l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2),lim)

    return coef * restco
end

function  CTwoCenterNucAttractABA(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = CTwoCenterOverlap(n1 - RF(1), l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2),lim)

    return coef * restco
end

function  VTwoCenterNucAttractABA(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = VTwoCenterOverlap(n1 - RF(1), l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2),lim)

    return coef * restco
end

function  STwoCenterNucAttractABA(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = STwoCenterOverlap(n1 - RF(1), l1, n2, l2, λ, (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2),lim)

    return coef * restco
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
################################################### NON-ALIGNED MOLECULAR COORDINATE SYSTEM
function  TwoCenterNucAttractABA(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb,
    x::Int, lim::Int)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = TwoCenterOverlap(
    n1 - RF(1), l1, m1,
    n2, l2, m2,
    (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), theta, phi, x, lim)

    return coef * restco
end

function  CTwoCenterNucAttractABA(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = CTwoCenterOverlap(
    n1 - RF(1), l1, m1,
    n2, l2, m2,
    (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), theta, phi, x, lim)

    return coef * restco
end

function  VTwoCenterNucAttractABA(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = VTwoCenterOverlap(
    n1 - RF(1), l1, m1,
    n2, l2, m2,
    (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), theta, phi, x, lim)

    return coef * restco
end

function  STwoCenterNucAttractABA(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb,
    x::Int, lim::Float64)

    coef = (2 * ζ1) // (Sqrt(2 * n1 * (2 * n1 - 1)))
    restco = STwoCenterOverlap(
    n1 - RF(1), l1, m1,
    n2, l2, m2,
    (R//2)*(ζ1+ζ2), (ζ1-ζ2)//(ζ1+ζ2), theta, phi, x, lim)

    return coef * restco
end
###########################################################################################################
function TwoCenterNucAttractAAB(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Int
)
    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)
    n2z = RF(0.000000000000000000000000000000000000000000000001)

    coeff = (R // 2) * SlaterON(n1, n2, ρ, τ)

    sresaab = RF(0)
    for s1 in abs(l1 - l2) : l1 + l2
        for s2 in (EvenQ(l1 + l2 - λ) ? -λ : -λ + 1) : 2 : l1 + l2
            for s3 in 0 : s2
                sresaab1 = Sqrt(2 * s1 + 1)
                sresaab2 = GauntGF(l1, λ, l2, λ, s1)
                sresaab3 = SphPCG(s3, s2, 0, s1, 0, 0, 0, 0)
                sresaab4 = AuxiliaryG(RF(10//10), s3, n1 + n2 - s2 - 1, n2z, RF(10//10), ρ, ρ, lim)
                sresaab += sresaab1 * sresaab2 * sresaab3 * sresaab4
            end
        end
    end

    resaab = coeff * sresaab
    return resaab
end

function CTwoCenterNucAttractAAB(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
)
    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)
    n2z = RF(0.000000000000000000000000000000000000000000000001)

    coeff = (R // 2) * SlaterON(n1, n2, ρ, τ)

    sresaab = RF(0)
    for s1 in abs(l1 - l2) : l1 + l2
        for s2 in (EvenQ(l1 + l2 - λ) ? -λ : -λ + 1) : 2 : l1 + l2
            for s3 in 0 : s2
                sresaab1 = Sqrt(2 * s1 + 1)
                sresaab2 = GauntGF(l1, λ, l2, λ, s1)
                sresaab3 = SphPCG(s3, s2, 0, s1, 0, 0, 0, 0)
                sresaab4 = CuhreAuxiliaryG(RF(10//10), s3, n1 + n2 - s2 - 1, n2z, RF(10//10), ρ, ρ, lim)[1]
                sresaab += sresaab1 * sresaab2 * sresaab3 * sresaab4
            end
        end
    end

    resaab = coeff * sresaab
    return resaab
end

function VTwoCenterNucAttractAAB(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
)
    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)
    n2z = RF(0.000000000000000000000000000000000000000000000001)

    coeff = (R // 2) * SlaterON(n1, n2, ρ, τ)

    sresaab = RF(0)
    for s1 in abs(l1 - l2) : l1 + l2
        for s2 in (EvenQ(l1 + l2 - λ) ? -λ : -λ + 1) : 2 : l1 + l2
            for s3 in 0 : s2
                sresaab1 = Sqrt(2 * s1 + 1)
                sresaab2 = GauntGF(l1, λ, l2, λ, s1)
                sresaab3 = SphPCG(s3, s2, 0, s1, 0, 0, 0, 0)
                sresaab4 = VegasAuxiliaryG(RF(10//10), s3, n1 + n2 - s2 - 1, n2z, RF(10//10), ρ, ρ, lim)[1]
                sresaab += sresaab1 * sresaab2 * sresaab3 * sresaab4
            end
        end
    end

    resaab = coeff * sresaab
    return resaab
end

function STwoCenterNucAttractAAB(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
)
    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)
    n2z = RF(0.000000000000000000000000000000000000000000000001)

    coeff = (R // 2) * SlaterON(n1, n2, ρ, τ)

    sresaab = RF(0)
    for s1 in abs(l1 - l2) : l1 + l2
        for s2 in (EvenQ(l1 + l2 - λ) ? -λ : -λ + 1) : 2 : l1 + l2
            for s3 in 0 : s2
                sresaab1 = Sqrt(2 * s1 + 1)
                sresaab2 = GauntGF(l1, λ, l2, λ, s1)
                sresaab3 = SphPCG(s3, s2, 0, s1, 0, 0, 0, 0)
                sresaab4 = SuaveAuxiliaryG(RF(10//10), s3, n1 + n2 - s2 - 1, n2z, RF(10//10), ρ, ρ, lim)[1]
                sresaab += sresaab1 * sresaab2 * sresaab3 * sresaab4
            end
        end
    end

    resaab = coeff * sresaab
    return resaab
end
###########################################################################################################
#function  TwoCenterNucAttractAAB(
#    n1::arb, l1::Int, m1::Int, ζ1::arb,
#    n2::arb, l2::Int, m2::Int, ζ2::arb,
#    R::arb, theta::arb, phi::arb,
#    x::Int, lim::Int
#    )

#    lowL = abs(l1-l2)
#    upL = l1 + l2
#    n2z = RF(0.000000000000000000000000000000000000000000000001)
#    ρ = (R//2)*(ζ1+ζ2)
#    τ = (ζ1-ζ2)//(ζ1+ζ2)

#    resaab = RF(0)
#    for s1 in lowL : 2 : upL
#        for s2 in -s1 : s1
#            resaab += ChargeDOneC(n1,l1,m1,ζ1,n2,l2,m2,ζ2,n1+n2-1,s1,s2) *
#            TwoCenterOverlap(n1+n2-1,s1,s2,n2z,0,0,ρ,τ,theta,phi,1,lim)
#        end
#    end

#    return resaab
#end

function  TwoCenterNucAttractAAB(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, lim::Int
    )

    lowL = abs(l1-l2)
    upL = l1 + l2

    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    resaab = RF(0)
    for s1 in lowL : 2 : upL
        for s2 in -s1 : s1
            resaab1 = Sqrt((4 * Pi(RF)) // (2 * s1 + 1))
            resaab += (
            GGauntG(l1, m1, l2, m2, s1, s2) * RSphCA(s2, m1, m2) *
            OneCenterP(s1, n1, n2, ζ1, ζ2, R) * resaab1 * SphericalHarmonicsS(s1, s2, theta, phi)
            )
        end
    end

    return resaab//2
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
################################# TWO CENTER KINETIC ENERGY INTEGRALS #####################################(43)
################################################### ALIGNED MOLECULAR COORDINATE SYSTEM
###########################################################################################################
###########################################################################################################
###########################################################################################################
################################# TWO CENTER KINETIC ENERGY INTEGRALS #####################################(43)
################################################### ALIGNED MOLECULAR COORDINATE SYSTEM
#function  TwoCenterKinEnergy(
#    n1::arb, l1::Int, ζ1::arb,
#    n2::arb, l2::Int, ζ2::arb,
#    λ::Int, R::arb, lim::Int
#    )

#    ρ = (R // 2) * (ζ1 + ζ2)
#    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

#    coeff1 = -Power(ζ2, 2) // 2
#    coeff2 = 4 * n2 * Sqrt(Gamma(2*n2 - 1) // Gamma(2*n2 + 1))

#    res = coeff1 * (
#        TwoCenterOverlap(n1, l1, ζ1, n2, l2, ζ2, λ, R, lim) -
#        coeff2 * TwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim) +
#     	if 2 * n2 - 3 < 0
#	0
#    else
#        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
#        overlap3 = TwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
#        coeff3 * overlap3
#    end
#    )

#    return res

#end

function  TwoCenterKinEnergy(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Int
    )

    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    coeff1 = -Power(ζ2, 2) // 2
    coeff2 = 4 * n2 * Sqrt(Gamma(2*n2 - 1) // Gamma(2*n2 + 1))

#    if 2 * n2 - 3 < 0 || l2 > n2 - 2
#        coeff3 = 0
#        overlap3 = 0
#    else
#        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
#        overlap3 = TwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
#    end

    if 2 * n2 - 3 < 0
        coeff3 = 0
        overlap3 = 0
    else
        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
        overlap3 = TwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    end

    overlap1 = TwoCenterOverlap(n1, l1, ζ1, n2, l2, ζ2, λ, R, lim)
    overlap2 = TwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #if l2 > n2-1
    #    overlap2 = 0
    #else
    #    overlap2 = TwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #end

    res = coeff1 * (
        overlap1 - coeff2 * overlap2 + coeff3 * overlap3
    )

    return res

end

function  CTwoCenterKinEnergy(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
    )

    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    coeff1 = -Power(ζ2, 2) // 2
    coeff2 = 4 * n2 * Sqrt(Gamma(2*n2 - 1) // Gamma(2*n2 + 1))

#    if 2 * n2 - 3 < 0 || l2 > n2 - 2
#        coeff3 = 0
#        overlap3 = 0
#    else
#        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
#        overlap3 = CTwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
#    end

    if 2 * n2 - 3 < 0
        coeff3 = 0
        overlap3 = 0
    else
        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
        overlap3 = CTwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    end

    overlap1 = CTwoCenterOverlap(n1, l1, ζ1, n2, l2, ζ2, λ, R, lim)
    overlap2 = CTwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #if l2 > n2-1
    #    overlap2 = 0
    #else
    #    overlap2 = CTwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #end

    res = coeff1 * (
        overlap1 - coeff2 * overlap2 + coeff3 * overlap3
    )

    return res

end

function  VTwoCenterKinEnergy(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
    )

    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    coeff1 = -Power(ζ2, 2) // 2
    coeff2 = 4 * n2 * Sqrt(Gamma(2*n2 - 1) // Gamma(2*n2 + 1))

    #    if 2 * n2 - 3 < 0 || l2 > n2 - 2
    #        coeff3 = 0
    #        overlap3 = 0
    #    else
    #        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
    #        overlap3 = VTwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    #    end

    if 2 * n2 - 3 < 0
        coeff3 = 0
        overlap3 = 0
    else
        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
        overlap3 = VTwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    end

    overlap1 = VTwoCenterOverlap(n1, l1, ζ1, n2, l2, ζ2, λ, R, lim)
    overlap2 = VTwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)

    #if l2 > n2-1
    #    overlap2 = 0
    #else
    #    overlap2 = VTwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #end

    res = coeff1 * (
        overlap1 - coeff2 * overlap2 + coeff3 * overlap3
    )

    return res

end

function  STwoCenterKinEnergy(
    n1::arb, l1::Int, ζ1::arb,
    n2::arb, l2::Int, ζ2::arb,
    λ::Int, R::arb, lim::Float64
    )

    ρ = (R // 2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    coeff1 = -Power(ζ2, 2) // 2
    coeff2 = 4 * n2 * Sqrt(Gamma(2*n2 - 1) // Gamma(2*n2 + 1))

    #    if 2 * n2 - 3 < 0 || l2 > n2 - 2
    #        coeff3 = 0
    #        overlap3 = 0
    #    else
    #        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
    #        overlap3 = STwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    #    end

    if 2 * n2 - 3 < 0
        coeff3 = 0
        overlap3 = 0
    else
        coeff3 = 4* (n2 + l2) * (n2 - l2 - 1) * Sqrt(Gamma(2*n2 - 3) // Gamma(2*n2 + 1))
        overlap3 = STwoCenterOverlap(n1, l1, ζ1, n2 - 2, l2, ζ2, λ, R, lim)
    end

    overlap1 = STwoCenterOverlap(n1, l1, ζ1, n2, l2, ζ2, λ, R, lim)
    overlap2 = STwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)

    #if l2 > n2-1
    #    overlap2 = 0
    #else
    #    overlap2 = STwoCenterOverlap(n1, l1, ζ1, n2 - 1, l2, ζ2, λ, R, lim)
    #end

    res = coeff1 * (
        overlap1 - coeff2 * overlap2 + coeff3 * overlap3
    )

    return res

end
###########################################################################################################
###########################################################################################################
###########################################################################################################
################################################### NON-ALIGNED MOLECULAR COORDINATE SYSTEM
function  TwoCenterKinEnergy(
    n1::arb, l1::Int, m1:: Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, x::Int, lim::Int
    )

    NKinEnergy = RF(0)

    for s in 0 : min(l1,l2)
        NKinEnergy += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        TwoCenterKinEnergy(n1, l1, ζ1, n2, l2, ζ2, s, R, lim))
    end
    NKinEnergy
end

function  CTwoCenterKinEnergy(
    n1::arb, l1::Int, m1:: Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, x::Int, lim::Float64
    )

    NKinEnergy = RF(0)

    for s in 0 : min(l1,l2)
        NKinEnergy += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        CTwoCenterKinEnergy(n1, l1, ζ1, n2, l2, ζ2, s, R, lim))
    end
    NKinEnergy
end

function  VTwoCenterKinEnergy(
    n1::arb, l1::Int, m1:: Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, x::Int, lim::Float64
    )

    NKinEnergy = RF(0)

    for s in 0 : min(l1,l2)
        NKinEnergy += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        VTwoCenterKinEnergy(n1, l1, ζ1, n2, l2, ζ2, s, R, lim))
    end
    NKinEnergy
end

function  STwoCenterKinEnergy(
    n1::arb, l1::Int, m1:: Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    R::arb, theta::arb, phi::arb, x::Int, lim::Float64
    )

    NKinEnergy = RF(0)

    for s in 0 : min(l1,l2)
        NKinEnergy += (RotaT(s, l1, m1, l2, m2, theta, phi, x) *
        STwoCenterKinEnergy(n1, l1, ζ1, n2, l2, ζ2, s, R, lim))
    end
    NKinEnergy
end


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
#using Revise
using ZigZagBoomerang
using DataStructures
using LinearAlgebra
using Random
using SparseArrays
using Test
using FileIO
using Statistics
using Makie, AbstractPlotting
using ForwardDiff
using StructArrays
const ρ0 = 0.0
Σ = [1 ρ0; ρ0 1]
const d = 2
const dist = 1.5

#phi(x, y, rho) =  1/(2*pi*sqrt(1-rho^2))*exp(-0.5*(x^2 + y^2 - 2x*y*rho)/(1-rho^2))
#logdensity(x, y) = log(phi(x - dist, y - dist, 0.0) + phi(x + dist, y + dist, 0.0))

logdensity(x, y, sigma=1/sqrt(2)) = -2log(2sigma) - abs(x)/sigma - abs(y)/sigma

r = -2.5:0.01:2.5

surf = surface(r, r, [3*exp(logdensity(x1, x2)) for x1 in r, x2 in r], show_axis=false)
#=
function partiali()
    ith = zeros(d)
    function (x,i)
        ith[i] = 1
        sa = StructArray{ForwardDiff.Dual{}}((x, ith))
        δ = density(sa).partials[]
        ith[i] = 0
        return δ
    end
end
=#
function gradϕ(x,i)
    push!(evals, copy(x))
    -ForwardDiff.partials(logdensity(ForwardDiff.Dual{}(x[1], 1.0*(i==1)), ForwardDiff.Dual{}(x[2], 1.0*(i==2))))[]
end

function gradϕ!(y, x)
    y[1] = gradϕ(x,1)
    y[2] = gradϕ(x,2)
    y
end

#κ = [0.1, 0.1]
κ = [100, 100]

t0 = 0.0
x0 = randn(d)
evals = [x0]

θ0 = [1.0, 1.0]
c = 20.0
T = 1000.0

trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ, t0, x0, θ0, T, c*ones(d), ZigZag(sparse(I(d)), 0*x0), κ; structured=true, adapt=false)
#trace, (tT, xT, θT), (acc, num) = spdmp(gradϕ, t0, x0, θ0, T, c*ones(d), ZigZag(sparse(I(d)), 0*x0); adapt=false)
ts, xs = ZigZagBoomerang.sep(collect(ZigZagBoomerang.discretize(trace, 0.05)))
p1scatter = scatter(first.(xs), last.(xs), color=ts, markersize = 3, strokewidth = 0.1, alpha = 0.5)

ts, xs = ZigZagBoomerang.sep(collect(trace))

p1 = lines(first.(xs), last.(xs), color=ts,  alpha = 0.2)
scatter!(first.(evals), last.(evals), color=ts, markersize = 3, strokewidth = 0.1, alpha = 0.5)
error("stop")
scene, layout = layoutscene(resolution = (1200, 900))
layout[1, 1] = ax1 = Axis(scene)
layout[2, 1] = ax2 = Axis(scene)

linkyaxes!(ax1, ax2)
linkxaxes!(ax1, ax2)

lines!(ax1, ts, first.(xs))
lines!(ax2, ts, last.(xs))
p1a = scene

c = 20.0

evals = [x0]
trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ!, t0, x0, θ0, T, c, BouncyParticle(sparse(I(d)), 0*x0, 0.1), κ; adapt=false)

#ts, xs = ZigZagBoomerang.sep(collect(discretize(trace, 0.05)))
ts, xs = ZigZagBoomerang.sep(collect(trace))

p2 = lines(first.(xs), last.(xs), color=ts)

scene, layout = layoutscene(resolution = (1200, 900))
layout[1, 1] = ax1 = Axis(scene)
layout[2, 1] = ax2 = Axis(scene)

linkyaxes!(ax1, ax2)
linkxaxes!(ax1, ax2)

lines!(ax1, ts, first.(xs))
lines!(ax2, ts, last.(xs))
p2a = scene

Random.seed!(2)
c = 2.0
evals = [x0]
BB = Boomerang(sparse(I(d)), 0*x0, 0.1, 0.9)
trace, (tT, xT, θT), (acc, num) = sspdmp(gradϕ!, t0, x0, θ0, T, c, BB, κ; structured=true, adapt=false)

tsxs = sdiscretise(trace.events, BB, 0.05)

ts, xs= tsxs.t, tsxs.x

p3 = scatter(first.(xs), last.(xs), color=ts, markersize = 1, strokewidth = 0.1, alpha = 0.5)
p3c = scatter(first.(xs), last.(xs), color=ts, strokewidth = 0.05, markersize = 7, alpha = 0.5)
scatter!(first.(evals), last.(evals), color=:black, markersize = 10, marker='x', strokewidth = 0.1, alpha = 0.5)
p3c

scene, layout = layoutscene(resolution = (1200, 900))
    layout[1, 1] = ax1 = Axis(scene)
    layout[2, 1] = ax2 = Axis(scene)

    linkyaxes!(ax1, ax2)
    linkxaxes!(ax1, ax2)

    lines!(ax1, ts, first.(xs))
    lines!(ax2, ts, last.(xs))
    p3a = scene

save("zigzagphase.png", p1)
save("zigzagtrace.png", p1a)

save("bouncyphase.png", p2)
save("bouncytrace.png", p2a)

save("boomyphase.png", p3)
save("boomytrace.png", p3a)

ext: jl


text:
using FileManager
using Visualization
using Common
using OrthographicProjection


text:
demand = [12, 14, 16, 18, 20, 22, 24, 26, 28, 30] # Demand of newspapers in each scenario
lengthS = length(demand)
vec_prob = [0.05, 0.10, 0.10, 0.10, 0.15, 0.15, 0.10, 0.10, 0.10, 0.05 ] # probability of scenario
c = 20 # purchase price
p = 70 # selling price
h = 10 # scrap value
y = 20
# 1. Solved by L-Shaped Benders Algorithm ------------------------------------------------------------------------------
n_x = 10
n_y = 1
vec_min_y = hcat([0])
vec_max_y = hcat([30])
vec_f = hcat(c - 10)
vec_pi = hcat(vec_prob)
mat_c = zeros(10, 1, 1)
for i = 1:10
    mat_c[i, :, :] = hcat([- 70 + 10])
end
mat_h = zeros(10, 2, 1)
for i = 1: 10
    mat_h[i, 1, :] = hcat(- demand[i])
    mat_h[i, 2, :] = hcat(0)
end
mat3_t = zeros(10, 2, 1)
for i = 1:10
    mat3_t[i, 1, :] = hcat(0)
    mat3_t[i, 2, :] = hcat(1)
end
mat3_w = - ones(10, 2, 1)
epsilon = 0.001
timesIterationMax = 500


text:
"""
Use git to clone the meshcat javascript assets repository for local development.

$(SIGNATURES)

You should only do this if you plan on editing the javascript or HTML components of
meshcat itself. To undo this operation, you will need to delete the `assets/meshcat`
folder and then run `using Pkg; Pkg.build("MeshCat")`
"""
function develop_meshcat_assets(skip_confirmation=false)
    meshcat_dir = abspath(joinpath(@__DIR__, "..", "assets", "meshcat"))
    if !skip_confirmation
        println("CAUTION: This will delete all downloaded meshcat assets and replace them with a git clone.")
        println("The following path will be overwritten:")
        println(meshcat_dir)
        println("To undo this operation, you will need to manually remove that directory and then run `Pkg.build(\"MeshCat\")`")
        print("Proceed? (y/n) ")
        choice = chomp(readline())
        if isempty(choice) || lowercase(choice[1]) != 'y'
            println("Canceled.")
            return
        end
    end
    println("Removing $meshcat_dir")
    rm(meshcat_dir, force=true, recursive=true)
    run(`git clone https://github.com/rdeits/meshcat $meshcat_dir`)
    rm(joinpath(meshcat_dir, "..", "meshcat.stamp"))
end


text:
# Givens representation of Q matrix, memory efficient compared to dense representation
struct GivensQ{T, QT<:AbstractVector{<:Givens{T}}} <: Factorization{T}
    rotations::QT
    n::Int
    m::Int
end

# r controls how many columns of Q we instantiate
# NOTE: this also works if rotations is empty (returns identity matrix)
function Base.Matrix(F::GivensQ, r::Int = size(F, 2))
    n, m = size(F)
    Q = Matrix(one(eltype(F))*I, n, r)
    return lmul!(F, Q)
end
Base.size(F::GivensQ) = (F.n, F.n) # could be square but also (n, m)?
Base.size(F::GivensQ, i::Int) = i > 2 ? 1 : size(F)[i]
Base.copy(F::GivensQ) = GivensQ(copy(F.rotations), F.n, F.m)
Base.adjoint(F::GivensQ) = Adjoint(F)
# NOTE: this also works if rotations is empty (is identity operator)
function LinearAlgebra.lmul!(F::GivensQ, X::AbstractVecOrMat)
    m = size(X, 1)
    for G in Iterators.reverse(F.rotations) # the adjoint of a GivensQ reverses the order of the rotations
        lmul!(G', X) # ... and takes their adjoint
    end
    return X
end
function LinearAlgebra.lmul!(A::Adjoint{<:Any, <:GivensQ}, X::AbstractVecOrMat)
    F = A.parent
    for G in F.rotations # the adjoint of a GivensQ reverses the order of the rotations
        lmul!(G, X) # ... and takes their adjoint
    end
    return X
end
LinearAlgebra.rmul!(X::AbstractVecOrMat, F::GivensQ) = lmul!(F', X')'
LinearAlgebra.rmul!(X::AbstractVecOrMat, F::Adjoint{<:Any, <:GivensQ}) = lmul!(F', X')'

lmul(F, X) = lmul!(F, copy(X))
rmul(X, F) = rmul!(copy(X), F)

Base.:*(F::GivensQ, X::AbstractVector) = lmul(F, X)
Base.:*(F::GivensQ, X::AbstractMatrix) = lmul(F, X)
Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractVector) = lmul(F, X)
Base.:*(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul(F, X)

Base.:*(X::AbstractMatrix, F::GivensQ) = rmul(X, F)
Base.:*(X::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}) = rmul(X, F)

function Base.:*(F::GivensQ, G::Adjoint{<:Any, <:GivensQ})
    T = promote_type(eltype(F), eltype(G))
    F === G' ? (one(T) * I)(F.n) : F * Matrix(G)
end
function Base.:*(F::Adjoint{<:Any, <:GivensQ}, G::GivensQ)
    T = promote_type(eltype(F), eltype(G))
    F' === G ? (one(T) * I)(G.n) : F * Matrix(G)
end
function Base.:(==)(F::GivensQ, G::GivensQ)
    F.rotations == G.rotations && F.n == G.n && F.m == G.m
end

LinearAlgebra.ldiv!(F::GivensQ, x::AbstractVector) = lmul!(F', x)
LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector) = lmul!(F', x)
LinearAlgebra.ldiv!(F::GivensQ, X::AbstractMatrix) = lmul!(F', X)
LinearAlgebra.ldiv!(F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix) = lmul!(F', X)

function LinearAlgebra.mul!(y::AbstractVector, F::GivensQ, x::AbstractVector)
    @. y = x
    lmul!(F, y)
end
function LinearAlgebra.mul!(Y::AbstractMatrix, F::GivensQ, X::AbstractMatrix)
    @. Y = X
    lmul!(F, Y)
end
function LinearAlgebra.mul!(y::AbstractVector, F::Adjoint{<:Any, <:GivensQ}, x::AbstractVector)
    @. y = x
    lmul!(F, y)
end
function LinearAlgebra.mul!(Y::AbstractMatrix, F::Adjoint{<:Any, <:GivensQ}, X::AbstractMatrix)
    @. Y = X
    lmul!(F, Y)
end


text:
include("preprocessing.jl")
include("postprocessing.jl")
include("inprocessing.jl")


text:
using LTC
using ModelingToolkit
using LTC: rand_uniform

@register NNlib.sigmoid(t)

function SigmoidSynapse(;name)
  vars = @variables I(t), v_pre(t), v_post(t)
  ps = @parameters begin
    μ = rand_uniform(Float32, 0.3, 0.8), [lower=0.1f0, upper=1.0f0]
    σ = rand_uniform(Float32, 3, 8), [lower=1f0, upper=10f0]
    G = rand_uniform(Float32, 0.001, 1), [lower=0.000001f0, upper=1.1f0]
    E = rand_uniform(Float32, -0.3, 0.3), [lower=-1.1f0, upper=1.1f0]
  end
  eqs = [
    I ~ G * NNlib.sigmoid((v_pre - μ) * σ) * (v_post - E)
  ]
  ODESystem(eqs, t, vars, ps; name)
end

function LeakChannel(;name)
  vars = @variables I(t), v(t)
  ps = @parameters begin
    G = rand_uniform(Float32, 0.001, 1), [lower = 0.000001f0, upper = 1.1f0]
    E = rand_uniform(Float32, -0.3, 0.3), [lower = -1.1f0, upper = 1.1f0]
  end
  eqs = [
    I ~ G * (v - E)
  ]
  ODESystem(eqs, t, vars, ps; name)
end

function Neuron(; name)
  @variables begin
    (v(t) = rand_uniform(Float32, 0.001, 0.2)), [lower=-2f0, upper=2f0]
    (I_comps(t))
  end
  ps = @parameters begin
    Cm = rand_uniform(Float32, 1, 3), [lower = 0.8f0, upper = 4f0]
  end
  @named leak = LeakChannel()
  eqs = [
    D(v) ~ -Cm * (leak.I + I_comps)
    leak.v ~ v
  ]
  systems = ODESystem[leak]
  ODESystem(eqs, t, [v, I_comps], ps; systems, name)
end

function Net(wiring; name)
  vars = Num[]
  ps = Num[]
  eqs = Equation[]
  systems = ODESystem[]

  N = wiring.n_total
  inputs, outputs = create_pins(wiring.n_in, wiring.n_out)
  push!(systems, inputs...)

  n = 1
  neurons = ODESystem[]
  for i in 1:wiring.n_sensory
    push!(neurons, Neuron(;name=Symbol("n$(n)_SensoryNeuron")))
    n += 1
  end
  for i in 1:wiring.n_inter
    push!(neurons, Neuron(;name=Symbol("n$(n)_InterNeuron")))
    n += 1
  end
  for i in 1:wiring.n_command
    push!(neurons, Neuron(;name=Symbol("n$(n)_CommandNeuron")))
    n += 1
  end
  for i in 1:wiring.n_motor
    push!(neurons, Neuron(;name=Symbol("n$(n)_MotorNeuron")))
    n += 1
  end
  push!(systems, neurons...)

  n = 1
  for dst in 1:length(neurons)
    I_comps_dst = 0
    for src in 1:length(wiring.sens_mask[:,dst])
      wiring.sens_mask[src,dst] == 0 && continue
      syn = SigmoidSynapse(;name=Symbol("s$(n)_x$(src)-->n$(dst)_SigmoidSynapse"))
      n += 1
      push!(eqs, syn.v_pre ~ inputs[src].x)
      push!(eqs, syn.v_post ~ neurons[dst].v)
      push!(systems, syn)
      I_comps_dst += syn.I
    end
    for src in 1:length(wiring.syn_mask[:,dst])
      wiring.syn_mask[src,dst] == 0 && continue
      syn = SigmoidSynapse(;name=Symbol("s$(n)_n$(src)-->n$(dst)_SigmoidSynapse"))
      n += 1
      push!(eqs, syn.v_pre ~ neurons[src].v)
      push!(eqs, syn.v_post ~ neurons[dst].v)
      push!(systems, syn)
      I_comps_dst += syn.I
    end
    push!(eqs, neurons[dst].I_comps ~ I_comps_dst)
  end

  # connect output pins with n_out last neurons
  for i in 1:wiring.n_out
    push!(eqs, outputs[i].x ~ neurons[end-wiring.n_out+i].v)
  end
  push!(systems, outputs...)

  ODESystem(eqs, t, vars, ps; systems, name=name)
end


text:
using Distributed: RemoteChannel

"""
    Connector{T, S}

A channel-like interface that accepts objects of type `T` (through `put!`) and
may provide objects of type `S` (through `take!`).
"""
struct Connector{T, S}
    a::RemoteChannel{Channel{T}}
    b::RemoteChannel{Channel{S}}
end

Connector(;proc = 1) = Connector{Any, Any}(;proc)
Connector{T, S}(;proc = 1) where {T, S} =
    Connector(
        RemoteChannel(() -> Channel{T}(0), proc),
        RemoteChannel(() -> Channel{S}(0), proc))

Base.show(io::IO, ::MIME"text/plain", con::Connector{T,S}) where {T,S} =
    print(io, "Connector: $T to $S, ", isopen(con) ? "open" : "closed")

# channel-like interface
Base.put!(con::Connector{T}, x::T) where T = (put!(con.a, x); con)
Base.take!(con::Connector) = take!(con.b)
Base.close(con::Connector) = (close(con.a); close(con.b); nothing)
Base.isopen(con::Connector) = isopen(con.a) && isopen(con.b)

"""
    reverse(::Connector)

Obtain the other matching end of the connector pair.
"""
Base.reverse(con::Connector) = Connector(con.b, con.a)


text:
using HTTP, JSON2, Sockets
using Dates
const ROUTER = HTTP.Router()

include("lib/Covid19Modelling.jl")

using .Covid19Modelling

include("models.jl")
include("router.jl")

println("[log] precompilation complete, starting server at port 8081")
HTTP.serve(ROUTER, "0.0.0.0", 8081)


text:
function action_space_size(act_space::ClosedInterval)
    return length(leftendpoint(act_space))
end
function action_space_size(act_space::Base.OneTo)
    return 1
end

function block_diagm(A::Vector{Float64}, rep_number::Int)
    return diagm(repeat(A, rep_number))
end

function block_diagm(A::Matrix{Float64}, rep_number::Int)
    r = size(A)[1]
    rm1 = r -1
    B = zeros(Float64, r*rep_number, r*rep_number)
    for ii in 1:r:(r*rep_number)
        B[ii:ii+rm1, ii:ii+rm1] = A
    end
    return B
end

function get_model_controls(action_space::ClosedInterval, V::Vector{Float64})
    return get_model_controls(action_space, V, 1)
end
function get_model_controls(action_space::Base.OneTo, V::Vector{Float64})
    return get_model_controls(action_space, V, 1)
end
function get_model_controls(action_space::ClosedInterval, V::Vector{Float64}, horizon::Int)
    as = action_space_size(action_space)
    min_controls = leftendpoint(action_space)
    max_controls = rightendpoint(action_space)
    control_mat = reshape(V, as, horizon)
    for rᵢ ∈ 1:size(control_mat)[1]
        control_mat[rᵢ,:] = clamp.(control_mat[rᵢ,:], min_controls[rᵢ], max_controls[rᵢ])
    end
    if as == 1
        control_mat = vec(control_mat)
    end
    return control_mat
end
function get_model_controls(action_space::Base.OneTo, V::Vector{Float64}, horizon::Int)
    as = 1
    min_controls = minimum(action_space)
    max_controls = maximum(action_space)

    control_mat = reshape(V, as, horizon)
    control_mat[1,:] = clamp.(control_mat[1,:], min_controls, max_controls)
    control_mat = round.(Int, vec(control_mat))
    return control_mat
end

function compute_weights(weight_method::Information_Theoretic, costs::Vector{Float64})
    λ = weight_method.λ
    ρ = minimum(costs)
    normalized_costs = -1/λ *(costs .- ρ)
    weights = exp.(normalized_costs)
    η = sum(weights)
    return weights ./ η
end

function get_controls_roll_U!(pol::AbstractPathIntegralPolicy, weighted_controls::Vector)
    as = pol.params.as
    # Get control (action set for the first time step)
    control = get_model_controls(action_space(pol.env), weighted_controls[1:as])

    # Roll the control policy so next interation we start with a mean of pol.U
    if pol.params.horizon > 1
        pol.U[1:(end-as)] = weighted_controls[(as+1):end]
        pol.U[(end-as):end] = pol.params.U₀[(end-as):end]
    else
        pol.U = weighted_controls
    end
    return control
end

function rollout_model(env::AbstractEnv, T::Int, model_controls::Vector,
    pol::AbstractPathIntegralPolicy, k::Int, n::Int)
    model_controls_mat = reshape(model_controls, size(model_controls, 1), 1)
    rollout_model(env, T, model_controls_mat, pol, k, n)
end

function rollout_model(env::AbstractEnv, T::Int, model_controls::Matrix,
    pol::AbstractPathIntegralPolicy, k::Int, n::Int,
)
    as = pol.params.as
    K = pol.params.num_samples
    traj_cost = 0.0
    for t ∈ 1:T
        controls = as == 1 ? model_controls[t] : model_controls[:,t]
        env(controls)
        traj_cost -= reward(env) # Subtracting based on "reward"
        if pol.params.log
            pol.logger.trajectories[k+(n-1)*K][t, :] = env.state
        end
    end
    return traj_cost
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
Base.isnull(x) = false
Base.get(x) = x

macro lift(call, T)
    esc(_lift(call, T))
end

function _lift(call, T)
    arg_cache = Dict{Union{Symbol, Expr}, Expr}()
    e_call = gen_calls(call, arg_cache, true)
    args = collect(keys(arg_cache))
    e_nullcheck = :(isnull($(args[1])))
    for i = 2:length(args)
        e_nullcheck = Expr(:||, e_nullcheck, :(isnull($(args[i]))))
    end

    res = quote
        if $e_nullcheck
            Nullable{$T}()
        else
            $e_call
        end
    end
    res
end

# base case for literals
gen_calls(e, arg_cache, top) = e

# base case for symbols
function gen_calls(e::Symbol, arg_cache, top)
    new_arg = get!(arg_cache, e, :(get($e)))
    return new_arg
end

function gen_calls(e::Expr, arg_cache, top)
    if e.head == :call
        if top == true
            res = Expr(:call, e.args[1], gen_calls(e.args[2:end], arg_cache, false)...)
            return :(Nullable($res))
        else
            return Expr(:call, e.args[1], gen_calls(e.args[2:end], arg_cache, false)...)
        end
    elseif e.head == :ref
        new_arg = get!(arg_cache, e, :(get($e)))
        return new_arg
    # `if` support
    elseif e.head == :if
        return Expr(:if, [ gen_calls(expr, arg_cache, top) for expr in e.args ]...)
    elseif e.head == :block
        return Expr(:block, [ gen_calls(expr, arg_cache, top) for expr in e.args ]...)
    elseif e.head == :comparison
        return Expr(:comparison, gen_calls(e.args[1], arg_cache, top), e.args[2], gen_calls(e.args[3], arg_cache, top))
    # /if support
    else
        return e
    end
end

function gen_calls(args::Array, arg_cache, top)
    return [ gen_calls(arg, arg_cache, top) for arg in args ]
end

max_stars_repo_name: davidagold/NullableUtils

max_issues_repo_name: davidagold/NullableUtils

max_forks_repo_name: davidagold/NullableUtils


text:
@testset "ConditionalCategorical" begin
    c = 3
    x = randn(2)
    y = 1
    xb = randn(2, 5)
    yb = [1,2,3,2,3]

    m = Chain(Dense(2, 3, swish), Dense(3,c), softmax)
    p = ConditionalCategorical(m)

    # test that the distribution is differentiable
    # and parameters are okay
    ps = Flux.params(p)
    @test !isempty(ps)
    @test length(ps) == length(Flux.params(m))
    fl() = logpdf(condition(p, x), y)
    @test_nowarn Flux.gradient(fl, ps)

    # test for vectors
    res = condition(p, x)
    @test res isa Categorical
    @test size(probs(res)) == (c, )

    # test for batches
    res = condition(p, xb)
    @test res isa ConditionalDists.BatchCategorical
    @test size(probs(res)) == (c, size(xb, 2))

    # test pdf, logpdf and entropy for single
    f(x) = [0.2,0.3,0.5]
    p = ConditionalCategorical(f)
    res = condition(p, x)

    @test entropy(res) == entropy([0.2,0.3,0.5])
    @test pdf.(res, [1,2,3]) == [0.2,0.3,0.5]
    @test logpdf.(res, [1,2,3]) == log.([0.2,0.3,0.5])

    # test pdf, logpdf, entropy for batch
    g(x) = mapreduce(xi -> f(xi), hcat, eachcol(x))
    p = ConditionalCategorical(g)
    res = condition(p, xb)

    @test length(unique(entropy(res))) == 1
    @test unique(entropy(res))[1] == entropy([0.2,0.3,0.5])
    @test length(pdf(res, yb)) == size(xb, 2)
    @test pdf(res, yb) == [0.2,0.3,0.5,0.3,0.5]
    @test logpdf(res, yb) == log.([0.2,0.3,0.5,0.3,0.5])

end


text:
using JuCUDA
if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

files = ["array", "arraymath"]
const T = Float32

for f in files
    include("$(f).jl")
end


text:
export @dual, @selfdual, @dualtype

"""
    @dual f invf

Define `f` and `invf` as a pair of dual instructions, i.e. reverse to each other.
"""
macro dual(f, invf)
    esc(quote
        if !$NiLangCore.isprimitive($f)
            $NiLangCore.isprimitive(::typeof($f)) = true
        end
        if !$NiLangCore.isprimitive($invf)
            $NiLangCore.isprimitive(::typeof($invf)) = true
        end
        if Base.:~($f) !== $invf
            Base.:~(::typeof($f)) = $invf;
        end
        if Base.:~($invf) !== $f
            Base.:~(::typeof($invf)) = $f;
        end
    end)
end

macro dualtype(t, invt)
    esc(quote
        $invtype($t) === $invt || begin
            $NiLangCore.invtype(::Type{$t}) = $invt
            $NiLangCore.invtype(::Type{T}) where T<:$t = $invt{T.parameters...}
        end
        $invtype($invt) === $t || begin
            $NiLangCore.invtype(::Type{$invt}) = $t
            $NiLangCore.invtype(::Type{T}) where T<:$invt = $t{T.parameters...}
        end
    end)
end
@dualtype PlusEq MinusEq
@dualtype DivEq MulEq
@dualtype XorEq XorEq

"""
    @selfdual f

Define `f` as a self-dual instructions.
"""
macro selfdual(f)
    esc(:(@dual $f $f))
end

export @const
@eval macro $(:const)(ex)
    esc(ex)
end

export @skip!
macro skip!(ex)
    esc(ex)
end

export @assignback

# TODO: include control flows.
"""
    @assignback f(args...) [invcheck]

Assign input variables with output values: `args... = f(args...)`, turn off invertibility error check if the second argument is false.
"""
macro assignback(ex, invcheck=true)
    ex = precom_ex(__module__, ex, PreInfo())
    esc(assignback_ex(ex, invcheck))
end

function assignback_ex(ex::Expr, invcheck::Bool)
    @smatch ex begin
        :($f($(args...))) => begin
            symres = gensym("results")
            ex = :($symres = $f($(args...)))
            res = assign_vars(seperate_kwargs(args)[1], symres, invcheck)
            pushfirst!(res.args, ex)
            return res
        end
        _ => error("assign back fail, got $ex")
    end
end

"""
    assign_vars(args, symres, invcheck)

Get the expression of assigning `symres` to `args`.
"""
function assign_vars(args, symres, invcheck)
    exprs = []
    for (i,arg) in enumerate(args)
        exi = @smatch arg begin
            :($ag...) => begin
                i!=length(args) && error("`args...` like arguments should only appear as the last argument!")
                assign_ex(ag, :($tailn($symres, Val($i-1), $ag)), invcheck)
            end
            _ => if length(args) == 1
                assign_ex(arg, symres, invcheck)
            else
                assign_ex(arg, :($symres[$i]), invcheck)
            end
        end
        exi !== nothing && push!(exprs, exi)
    end
    Expr(:block, exprs...)
end

error_message_fcall(arg) = """
function arguments should not contain function calls on variables, got `$arg`, try to decompose it into elementary statements, e.g. statement `z += f(g(x))` should be written as
    y += g(x)
    z += y

If `g` is a dataview (a function map an object to its field or a bijective function), one can also use the pipline like
    z += f(x |> g)
"""

assign_ex(arg, res, invcheck) = @smatch arg begin
    ::Number || ::String => _invcheck(invcheck, arg, res)
    ::Symbol || ::GlobalRef => _isconst(arg) ? _invcheck(invcheck, arg, res) : :($arg = $res)
    :(@skip! $line $x) => nothing
    :(@fields $line $x) => assign_ex(x, Expr(:call, default_constructor, :(typeof($x)), Expr(:..., res)), invcheck)
    :($x::∅) => assign_ex(x, res, invcheck)
    :($x::$T) => assign_ex(x, :($loaddata($T, $res)), invcheck)
    :($x.$k) => _isconst(x) ? _invcheck(invcheck, arg, res) : assign_ex(x, :(chfield($x, $(Val(k)), $res)), invcheck)
    # tuples must be index through (x |> 1)
    :($a |> tget($x)) => assign_ex(a, :($(TupleTools.insertat)($a, $x, ($res,))), invcheck)
    :($a |> subarray($(ranges...))) => :(($res===view($a, $(ranges...))) || (view($a, $(ranges...)) .= $res))
    :($x |> $f) => _isconst(x) ? _invcheck(invcheck, arg,res) : assign_ex(x, :(chfield($x, $f, $res)), invcheck)
    :($x .|> $f) => _isconst(x) ? _invcheck(invcheck, arg,res) : assign_ex(x, :(chfield.($x, Ref($f), $res)), invcheck)
    :($x') => _isconst(x) ? _invcheck(invcheck, arg, res) : assign_ex(x, :(chfield($x, adjoint, $res)), invcheck)
    :(-$x) => _isconst(x) ? _invcheck(invcheck, arg,res) : assign_ex(x, :(chfield($x, -, $res)), invcheck)
    :($t{$(p...)}($(args...))) => begin
        if length(args) == 1
            assign_ex(args[1], :($getfield($res, 1)), invcheck)
        else
            assign_vars(args, :($type2tuple($res)), invcheck)
        end
    end
    :($f($(args...))) => all(_isconst, args) || error(error_message_fcall(arg))
    :($f.($(args...))) => all(_isconst, args) || error(error_message_fcall(arg))
    :($a[$(x...)]) => begin
        :($a[$(x...)] = $res)
    end
    :(($(args...),)) => begin
        # TODO: avoid possible repeated evaluation (not here, in swap)
        Expr(:block, [assign_ex(args[i], :($res[$i]), invcheck) for i=1:length(args)]...)
    end
    _ => _invcheck(invcheck, arg, res)
end

# general
@inline tailn(t::Tuple, ::Val{n}) where n = tailn(TupleTools.tail(t), Val{n-1}())
@inline tailn(t::Tuple, ::Val{n}, input::Tuple) where n = tailn(t, Val{n}())
tailn(t::Tuple, ::Val{0}) = t
# single argument, the second field can be 0, 1
tailn(t, ::Val{0}, input::NTuple{1}) = (t,)
tailn(t, ::Val{1}, input::NTuple{0}) = ()
# avoid ambiguity errors
tailn(t::Tuple, ::Val{0}, input::NTuple{1}) = (t,)
tailn(t::Tuple, ::Val{1}, input::NTuple{0}) = ()

export @assign

"""
    @assign a b [invcheck]

Perform the assign `a = b` in a reversible program.
Turn off invertibility check if the `invcheck` is false.
"""
macro assign(a, b, invcheck=true)
    esc(assign_ex(a, b, invcheck))
end


text:
for x=1:100
    if mod(x,3) == 0 && mod(x,5) == 0
        println("FizzBuzz")
    elseif mod(x,3) == 0
        println("Fizz")
    elseif mod(x,5) == 0
        println("Buzz")
    else
        println(x)
    end
end

max_stars_repo_name: runandrerun/fizzbuzz

max_issues_repo_name: runandrerun/fizzbuzz

max_forks_repo_name: runandrerun/fizzbuzz


text:
module Jax

using PyCall
using ChainRulesCore
using DLPack
using Adapt

using ..PyCallChainRules: PyAdaptor, fmap

const inspect = PyNULL()
const jax = PyNULL()
const dlpack = PyNULL()
const stax = PyNULL()
const numpy = PyNULL()

const ispysetup = Ref{Bool}(false)

pyto_dlpack(x) = @pycall dlpack.to_dlpack(x)::PyObject
pyfrom_dlpack(x) = @pycall dlpack.from_dlpack(x)::PyObject

struct JaxFunctionWrapper
    jaxfn::PyObject
end

function (wrap::JaxFunctionWrapper)(args...; kwargs...)
    out = wrap.jaxfn(fmap(x->DLPack.share(x, PyObject, pyfrom_dlpack), args)...)
    return fmap(x->DLPack.wrap(x, pyto_dlpack), out)
end

function ChainRulesCore.rrule(wrap::JaxFunctionWrapper, args...; kwargs...)
    T = typeof(first(args))
    project = ProjectTo(args)
    pyargs = fmap(x->DLPack.share(x, PyObject, pyfrom_dlpack), args)
    jax_primal, jax_vjpfun = jax.vjp(wrap.jaxfn, pyargs...; kwargs...)
    function JaxFunctionWrapper_pullback(Δ)
        cΔ = fmap(x->Adapt.adapt(PyAdaptor{T}(), x), Δ)
        dlΔ = fmap(x->DLPack.share(x, PyObject, pyfrom_dlpack), cΔ)
        tangent_vals = fmap(x->DLPack.wrap(x, pyto_dlpack), jax_vjpfun(dlΔ))
        return (NoTangent(), project(tangent_vals)...)
    end
    return fmap(x->DLPack.wrap(x, pyto_dlpack), jax_primal), JaxFunctionWrapper_pullback
end

function __init__()
    try
        copy!(jax, pyimport("jax"))
        copy!(dlpack, pyimport("jax.dlpack"))
        copy!(numpy, pyimport("numpy"))
        copy!(stax, pyimport("jax.example_libraries.stax"))
        copy!(inspect, pyimport("inspect"))
        ispysetup[] = true
    catch err
        @warn """PyCallChainRules.jl has failed to import jax from Python.
                 Please make sure these are installed.
        """
        @debug err
        ispysetup[] = false
        #rethrow(err)
    end
end

end


text:
using GeometryBasics, RPRMakie
using Colors, FileIO
using Colors: N0f8

begin
    RPRMakie.activate!(plugin=RPR.Tahoe)
    fig = Figure(; resolution=(1000, 1000))
    ax = LScene(fig[1, 1])
    points = Point3f[]
    for i in 4:10
        n = i + 1
        y = LinRange(0, i, n)
        y2 = (y ./ 2) .- 2
        xyz = Point3f.((i - 5) ./ 2, y2, sin.(y) .+ 1)
        lines!(ax, xyz; linewidth=5, color=:red)
        append!(points, xyz)
    end
    meshscatter!(ax, points, color=:green)
    ax.scene
end


text:
mutable struct JacobianCache{CacheType1,CacheType2,CacheType3,ColorType,SparsityType,fdtype,returntype}
    x1  :: CacheType1
    fx  :: CacheType2
    fx1 :: CacheType3
    colorvec :: ColorType
    sparsity :: SparsityType
end

function JacobianCache(
    x,
    fdtype     :: Type{T1} = Val{:forward},
    returntype :: Type{T2} = eltype(x);
    inplace    :: Type{Val{T3}} = Val{true},
    colorvec = 1:length(x),
    sparsity = nothing) where {T1,T2,T3}

    if eltype(x) <: Real && fdtype==Val{:complex}
        x1  = false .* im .* x
        _fx = false .* im .* x
    else
        x1 = copy(x)
        _fx = copy(x)
    end

    if fdtype==Val{:complex}
        _fx1  = nothing
    else
        _fx1 = copy(x)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype;colorvec=colorvec,sparsity=sparsity)
end

function JacobianCache(
    x ,
    fx,
    fdtype     :: Type{T1} = Val{:forward},
    returntype :: Type{T2} = eltype(x);
    inplace    :: Type{Val{T3}} = Val{true},
    colorvec = 1:length(x),
    sparsity = nothing) where {T1,T2,T3}

    if eltype(x) <: Real && fdtype==Val{:complex}
        x1  = false .* im .* x
    else
        x1 = copy(x)
    end

    if eltype(fx) <: Real && fdtype==Val{:complex}
        _fx = false .* im .* fx
    else
        _fx = copy(fx)
    end

    if fdtype==Val{:complex}
        _fx1  = nothing
    else
        _fx1 = copy(fx)
    end

    JacobianCache(x1,_fx,_fx1,fdtype,returntype;colorvec=colorvec,sparsity=sparsity)
end

function JacobianCache(
    x1 ,
    fx ,
    fx1,
    fdtype     :: Type{T1} = Val{:forward},
    returntype :: Type{T2} = eltype(fx);
    inplace    :: Type{Val{T3}} = Val{true},
    colorvec = 1:length(x1),
    sparsity = nothing) where {T1,T2,T3}

    if fdtype==Val{:complex}
        !(returntype<:Real) && fdtype_error(returntype)

        if eltype(fx) <: Real
            _fx  = false .* im .* fx
        else
            _fx = fx
        end
        if eltype(x1) <: Real
            _x1  = false .* im .* x1
        else
            _x1 = x1
        end
    else
        _x1 = x1
        @assert eltype(fx) == T2
        @assert eltype(fx1) == T2
        _fx = fx
    end
    JacobianCache{typeof(_x1),typeof(_fx),typeof(fx1),typeof(colorvec),typeof(sparsity),fdtype,returntype}(_x1,_fx,fx1,colorvec,sparsity)
end

function _make_Ji(rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
    pick_inds = [i for i in 1:length(rows_index) if colorvec[cols_index[i]] == color_i]
    rows_index_c = rows_index[pick_inds]
    cols_index_c = cols_index[pick_inds]
    len_rows = length(pick_inds)
    unused_rows = setdiff(1:nrows,rows_index_c)
    perm_rows = sortperm(vcat(rows_index_c,unused_rows))
    cols_index_c = vcat(cols_index_c,zeros(Int,nrows-len_rows))[perm_rows]
    Ji = [j==cols_index_c[i] ? dx[i] : false for i in 1:nrows, j in 1:ncols]
    Ji
end

function _make_Ji(xtype, dx, color_i, nrows, ncols)
    Ji = mapreduce(i -> i==color_i ? dx : zeros(xtype,nrows), hcat, 1:ncols)
    size(Ji)!=(nrows, ncols) ? reshape(Ji,(nrows,ncols)) : Ji #branch when size(dx) == (1,) => size(Ji) == (1,) while size(J) == (1,1)
end

function Base.vec(x::Number)
    x
end

function finite_difference_jacobian(f, x::AbstractArray{<:Number},
    fdtype     :: Type{T1}=Val{:forward},
    returntype :: Type{T2}=eltype(x),
    f_in       :: Union{AbstractArray{<:T2},Nothing}=nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = 1:length(x),
    sparsity = nothing,
    jac_prototype = nothing,
    dir=true) where {T1,T2,T3}

    if f_in isa Nothing
        fx = f(x)
    else
        fx = f_in
    end
    cache = JacobianCache(x, fx, fdtype, returntype)
    finite_difference_jacobian(f, x, cache, fx; relstep=relstep, absstep=absstep, colorvec=colorvec, sparsity=sparsity, jac_prototype=jac_prototype, dir=dir)
end

function finite_difference_jacobian(
    f,
    x,
    cache::JacobianCache{T1,T2,T3,cType,sType,fdtype,returntype},
    f_in=nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = cache.colorvec,
    sparsity = cache.sparsity,
    jac_prototype = nothing,
    dir=true) where {T1,T2,T3,cType,sType,fdtype,returntype}

    x1, fx, fx1 = cache.x1, cache.fx, cache.fx1

    if !(f_in isa Nothing)
        vecfx = vec(f_in)
    elseif fdtype == Val{:forward}
        vecfx = vec(f(x))
    elseif fdtype == Val{:complex} && returntype <: Real
        vecfx = real(fx)
    else
        vecfx = vec(fx)
    end
    vecx = vec(x)
    vecx1 = vec(x1)
    J = jac_prototype isa Nothing ? (sparsity isa Nothing ? false.*vecfx.*vecx' : zeros(eltype(x),size(sparsity))) : zero(jac_prototype)
    nrows, ncols = size(J)

    if !(sparsity isa Nothing)
        rows_index, cols_index = ArrayInterface.findstructralnz(sparsity)
        rows_index = [rows_index[i] for i in 1:length(rows_index)]
        cols_index = [cols_index[i] for i in 1:length(cols_index)]
    end

    if fdtype == Val{:forward}
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x_save = vecx[color_i]
                epsilon = compute_epsilon(Val{:forward}, x_save, relstep, absstep, dir)
                _vecx1 = Base.setindex(vecx,x_save+epsilon,color_i)
                _x1 = reshape(_vecx1,size(x))
                vecfx1 = vec(f(_x1))
                dx = (vecfx1-vecfx)/epsilon
                J = J + _make_Ji(eltype(x), dx, color_i, nrows, ncols)
            else
                tmp = norm(vecx .* (colorvec .== color_i))
                epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir)
                _vecx = @. vecx + epsilon * (colorvec == color_i)
                _x = reshape(_vecx,size(x))
                vecfx1 = vec(f(_x))
                dx = (vecfx1-vecfx)/epsilon
                Ji = _make_Ji(rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                J = J + Ji
            end
        end
    elseif fdtype == Val{:central}
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x1_save = vecx1[color_i]
                x_save = vecx[color_i]
                epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir)
                _vecx1 = Base.setindex(vecx1,x1_save+epsilon,color_i)
                _vecx = Base.setindex(vecx,x_save-epsilon,color_i)
                _x1 = reshape(_vecx1,size(x))
                _x = reshape(_vecx,size(x))
                vecfx1 = vec(f(_x1))
                vecfx = vec(f(_x))
                dx = (vecfx1-vecfx)/(2epsilon)
                J = J + _make_Ji(eltype(x), dx, color_i, nrows, ncols)
            else
                tmp = norm(vecx1 .* (colorvec .== color_i))
                epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir)
                _vecx1 = @. vecx1 + epsilon * (colorvec == color_i)
                _vecx = @. vecx - epsilon * (colorvec == color_i)
                _x1 = reshape(_vecx1,size(x))
                _x = reshape(_vecx, size(x))
                vecfx1 = vec(f(_x1))
                vecfx = vec(f(_x))
                dx = (vecfx1-vecfx)/(2epsilon)
                Ji = _make_Ji(rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                J = J + Ji
            end
        end
    elseif fdtype == Val{:complex} && returntype <: Real
        epsilon = eps(eltype(x))
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x_save = vecx[color_i]
                _vecx = Base.setindex(complex.(vecx),x_save+im*epsilon,color_i)
                _x = reshape(_vecx,size(x))
                vecfx = vec(f(_x))
                dx = imag(vecfx)/epsilon
                J = J + _make_Ji(eltype(x), dx, color_i, nrows, ncols)
            else
                _vecx = @. vecx + im * epsilon * (colorvec == color_i)
                _x = reshape(_vecx,size(x))
                vecfx = vec(f(_x))
                dx = imag(vecfx)/epsilon
                Ji = _make_Ji(rows_index,cols_index,dx,colorvec,color_i,nrows,ncols)
                J = J + Ji
            end
        end
    else
        fdtype_error(returntype)
    end
    J
end

function finite_difference_jacobian!(J::AbstractMatrix,
    f,
    x::AbstractArray{<:Number},
    fdtype     :: Type{T1}=Val{:forward},
    returntype :: Type{T2}=eltype(x),
    f_in       :: Union{AbstractArray{<:T2},Nothing}=nothing;
    relstep=default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = 1:length(x),
    sparsity = ArrayInterface.has_sparsestruct(J) ? J : nothing) where {T1,T2}
    if f_in isa Nothing && fdtype == Val{:forward}
        if size(J,1) == length(x)
            fx = zero(x)
        else
            fx = zeros(returntype,size(J,1))
        end
        f(fx,x)
        cache = JacobianCache(x, fx, fdtype, returntype)
    elseif f_in isa Nothing
        cache = JacobianCache(x, fdtype, returntype)
    else
        cache = JacobianCache(x, f_in, fdtype, returntype)
    end
    finite_difference_jacobian!(J, f, x, cache, cache.fx; relstep=relstep, absstep=absstep, colorvec=colorvec, sparsity=sparsity)
end

function finite_difference_jacobian!(
    J::AbstractMatrix{<:Number},
    f,
    x::AbstractArray{<:Number},
    cache::JacobianCache{T1,T2,T3,cType,sType,fdtype,returntype},
    f_in::Union{T2,Nothing}=nothing;
    relstep = default_relstep(fdtype, eltype(x)),
    absstep=relstep,
    colorvec = cache.colorvec,
    sparsity::Union{AbstractArray,Nothing} = cache.sparsity,
    dir = true) where {T1,T2,T3,cType,sType,fdtype,returntype}

    m, n = size(J)
    _color = reshape(colorvec,size(x)...)

    x1, fx, fx1 = cache.x1, cache.fx, cache.fx1
    copyto!(x1, x)
    vfx = vec(fx)

    rows_index = nothing
    cols_index = nothing
    if _use_findstructralnz(sparsity)
        rows_index, cols_index = ArrayInterface.findstructralnz(sparsity)
    end

    if sparsity !== nothing
        fill!(J,false)
    end

    if fdtype == Val{:forward}
        vfx1 = vec(fx1)

        if f_in isa Nothing
            f(fx, x)
            vfx = vec(fx)
        else
            vfx = vec(f_in)
        end

        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
                epsilon = compute_epsilon(Val{:forward}, x1_save, relstep, absstep, dir)
                ArrayInterface.allowed_setindex!(x1,x1_save + epsilon,color_i)
                f(fx1, x1)
                # J is dense, so either it is truly dense or this is the
                # compressed form of the coloring, so write into it.
                @. J[:,color_i] = (vfx1 - vfx) / epsilon
                # Now return x1 back to its original value
                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
            else # Perturb along the colorvec vector
                @. fx1 = x1 * (_color == color_i)
                tmp = norm(fx1)
                epsilon = compute_epsilon(Val{:forward}, sqrt(tmp), relstep, absstep, dir)
                @. x1 = x1 + epsilon * (_color == color_i)
                f(fx1, x1)
                # J is a sparse matrix, so decompress on the fly
                @. vfx1 = (vfx1 - vfx) / epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                else
                    #=
                    J.nzval[rows_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    or
                    J[rows_index, cols_index] .+= (colorvec[cols_index] .== color_i) .* vfx1[rows_index]
                    += means requires a zero'd out start
                    =#
                    if J isa SparseMatrixCSC
                        @. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index)
                    else
                        @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index, cols_index)
                    end
                end
                # Now return x1 back to its original value
                @. x1 = x1 - epsilon * (_color == color_i)
            end
        end #for ends here
    elseif fdtype == Val{:central}
        vfx1 = vec(fx1)
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x_save = ArrayInterface.allowed_getindex(x,color_i)
                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
                epsilon = compute_epsilon(Val{:central}, x_save, relstep, absstep, dir)
                ArrayInterface.allowed_setindex!(x1,x1_save+epsilon,color_i)
                ArrayInterface.allowed_setindex!(x,x_save-epsilon,color_i)
                f(fx1, x1)
                f(fx, x)
                @. J[:,color_i] = (vfx1 - vfx) / 2epsilon
                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
                ArrayInterface.allowed_setindex!(x,x_save,color_i)
            else # Perturb along the colorvec vector
                @. fx1 = x1 * (_color == color_i)
                tmp = norm(fx1)
                epsilon = compute_epsilon(Val{:central}, sqrt(tmp), relstep, absstep, dir)
                @. x1 = x1 + epsilon * (_color == color_i)
                @. x  = x  - epsilon * (_color == color_i)
                f(fx1, x1)
                f(fx, x)
                @. vfx1 = (vfx1 - vfx) / 2epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx1,colorvec,color_i,n)
                else
                    if J isa SparseMatrixCSC
                        @. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index)
                    else
                        @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx1,),rows_index),rows_index, cols_index)
                    end
                end
                @. x1 = x1 - epsilon * (_color == color_i)
                @. x  = x  + epsilon * (_color == color_i)
            end
        end
    elseif fdtype==Val{:complex} && returntype<:Real
        epsilon = eps(eltype(x))
        @inbounds for color_i ∈ 1:maximum(colorvec)
            if sparsity isa Nothing
                x1_save = ArrayInterface.allowed_getindex(x1,color_i)
                ArrayInterface.allowed_setindex!(x1,x1_save + im*epsilon, color_i)
                f(fx,x1)
                @. J[:,color_i] = imag(vfx) / epsilon
                ArrayInterface.allowed_setindex!(x1,x1_save,color_i)
            else # Perturb along the colorvec vector
                @. x1 = x1 + im * epsilon * (_color == color_i)
                f(fx,x1)
                @. vfx = imag(vfx) / epsilon
                if ArrayInterface.fast_scalar_indexing(x1)
                    _colorediteration!(J,sparsity,rows_index,cols_index,vfx,colorvec,color_i,n)
                else
                   if J isa SparseMatrixCSC
                        @. setindex!((J.nzval,),getindex((J.nzval,),rows_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx,),rows_index),rows_index)
                    else
                        @. setindex!((J,),getindex((J,),rows_index, cols_index) + (getindex((_color,),cols_index) == color_i) * getindex((vfx,),rows_index),rows_index, cols_index)
                    end
                end
                @. x1 = x1 - im * epsilon * (_color == color_i)
            end
        end
    else
        fdtype_error(returntype)
    end
    nothing
end

function resize!(cache::JacobianCache, i::Int)
    resize!(cache.x1,  i)
    resize!(cache.fx,  i)
    cache.fx1 != nothing && resize!(cache.fx1, i)
    cache.colorvec = 1:i
    nothing
end


text:
function checksolver(sol::DeepCFRSolver)
    game = sol.game
    h0 = initialhist(game)
    k = vectorized(game, infokey(game, h0))

    K = infokeytype(game)
    VK = first(Base.return_types(vectorized, (typeof(game),K)))
    @assert VK <: AbstractVector "vectorized game keys must be vectors"

    A = actions(game, h0)
    in_size = length(k)
    out_size = length(A)

    for (i,net) in enumerate(sol.V)
        try
            output = net(k)
        catch e
            if e isa DimensionMismatch
                throw(DimensionMismatch("Value network $i input layer size should be $(length(k))"))
            else
                rethrow(e)
            end
        end
        output = net(k)
        @assert(
            length(output) == out_size,
            "Value network $i output dim : $(length(output)) \n
            Required output dim: $out_size"
        )
    end

    # TODO: is there some `Flux.nfan` for chain?
    try
        output = sol.Π(k)
    catch e
        if e isa DimensionMismatch
            throw(DimensionMismatch("Strategy input layer size should be $(length(k))"))
        else
            rethrow(e)
        end
    end
    output = sol.Π(k)
    @assert(
        length(output) == out_size,
        "Strategy network output dim : $(length(output)) \n
        Required output dim: $out_size"
    )

    @assert(
        sum(output) ≈ 1.0,
        "Strategy network output should be normalized to 1.0"
    )

    return "Tests Pass" # better way to do this?
end


text:
function Black76(CallPutFlag, F, X, T, r, v)

  d1 = (log(F / X) + (v^2 / 2) * T) / (v * sqrt(T))
  d2 = d1 - v * sqrt(T)

  if CallPutFlag == "c"
          exp(-r * T) * (F * CND(d1) - X * CND(d2))
  elseif CallPutFlag == "p"
          exp(-r * T) * (X * CND(-d2) - F * CND(-d1))
  else
          println(STDERR, "Falsche Eingabe! Abbruch")
  end

end


text:
###############################################################################
#
# Download ERA5 regenerated weather file
#
###############################################################################
const ERA5_LAND_ALL = [
            "10m_u_component_of_wind",
            "10m_v_component_of_wind",
            "2m_dewpoint_temperature",
            "2m_temperature",
            "evaporation_from_bare_soil",
            "evaporation_from_open_water_surfaces_excluding_oceans",
            "evaporation_from_the_top_of_canopy",
            "evaporation_from_vegetation_transpiration",
            "forecast_albedo",
            "lake_bottom_temperature",
            "lake_ice_depth",
            "lake_ice_temperature",
            "lake_mix_layer_depth",
            "lake_mix_layer_temperature",
            "lake_shape_factor",
            "lake_total_layer_temperature",
            "leaf_area_index_high_vegetation",
            "leaf_area_index_low_vegetation",
            "potential_evaporation",
            "runoff",
            "skin_reservoir_content",
            "skin_temperature",
            "snow_albedo",
            "snow_cover",
            "snow_density",
            "snow_depth",
            "snow_depth_water_equivalent",
            "snow_evaporation",
            "snowfall",
            "snowmelt",
            "soil_temperature_level_1",
            "soil_temperature_level_2",
            "soil_temperature_level_3",
            "soil_temperature_level_4",
            "sub_surface_runoff",
            "surface_latent_heat_flux",
            "surface_net_solar_radiation",
            "surface_net_thermal_radiation",
            "surface_pressure",
            "surface_runoff",
            "surface_sensible_heat_flux",
            "surface_solar_radiation_downwards",
            "surface_thermal_radiation_downwards",
            "temperature_of_snow_layer",
            "total_evaporation",
            "total_precipitation",
            "volumetric_soil_water_layer_1",
            "volumetric_soil_water_layer_2",
            "volumetric_soil_water_layer_3",
            "volumetric_soil_water_layer_4"];
const ERA5_LAND_SELECTION = [
            "evaporation_from_vegetation_transpiration"];
const ERA5_LAND_TIMES = [
            "00:00", "01:00", "02:00", "03:00", "04:00", "05:00",
            "06:00", "07:00", "08:00", "09:00", "10:00", "11:00",
            "12:00", "13:00", "14:00", "15:00", "16:00", "17:00",
            "18:00", "19:00", "20:00", "21:00", "22:00", "23:00"];
const ERA5_SINGLE_LEVELS_ALL =[
            "100m_u_component_of_wind",
            "100m_v_component_of_wind",
            "10m_u_component_of_neutral_wind",
            "10m_u_component_of_wind",
            "10m_v_component_of_neutral_wind",
            "10m_v_component_of_wind",
            "10m_wind_gust_since_previous_post_processing",
            "2m_dewpoint_temperature",
            "2m_temperature",
            "air_density_over_the_oceans",
            "altimeter_corrected_wave_height",
            "altimeter_range_relative_correction",
            "altimeter_wave_height",
            "angle_of_sub_gridscale_orography",
            "anisotropy_of_sub_gridscale_orography",
            "benjamin_feir_index",
            "boundary_layer_dissipation",
            "boundary_layer_height",
            "charnock",
            "clear_sky_direct_solar_radiation_at_surface",
            "cloud_base_height",
            "coefficient_of_drag_with_waves",
            "convective_available_potential_energy",
            "convective_inhibition",
            "convective_precipitation",
            "convective_rain_rate",
            "convective_snowfall",
            "convective_snowfall_rate_water_equivalent",
            "downward_uv_radiation_at_the_surface",
            "duct_base_height",
            "eastward_gravity_wave_surface_stress",
            "eastward_turbulent_surface_stress",
            "evaporation",
            "forecast_albedo",
            "forecast_logarithm_of_surface_roughness_for_heat",
            "forecast_surface_roughness",
            "free_convective_velocity_over_the_oceans",
            "friction_velocity",
            "gravity_wave_dissipation",
            "high_cloud_cover",
            "high_vegetation_cover",
            "ice_temperature_layer_1",
            "ice_temperature_layer_2",
            "ice_temperature_layer_3",
            "ice_temperature_layer_4",
            "instantaneous_10m_wind_gust",
            "instantaneous_eastward_turbulent_surface_stress",
            "instantaneous_large_scale_surface_precipitation_fraction",
            "instantaneous_moisture_flux",
            "instantaneous_northward_turbulent_surface_stress",
            "instantaneous_surface_sensible_heat_flux",
            "k_index",
            "lake_bottom_temperature",
            "lake_cover",
            "lake_depth",
            "lake_ice_depth",
            "lake_ice_temperature",
            "lake_mix_layer_depth",
            "lake_mix_layer_temperature",
            "lake_shape_factor",
            "lake_total_layer_temperature",
            "land_sea_mask",
            "large_scale_precipitation",
            "large_scale_precipitation_fraction",
            "large_scale_rain_rate",
            "large_scale_snowfall",
            "large_scale_snowfall_rate_water_equivalent",
            "leaf_area_index_high_vegetation",
            "leaf_area_index_low_vegetation",
            "low_cloud_cover",
            "low_vegetation_cover",
            "maximum_2m_temperature_since_previous_post_processing",
            "maximum_individual_wave_height",
            "maximum_total_precipitation_rate_since_previous_post_processing",
            "mean_boundary_layer_dissipation",
            "mean_convective_precipitation_rate",
            "mean_convective_snowfall_rate",
            "mean_direction_of_total_swell",
            "mean_direction_of_wind_waves",
            "mean_eastward_gravity_wave_surface_stress",
            "mean_eastward_turbulent_surface_stress",
            "mean_evaporation_rate",
            "mean_gravity_wave_dissipation",
            "mean_large_scale_precipitation_fraction",
            "mean_large_scale_precipitation_rate",
            "mean_large_scale_snowfall_rate",
            "mean_northward_gravity_wave_surface_stress",
            "mean_northward_turbulent_surface_stress",
            "mean_period_of_total_swell",
            "mean_period_of_wind_waves",
            "mean_potential_evaporation_rate",
            "mean_runoff_rate",
            "mean_sea_level_pressure",
            "mean_snow_evaporation_rate",
            "mean_snowfall_rate",
            "mean_snowmelt_rate",
            "mean_square_slope_of_waves",
            "mean_sub_surface_runoff_rate",
            "mean_surface_direct_short_wave_radiation_flux",
            "mean_surface_direct_short_wave_radiation_flux_clear_sky",
            "mean_surface_downward_long_wave_radiation_flux",
            "mean_surface_downward_long_wave_radiation_flux_clear_sky",
            "mean_surface_downward_short_wave_radiation_flux",
            "mean_surface_downward_short_wave_radiation_flux_clear_sky",
            "mean_surface_downward_uv_radiation_flux",
            "mean_surface_latent_heat_flux",
            "mean_surface_net_long_wave_radiation_flux",
            "mean_surface_net_long_wave_radiation_flux_clear_sky",
            "mean_surface_net_short_wave_radiation_flux",
            "mean_surface_net_short_wave_radiation_flux_clear_sky",
            "mean_surface_runoff_rate",
            "mean_surface_sensible_heat_flux",
            "mean_top_downward_short_wave_radiation_flux",
            "mean_top_net_long_wave_radiation_flux",
            "mean_top_net_long_wave_radiation_flux_clear_sky",
            "mean_top_net_short_wave_radiation_flux",
            "mean_top_net_short_wave_radiation_flux_clear_sky",
            "mean_total_precipitation_rate",
            "mean_vertical_gradient_of_refractivity_inside_trapping_layer",
            "mean_vertically_integrated_moisture_divergence",
            "mean_wave_direction",
            "mean_wave_direction_of_first_swell_partition",
            "mean_wave_direction_of_second_swell_partition",
            "mean_wave_direction_of_third_swell_partition",
            "mean_wave_period",
            "mean_wave_period_based_on_first_moment",
            "mean_wave_period_based_on_first_moment_for_swell",
            "mean_wave_period_based_on_first_moment_for_wind_waves",
            "mean_wave_period_based_on_second_moment_for_swell",
            "mean_wave_period_based_on_second_moment_for_wind_waves",
            "mean_wave_period_of_first_swell_partition",
            "mean_wave_period_of_second_swell_partition",
            "mean_wave_period_of_third_swell_partition",
            "mean_zero_crossing_wave_period",
            "medium_cloud_cover",
            "minimum_2m_temperature_since_previous_post_processing",
            "minimum_total_precipitation_rate_since_previous_post_processing",
            "minimum_vertical_gradient_of_refractivity_inside_trapping_layer",
            "model_bathymetry",
            "near_ir_albedo_for_diffuse_radiation",
            "near_ir_albedo_for_direct_radiation",
            "normalized_energy_flux_into_ocean",
            "normalized_energy_flux_into_waves",
            "normalized_stress_into_ocean",
            "northward_gravity_wave_surface_stress",
            "northward_turbulent_surface_stress",
            "ocean_surface_stress_equivalent_10m_neutral_wind_direction",
            "ocean_surface_stress_equivalent_10m_neutral_wind_speed",
            "orography",
            "peak_wave_period",
            "period_corresponding_to_maximum_individual_wave_height",
            "potential_evaporation",
            "precipitation_type",
            "runoff",
            "sea_ice_cover",
            "sea_surface_temperature",
            "significant_height_of_combined_wind_waves_and_swell",
            "significant_height_of_total_swell",
            "significant_height_of_wind_waves",
            "significant_wave_height_of_first_swell_partition",
            "significant_wave_height_of_second_swell_partition",
            "significant_wave_height_of_third_swell_partition",
            "skin_reservoir_content",
            "skin_temperature",
            "slope_of_sub_gridscale_orography",
            "snow_albedo",
            "snow_density",
            "snow_depth",
            "snow_evaporation",
            "snowfall",
            "snowmelt",
            "soil_temperature_level_1",
            "soil_temperature_level_2",
            "soil_temperature_level_3",
            "soil_temperature_level_4",
            "soil_type",
            "standard_deviation_of_filtered_subgrid_orography",
            "standard_deviation_of_orography",
            "sub_surface_runoff",
            "surface_latent_heat_flux",
            "surface_net_solar_radiation",
            "surface_net_solar_radiation_clear_sky",
            "surface_net_thermal_radiation",
            "surface_net_thermal_radiation_clear_sky",
            "surface_pressure",
            "surface_runoff",
            "surface_sensible_heat_flux",
            "surface_solar_radiation_downward_clear_sky",
            "surface_solar_radiation_downwards",
            "surface_thermal_radiation_downward_clear_sky",
            "surface_thermal_radiation_downwards",
            "temperature_of_snow_layer",
            "toa_incident_solar_radiation",
            "top_net_solar_radiation",
            "top_net_solar_radiation_clear_sky",
            "top_net_thermal_radiation",
            "top_net_thermal_radiation_clear_sky",
            "total_cloud_cover",
            "total_column_cloud_ice_water",
            "total_column_cloud_liquid_water",
            "total_column_ozone",
            "total_column_rain_water",
            "total_column_snow_water",
            "total_column_supercooled_liquid_water",
            "total_column_water",
            "total_column_water_vapour",
            "total_precipitation",
            "total_sky_direct_solar_radiation_at_surface",
            "total_totals_index",
            "trapping_layer_base_height",
            "trapping_layer_top_height",
            "type_of_high_vegetation",
            "type_of_low_vegetation",
            "u_component_stokes_drift",
            "uv_visible_albedo_for_diffuse_radiation",
            "uv_visible_albedo_for_direct_radiation",
            "v_component_stokes_drift",
            "vertical_integral_of_divergence_of_cloud_frozen_water_flux",
            "vertical_integral_of_divergence_of_cloud_liquid_water_flux",
            "vertical_integral_of_divergence_of_geopotential_flux",
            "vertical_integral_of_divergence_of_kinetic_energy_flux",
            "vertical_integral_of_divergence_of_mass_flux",
            "vertical_integral_of_divergence_of_moisture_flux",
            "vertical_integral_of_divergence_of_ozone_flux",
            "vertical_integral_of_divergence_of_thermal_energy_flux",
            "vertical_integral_of_divergence_of_total_energy_flux",
            "vertical_integral_of_eastward_cloud_frozen_water_flux",
            "vertical_integral_of_eastward_cloud_liquid_water_flux",
            "vertical_integral_of_eastward_geopotential_flux",
            "vertical_integral_of_eastward_heat_flux",
            "vertical_integral_of_eastward_kinetic_energy_flux",
            "vertical_integral_of_eastward_mass_flux",
            "vertical_integral_of_eastward_ozone_flux",
            "vertical_integral_of_eastward_total_energy_flux",
            "vertical_integral_of_eastward_water_vapour_flux",
            "vertical_integral_of_energy_conversion",
            "vertical_integral_of_kinetic_energy",
            "vertical_integral_of_mass_of_atmosphere",
            "vertical_integral_of_mass_tendency",
            "vertical_integral_of_northward_cloud_frozen_water_flux",
            "vertical_integral_of_northward_cloud_liquid_water_flux",
            "vertical_integral_of_northward_geopotential_flux",
            "vertical_integral_of_northward_heat_flux",
            "vertical_integral_of_northward_kinetic_energy_flux",
            "vertical_integral_of_northward_mass_flux",
            "vertical_integral_of_northward_ozone_flux",
            "vertical_integral_of_northward_total_energy_flux",
            "vertical_integral_of_northward_water_vapour_flux",
            "vertical_integral_of_potential_and_internal_energy",
            "vertical_integral_of_potential_internal_and_latent_energy",
            "vertical_integral_of_temperature",
            "vertical_integral_of_thermal_energy",
            "vertical_integral_of_total_energy",
            "vertically_integrated_moisture_divergence",
            "volumetric_soil_water_layer_1",
            "volumetric_soil_water_layer_2",
            "volumetric_soil_water_layer_3",
            "volumetric_soil_water_layer_4",
            "wave_spectral_directional_width",
            "wave_spectral_directional_width_for_swell",
            "wave_spectral_directional_width_for_wind_waves",
            "wave_spectral_kurtosis",
            "wave_spectral_peakedness",
            "wave_spectral_skewness",
            "zero_degree_level",];
const ERA5_SINGLE_LEVELS_SELECTION =[
            "10m_u_component_of_wind",
            "10m_v_component_of_wind",
            "2m_dewpoint_temperature",
            "2m_temperature",
            "mean_surface_direct_short_wave_radiation_flux",
            "mean_surface_downward_short_wave_radiation_flux",
            "skin_temperature",
            "soil_temperature_level_1",
            "soil_temperature_level_2",
            "soil_temperature_level_3",
            "soil_temperature_level_4",
            "surface_pressure",
            "volumetric_soil_water_layer_1",
            "volumetric_soil_water_layer_2",
            "volumetric_soil_water_layer_3",
            "volumetric_soil_water_layer_4"];

"""
    fetch_ERA5!(dt::ERA5LandHourly,
                lat_ind::Int,
                lon_ind::Int,
                year::Int = 2020,
                month::Int = 6;
                vars::Array{String,1} = ERA5_LAND_ALL,
                times::Array{String,1} = ERA5_LAND_TIMES,
                folder::String = ""
    )
    fetch_ERA5!(dt::ERA5LandHourly,
                lat_ind::Int,
                lon_ind::Int,
                year::Int;
                vars::Array{String,1} = ERA5_LAND_ALL,
                times::Array{String,1} = ERA5_LAND_TIMES,
                folder::String = ""
    )

Download dataset per 1 by 1 degree, given
- `dt` Dataset type
- `lat_ind` Latitude index (0 for -90~-89, 90 for -1~0)
- `lon_ind` Longitude index (0 for -180~-179, 180 for -1~0)
- `year` Which year of data to download, default = 2020
- `month` Which month of data to download, default = 6 (June)
- `vars` Variables to download, default is `ERA5_LAND_ALL`
- `times` Time list to download, default is all 24 hours
- `folder` Folder to save the downloaded dataset, default is "" (current dir)
"""
function fetch_ERA5!(
            dt::ERA5LandHourly,
            lat_ind::Int,
            lon_ind::Int,
            year::Int,
            month::Int;
            vars::Array{String,1} = ERA5_LAND_ALL,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading land data for year $(year) month $(month)...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_LAND_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_LAND_ALL;
    end

    # set lat and lon bounds
    lat1 = lat_ind - 91;
    lat2 = lat_ind - 90;
    lon1 = lon_ind - 181;
    lon2 = lon_ind - 180;

    # file name to save the downloaded dataset
    file_name = "$(folder)ERA5_LAND_$(year)_$(month)_$(lat_ind)_$(lon_ind).nc";

    # retrieve data if file does not exist
    if !isfile(file_name)
        CDSAPI_CLIENT.retrieve(
            "reanalysis-era5-land",
            Dict(
                "format"   => "netcdf",
                "variable" => vars,
                "year"     => "$(year)",
                "month"    => "$(month)",
                "day"      => ["$(i)" for i in 1:month_days(year,month)],
                "time"     => times,
                "area"     => ["$(lat2)", "$(lon1)", "$(lat1)", "$(lon2)"],
            ),
            file_name
        );
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5LandHourly,
            lat_ind::Int,
            lon_ind::Int,
            year::Int;
            vars::Array{String,1} = ERA5_LAND_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading land data for year $(year)...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_LAND_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_LAND_ALL;
    end

    # set lat and lon bounds
    lat1 = lat_ind - 91;
    lat2 = lat_ind - 90;
    lon1 = lon_ind - 181;
    lon2 = lon_ind - 180;

    # file name to save the downloaded dataset
    file_name = "$(folder)ERA5_LAND_$(year)_$(lat_ind)_$(lon_ind).nc";

    # retrieve data if file does not exist
    if !isfile(file_name)
        CDSAPI_CLIENT.retrieve(
            "reanalysis-era5-land",
            Dict(
                "format"   => "netcdf",
                "variable" => vars,
                "year"     => "$(year)",
                "month"    => ["$(i)" for i in 1:12],
                "day"      => ["$(i)" for i in 1:31],
                "time"     => times,
                "area"     => ["$(lat2)", "$(lon1)", "$(lat1)", "$(lon2)"],
            ),
            file_name
        );
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5LandHourly,
            year::Int;
            vars::Array{String,1} = ERA5_LAND_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading land data for year $(year)...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_LAND_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_LAND_ALL;
    end

    # file name to save the downloaded dataset
    file_name = "$(folder)ERA5_LAND_$(year).nc";

    # retrieve data if file does not exist
    if !isfile(file_name)
        CDSAPI_CLIENT.retrieve(
            "reanalysis-era5-land",
            Dict(
                "format"   => "netcdf",
                "variable" => vars,
                "year"     => "$(year)",
                "month"    => ["$(i)" for i in 1:12],
                "day"      => ["$(i)" for i in 1:31],
                "time"     => times,
            ),
            file_name
        );
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5LandHourly,
            year::Int,
            per_item::Bool;
            vars::Array{String,1} = ERA5_LAND_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading land data for year $(year) per item...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_LAND_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_LAND_ALL;
    end

    for _item in vars
        @info "Querying $(_item)";

        # file name to save the downloaded dataset
        file_name = "$(folder*_item)_LAND_$(year).nc";

        # retrieve data if file does not exist
        if !isfile(file_name)
            CDSAPI_CLIENT.retrieve(
                "reanalysis-era5-land",
                Dict(
                    "format"   => "netcdf",
                    "variable" => _item,
                    "year"     => "$(year)",
                    "month"    => ["$(i)" for i in 1:12],
                    "day"      => ["$(i)" for i in 1:31],
                    "time"     => times,
                ),
                file_name
            );
        end
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5SingleLevelsHourly,
            lat_ind::Int,
            lon_ind::Int,
            year::Int;
            vars::Array{String,1} = ERA5_SINGLE_LEVELS_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading single levels data for year $(year)...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_SINGLE_LEVELS_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_SINGLE_LEVELS_ALL;
    end

    # set lat and lon bounds
    lat1 = lat_ind - 91;
    lat2 = lat_ind - 90;
    lon1 = lon_ind - 181;
    lon2 = lon_ind - 180;

    # file name to save the downloaded dataset
    file_name = "$(folder)ERA5_SL_$(year)_$(lat_ind)_$(lon_ind).nc";

    # retrieve data if file does not exist
    if !isfile(file_name)
        CDSAPI_CLIENT.retrieve(
            "reanalysis-era5-single-levels",
            Dict(
                "product_type" => "reanalysis",
                "format"       => "netcdf",
                "variable"     => vars,
                "year"         => "$(year)",
                "month"        => ["$(i)" for i in 1:12],
                "day"          => ["$(i)" for i in 1:31],
                "time"         => times,
                "area"         => ["$(lat2)", "$(lon1)", "$(lat1)", "$(lon2)"],
            ),
            file_name
        );
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5SingleLevelsHourly,
            year::Int;
            vars::Array{String,1} = ERA5_SINGLE_LEVELS_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading single levels data for year $(year)...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_SINGLE_LEVELS_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_SINGLE_LEVELS_ALL;
    end

    # file name to save the downloaded dataset
    file_name = "$(folder)ERA5_SL_$(year).nc";

    # retrieve data if file does not exist
    if !isfile(file_name)
        CDSAPI_CLIENT.retrieve(
            "reanalysis-era5-single-levels",
            Dict(
                "product_type" => "reanalysis",
                "format"       => "netcdf",
                "variable"     => vars,
                "year"         => "$(year)",
                "month"        => ["$(i)" for i in 1:12],
                "day"          => ["$(i)" for i in 1:31],
                "time"         => times,
            ),
            file_name
        );
    end

    return nothing
end

function fetch_ERA5!(
            dt::ERA5SingleLevelsHourly,
            year::Int,
            per_item::Bool;
            vars::Array{String,1} = ERA5_SINGLE_LEVELS_SELECTION,
            times::Array{String,1} = ERA5_LAND_TIMES,
            folder::String = ""
)
    @info "Downloading single levels data for year $(year) per item...";

    # set up CDSAPI_CLIENT
    global CDSAPI_CLIENT, ERA5_SINGLE_LEVELS_ALL;
    update_CDSAPI_password!();

    # make sure data label is supported
    for _var in vars
        @assert _var in ERA5_SINGLE_LEVELS_ALL;
    end

    for _item in vars
        @info "Querying $(_item)";

        # file name to save the downloaded dataset
        file_name = "$(folder*_item)_SL_$(year).nc";

        # retrieve data if file does not exist
        if !isfile(file_name)
            CDSAPI_CLIENT.retrieve(
                "reanalysis-era5-single-levels",
                Dict(
                    "product_type" => "reanalysis",
                    "format"       => "netcdf",
                    "variable"     => _item,
                    "year"         => "$(year)",
                    "month"        => ["$(i)" for i in 1:12],
                    "day"          => ["$(i)" for i in 1:31],
                    "time"         => times,
                ),
                file_name
            );
        end
    end

    return nothing
end


text:
using Requires
@init @require CUDAnative = "be33ccc6-a3ff-5ff2-a52e-74243cff1e17" begin
  using .CUDAnative
end

function update!(fast_dQ, slow_dQ, δ, slow_rka = nothing)
  @inbounds @loop for i = (1:length(fast_dQ);
                           (blockIdx().x - 1) * blockDim().x + threadIdx().x)
    fast_dQ[i] += δ * slow_dQ[i]
    if slow_rka !== nothing
      slow_dQ[i] *= slow_rka
    end
  end
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA