text:
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                           SPECİAL FUNCTIONS                                             #
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
###########################################################################################################
export LegendreP, NLegendreP, LegendrePG, LegendreQ, NLegendreQ, LaguerreL, LaguerreG, SphericalHarmonicsY,
SphericalHarmonicsYG, SphericalHarmonicsS, SphericalHarmonicsSG, SphericalSpinors, RadoslawSSpinors, RotaD,
Rotad, RotaT, SlaterORadial, SlaterOrbitalY, SlaterOrbitalS, ChargeDOneC, OneCenterP, SlaterSRadial,
Hypergeometric0F0, Hypergeometric1F0, Hypergeometric0F1, Hypergeometric1F1, HypergeometricU, WhittakerM,
WhittakerW, Hypergeometric2F1, AuxiliaryGk, AuxiliaryGm, AuxiliaryGl, AuxiliaryGh, AuxiliaryG,
AuxiliaryGrlm, AuxiliaryGrlp, AuxiliaryGrh
###########################################################################################################
########################################## LEGENDRE POLYNOMIALS ###########################################(43)
################# Condon-Shortley phases.
################# First kind associated Legendre polynomials.
################# Type 0 and type 1 respectively correspond to type 2 and type 3 in Mathematica and mpmath.
################# Mathematica formoat: LegendreP[n,m,type,x]
function LegendreP(n::acb, m::acb, x::acb, t::UInt)
    z = parent(x)()
    ccall((:acb_hypgeom_legendre_p, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, UInt, Int),
    z, n, m, x, t, parent(n).prec)
    return z
end
function LegendreP(n::acb, m::acb, x::acb, t::Int)
    LegendreP(n, m, x, UInt(t))
end
function LegendreP(n::Complex, m::Complex, x::Complex, t::Int)
    LegendreP(CF(real(n), imag(n)), CF(real(m), imag(m)), CF(real(x), imag(x)), UInt(t))
end
function LegendreP(n::arb, m::arb, x::arb, t::Int)
    LegendreP(CF(n), CF(m), CF(x), t)
end
function LegendreP(n::Real, m::Real, x::Real, t::Int)
    LegendreP(CF(n), CF(m), CF(x), t)
end
###########################################################################################################
################# First kind normalized associated Legendre polynomials
function NLegendreP(n::acb, m::acb, x::acb, t::UInt)
    Sqrt(((2n+1)//2) * (Factorial(n-m)//Factorial(n+m)))*LegendreP(n, m, x, t)
end
function NLegendreP(n::acb, m::acb, x::acb, t::Int)
    NLegendreP(n, m, x, UInt(t))
end
function NLegendreP(n::Complex, m::Complex, x::Complex, t::Int)
    NLegendreP(CF(real(n), imag(n)), CF(real(m), imag(m)), CF(real(x), imag(x)), UInt(t))
end
function NLegendreP(n::arb, m::arb, x::arb, t::Int)
    NLegendreP(CF(n), CF(m), CF(x), t)
end
function NLegendreP(n::Real, m::Real, x::Real, t::Int)
    NLegendreP(CF(n), CF(m), CF(x), t)
end
###########################################################################################################
################# LegendrePG is a Guseinov's formulations for normalized associated Legendre polynomials
################# ((-1)^m) for Condon-Shortley phase
function LegendrePG(n::Int, m::Int, x::acb)
    lowl = Zero(0)
    upl = Int(floor((n-abs(m))/2))

    res1 = Sqrt((2*n+1)/2)*((1-x^2)^(Abs(m)/2))
    res2 = zeros(CF, upl - lowl + 1)

    for k in lowl : upl
        res2[k - lowl + 1] = SphCB(k, n, abs(m))*(x^(n-abs(m)-2*k))
    end
    cs = m < 0 ? RF(1) : Power(-1,m) #For Condon-Shortley Phases#
    return cs * res1 * sum(res2)
end
###########################################################################################################
################# Second kind associated Legendre polynomials
function LegendreQ(n::acb, m::acb, x::acb, t::UInt)
    z = parent(x)()
    ccall((:acb_hypgeom_legendre_q, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, UInt, Int),
    z, n, m, x, t, parent(n).prec)
    return z
end

function LegendreQ(n::acb, m::acb, x::acb, t::Int)
    LegendreQ(n, m, x, UInt(t))
end
function LegendreQ(n::Complex, m::Complex, x::Complex, t::Int)
    LegendreQ(CF(real(n), imag(n)), CF(real(m), imag(m)), CF(real(x), imag(x)), UInt(t))
end
function LegendreQ(n::arb, m::arb, x::arb, t::Int)
    LegendreQ(CF(n), CF(m), CF(x), t)
end
function LegendreQ(n::Real, m::Real, x::Real, t::Int)
    LegendreQ(CF(n), CF(m), CF(x), t)
end
###########################################################################################################
################# Second kind normalized associated Legendre polynomials
function NLegendreQ(n::acb, m::acb, x::acb, t::UInt)
    Sqrt(((2n+1)//2) * (Factorial(n-m)//Factorial(n+m)))*LegendreP(n, m, x, t)
end
function NLegendreQ(n::acb, m::acb, x::acb, t::Int)
    NLegendreQ(n, m, x, UInt(t))
end
function NLegendreQ(n::Complex, m::Complex, x::Complex, t::Int)
    NLegendreQ(CF(real(n), imag(n)), CF(real(m), imag(m)), CF(real(x), imag(x)), UInt(t))
end
function NLegendreQ(n::arb, m::arb, x::arb, t::Int)
    NLegendreQ(CF(n), CF(m), CF(x), t)
end
function NLegendreQ(n::Real, m::Real, x::Real, t::Int)
    NLegendreQ(CF(n), CF(m), CF(x), t)
end
###########################################################################################################
########################################## LAGUERRE POLYNOMIALS ###########################################(43)
################# Associated Laguerre polynomials L_{p}^{q}(x)
function LaguerreL(q::acb, p::acb, x::acb)
    z = parent(x)()
    ccall((:acb_hypgeom_laguerre_l, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, Int),
    z, q, p, x, parent(x).prec)
    return z
end
function LaguerreL(q::Complex, p::Complex, x::Complex)
    LaguerreL(CF(real(q), imag(q)), CF(real(p), imag(p)), CF(real(x), imag(x)))
end
function LaguerreL(q::arb, p::arb, x::arb)
    LaguerreL(CF(q), CF(p), CF(x))
end
function LaguerreL(q::Real, p::Real, x::Real)
    LaguerreL(CF(q), CF(p), CF(x))
end
################# Guseinoc Laguerre polynomials L_{p}^{q}(x)
function LaguerreG(q::acb, p::acb, x::acb)
    Power(CF(-1), p) * Gamma(q + CF(1)) * LaguerreL(q-p, p, x)
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                          SPHERICAL HARMONICS                                            #
###########################################################################################################
###########################################################################################################
###########################################################################################################
########################################### SPHERICAL HARMONICS ###########################################(43)
################# Condon-Shortley phases
################# Complex Spherical Harmonics
function SphericalHarmonicsY(n::UInt, m::UInt, theta::acb, phi::acb)
    z = parent(theta)()
    ccall((:acb_hypgeom_spherical_y, Nemo.libarb), Nothing, (Ref{acb}, UInt, UInt, Ref{acb}, Ref{acb}, Int),
    z, n, m, theta, phi, parent(theta).prec)
    return z
end
function SphericalHarmonicsY(n::Int, m::Int, theta::acb, phi::acb)
    if m >= 0
        SphericalHarmonicsY(UInt(n), UInt(m), theta, phi)
    else
        m = abs(m)
        SphericalHarmonicsY(UInt(n), -UInt(m), theta, phi)
    end
end
function SphericalHarmonicsY(n::Int, m::Int, theta::arb, phi::arb)
    SphericalHarmonicsY(n, m, CF(theta), CF(phi))
end
function SphericalHarmonicsY(n::Int, m::Int, theta::Real, phi::Real)
    SphericalHarmonicsY(n, m, CF(theta), CF(phi))
end

function SphericalHarmonicsY(n::UInt, m::UInt, x::arb, y::arb, z::arb)
    r = Sqrt(x^2 + y^2 + z^2)
    SphericalHarmonicsY(n, m, r * Sin(theta) * Cos(phi), r * Sin(theta) * Sin(phi), r * Cos(theta))
end
################# Complex Spherical Harmonics Guseinov formulae
function SphericalHarmonicsYG(n::Int, m::Int, x::arb, y::arb, z::arb)

    lowl = Zero(0)
    upl = Int(floor((n-abs(m))/2))
    res2 = zeros(RF, upl - lowl + 1)
    r = Sqrt(x^2 + y^2 + z^2)

    coef = Sqrt((2*n + 1) // (4*Pi(RF)))
    res1 = Power((x // r) + CF(0,1)* Epsilon(m,0) * ( y // r), Abs(m))

    for k in lowl : upl
        res2[k - lowl + 1] = SphCB(k, n, abs(m))*((z // r)^(n-abs(m)-2*k))
    end
    coef * res1 * sum(res2)

end

function SphericalHarmonicsYG(n::Int, m::Int, theta::arb, phi::arb)
    lowl = Zero(0)
    upl = Int(floor((n-abs(m))/2))
    res2 = zeros(RF, upl - lowl + 1)

    coef = Sqrt((2*n + 1) // (4*Pi(RF)))
    res1 = Power((Sin(theta) * Cos(phi)) + CF(0,1)* Epsilon(m,0) * (Sin(theta) * Sin(phi)), Abs(m))

    for k in lowl : upl
        res2[k - lowl + 1] = SphCB(k, n, abs(m))*((Cos(theta))^(n-abs(m)-2*k))
    end
    coef * res1 * sum(res2)
end
###########################################################################################################
################# SphericalHarmonicsS already in Condon-Shortley phase
function SphericalHarmonicsS(n::Int, m::Int, theta::acb, phi::acb)
    res1 = RSqrt(Pi(CF) * (One(m) + KroneckerDelta(m, Zero(m))))
    res2 = NLegendreP(CF(n), Abs(CF(m)), Cos(theta), UInt(0))
    res3 = ifelse(m >= 0 , Cos(Abs(m)*phi), Sin(Abs(m)*phi))
    res1 * res2 * res3
end
function SphericalHarmonicsS(n::Int, m::Int, theta::arb, phi::arb)
    SphericalHarmonicsS(n, m, CF(theta), CF(phi))
end
function SphericalHarmonicsS(n::Int, m::Int, theta::Real, phi::Real)
    SphericalHarmonicsS(n, m, CF(theta), CF(phi))
end

function SphericalHarmonicsS(n::Int, m::Int, x::arb, y::arb, z::arb)
    r = Sqrt(x^2 + y^2 + z^2)
    theta = ArcTan(y // x)
    phi = ArcCos(z // r)
    SphericalHarmonicsS(n, m, theta, phi)
end
################# Real Spherical Harmonics Guseinov formulae
################# ((-1)^m) for Condon-Shortley phase
function SphericalHarmonicsSG(n::Int, m::Int, theta::arb, phi::arb)

    lowl = Zero(0)
    upl = Int(floor((n-abs(m))/2))
    res2 = zeros(RF, upl - lowl + 1)

    if m > 0
        coef = (1 // 2) * Sqrt((2*n + 1) // (2 * Pi(RF)))
    elseif m < 0
        coef = -CF(0,1) * (1 // 2) * Sqrt((2*n + 1) // (2 * Pi(RF)))
    else
        coef = (1 // 2) * Sqrt((2*n + 1) // (1 * Pi(RF)))
    end

    xdr = Sin(theta) * Cos(phi)
    ydr = Sin(theta) * Sin(phi)
    zdr = Cos(theta)
    resxyp = Power(xdr + CF(0,1) * (ydr), Abs(m))
    resxyn = Power(xdr - CF(0,1) * (ydr), Abs(m))

    if m > 0
        res1 = resxyp + resxyn
    elseif m < 0
        res1 = resxyp - resxyn
    else
        res1 = RF(1)
    end

    for k in lowl : upl
        res2[k - lowl + 1] = SphCB(k, n, abs(m))*((zdr)^(n-abs(m)-2*k))
    end

    coef * res1 * sum(res2)
end

function SphericalHarmonicsSG(n::Int, m::Int, x::arb, y::arb, z::arb)

    lowl = Zero(0)
    upl = Int(floor((n-abs(m))/2))
    res2 = zeros(RF, upl - lowl + 1)

    if m > 0
        coef = (1 // 2) * Sqrt((2*n + 1) // (2 * Pi(RF)))
    elseif m < 0
        coef = -CF(0,1) * (1 // 2) * Sqrt((2*n + 1) // (2 * Pi(RF)))
    else
        coef = (1 // 2) * Sqrt((2*n + 1) // (1 * Pi(RF)))
    end

    r = Sqrt(x^2 + y^2 + z^2)
    xdr = x // r
    ydr = y // r
    zdr = z // r
    resxyp = Power(xdr + CF(0,1) * (ydr), Abs(m))
    resxyn = Power(xdr - CF(0,1) * (ydr), Abs(m))

    if m > 0
        res1 = resxyp + resxyn
    elseif m < 0
        res1 = resxyp - resxyn
    else
        res1 = RF(1)
    end

    for k in lowl : upl
        res2[k - lowl + 1] = SphCB(k, n, abs(m))*((zdr)^(n-abs(m)-2*k))
    end
    coef * res1 * sum(res2)

end

function SphericalHarmonicsSG(n::Int, m::Int, x::Real, y::Real, z::Real)
    SphericalHarmonicsSG(n, m, RF(x), RF(y), RF(z))
end
###########################################################################################################
############################################ SPHERICAL SPINORS ############################################(43)
function SphericalSpinors(β::Int, κ::Int, m::arb, theta::arb, phi::arb)
    lᵦ = ifelse(β*κ > 0, β*κ - 1, -β*κ)
    lᵦ = Int(convert(Float64, lᵦ))
    ml0 = Int(convert(Float64,mλ(m,0)))
    ml1 = Int(convert(Float64,mλ(m,1)))
    ucg = Sign(β * κ) * Sqrt((-β*κ + 1//2 - m) // (-2*β*κ + 1))
    lcg = Sqrt((-β*κ + 1//2 + m) // (-2*β*κ + 1))
    uc = ucg * Eta(mλ(m, 0)) * SphericalHarmonicsY(lᵦ, ml0, theta, phi)
    lc = lcg * Eta(mλ(m, 1)) * SphericalHarmonicsY(lᵦ, ml1, theta, phi)
    [uc, lc]
end

function SphericalSpinors(β::Int, κ::Int, m::arb, x::arb, y::arb, z::arb)
    r = Sqrt(x^2 + y^2 + z^2)
    theta = ArcTan(y // x)
    phi = ArcCos(z // r)

    lᵦ = ifelse(β*κ > 0, β*κ - 1, -β*κ)
    lᵦ = Int(convert(Float64, lᵦ))
    ml0 = Int(convert(Float64,mλ(m,0)))
    ml1 = Int(convert(Float64,mλ(m,1)))
    ucg = Sign(β * κ) * Sqrt((-β*κ + 1//2 - m) // (-2*β*κ + 1))
    lcg = Sqrt((-β*κ + 1//2 + m) // (-2*β*κ + 1))
    uc = ucg * Eta(mλ(m, 0)) * SphericalHarmonicsY(lᵦ, ml0, theta, phi)
    lc = lcg * Eta(mλ(m, 1)) * SphericalHarmonicsY(lᵦ, ml1, theta, phi)
    [uc, lc]
end

function RadoslawSSpinors(κ::Int, m::arb, theta::arb, phi::arb)

    l = ifelse(κ > 0, κ, -κ - 1)
    ml0 = Int(convert(Float64,mλ(m,0)))
    ml1 = Int(convert(Float64,mλ(m,1)))
    ucg = Sign(-κ) * Sqrt((κ + 1//2 -m) // (2 * κ + 1))
    lcg = Sqrt((κ + 1//2 + m) // (2 * κ + 1))
    uc = ucg * Eta(mλ(m, 0)) * SphericalHarmonicsY(l, ml0, theta, phi)
    lc = lcg * Eta(mλ(m, 1)) * SphericalHarmonicsY(l, ml1, theta, phi)
    [uc, lc]
end

function SphericalSpinors(l::Int, j::arb, m::arb, theta::arb, phi::arb)
    κ = 2*(l-j)*(j + 1/2)
    κ = Int(convert(Float64, κ))
    uc = SphericalSpinors(-1, κ, m, theta, phi)
    lc = SphericalSpinors(1, κ, m, theta, phi)
    reshape([uc lc], 4,1)
end

function SphericalSpinors(l::Int, j::arb, m::arb, x::arb, y::arb, z::arb)
    κ = 2*(l-j)*(j + 1/2)
    κ = Int(convert(Float64, κ))
    uc = SphericalSpinors(-1, κ, m, x, y, z)
    lc = SphericalSpinors(1, κ, m, x, y, z)
    reshape([uc lc], 4,1)
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                     ROTATED ANGULAR FUNCTIONS                                           #
###########################################################################################################
###########################################################################################################
###########################################################################################################(43)
################# x = 0 for the complex spherical harmonics, x = 1 for the real ones
function RotaD(λ::Int, l1::Int, m1::Int, l2::Int, m2::Int, theta::arb, phi::arb)
    coef = (RF(2) // (RF(1) + KroneckerDelta(λ,0)))
    resr = RF(0)
    for L in abs(l1-l2) : 2 : l1 + l2
        res1 = ClebschGordanG(l1, -m1, l2, m2, L, (-m1 + m2))
        res2 = ClebschGordanG(l1, -λ, l2, λ, L, 0)
        res3 = Sqrt((RF(4) * Pi(RF)) // (2 * L + 1))
        res4 = Eta(-m1 + m2) * SphericalHarmonicsY(L, -m1 + m2, theta, phi)
        resr += res1 * res2 * res3 * res4
    end
    res = coef * resr
    return res
end

function Rotad(λ::Int, l1::Int, m1::Int, l2::Int, m2::Int, theta::arb, phi::arb)
    if (theta == 0) && (phi == 0)
        res = KroneckerDelta(m1,m2) * KroneckerDelta(λ, abs(m1))
    elseif (theta == Pi(RF)) && (phi == 0)
        res = KroneckerDelta(m1,m2) * KroneckerDelta(λ, abs(m1)) * ((-1)^(l1+l2))
    else
        coef1 = (
            Power(CF(0,1), m1 >= 0 ? 1 : 0) *
            Power(CF(0,-1),m2 >= 0 ? 1 : 0)
        )
        coef2 = Sqrt(
            (RF(1) + KroneckerDelta(m1,0)) *
            (RF(1) + KroneckerDelta(m2,0))
        )
        coef = coef1//(2*coef2)
        RotaD1 = RotaD(λ, l1, abs(m1), l2, abs(m2), theta, phi)
        RotaD2 = Epsilon(m1, 0) * RotaD(λ, l1, -abs(m1), l2, abs(m2), theta, phi)
        RotaD3 = Epsilon(0, m2) * RotaD(λ, l1, abs(m1), l2, -abs(m2), theta, phi)
        RotaD4 = Epsilon(m1, m2) * RotaD(λ, l1, -abs(m1), l2, -abs(m2), theta, phi)
        res = coef * (RotaD1 + RotaD2 + RotaD3 + RotaD4)
    end
    return res
end

function RotaT(λ::Int, l1::Int, m1::Int, l2::Int, m2::Int, theta::arb, phi::arb, x::Int)
    if x == 0
        RotaD(λ::Int, l1::Int, m1::Int, l2::Int, m2::Int, theta::arb, phi::arb)
    else
        Rotad(λ::Int, l1::Int, m1::Int, l2::Int, m2::Int, theta::arb, phi::arb)
    end
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                          SLATER-TYPE FUNCTIONS                                          #
###########################################################################################################
###########################################################################################################
###########################################################################################################
function SlaterORadial(n::arb, ζ::arb, r::arb)
    (r^(n - 1)) * Exp(-ζ * r)
end
function SlaterORadial(n::Real, ζ::Real, r::Real)
    SlaterORadial(RF(n), RF(ζ), RF(r))
end

function SlaterOrbitalY(n, l, m, ζ, r, theta, phi)
    SlaterON(n, ζ) * SlaterORadial(n, ζ, r) * SphericalHarmonicsY(l, m, theta, phi)
end
function SlaterOrbitalY(n, l, m, ζ, r, x, y, z)
    r = Sqrt(x^2 + y^2 + z^2)
    theta = ArcTan(y // x)
    phi = ArcCos(z // r)

    SlaterON(n, ζ) * SlaterORadial(n, ζ, r) * SphericalHarmonicsY(l, m, theta, phi)
end

function SlaterOrbitalS(n, l, m, ζ, r, theta, phi)
    SlaterON(n, ζ) * SlaterORadial(n, ζ, r) * SphericalHarmonicsS(l, m, theta, phi)
end
function SlaterOrbitalS(n, l, m, ζ, r, x, y, z)
    r = Sqrt(x^2 + y^2 + z^2)
    theta = ArcTan(y // x)
    phi = ArcCos(z // r)

    SlaterON(n, ζ) * SlaterORadial(n, ζ, r) * SphericalHarmonicsS(l, m, theta, phi)
end
###########################################################################################################
###################################### CHARGE DENSITY DISTRIBUTION ########################################
function ChargeDOneC(
    n1::arb, l1::Int, m1::Int, ζ1::arb,
    n2::arb, l2::Int, m2::Int, ζ2::arb,
    N::arb, L::Int, M::Int
    )

    z = ζ1 + ζ2
    τ = ((ζ1 - ζ2) // (ζ1 + ζ2))

    res1 = (Power(z, 3//2) // Power(2, N))
    res2 =  ((2*L +1) // 2)
    res3 = ((Gamma(2*N + 1)) // (Gamma(2*n1 + 1) * Gamma(2*n2 + 1)))
    res4 = (1 + τ)
    res5 = (1 - τ)
    res6 = GGauntG(l1, m1, l2, m2, L, abs(M))
    res7 = RSphCA(M, m1, m2)

    res = res1 * Sqrt(res2 * res3) * res4 * res5 * res6 * res7
    return res
end
###########################################################################################################
######################################### ONE-CENTER POTENTIAL ############################################
function OneCenterP(L::Int, n1::arb, n2::arb, ζ1::arb, ζ2::arb, R::arb)
    ρ = (R//2) * (ζ1 + ζ2)
    τ = (ζ1 - ζ2) // (ζ1 + ζ2)

    coeff1 = SlaterON(n1, n2, RF(1), τ) * (2 * (ζ1 + ζ2))
    coeff2 = Gamma(n1 + n2 + L + 1)
    coeff3 = Power(2 * ρ, L + 1)
    coeff = coeff1 * coeff2 * (RF(1) // coeff3)

    res1 = PGamma(n1 + n2 + L + 1, 2 * ρ)
    res2 = Power(2 * ρ, 2 * L + 1)
    res3 = Pochhammer(n1 + n2 - L, 2 * L + 1)
    res4 = QGamma(n1 + n2 - L, 2 * ρ)
    res5 = res1 + (res2 // res3) * res4

    res = coeff * res5

    return res
end
###########################################################################################################
########################################### SLATER TYPE SPINORS ###########################################
################# Two-component form radial
function SlaterSRadial(β::Int, n::arb, κ::Int, ζ::arb, r::arb)
    ((SlaterSA(β, n, κ) * (r^n)) + (SlaterSB(β, n, κ) * (r^(n-1)))) * Exp(-ζ * r)
end

function SlaterSRadial(β::Int, n::arb, l::Int, j::arb, ζ::arb, r::arb)
    κ = 2*(l-j)*(j + 1/2)
    κ = Int(convert(Float64, κ))
    ((SlaterSA(β, n, κ) * (r^n)) + (SlaterSB(β, n, κ) * (r^(n-1)))) * Exp(-ζ * r)
end
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                          HYPERGEOMETRIC FUNCTIONS                                       #
###########################################################################################################
###########################################################################################################
############################################ HYPERGEOMETRIC0F0 ############################################(44)
Hypergeometric0F0(x::acb) = Exp(x)
Hypergeometric0F0(x::Complex) = Exp(x)
Hypergeometric0F0(x::arb) = Exp(x)
Hypergeometric0F0(x::Real) = Exp(x)
###########################################################################################################
############################################ HYPERGEOMETRIC1F0 ############################################(44)
Hypergeometric1F0(a::acb, x::acb) = Surd(1-x, a)
Hypergeometric1F0(a::Complex, x::Complex) = Surd(1-x, a)
Hypergeometric1F0(a::arb, x::arb) = Surd(1-x, a)
Hypergeometric1F0(a::Real, x::Real) = Surd(1-x, a)
###########################################################################################################
############################################ HYPERGEOMETRIC0F1 ############################################(44)
################# Confluent hypergeometric limit functions
################# Confluent hypergeometric function 0f1. For r = (0,1) = (irregularized, regularized)
################# r = 1 -> 0F1 = 0F1/Gamma(a)
function Hypergeometric0F1(a::acb, x::acb, r::Int)
    z = parent(x)()
    ccall((:acb_hypgeom_0f1, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Int, Int),
    z, a, x, r, parent(x).prec)
    return z
end
Hypergeometric0F1(a::Complex, x::Complex, r::Int) =
Hypergeometric0F1(CF(real(a), imag(a)), CF(real(x), imag(x)), r)

Hypergeometric0F1(a::arb, x::arb, r::Int) = Hypergeometric0F1(CF(a), CF(x), r)
Hypergeometric0F1(a::Real, x::Real, r::Int) = Hypergeometric0F1(CF(a), CF(x), r)
###########################################################################################################
############################################ HYPERGEOMETRIC1F1 ############################################(44)
################# Confluent hypergeometric functions of the first kind also writtten as M(a;b;x)
################# Kummer confluent hypergeometric functions
################# (0,1) = (irregularized, regularized)
function Hypergeometric1F1(a::acb, b::acb, x::acb, r::Int)
    z = parent(x)()
    ccall((:acb_hypgeom_m, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, Int, Int),
    z, a, b, x, r, parent(x).prec)
    return z
end
Hypergeometric1F1(a::Complex, b::Complex, x::Complex, r::Int) =
Hypergeometric1F1(CF(real(a), imag(a)), CF(real(b), imag(b)), CF(real(x), imag(x)), r)

Hypergeometric1F1(a::arb, b::arb, x::arb, r::Int) = Hypergeometric1F1(CF(a), CF(b), CF(x), r)
Hypergeometric1F1(a::Real, b::Real, x::Real, r::Int) = Hypergeometric1F1(CF(a), CF(b), CF(x), r)
###########################################################################################################
############################################# HYPERGEOMETRICU #############################################(45)
function HypergeometricU(a::acb, b::acb, x::acb)
    z = parent(x)()
    ccall((:acb_hypgeom_u, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, Int),
    z, a, b, x, parent(x).prec)
    return z
end
HypergeometricU(a::Complex, b::Complex, x::Complex) =
HypergeometricU(CF(real(a), imag(a)), CF(real(b), imag(b)), CF(real(x), imag(x)))

HypergeometricU(a::arb, b::arb, x::arb) = HypergeometricU(CF(a), CF(b), CF(x))
HypergeometricU(a::Real, b::Real, x::Real) = HypergeometricU(CF(a), CF(b), CF(x))
###########################################################################################################
################################################ WHITTAKER ################################################(48)
################# WhittakerM(k,m,z) = Exp(-z/2)*Power(z, m+1/2)*1F1(1/2+m-k,1+2m;z)
function WhittakerM(k, m, x)
    Exp(-x//2) * Power(x, m + (1//2)) * Hypergeometric1F1(1//2 + m - k, One(m) + 2*m, x, sprec)
end
################# WhittakerW(k,m,z) = Exp(-z/2)*Power(z, m+1/2)*U(1/2+m-k,1+2m;z)
function WhittakerW(k, m, x)
    Exp(-x//2) * Power(x, m + (1//2)) * HypergeometricU(1//2 + m - k, One(m) + 2*m, x)
end
###########################################################################################################
############################################ HYPERGEOMETRIC2F1 ############################################(44)
function Hypergeometric2F1(a::acb, b::acb, c::acb, x::acb, r::Int)
    z = parent(x)()
    ccall((:acb_hypgeom_2f1, Nemo.libarb), Nothing, (Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, Ref{acb}, Int, Int),
    z, a, b, c, x, r, parent(x).prec)
    return z
end
Hypergeometric2F1(a::Complex, b::Complex, c::Complex, x::Complex, r::Int) = Hypergeometric2F1(
    CF(real(a), imag(a)),
    CF(real(b), imag(b)),
    CF(real(c), imag(c)),
    CF(real(x), imag(x)),
    r
    )
Hypergeometric2F1(a::arb, b::arb, c::arb, x::arb, r::Int) = Hypergeometric2F1(
    CF(a), CF(b), CF(c), CF(x), r
    )
Hypergeometric2F1(a::Real, b::Real, c::Real, x::Real, r::Int) = Hypergeometric2F1(
    CF(a), CF(b), CF(c), CF(x), r
    )
###########################################################################################################
###########################################################################################################
###########################################################################################################
#                                     MOLECULAR AUXILIARY FUNCTIONS                                       #
###########################################################################################################
###########################################################################################################
############################################## OVERLAP-LİKE G #############################################(45)
################# Overlap-like G used for calculationg the two-center one-electron integrals
###########################################################################################################
###########################################################################################################
################# Case 1 The parameter p_{3} = 0
function AuxiliaryGk(n1::arb, q::Int, n2::arb, n3::arb, p1::arb, p2::arb)
    two = parent(n1)(2)
    res1 = ((p1^n1) // (Gamma(n1+1))) * (two ^ (n2+n3+q+One(n2)))
    res2 = Beta(n2+1, n3+1, parent(n2)(1/2)) * ExpIntegralE(-(n2+n3+q+1), p2)
    res = res1 * res2
end

function AuxiliaryGm(n1::arb, n2::arb, p::arb)
    two = parent(n1)(2)
    (two ^ n1) * HypergeometricU(n2+One(n2), n1+n2+two, p) * Gamma(n2+One(n2)) * Exp(-p)
end

function AuxiliaryGl(n1::arb, q::Int, n2::arb, n3::arb, p1::arb, p2::arb, lim::Int)
    two = parent(n1)(2)
    lowl = Zero(lim)
    upl = lim

    a = (n3) // (n2 + q + One(n2))
    b = (n2-n3+q+One(n2)-p2) // (n2+q+One(n2))

    res = zeros(RF, upl-lowl+1)
    p = zeros(RF, upl-lowl+1)
    f = zeros(RF, upl-lowl+1)
    m1 = zeros(RF, upl-lowl+1)
    m2 = zeros(RF, upl-lowl+1)
    m = zeros(RF, upl-lowl+1)

    p[1] = Pochhammer(-n2,0)
    f[1] = RF(1)
    m1[1] = real(AuxiliaryGm(n2 + q + two, n3 - One(n3), p2))
    m2[1] = real(AuxiliaryGm(n2 + q + One(n2), n3, p2))
    m[1] = (1//4) * a * m1[1] + (1//2) * b * m2[1]

    for s in lowl : upl-1
        s1 = s-lowl+1
        p[s1+1] = (-n2 + s1 -One(n2)) * p[s1]
        f[s1+1] = (s+1) * f[s1]
        a = (n3+s1) // (n2 + q + One(n2) - s1)
        b = (n2-n3-2*s1+q+One(n2)-p2) // (n2+q+One(n2)-s1)
        m[s1+1] = (1//4) * a * m2[s1] + (1//2) * b * m[s1]
        m2[s1+1] = m[s1]
    end
    for s in lowl : upl
        s1 = s-lowl+1
        res[s1] = (p[s1] // ((n3 + s + One(n3)) * f[s1])) * m[s1]
    end
    ((p1^n1) // Gamma(n1 + One(n1))) * sum(res)
end

function AuxiliaryGh(n1::arb, q::Int, n2::arb, n3::arb, p1::arb, p2::arb, lim::Int)
    two = parent(n1)(2)
    res1 = ((p1^n1) // Gamma(n1 + One(n1))) * (two^(n2+n3+q+One(n1))) * Beta(n2+One(n2), n3+One(n3))
    res2 = ExpIntegralE(-(n2 + n3 + q + One(n1)), p2)
    res3 = AuxiliaryGl(n1, q, n2, n3, p1, p2, lim)
    res1 * res2 - res3
end

function AuxiliaryG(n1::arb, n2::arb, n3::arb, p1::arb, p2::arb, lim::Int)
    two = parent(n1)(2)
    res1 = ((p1^n1) // Gamma(n1 + One(n1))) * (two^(n2+n3+One(n1))) * Beta(n2+One(n2), n3+One(n3))
    res2 = ExpIntegralE(-(n2 + n3 + One(n1)), p2)
    res3 = AuxiliaryGl(n1, 0, n2, n3, p1, p2, lim)
    res4 = AuxiliaryGl(n1, 0, n3, n2, p1, p2, lim)
    res1 * res2 - res3 - res4
end

function AuxiliaryGrlm(n1::arb, n2::arb, n3::arb, p1::arb, p2::arb, p3::arb)
    two = parent(n1)(2)
    res1 = ((p1^n1) // Gamma(n1 + One(n1))) * (two^(n3 + One(n1))) * Exp(-p2) * Exp(p3)
    res2 = AuxiliaryGm(n2, n3, two*p3)
    res1*res2
end

function AuxiliaryGrlp(n1::arb, n2::arb, n3::arb, p1::arb, p2::arb, p3::arb)
    two = parent(n1)(2)
    res1 = ((p1^n1) // Gamma(n1 + One(n1))) * (two^(n3+One(n1))) * Exp(p2) * Exp(p3)
    res2 = AuxiliaryGm(n2, n3, two*p3)
    res1*res2
end

function AuxiliaryGrh(n1::arb, n2::arb, n3::arb, p1::arb, p2::arb, p3::arb, prec::Int)
    two = parent(n1)(2)
    res1 = ((p1^n1) // Gamma(n1 + One(n1))) * (two^(n2 + n3 + One(n2 + n3)))
    res2 = Gamma(n2 + One(n2)) * Gamma(n3 + One(n3))
    res3 = Hypergeometric1F1(n2 + One(n2), n2 + n3 + two, -two * p3, prec)
    res = Exp(p3 - p2) * res1 * res2 * res3
end
###########################################################################################################


text:
function Base.show(io::IO, est::EstimationStats{D, T}) where {D, T}
    println(io, "--------------------------------------------------------")
    println(io, "Distribution:                 ", D)
    println(io, "Number of observations:       ", Int(est.num_obs))
    println(io, "Number of unknown parameters: ", Int(est.np))
    println(io, "Log-likelihood:               ", @sprintf("%.4f", est.loglikelihood))
    println(io, "AIC:                          ", @sprintf("%.4f", est.aic))
    println(io, "BIC:                          ", @sprintf("%.4f", est.bic))
    print_coefs_stats(est.coefs_stats)
    return nothing
end

function print_coefs_stats(coefs_stats::CoefsStats{T}) where T
    println("--------------------------------------------------------")
    println("Parameter      Estimate   Std.Error     t stat   p-value")
    offset = 1
    for i in coefs_stats.unknowns.ω
        p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
        p = build_print(p_c, p_std, p_t_stat, p_p_val, "omega_$i")
        println(p)
        offset += 1
    end
    sorted_keys_A = sort(collect(keys(coefs_stats.unknowns.A)))
    for k in sorted_keys_A
        for i in coefs_stats.unknowns.A[k]
            p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
            ind = round(Int, sqrt(i))
            p = build_print(p_c, p_std, p_t_stat, p_p_val, "A_$(k)_$ind$ind")
            println(p)
            offset += 1
        end
    end
    sorted_keys_B = sort(collect(keys(coefs_stats.unknowns.B)))
    for k in sorted_keys_B
        for i in coefs_stats.unknowns.B[k]
            p_c, p_std, p_t_stat, p_p_val = print_coefs_sta(coefs_stats, offset)
            ind = round(Int, sqrt(i))
            p = build_print(p_c, p_std, p_t_stat, p_p_val, "B_$(k)_$ind$ind")
            println(p)
            offset += 1
        end
    end
    return nothing
end

function print_coefs_sta(sta, offset::Int)
    p_c      = @sprintf("%.4f", sta.coefs[offset])
    p_std    = @sprintf("%.4f", sta.std_errors[offset])
    p_t_stat = @sprintf("%.4f", sta.t_stat[offset])
    p_p_val  = @sprintf("%.4f", sta.p_values[offset])
    return p_c, p_std, p_t_stat, p_p_val
end

function build_print(p_c, p_std, p_t_stat, p_p_val, param)
    p = ""
    p *= param
    p *= " "^max(0, 23 - length(p_c) - length(param))
    p *= p_c
    p *= " "^max(0, 12 - length(p_std))
    p *= p_std
    p *= " "^max(0, 11 - length(p_t_stat))
    p *= p_t_stat
    p *= " "^max(0, 10 - length(p_p_val))
    p *= p_p_val
    return p
end

ext: jl


text:
@benchgroup "get_n_words" ["IO", "Pfam"] begin

    ascii_str = "#=GR O31698/18-71 SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH\n"
    utf8_str  = "#=GF CC   (Römling U.  and Galperin M.Y. “Bacterial cellulose\n"

    @bench "ascii" get_n_words($ascii_str, 3)
    @bench "utf8"  get_n_words($utf8_str,  3)
end


text:
module ThreadingTools
include("mapreduce.jl")
end # module


text:
using PSCOPF

using .PSCOPFFixtures

using Test
using JuMP
using Dates
using DataStructures
using Printf

@testset verbose=true "test_bilevel_pilotables_start" begin

    TS = [DateTime("2015-01-01T11:00:00")]
    ech = DateTime("2015-01-01T07:00:00")
    next_ech = DateTime("2015-01-01T07:30:00")
    function create_instance(load_1, load_2,
                            market_decision_prod1, market_decision_prod2,
                            limit::Float64=35.,
                            logs=nothing)
        network = PSCOPFFixtures.network_2buses(limit=limit)
        # Pilotables
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_1", "prod_1_1", PSCOPF.Networks.PILOTABLE,
                                                20., 200.,
                                                1000., 10.,
                                                Dates.Second(10), Dates.Second(0))
        PSCOPF.Networks.add_new_generator_to_bus!(network, "bus_2", "prod_2_1", PSCOPF.Networks.PILOTABLE,
                                                20., 200.,
                                                5000., 50.,
                                                Dates.Second(10), Dates.Second(0))

        # Uncertainties
        uncertainties = PSCOPF.Uncertainties()
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_1", TS[1], "S1", load_1)
        PSCOPF.add_uncertainty!(uncertainties, ech, "bus_2", TS[1], "S1", load_2)

        gen_initial_state = SortedDict{String,PSCOPF.GeneratorState}(
                                "prod_1_1" => PSCOPF.ON,
                                "prod_2_1" => PSCOPF.ON,
                            )

        context = PSCOPF.PSCOPFContext(network, TS, PSCOPF.PSCOPF_MODE_1,
                                    gen_initial_state,
                                    uncertainties, nothing, logs)
        if !ismissing(market_decision_prod1)
            PSCOPF.set_commitment_value!(context.market_schedule, "prod_1_1", TS[1], "S1", market_decision_prod1)
        end
        if !ismissing(market_decision_prod2)
            PSCOPF.set_commitment_value!(context.market_schedule, "prod_2_1", TS[1], "S1", market_decision_prod2)
        end

        return context
    end

    @testset "importance_of_input_unit_states_for_tso_bilevel" begin
        #=
        These tests show the importance of the preceding ech decisions as it will affect the TSO decisions
        Here, we illustrate different starting cases.
        Since the preceding modules may have been launched on different uncertainties,
         we can consider arbitrary starting scenarios.
        This proves the importance of the hypothesis that a TSOBilevel model considers a balanced state as input.

        The current model does not allow the TSO to make efficient starting decisions.
         In addition, its limitation decisions affect the performance of the markets.

        The only thing that changes in the testcases is the preceding decision on the unit's start state.
        (these are the decisions at ech-1 for ts and not the generator initial state at ts-1)
        =#

        #=
        TS: [11h]
        S: [S1]
                            bus 1                   bus 2
                            |                      |
        (limitable) prod_1_1|       "1_2"          |prod_2_1
        Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
        Csta=1k, Cprop=10   |                      |Csta=5k, Cprop=50
                            |----------------------|
                            |         35           |
                            |                      |
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 30              |                      | S1: 30

        Suppose market started prod_1_1 and prod_2_1 at the preceding ech.

        There are no RSO constraints
        => TSO does not have to undertake an action
        => TSO keeps the started units at their bounds
        => prod_1_1 : 20 -> 200
           prod_2_1 : 20 -> 200

        Bilevel Market decision would be :
            prod_1_1 : 40 (economically preferred)
            prod_2_1 : 20 (to respect imposed minimum)

        However, the next market module decision without the limitation constraints would have been:
            prod_1_1 : 60 (economically preferred)
            prod_2_1 : 0 (shut down)
            => the limitation constraints lead to higher costs for no RSO reason

        FIXME? we may want the TSO to shutdown one of the units
        =#
        @testset "tso_avoids_shutting_down_units" begin

            context = create_instance(30., 30.,
                                    PSCOPF.ON, PSCOPF.ON,
                                    35.,"start_imp_test")

            tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
            firmness = PSCOPF.compute_firmness(tso,
                                                ech, next_ech,
                                                TS, context)
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            #TSO RSO constraints are OK
            @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #Market EOD constraints are OK
            @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #TSO does not need to bound pilotables because there is no risk of breaking RSO constraints
            #prod_1_1 is not bound
            @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
            #prod_2_1 is not bound
            @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

            #Market chooses the levels respecting the bounds for pilotable production
            @test 40. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
            @test 20. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

            @test - 1e-09 < objective_value(result.upper.model) < 1e-09

        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1                   bus 2
                            |                      |
        (limitable) prod_1_1|                      |prod_2_1
        Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
        Csta=1k, Cprop=10   |        "1_2"         |Csta=5k, Cprop=50
                            |----------------------|
                            |         35           |
                            |                      |
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 30              |                      | S1: 30

        Suppose market did not start any of prod_1_1 nor prod_2_1 (first TSO launch ?)
        The EOD constraint would be violated
        but there are no RSO constraints
        => TSO does not have to undertake any actions cause EOD is not his responsibility
        => TSO keeps the units shut down

        COST for TSO : 0. (no limitations cause units were already shutdown, no LoL for RSO reasons)
        However, if the TSO had started any of the units he would have fixed its minimum production,
        Thus, he would have paid for it.

        PROBLEM : the following market will not be able to start the units due to the impositions
        => prod_1_1 : 0
           prod_2_1 : 0
           But the market will need to cut all the conso
           market cost : LoL(60MW)

        However, the next market module decision without the limitation constraints would have been:
            prod_1_1 : 60 (economically preferred)
            prod_2_1 : 0 (shut down)
            => the limitation constraints lead to LoL for no reason

        FIXME?
        we want TSO to start prod_1_1 without limiations (i.e 20-200) => costs 20. for tso
        =#
        @testset "tso_bilevel_does_not_start_units_for_eod_reasons_from_off" begin
            #no_risk_of_breaking_rso_constraint_off_off

            context = create_instance(30., 30.,
                                    PSCOPF.OFF, PSCOPF.OFF,
                                    35.,"start_imp_test")

            tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=true))
            firmness = PSCOPF.compute_firmness(tso,
                                                ech, next_ech,
                                                TS, context)
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)

            # Solution status
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            #TSO RSO constraints are OK
            @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #Market
            @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #TSO impositions
            #prod_1_1
            @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
            #prod_2_1
            @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
            @test value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"]) < 1e-09

            #Market chooses the levels respecting the bounds for pilotable production
            @test 60. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
            @test value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"]) < 1e-09

            @test value(PSCOPF.get_upper_obj_expr(result)) ≈ 20. * 10. #imposition cost of prod_1_1
            @test value(PSCOPF.get_lower_obj_expr(result)) ≈ (60 * 10.)

        end

        #=
        same as above :

        TS: [11h]
        S: [S1]
                            bus 1                   bus 2
                            |                      |
        (limitable) prod_1_1|       "1_2"          |prod_2_1
        Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
        Csta=1k, Cprop=10   |                      |Csta=5k, Cprop=50
                            |----------------------|
                            |         35           |
                            |                      |
                            |                      |
               load(bus_1)  |                      |load(bus_2)
        S1: 30              |                      | S1: 30

        Suppose we still haven't decided the state of prod_1_1 nor prod_2_1 (first TSO launch ?)
        This is similar to the preceding since we have the same cost expressions for missing or OFF units

        The EOD constraint would be violated
        but there are no RSO constraints
        => TSO does not have to undertake an action cause EOD is not his responsibility
        => he keeps the units shut down
        PROBLEM : the following market will not be able to start the units due to the impositions
        => prod_1_1 : 0
           prod_2_1 : 0
           These cost 0.
           But the market will need to cut all the conso

        If the TSO had started any of the units he would have fixed its minimum production,
        Thus, he would have paid for it.

        FIXME? issue if first tso bilevel (mode 3)
        =#
        @testset "tso_bilevel_does_not_start_units_for_eod_reasons_from_missing" begin
            # no_risk_of_breaking_rso_constraint_missing_missing
            context = create_instance(30., 30.,
                                    missing, missing,
                                    35.,"start_imp_test")

            tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=true))
            firmness = PSCOPF.compute_firmness(tso,
                                                ech, next_ech,
                                                TS, context)
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)

            # Solution status
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            #TSO RSO constraints are OK
            @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #Market
            @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #TSO impositions
            #prod_1_1
            @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
            #prod_2_1
            @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
            @test value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"]) < 1e-09

            #Market chooses the levels respecting the bounds for pilotable production
            @test 60. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
            @test value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"]) < 1e-09

            @test value(PSCOPF.get_upper_obj_expr(result)) ≈ 20. *10. #imposition of prod_1_1
            @test value(PSCOPF.get_lower_obj_expr(result)) ≈ (60 * 10.)

        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1                   bus 2
                            |                      |
        (limitable) prod_1_1|       "1_2"          |prod_2_1
        Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
        Csta=1k, Cprop=10   |                      |Csta=5k, Cprop=50
                            |----------------------|
                            |         35           |
                            |                      |
                            |                      |
            load(bus_1)     |                      |load(bus_2)
        S1: 30              |                      | S1: 30

        Suppose market only started prod_1_1
        There are no RSO constraints
        => TSO does not have to undertake an action
        => he keeps the unit prod_2_1 shutdown
        => prod_1_1 : 20 -> 200 : cost 0.
           prod_2_1 : 0 : cost 0.
           TSO cost : 0.

        The Bilevel market will use prod_1_1 to produce 60MW

        The same behaviour will be observed by the following market module
        =#
        @testset "no_risk_of_breaking_rso_constraint_on_off" begin

            context = create_instance(30., 30.,
                                    PSCOPF.ON, PSCOPF.OFF,
                                    35.,"start_imp_test")

            tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
            firmness = PSCOPF.compute_firmness(tso,
                                                ech, next_ech,
                                                TS, context)
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            #TSO RSO constraints are OK
            @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #Market EOD constraints are OK
            @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #TSO impositions
            #prod_1_1
            @test 20. ≈value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
            #prod_2_1
            @test value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"]) < 1e-09
            @test value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"]) < 1e-09

            #Market chooses the levels respecting the bounds for pilotable production
            @test 60. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
            @test value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"]) < 1e-09

            @test - 1e-09 < objective_value(result.upper.model) < 1e-09

        end

        #=
        TS: [11h]
        S: [S1]
                            bus 1                   bus 2
                            |                      |
        (limitable) prod_1_1|       "1_2"          |prod_2_1
        Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
        Csta=1k, Cprop=10   |                      |Csta=5k, Cprop=50
                            |----------------------|
                            |         35           |
                            |                      |
                            |                      |
            load(bus_1)     |                      |load(bus_2)
        S1: 30              |                      | S1: 30

        Suppose market only started prod_2_1
        There are no RSO constraints
        => TSO does not have to undertake an action
        => TSO keeps the unit prod_1_1 shutdown
        => prod_1_1 : 0 : cost 0.
           prod_2_1 : 20 -> 200 : cost 0.
           These cost 0.

        If TSO had imposed prod_1_1's starting, TSO would have paid : 20
        =#
        @testset "tso_follows_preceding_market_starts" begin

            context = create_instance(30., 30.,
                                    PSCOPF.OFF, PSCOPF.ON,
                                    35.,"start_imp_test")

            tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
            firmness = PSCOPF.compute_firmness(tso,
                                                ech, next_ech,
                                                TS, context)
            result = PSCOPF.run(tso, ech, firmness,
                        PSCOPF.get_target_timepoints(context),
                        context)

            # Solution is optimal
            @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

            #TSO RSO constraints are OK
            @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #Market EOD constraints are OK
            @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
            @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

            #TSO does not need to bound pilotables
            #prod_1_1 is not bound
            @test value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"]) < 1e-09
            @test value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"]) < 1e-09
            #prod_2_1 is not bound
            @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
            @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

            #Market chooses the levels respecting the bounds for pilotable production
            @test value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"]) < 1e-09
            @test 60. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

            @test objective_value(result.upper.model)  < 1e-09

        end

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 100

    Market chooses prod_1_1 cause cheaper than prod_2_1
    => a flow of 100. => RSO problem

    TSO needs to impose using prod_2_1. two options:

    option 1 : impose prod_1_1 to produce at most 10+35=45 => P_1_1 in [0,45]
        If the TSO chooses to shut down prod_1_1 entirely i.e. [O, O] it will pay :
            => cost (unit was ON): (0-20) + (200-0) = -20 + 200 = 180
        If the TSO chooses to allow using prod_1_1 i.e. [20, 45] it will pay :
            => cost (unit was ON): (20-20) + (200-45) = 155
    (*) If the TSO chooses to allow using prod_1_1 freely i.e. [20, 200] and use option 2 it will pay :
            => cost (unit was ON): (20-20) + (200-200) = 0.    + cost(option 2)

    option 2 : impose prod_2_1 to produce at least 100-35=65 => P_2_1 in [65,200] => cost (unit was OFF): 65

    Imposing one of the preceding conditions guarantees RSO constraint
        => TSO will impose the cheaper one => option 1

    FIXME Note: this might not be what we want, the TSO could have shut bus1 completely and only used bus2
    =#
    @testset "rso_constraint_requires_setting_pilotables_bounds_min" begin

        context = create_instance(10., 100.,
                                PSCOPF.ON, PSCOPF.OFF,
                                35.,"pilotable_test")

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_OPTIMAL

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test 65. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels sets the bounds for pilotable production
        @test 45. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 65. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test 65. ≈ objective_value(result.upper.model)

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         35           |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 10            |                      | S1: 300

    Suppose market starts prod_1_1 and prod_2_1
        because prod_1_1 is cheaper, it will be used at max capa i.e. 200.
        This will cause an RSO constraint on branch 1_2

    TSO needs to solve the RSO constraint by
    option 1 : limit using prod_1_1 to 10+35=45
        => limitng prod_1_1 at [20, 45]
        prod_1_1 : 20 -> 45, cost : 200-45 = 155
        prod_2_1 : 20 -> 200, cost : 0.
        => costs 155
    option 2 : cut conso on bus 2 (expensive)

    option 1 is adopted.
    This however will cause an EOD disbalance in market
     cause demand=310, capacity=245
     => market needs to loss_of_load by 65
    =#
    @testset "rso_constraint_requires_setting_pilotables_bounds_max" begin

        context = create_instance(10., 300.,
                                PSCOPF.ON, PSCOPF.ON,
                                35.,"pilotable_test")

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints are OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market EOD constraints are OK
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 65. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        # RSO locates cut conso on bus 2 to assure flow limit
        @test value(result.upper.lol_model.p_loss_of_load["bus_1", TS[1],"S1"]) < 1e-09
        @test 65. ≈ value(result.upper.lol_model.p_loss_of_load["bus_2", TS[1],"S1"])

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
        @test 45. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels sets the bounds for pilotable production
        @test 45. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

        @test objective_value(result.upper.model) ≈ ( 155.
                                                        + 65. *tso.configs.TSO_LOL_PENALTY #market LoL
                                                    )

    end

    #=
    TS: [11h]
    S: [S1]
                        bus 1                   bus 2
                        |                      |
    (limitable) prod_1_1|       "1_2"          |prod_2_1
    Pmin=20, Pmax=200   |                      |Pmin=20, Pmax=200
    Csta=0, Cprop=10    |                      |Csta=0, Cprop=50
                        |----------------------|
                        |         500          |
                        |                      |
                        |                      |
           load(bus_1)  |                      |load(bus_2)
      S1: 400           |                      | S1: 20

    No risk of RSO constraint breaking, but EOD cannot be satisfied
        => the market needs to cut conso
        => the TSO will locate the bus where we cut
    =#
    @testset "eod_problem_requires_cutting_conso" begin

        context = create_instance(400., 20.,
                                PSCOPF.ON, PSCOPF.ON,
                                500.)

        tso = PSCOPF.TSOBilevel(PSCOPF.TSOBilevelConfigs(USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST=false))
        firmness = PSCOPF.compute_firmness(tso,
                                            ech, next_ech,
                                            TS, context)
        result = PSCOPF.run(tso, ech, firmness,
                    PSCOPF.get_target_timepoints(context),
                    context)

        # Solution is optimal
        @test PSCOPF.get_status(result) == PSCOPF.pscopf_HAS_SLACK

        #TSO RSO constraints can be respected without capping or losing load OK
        @test value(result.upper.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test value(result.upper.lol_model.p_global_loss_of_load[TS[1],"S1"]) < 1e-09

        #Market needs to cut conso to assure EOD constraint
        @test value(result.lower.limitable_model.p_global_capping[TS[1],"S1"]) < 1e-09
        @test 20. ≈ value(result.lower.lol_model.p_global_loss_of_load[TS[1],"S1"])

        # RSO locates cut conso : here any combination will do
        @test 20. ≈ ( value(result.upper.lol_model.p_loss_of_load["bus_1", TS[1],"S1"])
                    + value(result.upper.lol_model.p_loss_of_load["bus_2", TS[1],"S1"]) )

        #TSO sets the bounds for pilotable production respecting units' pmin and pmax
        #prod_1_1 is not bound
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_1_1",TS[1],"S1"])
        #prod_2_1 bounding + EOD, guarantees RSO constraints
        @test 20. ≈ value(result.upper.pilotable_model.p_imposition_min["prod_2_1",TS[1],"S1"])
        @test 200. ≈ value(result.upper.pilotable_model.p_imposition_max["prod_2_1",TS[1],"S1"])

        #Market chooses the levels sets the bounds for pilotable production
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_1_1",TS[1],"S1"])
        @test 200. ≈ value(result.lower.pilotable_model.p_injected["prod_2_1",TS[1],"S1"])

    end

end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using Revise
using DungBase, UUIDs
import DungBase:AbstractTimeLine, Temporal, Board, Calibration, Run, Experiment
dbpath = "/home/yakir/dungProject/database"
cofeesource = "/home/yakir/mnt/coffee source"
using DataDeps, JuliaDB, Dates
register(DataDep("coffeebeetle", "the coffee beetle database", "https://s3.eu-central-1.amazonaws.com/vision-group-file-sharing/Data%20backup%20and%20storage/Yakir/coffee%20beetles/database.zip", "557ec5027d34d2641db4f2cc61c993edd2ce2bee1530254d45d26cc41c5adf7a", post_fetch_method = unpack))
video = loadtable(datadep"coffeebeetle/video.csv", indexcols = :video)
videofile = loadtable(datadep"coffeebeetle/videofile.csv", indexcols = :file_name)
x = join(video, videofile, rkey = :video)
function AbstractTimeLine(x)
    n = length(x)
    if n == 1
        y = x[]
        comment = ismissing(y.comment) ? "" : y.comment
        WholeVideo(VideoFile(y.file_name, y.date_time, Nanosecond(y.duration)), comment)
    else
        files = Vector{VideoFile}(undef, length(x))
        for i in rows(x)
            files[i.index] = VideoFile(i.file_name, i.date_time, Nanosecond(i.duration))
        end
        y = x[1]
        comment = ismissing(y.comment) ? "" : y.comment
        FragmentedVideo(files, comment)
    end
end
timeline = groupby(AbstractTimeLine,  x)
interval = dropmissing(loadtable(datadep"coffeebeetle/interval.csv", indexcols = :interval), :start)
x = join(interval, timeline, lkey = :video)
function Temporal(x)
    video = x.AbstractTimeLine
    time = AbstractPeriod(Nanosecond(x.start), x.stop)
    Temporal(video, time, ismissing(x.comment) ? "" : x.comment)
end
temporal = Dict(i.interval => Temporal(i) for i in x)

board = loadtable(datadep"coffeebeetle/board.csv", indexcols = :designation)
calibration = loadtable(datadep"coffeebeetle/calibration.csv", indexcols = :calibration)
x = join(calibration, board, lkey = :board)
function Calibration(x)
    intrinsic = get(temporal, x.intrinsic, missing)
    extrinsic = get(temporal, x.extrinsic, missing)
    board = Board(x.board, x.checker_width_cm, (x.checker_per_width, x.checker_per_height), x.board_description)
    comment = ismissing(x.comment) ? "" : x.comment
    Calibration(UUID(x.calibration), intrinsic, extrinsic, board, comment)
end
calibration = Dict(i.calibration => Calibration(i) for i in x)

connector = filter(!isempty, loadtable(datadep"coffeebeetle/poi.csv", indexcols = :poi), select = :calibration)
using CSV, Tables
x = CSV.File(datadep"coffeebeetle/run.csv", types = Dict(:date => Date)) |> rowtable
fname = setdiff(keys(x[1]), (:run, :date, :experiment, :comment))
x = map(x) do i
    (run = i.run, date = i.date, experiment = i.experiment, comment = i.comment, factors = Dict(s => getfield(i, s) for s in fname))
end
x = table(x, pkey = :run)
experiment = loadtable(datadep"coffeebeetle/experiment.csv", indexcols = :experiment)
x = join(x, experiment, lkey = :experiment)
x = join(connector, x, lkey = :run, rkey = :run)
function Run(x)
    pois = Dict(Symbol(i.type) => POI(UUID(i.poi), calibration[i.calibration], temporal[i.interval]) for i in x)
    setup = x[1].factors
    date = x[1].date
    comment = ismissing(x[1].comment) ? "" : x[1].comment
    Run(pois, setup, date, comment)
end
function Experiment(x)
    runs = select(groupby(Run, table(x), :run), :Run)
    Experiment(runs, x[1].experiment_description)
end
t = groupby(Experiment, x, :experiment)

data = Dict(i.experiment => i.Experiment for i in t)

using Serialization
serialize(joinpath(datadep"coffeebeetle", "data"), data)

using Revise
using Serialization
using DungBase, Dates, UUIDs
using DataDeps
register(DataDep("coffeebeetle", "the coffee beetle database", "https://s3.eu-central-1.amazonaws.com/vision-group-file-sharing/Data%20backup%20and%20storage/Yakir/coffee%20beetles/database.zip", "557ec5027d34d2641db4f2cc61c993edd2ce2bee1530254d45d26cc41c5adf7a", post_fetch_method = unpack))
data = deserialize(joinpath(datadep"coffeebeetle", "data"))

dbpath = "/home/yakir/dungProject/database"
for (k,v) in data, (i, run) in enumerate(v.runs)
    path = joinpath(dbpath, "experiments", k, string(i))
    mkpath(path)
end

t = data["nest#closed person#therese"].runs[1].pois[:nest].calib.intrinsic
imgs, sar = DungBase.fetchimages(t)

using Images, CoordinateTransformations, StaticArrays

tfm = LinearMap(SDiagonal(SVector(sar, 1)))
img = imgs[1]
imgw = warp(img, tfm);
imshow(imgw)

pixel = Dict{
for file in readdir(datadep"coffeebeetle/pixel")

ENV["FFREPORT"] = "8"
using VideoIO

file = download("https://s3.eu-central-1.amazonaws.com/vision-group-file-sharing/Fun%20Stuff/a.MTS")

file = "/home/yakir/videos/scallop/Swimming Scallops-KQHg6M7-Du0.mp4"
f = openvideo(file)
seek(f, 3.0)

read(f)

file = "/home/yakir/Brooklyn Nine-Nine Season 5/Brooklyn Nine-Nine S05E01.mkv"
f = openvideo(file)
sar = f.aspect_ratio
if sar ≠ 0 && sar ≠ 1
    <deal with a freaky SAR>
end

playvideo(f, pixelaspectratio=1.0)

While the Data Aspect Ratio (DAR) or Pixel Aspect Ration (PAR) can be safely assumed to be equal to one (in modern times this means the shape of the pixels is square), in some cases this aspect ratio may be smaller or larger than one. In such cases, an image will look unnaturally wide/tall. To retrieve the reported aspect ratio form a file use:
```
f = openvideo(file)
sar = f.aspect_ratio
```
Note that not all videos store a correct apsect ratio. In such umbigous cases, `ffmpeg` will report a ratio of `0//1`. To help identify a truly ≠1 ratio use:
```
if sar ≠ 1 && sar ≠ 0
    <...>
end
```


text:
"""
Find all types that may be considered ansestors of this type.
For number types it is the same as Julia type hierarchy.
For array types it includes unparametrized versions of types, i.e.:

    type_ansestors(Vector{Float64})
    # ==> [Array{Float64,1}, Array{T,1}, DenseArray{T,1},
           AbstractArray{T,1}, AbstractArray{T,N}, Any]
"""
function type_ansestors(t::Type{T}) where T<:Number
    types = Type[]
    while t != Any
        push!(types, t)
        t = supertype(t)
    end
    push!(types, Any)
    return types
end

type_ansestors(t::Type{Vector{T}}) where {T} =
    [t, Vector, DenseVector, AbstractVector, AbstractArray, Any]
type_ansestors(t::Type{Matrix{T}}) where {T} =
    [t, Matrix, DenseMatrix, AbstractMatrix, AbstractArray, Any]

function bcast_to_call(pex::Expr)
    @assert pex.head == :(.)
    return Expr(:call, pex.args[1], pex.args[2].args...)
end

deriv_name(z::Symbol, x::Symbol) = Symbol("d$(z)!d$(x)")
split_deriv_name(vname) = Symbol.(split(String(vname), "!"))

function replace_node(g::AbstractExGraph, vname::Symbol, nds::Vector{ExNode})
    i = indexof(g, vname)
    delete!(g, vname)
    insert!(g, nds)
end

function find_related(g::AbstractExGraph, dydx_v::Symbol)
    subderivs = Symbol[]
    i = 1
    name = Symbol("$(dydx_v)__$(i)")
    while haskey(g, name)
        push!(subderivs, name)
        i += 1
        name = Symbol("$(dydx_v)__$(i)")
    end
    return subderivs
end

# (symbolic) derivative size propagation

const DERIV_NAME_PATTERN = r"(d.+)!(d.+)"

function propagate_deriv_size!(g::AbstractExGraph, dd_name::Symbol)
    sizes = @get_or_create(g.ctx, :sizes, Dict())
    rg = match(DERIV_NAME_PATTERN, String(dd_name))
    @assert length(rg.captures) == 2
    str_dnames = rg.captures
    zname = Symbol(str_dnames[1][2:end])
    xname = Symbol(split(str_dnames[2][2:end], "__")[1]) # cut down `__$(i)` part if any
    zsize, xsize = (sizes[zname], sizes[xname])
    if zsize == :(())
        # output var is constant
        sizes[dd_name] = xsize
    else
        sizes[dd_name] = :(($zsize..., $xsize...)) |> simplify
    end
end

function propagate_deriv_size!(g::AbstractExGraph)
    for nd in g.tape
        vname = varname(nd)
        if match(DERIV_NAME_PATTERN, String(vname)) != nothing
            propagate_deriv_size!(g, vname)
        end
    end
end

# (numeric) derivative size propagation

function infer_deriv_size!(g::AbstractExGraph, dd_name::Symbol)
    rg = match(DERIV_NAME_PATTERN, String(dd_name))
    @assert length(rg.captures) == 2
    str_dnames = rg.captures
    zname = Symbol(str_dnames[1][2:end])
    xname = Symbol(split(str_dnames[2][2:end], "__")[1]) # cut down `__$(i)` part if any
    # in case z or x haven't been evaluated and their size isn't known yet
    evaluate!(g, zname)
    evaluate!(g, xname)
    sizes = @get_or_create(g.ctx, :rsizes, Dict())
    zsize, xsize = (sizes[zname], sizes[xname])
    sizes[dd_name] = (zsize..., xsize...)
end

function infer_deriv_size!(g::AbstractExGraph)
    for nd in g.tape
        vname = varname(nd)
        if match(DERIV_NAME_PATTERN, String(vname)) != nothing
            infer_deriv_size!(g, vname)
        end
    end
end

# top type

"The top type describing given data"
top_type(x::AbstractArray{T,N}) where {T,N} = AbstractArray{T,N}
top_type(x::Number) = Number

top_type(::Type{AT}) where {AT <: AbstractArray{T,N}} where {T,N} = AbstractArray{T,N}
top_type(::Type{T}) where {T <: Number} = Number


text:
using Plots

function show!(p::Plots.Plot, boundary::StraightBoundary, meshdef::MeshDef)
    vert_start = boundary.vertice_start
    vert_end = boundary.vertice_end
    coords = meshdef.vertices[[vert_start, vert_end]]

    plot!(p, [coords[1][1], coords[2][1]], [coords[1][2], coords[2][2]], linecolor="blue")
end

function show!(p::Plots.Plot, boundary::CircularBoundary, meshdef::MeshDef)
    vert_start = boundary.vertice_start
    vert_end = boundary.vertice_end
    coords = meshdef.vertices[[vert_start, vert_end]]

    center = circlecenter(coords..., boundary.radius)

    leg1 = coords[1] .- center
    leg2 = coords[2] .- center

    γ = angle(vec(leg1), vec(leg2))

    x(α) = leg1[1] * cos(α) - leg1[2] * sin(α) + center[1]
    y(α) = leg1[1] * sin(α) + leg1[2] * cos(α) + center[2]

    plot!(p, x, y, 0, γ, leg=false, linecolor="blue")
end

function show!(p::Plots.Plot, vertice::Vector{Float64})
    scatter!(p, [vertice[1]], [vertice[2]], marker=(:circle, :red, 6))
end

function show!(p::Plots.Plot, element::Vector{Integer}, mesh::Mesh)
    start_node = mesh.nodes[element[1]]

    for node_id in element[2:end]

        cur_node = mesh.nodes[node_id]

        plot!(p, [start_node[1], cur_node[1]], [start_node[2], cur_node[2]], linecolor="grey")

        start_node = cur_node
    end

    last_node = mesh.nodes[element[1]]
    plot!(p, [start_node[1], last_node[1]], [start_node[2], last_node[2]], linecolor="grey")

end

function show!(p::Plots.Plot, mesh::Mesh)

    # TODO vectorize!
    for node in mesh.nodes
        scatter!(p, [node[1]], [node[2]], marker=(:circle, :blue, 3))
    end

    for element in mesh.elements
        show!(p, element, mesh)
    end
end

function show(meshdef::MeshDef; hideVertices::Bool=false)
    p = plot(legend=false)

    for boundary in meshdef.bounds
        show!(p, boundary, meshdef)
    end

    mesh = emptyMesh()
    meshBoundaries!(mesh, meshdef)

    show!(p, mesh)

    if !hideVertices
        for vertice in meshdef.vertices
            show!(p, vertice)
        end
    end

    p
end

function show(mesh::Mesh)
    p = plot(legend=false)

    show!(p, mesh)

    p
end


text:
"""`newton1d` performs the interval Newton method on the given function `f`
with its derivative `f′` and initial interval `x`.
Optional keyword arguments give the tolerances `reltol` and `abstol`.
`reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|,
and a `debug` boolean argument that prints out diagnostic information."""
function newton1d(f::Function, f′::Function, x::Interval{T};
                    reltol=eps(T), abstol=eps(T), debug=false, debugroot=false) where {T}

    L = Interval{T}[] # Array to hold the intervals still to be processed

    R = Root{Interval{T}}[] # Array to hold the `root` objects obtained
    reps = reps1 = 0

    push!(L, x) # Initialize
    initial_width = ∞
    X = emptyinterval(T) # Initialize
    while !isempty(L)   # Until all intervals have been processed
        X = pop!(L)     # Process next interval

        debug && (print("Current interval popped: "); @show X)

        m = mid(X)
        if (isempty(X))
            continue
        end

        if 0 ∉ f′(X)    # if 0 ∉ f′(X), Newton steps can be made normally until either X becomes empty, or is known to contain a unique root

            debug && println("0 ∉ f′(X)")

            while true

                m = mid(X)
                N = m - (f(interval(m)) / f′(X))    # Newton step

                debug && (print("Newton step1: "); @show (X, X ∩ N))
                if X == X ∩ N   # Checking if Newton step was redundant
                    reps1 += 1
                    if reps1 > 20
                        reps1 = 0
                        break
                    end
                end
                X = X ∩ N

                if (isempty(X))         # No root in X
                    break

                elseif 0 ∈ f(wideinterval(mid(X)))  # Root guaranteed to be in X
                    n = fa = fb = 0
                    root_exist = true
                    while (n < 4 && (fa == 0 || fb == 0))   # Narrowing the interval further
                        if fa == 0
                            if 0 ∈ f(wideinterval(X.lo))
                                fa = 1
                            else
                                N = X.lo - (f(interval(X.lo)) / f′(X))
                                X = X ∩ N
                                if (isempty(X))
                                    root_exist = false
                                    break
                                end
                            end
                        end
                        if fb == 0
                            if 0 ∈ f(wideinterval(X.hi))
                                fb = 1
                            else
                                if 0 ∈ f(wideinterval(mid(X)))
                                    N = X.hi - (f(interval(X.hi)) / f′(X))
                                else
                                    N = mid(X) - (f(interval(mid(X))) / f′(X))
                                end
                                X = X ∩ N
                                if (isempty(X))
                                    root_exist = false
                                    break
                                end
                            end
                        end
                        N = mid(X) - (f(interval(mid(X))) / f′(X))
                        X = X ∩ N
                        if (isempty(X))
                            root_exist = false
                            break
                        end
                        n += 1
                    end
                    if root_exist
                        push!(R, Root(X, :unique))
                        debugroot && @show "Root found", X  # Storing determined unique root
                    end

                    break
                end
            end

        else                                # 0 ∈ f′(X)
            if diam(X) == initial_width     # if no improvement occuring for a number of iterations
                reps += 1
                if reps > 10
                    push!(R, Root(X, :unknown))
                    debugroot && @show "Repeated root found", X
                    reps = 0
                    continue
                end
            end
            initial_width = diam(X)
            debug && println("0 ∈ f'(X)")

            expansion_pt = Inf
            # expansion point for the newton step might be m, X.lo or X.hi according to some conditions

            if 0 ∈ f(wideinterval(mid(X)))  # if root in X, narrow interval further
                # 0 ∈ fⁱ(x)

                debug && println("0 ∈ fⁱ(x)")

                if 0 ∉ f(wideinterval(X.lo))
                    expansion_pt = X.lo

                elseif 0 ∉ f(wideinterval(X.hi))
                    expansion_pt = X.hi

                else
                    x1 = mid(interval(X.lo, mid(X)))
                    x2 = mid(interval(mid(X), X.hi))
                    if 0 ∉ f(wideinterval(x1)) || 0 ∉ f(wideinterval(x2))
                        push!(L, interval(X.lo, m))
                        push!(L, interval(m, X.hi))
                        continue

                    else
                        push!(R, Root(X, :unknown))

                        debugroot && @show "Multiple root found", X

                        continue
                    end
                end

            else
                # 0 ∉ fⁱ(x)

                debug && println("0 ∉ fⁱ(x)")

                if (diam(X)/mag(X)) < reltol && diam(f(X)) < abstol     # checking if X is still within tolerances
                    push!(R, Root(X, :unknown))

                    debugroot && @show "Tolerance root found", X

                    continue
                end
            end
            # Step 8

            if isinf(expansion_pt)
                expansion_pt = mid(X)
            end

            a = f(interval(expansion_pt))
            b = f′(X)
            # Newton steps with extended division creating two intervals
            if 0 < b.hi && 0 > b.lo && 0 ∉ a
                if a.hi < 0
                    push!(L, X ∩ (expansion_pt - interval(-Inf, a.hi / b.hi)))
                    push!(L, X ∩ (expansion_pt - interval(a.hi / b.lo, Inf)))

                elseif a.lo > 0
                    push!(L, X ∩ (expansion_pt - interval(-Inf, a.lo / b.lo)))
                    push!(L, X ∩ (expansion_pt - interval(a.lo / b.hi, Inf)))

                end

                continue

            else
                N = expansion_pt - (f(interval(expansion_pt))/f′(X))

                debug && (print("Newton step2: "); @show (X, X ∩ N))

                X = X ∩ N
                m = mid(X)

                if isempty(X)
                    continue
                end
            end

            push!(L, X) # Pushing X into L to be processed again
        end
    end

    return R
end

"""`newton1d` performs the interval Newton method on the given function `f` and initial interval `x`.
Optional keyword arguments give the tolerances `reltol` and `abstol`.
`reltol` is the tolerance on the relative error whereas `abstol` is the tolerance on |f(X)|,
and a `debug` boolean argument that prints out diagnostic information."""
newton1d(f::Function, x::Interval{T}; args...) where {T} =
    newton1d(f, x->D(f,x), x; args...)


text:
"""Merge vector of UnitRanges into Vector of UnitRanges that do not overlap."""
function merge(ranges::Vector{UnitRange{Integer}})
    # prevent unwanted mutation
    ranges = copy(ranges)

    # Skip if there's nothing to merge
    length(ranges) <= 1 && (return ranges)

    merged_ranges = Vector{UnitRange{Integer}}()

    # smallest start at the end
    sort!(ranges, by=x->(x.start, x.stop), rev=true)
    while length(ranges) > 1
        r1 = pop!(ranges)

        # Overlap checking
        overlaps = r1.stop >= ranges[end].start - 1
        exceeds = r1.stop > ranges[end].stop

        # No overlap => push to merged_ranges
        if !overlaps
            push!(merged_ranges, r1)
        elseif exceeds
            ranges[end] = r1.start:r1.stop
        else
            ranges[end] = r1.start:ranges[end].stop
        end
    end
    push!(merged_ranges, ranges[end])
    merged_ranges
end


text:
@testset "318.maximum-product-of-word-lengths.jl" begin
    @test max_product(["abcw", "baz", "foo", "bar", "xtfn", "abcdef"]) == 16
    @test max_product(["a", "ab", "abc", "d", "cd", "bcd", "abcd"]) == 4
    @test max_product(["a", "aa", "aaa", "aaaa"]) == 0
end


text:
# TODO: Write documentation and comments
struct MutableNamedTuple{K,T}
    data::NamedTuple{K,T}
    function MutableNamedTuple(nt::NamedTuple)
        data = []
        for (key, val) in zip(keys(nt), values(nt))
            if val isa NamedTuple
                val = MutableNamedTuple(val)
            end
            type = val isa Function ? Function : typeof(val)
            # type = typeof(val)
            push!(data, key => type[val])
        end
        nt = (; data...)
        K,T = typeof(nt).parameters
        return new{K,T}(nt)
    end
end
MutableNamedTuple(; kwargs...) = MutableNamedTuple((;kwargs...))

Base.getproperty(mnt::MutableNamedTuple, key::Symbol) = getfield(mnt, :data)[key][1]

Base.propertynames(mnt::MutableNamedTuple) = propertynames(getfield(mnt, :data))

function Base.setproperty!(mnt::MutableNamedTuple, key::Symbol, val)
    d = getfield(mnt, :data)
    prop = getproperty(d, key)
    _setvalue!(prop, val)
    return nothing
end

# Need this function barrier for type stability
function _setvalue!(arr, val)
    arr[1] = val
    return nothing
end

Base.getindex(mnt::MutableNamedTuple, key::Symbol) = getproperty(mnt, key)

function Base.setindex!(mnt::MutableNamedTuple, val, key::Symbol)
    set_property!(mnt, key, val)
    # getfield(mnt, :data)[key][1] = val
    return nothing
end

function namedtuple(mnt::MutableNamedTuple)
    data = []
    for (key, val) in zip(keys(mnt), values(mnt))
        val = val[1]
        if val isa MutableNamedTuple
            val = namedtuple(val)
        end
        push!(data, key => val)
    end
    return (; data...)
end

function Base.show(io::IO, mnt::MutableNamedTuple)
    show(io, namedtuple(mnt))
    return nothing
end
function Base.show(io::IO, ::MIME"text/plain", mnt::MutableNamedTuple)
    print(io, "MutableNamedTuple")
    show(io, namedtuple(mnt))
    return nothing
end

Base.keys(mnt::MutableNamedTuple) = keys(getfield(mnt, :data))
Base.values(mnt::MutableNamedTuple) = values(getfield(mnt, :data))

Base.length(mnt::MutableNamedTuple) = length(getfield(mnt, :data))
Base.size(mnt::MutableNamedTuple) = size(getfield(mnt, :data))

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
# # copied from autogenerated SnoopCompile files
function _precompile_()
  nothing
end


text:
struct Layer
    name::String
    inputs::Int
    outputs::Int
    connect::Vector{Tuple{Int,Int}}
end

function Layer(layer::Any)
    if !hasproperty(layer, :W) ||!hasproperty(layer, :σ)
        error("type $(typeof(layer)) is not supported")
    end
    name = string(Base.typename(typeof(layer)))
    activation = layer.σ
    if activation !== identity
        name *= " ($activation)"
    end
    outputs, inputs = size(layer.W)
    connect = [(i, j) for i = 1:inputs for j = 1:outputs]
    Layer(name, inputs, outputs, connect)
end

Layer(layer::Function) = Layer(string(Symbol(layer)), 0, 0, [])

Layer(name::String, neurons::Int) =
    Layer(name, neurons, neurons, collect(zip(1:neurons, 1:neurons)))


text:
"""
    LearningPhase

Might be changed to SolvingPhase or SearchingPhase.

They are four phases encountered during Learning (or rather Solving).
- The InitializingPhase: before starting the search, some elements have to be initialised
in accordance to the current model.
- The StepPhase: a step is the simple fact of going from a CP State to another eather by assigning a value, removing a value.
- The DecisionPhase: the moment when a value selection heuristic takes a decision of assigning a value to a variable.
- The EndingPhase: when all nodes have been visited or a limit has been reached, this is just after the end of the search.
"""
abstract type LearningPhase end

struct InitializingPhase <: LearningPhase end
struct StepPhase <: LearningPhase end
struct DecisionPhase <: LearningPhase end
struct EndingPhase <: LearningPhase end

include("../../RL/RL.jl")
include("searchmetrics.jl")

abstract type ValueSelection end

include("classic/basicheuristic.jl")

abstract type ActionOutput end
struct VariableOutput <: ActionOutput end
struct FixedOutput <: ActionOutput end

include("learning/learning.jl")


text:
"""
    runExperiments(ns::Array{Integer}, ms::Array{Integer};  αs = collect(0.1:0.1:0.9), numberOfInstances = 5)

Entry point of experiments. This function runs experiments for minimum knapsack problem
with problem size `n` specified by the list `ns` and minimum assignment problem with problem size `m` specified
by the list `ms`. Optional argument `αs` specify a list of parameters defining
neighbourhood of some solution ``x`` and optional argument `numberOfInstances`
specify number of problem instances to be generated for each value of `alpha`.

# Examples:
```julia-repl
julia> using RobRecSolver.Experiments
julia> runExperiments([100, 400, 1000], [10, 25, 100])
```
"""
function runExperiments(ns::Array{Int}, ms::Array{Int}; αs = collect(0.1:0.1:0.9), numberOfInstances = 5)

    @info "Starting experiments for minimum knapsack problem for n in $ns"
    Δt = @elapsed runKnapsackExperiments(ns, αs = αs, numberOfInstances = numberOfInstances)
    @info "Experiments for minimum knapsack problem have finished in $(Δt)sec"

    @info "Starting experiments for minimum assignment problem for m in $ms"
    Δt = @elapsed runAssignmentExperiments(ms, αs = αs, numberOfInstances = numberOfInstances)
    @info "Experiments for minimum assignment problem have finished in $(Δt)sec"
end

"""
    runKnapsackExperiments(ns; αs = collect(0.1:0.1:0.9), numberOfInstances = 5)

Runs experiments for minimum knapsack problem.
"""
function runKnapsackExperiments(ns; αs = collect(0.1:0.1:0.9), numberOfInstances = 5)
    for n in ns
        @info "Starting minimum knapsack problem experiment for n = $(n)..."

        srand(310923)

        problemDescriptor = KnapsackProblemDescriptor(n)
        resultss = []

        for α in αs

            reducer(m₁, m₂) = vcat(m₁, m₂)
            reduced = @parallel (reducer) for i = 1:numberOfInstances
                generateInstanceAndCalculateRatios(α, problemDescriptor, i)
            end
            results = squeeze(mean(convert(DataArray, reduced), 1; skipmissing = true), 1)

            @info "Average recoverable ratio for α=$α is $(mean(results[1, 1])) was computed in $(@sprintf("%.2f", mean(results[1, 2])))sec on average"
            if getProblemSize(problemDescriptor) ≤ getSaneComputationLimit(problemDescriptor)
                @info "Average adversarial lower bound for α=$α is $(mean(results[2, 1])) was computed in $(@sprintf("%.2f", mean(results[2, 2])))sec on average"
                @info "Average recoverable lower bound for α=$α is $(mean(results[3, 1])) was computed in $(@sprintf("%.2f", mean(results[3, 2])))sec on average"
                @info "Average selection lower bound for α=$α is $(mean(results[4, 1])) was computed in $(@sprintf("%.2f", mean(results[4, 2])))sec on average"
            end

            push!(resultss, results)
        end

        resultss = cat(3, resultss...)
        exportKnapsackResults(problemDescriptor, αs, resultss)
    end
end

"""
    runAssignmentExperiments(ms; αs = collect(0.1:0.1:0.9), numberOfInstances = 5)

Runs experiments for minimum assignment problem.
"""
function runAssignmentExperiments(ms; αs = collect(0.1:0.1:0.9), numberOfInstances = 5)
    for m in ms
        @info "Starting minimum assignment problem experiment for m = $(m)..."

        srand(310923)

        problemDescriptor = AssignmentProblemDescriptor(m)
        resultss = []

        for α in αs

            reducer(m₁, m₂) = vcat(m₁, m₂)
            reduced = @parallel (reducer) for i = 1:numberOfInstances
                generateInstanceAndCalculateRatios(α, problemDescriptor, i)
            end
            results = squeeze(mean(convert(DataArray, reduced), 1; skipmissing = true), 1)

            @info "Average recoverable ratio for α=$α is $(mean(results[1, 1])) was computed in $(@sprintf("%.2f", mean(results[1, 2])))sec on average"
            if getProblemSize(problemDescriptor) ≤ getSaneComputationLimit(problemDescriptor)
                @info "Average adversarial lower bound for α=$α is $(mean(results[2, 1])) was computed in $(@sprintf("%.2f", mean(results[2, 2])))sec on average"
                @info "Average recoverable lower bound for α=$α is $(mean(results[3, 1])) was computed in $(@sprintf("%.2f", mean(results[3, 2])))sec on average"
                @info "Average selection lower bound for α=$α is $(mean(results[4, 1])) was computed in $(@sprintf("%.2f", mean(results[4, 2])))sec on average"
                @info "Average lagrangian lower bound for α=$α is $(mean(results[5, 1])) was computed in $(@sprintf("%.2f", mean(results[5, 2])))sec on average"

                solvedLagrangianPoblems = count(i -> !ismissing(i), reduced[:, 5, 1])
                @warn "$solvedLagrangianPoblems out of $numberOfInstances lagrangian lower bound problems was solved optimaly for α=$α"
            end

            push!(resultss, results)
        end

        resultss = cat(3, resultss...)
        exportAssignmentResults(problemDescriptor, αs, resultss)
    end
end

function generateInstanceAndCalculateRatios(α, problemDescriptor::ProblemDescriptor, i)

    @info "Generating data for instance #$(i) with α=$(α)"
    (C, c, d, Γ, X) = generateData(problemDescriptor)
    c₀ = initialScenario(c, d, Γ)

    @info "Computing recoverable ratio for instance #$(i) with α=$(α)"
    Δt₀ = @elapsed (ρ₀, x̲, x̅) = computeRecoverableRatio(C, c, d, Γ, X, α, c₀, problemDescriptor)
    @info "Computation of recoverable ratio for instance #$(i) with α=$(α) has finished in $(Δt₀)sec. with result $(ρ₀)"

    if getProblemSize(problemDescriptor) > getSaneComputationLimit(problemDescriptor)
        return cat(3, [ρ₀], [Δt₀])
    end

    Δtₙ = @elapsed numerator = computeRatioNumerator(C, c, d, Γ, X, α, x̲, x̅, problemDescriptor)
    @debug "Ratio numerator for for instance #$(i) with α=$(α) was calculated in $(Δtₙ)sec."

    @info "Computing adversarial lower bound for instance #$(i) with α=$(α)"
    Δtₐ = @elapsed ρₐ = computeAdversarialLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor)
    @info "Computation of adversarial lower bound for instance #$(i) with α=$(α) has finished in $(Δtₐ + Δtₙ)sec. with result $(ρₐ)"

    @info "Computing recoverable lower bound for instance #$(i) with α=$(α)"
    Δtₕ = @elapsed ρₕ = computeRecoverableLowerBound(C, X, α, c₀, numerator, problemDescriptor)
    @info "Computation of recoverable lower bound for instance #$(i) with α=$(α) has finished in $(Δtₕ + Δtₙ)sec. with result $(ρₕ)"

    @info "Computing selection lower bound for instance #$(i) with α=$(α)"
    Δtₛ = @elapsed ρₛ = computeSelectionLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor)
    @info "Computation of selection lower bound for instance #$(i) with α=$(α) has finished in $(Δtₛ + Δtₙ)sec. with result $(ρₛ)"

    if hasEqualCardinalityProperty(problemDescriptor)
        @assert hasEqualCardinalityProperty(problemDescriptor) "The problem is expected to has equal cardinality property"

        @info "Computing Lagrangian lower bound for instance #$(i) with α=$(α)"
        Δtₗ = @elapsed ρₗ = computeLagrangianLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor)
        @info "Computation of Lagrangian lower bound for instance #$(i) with α=$(α) has finished in $(Δtₗ + Δtₙ)sec. with result $(ρₗ)"

        cat(3, [ρ₀ ρₐ ρₕ ρₛ ρₗ numerator], [Δt₀ (Δtₐ + Δtₙ) (Δtₕ + Δtₙ) (Δtₛ + Δtₙ) (Δtₗ + Δtₙ) Δtₙ])
    else
        cat(3, [ρ₀ ρₐ ρₕ ρₛ numerator], [Δt₀ (Δtₐ + Δtₙ) (Δtₕ + Δtₙ) (Δtₛ + Δtₙ) Δtₙ])
    end
end

function computeRatioNumerator(C, c, d, Γ, X, α, x̲, x̅, problemDescriptor)
    obj₁ = evaluationProblem(C, c, d, Γ, α, x̲, X, problemDescriptor)
    obj₂ = evaluationProblem(C, c, d, Γ, α, x̅, X, problemDescriptor)
    min(obj₁, obj₂)
end

function computeRecoverableRatio(C, c, d, Γ, X, α, c₀, problemDescriptor)
    (x̲, y̲, t̲) = recoverableProblem(C, c, X, α, problemDescriptor)
    (x̅, y̅, t̅) = recoverableProblem(C, c + d, X, α, problemDescriptor)
    (x₀, y₀, t₀) = recoverableProblem(C, c₀, X, α, problemDescriptor)
    (min(t̲ + Γ, t̅) / t₀, x̲, x̅)
end

function computeAdversarialLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor)
    lb = adversarialProblem(C, c, d, Γ, X, α, problemDescriptor)
    numerator / lb
end

function computeRecoverableLowerBound(C, X, α, c₀, numerator, problemDescriptor)
    (x₀, y₀, t₀) = recoverableProblem(C, c₀, X, α, problemDescriptor)
    numerator / t₀
end

function computeSelectionLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor::ProblemDescriptor)
    t₀ = selectionLowerBound(C, c, d, Γ, X, α, problemDescriptor)
    numerator / t₀
end

function computeLagrangianLowerBound(C, c, d, Γ, X, α, numerator, problemDescriptor::ProblemDescriptor)
    l = ceil(size(c, 1) * (1 - α))
    t₀ = lagrangianLowerBound(C, c, d, Γ, X, l, problemDescriptor)
    numerator / t₀
end