text:
"""
Description createSystem
"""
function createSystem(sys::System,G=6.673e-11::Float64)
    if      typeof(sys) == EMsystem
        return SystemInfo(EM.mu,EM.Rp,EM.Rs,EM.r12,EM.m1,Float64(G))
    elseif  typeof(sys) == SEsystem
        return SystemInfo(SE.mu,SE.Rp,SE.Rs,SE.r12,SE.m1,Float64(G))
    end
end
"""
Description createConstants
"""
function createConstants(sysData::SystemInfo)
    G       = sysData.G # Gravity constant      [m3/Kg/s2]
    ### Make sure that the arguments are Floats
    mu      = sysData.mu
    if mu > 1 || mu < 0
        error("mu parameter has to be 0<mu<1")
    end
    Rp      = sysData.Rp
    Rs      = sysData.Rs
    r12     = sysData.r12
    m1      = sysData.m1
    ###
    m2      = m1*mu/(1-mu)
    # Rotation of reference system [rad/s]
    w       =  sqrt(G*(m1+m2)/r12^3)
    # Gravity constant of Earth [m^3/s^2]
    mup     = G*m1
    return ProblemInfo(mu,Rp/r12,Rs/r12,r12,w,mup/r12^3/w^2,-mu,1-mu)
end
"""
Description createPar
"""
function createPar(prob_info;   tspan = -2.0, hgoal = 300, hmmin = 100 ,
                        hmmax = 20e3, sType = "global", points = 20, x0 = [],
                        maxt0 = -1.0, rest = "basic", lim = 2000,
                        cTol = 1e-4, tRest = "ineq", alg = "lowOrder",
                        atol = 1e-7, rtol = 1e-6, retorno = true, detMin = true, y = [], optTol = 1e-8, neval = 400)

    if length(x0) == 0
        println("There has to be an initial guess ('x0' argument)")
        return nothing
    elseif length(x0) == 3
        stop = true
    elseif length(x0) == 4 || length(x0) == 6
        detMin  = 1.0
        stop    = true
    elseif length(x0) == 5 || length(x0) == 7
        stop = false
    end

    if length(x0) == 4 || length(x0) == 5
        col = true
    else
        col = false
    end

    if sType == "global"
        type = true
    elseif sType == "local"
        type = false
    else
        println("The search type has to be either 'global' or 'local' ('sType' argument)")
        return nothing
    end

    if tRest == "ineq"
        restType = true
    elseif tRest == "mixed"
        restType = false
    else
        println("The restriction type has to be 'ineq' or 'mixed' ('tRest' argument)")
        return nothing
    end
    restriction = dictRest[rest]

    if alg == "lowOrder"
        odeAlg = true
    elseif alg == "highOrder"
        odeAlg = false
    else
        println("Wrong ODE algorithm, choose between 'lowOrder' or 'highOrder'")
        return nothing
    end

    if y == []
        println("There has to be an insertion point ('y' argument)")
        return nothing
    elseif length(y) > 6 || length(y) < 6
        println("The insertion point ('y' argument) has to be of length 6")
        return nothing
    end

    println("--- Problem parameters created correctly ---")
    return ProblemParameters(Float64(tspan),Float64(hgoal*1e3/prob_info.r12),Float64(hmmin*1e3/prob_info.r12),Float64(hmmax*1e3/prob_info.r12),type,Int64(points),stop,Float64(maxt0),
                            restriction,Float64(lim/prob_info.r12/prob_info.w),Float64(cTol),restType,odeAlg,Float64(atol),Float64(rtol),
                            retorno,Float64(detMin),convert.(Float64,x0),convert.(Float64,y),optTol,Int64(neval),col)
end
"""
Description dineq!
"""
function dineq!(dy,y,p,t,mu::Float64)
    #mu  = p[1]
    xp, yp, zp, vx, vy, vz = y
    ro1 = sqrt((xp+mu)^2      +yp^2 +zp^2)
    ro2 = sqrt((xp-(1-mu))^2  +yp^2 +zp^2)
    dy[1] = vx
    dy[2] = vy
    dy[3] = vz
    dy[4] = 2*vy+xp-(1-mu)*(xp+mu)/ro1^3-mu*(xp-(1-mu))/ro2^3
    dy[5] = -2*vx+yp*(1 - (1-mu)/ro1^3-mu/ro2^3)
    dy[6] = -zp*((1-mu)/ro1^3+mu/ro2^3)
end
"""
Description calc_data
"""
function calc_data(x::Array{Float64,1},par::ProblemParameters, prob_info::ProblemInfo, PROB::Union{Tuple{ODEProblem,ODEProblem},Tuple{ODEProblem}})
    x1      = prob_info.x1      ::Float64
    x2      = prob_info.x2      ::Float64
    rearth  = prob_info.rp      ::Float64
    rmoon   = prob_info.rs      ::Float64
    prob    = PROB[1]
    y       = prob.u0           ::Array{Float64,1}

    ny0 = y - vcat(zeros(3),x[1:3])

    if par.stop
        if length(x) == 3
            intProb          = remake(prob, u0 = ny0 )
        elseif length(x) == 4
            intProb          = remake(prob, u0 = ny0 , p = [x[4], 0.0, 0.0, 0.0, par.detMin])
        else
            intProb          = remake(prob, u0 = ny0 , p = [x[4], x[5], x[6], 0.0, par.detMin])
        end
        tstop = [0.0]
    else
        if length(x) == 5
            intProb          = remake(prob, u0 = ny0 , p = [x[4], 0.0, 0.0, x[5], par.detMin])
        elseif length(x) == 7
            intProb          = remake(prob, u0 = ny0 , p = [x[4], x[5], x[6], x[7], par.detMin])
        end
        tstop = [x[end]]
    end

    if par.odeAlg
        sol     = solve(intProb, DP5(), abstol = par.atol, reltol= par.rtol, tstops= tstop, dense=false, save_everystep=false)
    else
        sol     = solve(intProb, Vern7(), abstol = par.atol, reltol= par.rtol, tstops= tstop, dense=false, save_everystep=false)
    end

    ts, tp, ps, pp, ImpF, ImpS, ImpI, tImpF, tImpS, tImpI = calcOuts(x,sol,par,prob_info)

    if par.retorno
        probR = PROB[2]
        if length(x) == 3
            Y0 = ImpF
        elseif length(sol.u) > 2
            Y0 = ImpS
        else
            Y0 = ImpF
        end
        prob2   = remake(probR, u0 = Y0 , tspan = (0,2.5) )
        if par.odeAlg
            solr    = solve(prob2, DP5(), abstol = par.atol, reltol= par.rtol, dense=false, save_everystep=false)
        else
            solr    = solve(prob2, Vern7(), abstol = par.atol, reltol= par.rtol, dense=false, save_everystep=false)
        end
        eRet       = calcDists(solr.u[end], x1,rearth)
    else
        eRet       = 0.0
    end

    ds, dp   = calcDists(ps, x2, rmoon), calcDists(pp, x1, rearth) # distances to primaries in the minimum distance points

    fimp = fImp(pp,prob_info)

    return ds, dp, ts, tp, eRet, ImpF, ImpS, ImpI, tImpF, tImpS, tImpI, fimp
end
# Add method to calc_data
function calc_data(x::Vector{T},par::ProblemParameters, prob_info::ProblemInfo, PROB::Tuple{ODEProblem}) where {T<:ForwardDiff.Dual}
    x1      = prob_info.x1      ::Float64
    x2      = prob_info.x2      ::Float64
    rearth  = prob_info.rp      ::Float64
    rmoon   = prob_info.rs      ::Float64
    prob    = PROB[1]

    xType = eltype(x)
    if par.stop
        if length(x) == 3
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x), tspan = convert.(xType, prob.tspan) )
        elseif length(x) == 4
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], 0.0, 0.0, 0.0, par.detMin]) )
        else
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], x[5], x[6], 0.0, par.detMin]) )
        end
        tstop = convert.(xType,[0.0])
    else
        if length(x) == 5
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], 0.0, 0.0, x[5], par.detMin]) )
        elseif length(x) == 7
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], x[5], x[6], x[7], par.detMin]) )
        end
        tstop = x[end]
    end

    if par.odeAlg
        sol     = solve(intProb, DP5(), abstol = par.atol, reltol= par.rtol, tstops=tstop, dense=false, save_everystep=false)
    else
        sol     = solve(intProb, Vern7(), abstol = par.atol, reltol= par.rtol, tstops=tstop, dense=false, save_everystep=false)
    end

    ts, tp, ps, pp = calcOuts(x,sol,par,prob_info,PROB)

    #ds, dp   = calcDists(ps, x2, rmoon), calcDists(pp, x1, rearth) # distances to primaries in the minimum distance points
    ds = calcDists(ps, x2, rmoon)
    dp = calcDists(pp, x1, rearth)

    return ds, dp, ts, tp
end
function calc_data(x::Vector{T},par::ProblemParameters, prob_info::ProblemInfo, PROB::Tuple{ODEProblem,ODEProblem}) where {T<:ForwardDiff.Dual}
    x1      = prob_info.x1      ::Float64
    x2      = prob_info.x2      ::Float64
    rearth  = prob_info.rp      ::Float64
    rmoon   = prob_info.rs      ::Float64
    prob    = PROB[1]
    probR   = PROB[2]

    xType = eltype(x)
    if par.stop
        if length(x) == 3
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x), tspan = convert.(xType, prob.tspan) )
        elseif length(x) == 4
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], 0.0, 0.0, 0.0, par.detMin]) )
        else
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], x[5], x[6], 0.0, par.detMin]) )
        end
        tstop = convert.(xType,[0.0])
    else
        if length(x) == 5
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], 0.0, 0.0, x[5], par.detMin]) )
        elseif length(x) == 7
            intProb          = remake(prob, u0 = prob.u0 - vcat(zeros(3),x[1:3]), tspan = convert.(xType, prob.tspan), p = convert.(xType,[x[4], x[5], x[6], x[7], par.detMin]) )
        end
        tstop = x[end]
    end

    if par.odeAlg
        sol     = solve(intProb, DP5(), abstol = par.atol, reltol= par.rtol, tstops=tstop, dense=false, save_everystep=false)
    else
        sol     = solve(intProb, Vern7(), abstol = par.atol, reltol= par.rtol, tstops=tstop, dense=false, save_everystep=false)
    end

    ts, tp, ps, pp, IMP = calcOuts(x,sol,par,prob_info,PROB)

    prob2   = remake(probR, u0 = IMP , tspan = convert.(xType,(0.0,2.5)) )
    if par.odeAlg
        solr    = solve(prob2, DP5(), abstol = par.atol, reltol= par.rtol, dense=false, save_everystep=false, callback = callRet )
    else
        solr    = solve(prob2, Vern7(), abstol = par.atol, reltol= par.rtol, dense=false, save_everystep=false, callback = callRet )
    end
    eRet       = calcDists(solr.u[end], x1,rearth)

    #ds, dp   = calcDists(ps, x2, rmoon), calcDists(pp, x1, rearth) # distances to primaries in the minimum distance points
    ds = calcDists(ps, x2, rmoon)
    dp = calcDists(pp, x1, rearth)

    return ds, dp, ts, tp, eRet
end
"""
Description calcOuts
"""
function calcOuts(x::Array{Float64,1},sol::ODESolution,par::ProblemParameters,prob_info::ProblemInfo)
    ts = sol.t[1] # Instant of close approach to secondary body
    tp = sol.t[end] # Instant of close approach to primary body
    ps = sol.u[1] # State of close approach to secondary body
    pp = sol.u[end] # State of close approach to primary body

    ImpF    = sol.u[1]
    tImpF   = sol.t[1]
    ImpI    = sol.u[end]
    tImpI   = sol.t[end]

    ImpS    = sol.u[2]
    tImpS   = sol.t[2]
    if length(x) == 3
        if length(sol.t) == 3
            ts = sol.t[2]
            ps = sol.u[2]
        end
    else
        if par.stop
            if length(sol.t) == 3
                ts = sol.t[2]
                ps = sol.u[2]
            end
        else
            if length(sol.t) >= 4
                tImp        = x[end]

                indexImp = findfirst( x -> sol.t[x] == tImp, sol.t)
                ImpS  = sol.u[indexImp]
                tImpS   = sol.t[indexImp]

                Mclose = [ (calcDists(sol.u[i],prob_info.x2,prob_info.rs),i) for i in 2:(length(sol.t)-1) if i != indexImp ]
                sort!(Mclose,by = x -> x[1])
                ts = sol.t[Mclose[1][2]]
                ps = sol.u[Mclose[1][2]]
            end
        end
    end
    return ts, tp, ps, pp, ImpF, ImpS, ImpI, tImpF, tImpS, tImpI
end
function calcOuts(x::Vector{T},sol::ODESolution,par::ProblemParameters,prob_info::ProblemInfo,PROB::Tuple{ODEProblem}) where {T<:ForwardDiff.Dual}
    ts = sol.t[1] # Instant of close approach to secondary body
    tp = sol.t[end] # Instant of close approach to primary body
    ps = sol.u[1] # State of close approach to secondary body
    pp = sol.u[end] # State of close approach to primary body

    if length(x) == 3
        if length(sol.t) == 3
            ts = sol.t[2]
            ps = sol.u[2]
        end
    else
        if par.stop
            if length(sol.t) == 3
                ts = sol.t[2]
                ps = sol.u[2]
            end
        else
            if length(sol.t) >= 4
                tImp        = x[end]

                indexImp = findfirst( x -> sol.t[x] == tImp, sol.t)

                Mclose = [ (calcDists(sol.u[i],prob_info.x2,prob_info.rs),i) for i in 2:(length(sol.t)-1) if i != indexImp ]
                sort!(Mclose,by = x -> x[1])
                ts = sol.t[Mclose[1][2]]
                ps = sol.u[Mclose[1][2]]
            end
        end
    end
    return ts, tp, ps, pp
end
function calcOuts(x::Vector{T},sol::ODESolution,par::ProblemParameters,prob_info::ProblemInfo,PROB::Tuple{ODEProblem,ODEProblem}) where {T<:ForwardDiff.Dual}
    ts = sol.t[1] # Instant of close approach to secondary body
    tp = sol.t[end] # Instant of close approach to primary body
    ps = sol.u[1] # State of close approach to secondary body
    pp = sol.u[end] # State of close approach to primary body

    if length(sol.u) > 2
        IMP     = sol.u[2]
    else
        IMP     = sol.u[1]
    end

    if length(x) == 3
        IMP     = sol.u[1]
        if length(sol.t) == 3
            ts = sol.t[2]
            ps = sol.u[2]
        end
    else
        if par.stop
            if length(sol.t) == 3
                ts = sol.t[2]
                ps = sol.u[2]
            end
        else
            if length(sol.t) >= 4
                tImp        = x[end]

                indexImp = findfirst( x -> sol.t[x] == tImp, sol.t)
                IMP  = sol.u[indexImp]

                Mclose = [ (calcDists(sol.u[i],prob_info.x2,prob_info.rs),i) for i in 2:(length(sol.t)-1) if i != indexImp ]
                sort!(Mclose,by = x -> x[1])
                ts = sol.t[Mclose[1][2]]
                ps = sol.u[Mclose[1][2]]
            end
        end
    end
    return ts, tp, ps, pp, IMP
end
"""
Description myObj
"""
function myObj(X::Vector,grad::Vector,par::ProblemParameters)
    if par.stop == 0
        x = X[1:end-1]
    else
        x = X
    end
    #println(x)
    if length(grad) > 0
        if par.stop == 0
            grad[:] = vcat(2x,[0])
        else
            grad[:] = 2x
        end
    end
    return dot(x,x)
end
"""
Description fConstraint
"""
function fConstraint(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem}) where {T<:ForwardDiff.Dual}

    ds, dp, ts, tp  = calc_data(x, par, prob_info, PROB)

    n = par.rest[1]
    if      n == 1
        res = [ dp - par.hgoal ]
    elseif n == 2
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9 ]
    elseif n == 3
        res = [ dp - par.hgoal,
                -tp+par.tspan*0.95 ]
    elseif  n == 4
        res = [ dp - par.hgoal,
                -ds + par.hmmin ]
    elseif n == 8
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9,
                -ds + par.hmmin ]
    end
    return  res
end
function fConstraint(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem,ODEProblem}) where {T<:ForwardDiff.Dual}

    ds, dp, ts, tp, eRet  = calc_data(x, par, prob_info, PROB)

    n = par.rest[1]
    if      n == 1
        res = [ dp - par.hgoal ]
    elseif n == 2
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9 ]
    elseif n == 3
        res = [ dp - par.hgoal,
                -tp+par.tspan*0.95 ]
    elseif  n == 4
        res = [ dp - par.hgoal,
                -ds + par.hmmin ]
    elseif n == 5
        res = [ dp - par.hgoal,
                eRet-1.1*par.hgoal ]
    elseif n == 6
        res = [ dp - par.hgoal,
                -ds+par.hmmin,
                eRet-1.1*par.hgoal ]
    elseif n == 7
        res = [ dp - par.hgoal,
                -eRet,
                eRet-1.1*par.hgoal ]
    elseif n == 8
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9,
                -ds + par.hmmin ]
    end
    return  res
end
function fConstraint(x::Array{Float64,1}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem})

    ds, dp, ts, tp, eRet, ImpF, ImpS, ImpI, tImpF, tImpS, tImpI, fimp  = calc_data(x, par, prob_info, PROB)

    n = par.rest[1]
    if      n == 1
        res = [ dp - par.hgoal ]
    elseif n == 2
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9 ]
    elseif n == 3
        res = [ dp - par.hgoal,
                -tp+par.tspan*0.95 ]
    elseif  n == 4
        res = [ dp - par.hgoal,
                -ds + par.hmmin ]
    elseif n == 5
        res = [ dp - par.hgoal,
                eRet-1.1*par.hgoal ]
    elseif n == 6
        res = [ dp - par.hgoal,
                -ds+par.hmmin,
                eRet-1.1*par.hgoal ]
    elseif n == 7
        res = [ dp - par.hgoal,
                -eRet,
                eRet-1.1*par.hgoal ]
    elseif n == 8
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9,
                -ds + par.hmmin ]
    end
    return  res
end
"""
Description fConstraint
"""
function fAltitude(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem}) where {T<:ForwardDiff.Dual}
    ds, dp, ts, tp  = calc_data(x, par, prob_info, PROB)
    return  dp - par.hgoal
end
function fAltitude(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem,ODEProblem}) where {T<:ForwardDiff.Dual}

    ds, dp, ts, tp, eRet  = calc_data(x, par, prob_info, PROB)

    return  dp - par.hgoal
end
"""
Description fRest
"""
function fRest(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem}) where {T<:ForwardDiff.Dual}

    ds, dp, ts, tp  = calc_data(x, par, prob_info, PROB)

    n = par.rest[1]
    if n == 2
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9 ]
    elseif n == 3
        res = [ dp - par.hgoal,
                -tp+par.tspan*0.95 ]
    elseif  n == 4
        res = [ dp - par.hgoal,
                -ds + par.hmmin ]
    elseif n == 8
        res = [ dp - par.hgoal,
                -dp +  par.hgoal*0.9,
                -ds + par.hmmin ]
    end
    return  res
end
function fRest(x::Vector{T}, par::ProblemParameters,prob_info::ProblemInfo, PROB::Tuple{ODEProblem,ODEProblem}) where {T<:ForwardDiff.Dual}

    ds, dp, ts, tp, eRet  = calc_data(x, par, prob_info, PROB)

    if n == 2
        res = [ -dp +  par.hgoal*0.9 ]
    elseif n == 3
        res = [ -tp+par.tspan*0.95 ]
    elseif  n == 4
        res = [ -ds + par.hmmin ]
    elseif n == 5
        res = [ eRet-1.1*par.hgoal ]
    elseif n == 6
        res = [ -ds+par.hmmin,
                eRet-1.1*par.hgoal ]
    elseif n == 7
        res = [ -eRet,
                eRet-1.1*par.hgoal ]
    elseif n == 8
        res = [ -dp +  par.hgoal*0.9,
                -ds + par.hmmin ]
    end
    return  res
end
"""
Description myConstraint, myConstAlt, myConstRest
"""
function myConstraint(result::Vector,x::Vector,grad::Matrix, par::ProblemParameters,prob_info::ProblemInfo, PROB, results)
    ForwardDiff.jacobian!(results, (x0::Vector) -> fConstraint(x0, par,prob_info, PROB) , x)
    if length(grad) > 0
        grad[:] = DiffResults.jacobian(results)'
    end
    result[:] = DiffResults.value(results)
end
function myConstAlt(result::Vector,x::Vector,grad::Matrix, par::ProblemParameters,prob_info::ProblemInfo, PROB, results)
    ForwardDiff.gradient!(results, (x0::Vector) -> fAltitude(x0, par,prob_info, PROB) , x)
    if length(grad) > 0
        grad[:] = reshape( DiffResults.gradient(results) ,length(x),1)
    end
    result[:] = [DiffResults.value(results)]
end
function myConstRest(result::Vector,x::Vector,grad::Matrix, par::ProblemParameters,prob_info::ProblemInfo, PROB, results)
    ForwardDiff.jacobian!(results, (x0::Vector) -> fRest(x0, par,prob_info, PROB) , x)
    if length(grad) > 0
        grad[:] = DiffResults.jacobian(results)'
    end
    result[:] = DiffResults.value(results)
end
"""
Description fImp
"""
function fImp(p::Vector{<:Real},prob_info::ProblemInfo) # Computes the value of the collinear impulse to perform in p to circularize (adimensional velocity)
    xp, yp, zp, vx, vy, vz = p
    mu  = prob_info.mu  ::Float64
    mup = prob_info.mup ::Float64
    return sqrt((vx-yp)^2+(vy+xp+mu)^2+vz^2) - sqrt( mup / sqrt( (xp+mu)^2 + yp^2 + zp^2 ) )
end
"""
Description fImp
"""
calcDists(y::Vector{<:Real},x::Float64,r::Float64)       = sqrt((y[1]-x)^2+y[2]^2+y[3]^2)-r
"""
Description plotCC
"""
function plotCC(body::Int, n = 16 ::Int, d=prob_info::ProblemInfo) # Function to plot Earth and Moon
    if body == 1 # plot prymary body
        pos             = -d.mu
        rad             = d.rp
        colorSurface    = :blues
    elseif body == 0 # plot secondary
        pos             = 1-d.mu
        rad             = d.rs
        colorSurface    = :Greys
    end
    u = collect(range(0,stop=2,length=n));
    v = collect(range(0,stop=1,length=n));
    x = rad * cospi.(u) * sinpi.(v)' .+ pos;
    y = rad * sinpi.(u) * sinpi.(v)';
    z = rad * repeat(cospi.(v)',outer=[n, 1])
    surface!(x,y,z, color=colorSurface)
end
"""
Description fsObj
"""
function fsObj(x::Vector,grad::Vector, intF)
    if length(grad) > 0
        grad[:] = ForwardDiff.jacobian(intF, x)'
    end
    return intF(x)[1]
end
"""
Description point
"""
function point(i::Int,orbit::Array{Float64,2}) # Uses interpolation to extract the relevant points
    ORB     = interpolate(orbit, BSpline(Quadratic(Reflect(OnCell()))))
    Y(x)    = ORB(x,2)

    if i == 1 # Point of maximum Z coordinate
        imax = argmax(orbit[:,3])
        if (imax==0) || (imax==length(orbit[:,1]))
            ind = imax
        else
            ind = find_zero(Y, imax, Order2())
        end
    elseif i == 4 # Point of minimum Z coordinate
        imax = argmax(-orbit[:,3])
        if (imax==0) || (imax==length(orbit[:,1]))
            ind = imax
        else
            ind = find_zero(Y, imax, Order2())
        end
    elseif i == 6 # Point of most negative Y coordinate
        imax = argmax(-orbit[:,2])
        if (imax==0) || (imax==length(orbit[:,1]))
            ind = imax
        else
            OPT = Opt(:LD_SLSQP, 1)
            #LN_COBYLA
            #LD_SLSQP
            lower_bounds!(OPT, [1.1])
            upper_bounds!(OPT, [200.9])
            ftol_rel!(OPT::Opt, 1e-14::Float64)
            ftol_abs!(OPT::Opt, 1e-14::Float64)
            xtol_rel!(OPT::Opt, 1e-14::Float64)
            xtol_abs!(OPT::Opt, 1e-14::Float64)

            min_objective!(OPT::Opt, (x,grad) -> fsObj(x,grad,Y))
            (optf,optx,ret) = optimize(OPT, [imax])
            ind = optx[1]
        end
    elseif i == 2 # Point of most positive Y coordinate
        imax = argmax(orbit[:,2])
        if (imax==0) || (imax==length(orbit[:,1]))
            ind = imax
        else
            OPT = Opt(:LD_SLSQP, 1)
            lower_bounds!(OPT, [1.1])
            upper_bounds!(OPT, [200.9])
            ftol_rel!(OPT::Opt, 1e-14::Float64)
            ftol_abs!(OPT::Opt, 1e-14::Float64)
            xtol_rel!(OPT::Opt, 1e-14::Float64)
            xtol_abs!(OPT::Opt, 1e-14::Float64)

            max_objective!(OPT::Opt, (x::Vector,grad::Vector) -> fsObj(x,grad,Y))
            (optf,optx,ret) = optimize(OPT, [imax])
            ind = optx[1]

        end
    elseif (i == 3) || (i == 5) # Points of zero Z coordinate
        Z(x)    = ORB(x,3)
        indexes = [i for i in 1:length(orbit[:,1])-1 if orbit[i,3].*orbit[i+1,3] < 0 ]
        if      i == 5
            iz  = [indexes[i] for i in 1:2 if Y(indexes[i]) < 0 ]
            ind = find_zero(Z, iz[1], Order2())
        elseif  i == 3
            iz  = [indexes[i] for i in 1:2 if Y(indexes[i]) > 0 ]
            ind = find_zero(Z, iz[1], Order2())
        end
    end
    return ORB(ind,1:6), ind
end
"""
Description Pset
"""
function Pset(lbounds::Array{Float64,1},ubounds::Array{Float64,1}, x0::Array{Float64,1}, n = 50 ::Int64)::Array{Array{Float64,1},1} # Creates the starting points of the GlobalSearch like MultiStart
    cuts = Int(floor(n^(1/length(x0)))) # number of cuts in the grid
    C = ones(Int, length(x0))*(cuts + 1)
    p = prod(C); i = 1
    while p > n
        C[i]    -= 1
        p       = prod(C)
        i       += 1
        if i > length(x0)
            i = 1
        end
    end

    for i in 1:length(C)
        if C[i] == 1
            C[i] += 1
        end
    end

    rangeVar = [collect(range(lbounds[i],stop = ubounds[i], length = C[i])) for i in 1:length(x0)]
    if length(x0) == 3
        setOfPoints = [ [i, j, k] for i=rangeVar[1] for j=rangeVar[2] for k=rangeVar[3]  ]
    elseif length(x0) == 4
        setOfPoints = [ [i,j,k,l] for i=rangeVar[1] for j=rangeVar[2] for k=rangeVar[3] for l=rangeVar[4]  ]
    elseif length(x0) == 5
        setOfPoints = [ [i,j,k,l,m] for i=rangeVar[1] for j=rangeVar[2] for k=rangeVar[3] for l=rangeVar[4] for m=rangeVar[5]  ]
    elseif length(x0) == 6
        setOfPoints = [ [i,j,k,l,m,n] for i=rangeVar[1] for j=rangeVar[2] for k=rangeVar[3] for l=rangeVar[4] for m=rangeVar[5] for n=rangeVar[6]  ]
    elseif length(x0) == 7
        setOfPoints = [ [i,j,k,l,m,n,p] for i=rangeVar[1] for j=rangeVar[2] for k=rangeVar[3] for l=rangeVar[4] for m=rangeVar[5] for n=rangeVar[6] for p=rangeVar[7]  ]
    else
        println("Wrong input, x0 dimension greater than 7")
        return [[0.0]]
    end

    while length(setOfPoints) < n-1
        point = ubounds.*rand(length(x0))+lbounds.*rand(length(x0))
        setOfPoints = cat(setOfPoints, [point], dims = 1)
    end
    return cat(setOfPoints, [x0], dims = 1)
end
"""
Description setOptProblem
"""
function setOptProblem(par::ProblemParameters,prob_info::ProblemInfo,PROB)

    if par.stop == 1
        lbounds = -ones(length(par.x0))*par.lim
        ubounds =  ones(length(par.x0))*par.lim
    else
        if par.maxt0 < 0
            lbounds = vcat(-ones(length(par.x0)-1)*par.lim, [par.maxt0] )
            ubounds = vcat( ones(length(par.x0)-1)*par.lim, [-1e-3] )
        else
            lbounds = vcat(-ones(length(par.x0)-1)*par.lim, [-1e-3] )
            ubounds = vcat( ones(length(par.x0)-1)*par.lim, [par.maxt0] )
        end
    end

    # Optimization
    Nopt = length(par.x0) # Number of optimization Variables
    alg = :LD_SLSQP
    OPT = Opt(alg, Nopt)
    lower_bounds!(OPT::Opt, lbounds)
    upper_bounds!(OPT::Opt, ubounds)
    min_objective!(OPT::Opt, (X::Vector,grad::Vector) -> myObj(X,grad,par))
    ftol_rel!(OPT::Opt, par.optTol)
    ftol_abs!(OPT::Opt, par.optTol)
    xtol_rel!(OPT::Opt, par.optTol)
    xtol_abs!(OPT::Opt, par.optTol)
    maxeval!(OPT::Opt, par.neval)

    if par.tRest
        # Create DiffResults object
        results = DiffResults.JacobianResult(ones(par.rest[2]),par.x0)
        inequality_constraint!(OPT, (result::Vector,x::Vector,grad::Matrix) -> myConstraint(result,x,grad, par,prob_info, PROB, results), par.cTol*ones(par.rest[2]) )
    else
        if par.rest[2] > 1
            results = DiffResults.JacobianResult(ones(par.rest[2]-1),par.x0)
            inequality_constraint!(OPT, (result::Vector,x::Vector,grad::Matrix) -> myConstRest(result,x,grad, par,prob_info, PROB, results), par.cTol*ones(par.rest[2]-1) )
        end
        Eresults = DiffResults.GradientResult(par.x0)
        equality_constraint!(OPT, (result::Vector,x::Vector,grad::Matrix) -> myConstAlt(result,x,grad, par,prob_info, PROB, Eresults), par.cTol*ones(1) )
    end

    if par.sType
        return OPT, Pset(lbounds,ubounds,par.x0,par.points)
    else
        return OPT, [par.x0]
    end
end
"""
Description createPROB
"""
function createPROB(par::ProblemParameters, prob_info::ProblemInfo)
    call, cRet          = createCallback(par,prob_info)
    # create ODE problem
    if par.retorno
        prob    = ODEProblem((dy,y,p,t) -> dineq!(dy,y,p,t,prob_info.mu), par.y, par.tspan, [0.0, 0.0, 0.0, 0.0, 1.0], callback = call)
        probRet = ODEProblem((dy,y,p,t) -> dineq!(dy,y,p,t,prob_info.mu), par.y, par.tspan, [0.0, 0.0, 0.0, 0.0, 1.0], callback = cRet)
        PROB    = tuple(prob,probRet)
    else
        prob    = ODEProblem((dy,y,p,t) -> dineq!(dy,y,p,t,prob_info.mu), par.y, par.tspan, [0.0, 0.0, 0.0, 0.0, 1.0], callback = call)
        PROB    = tuple(prob)
    end
    return PROB
end
"""
Description createPROB
"""
function search(par::ProblemParameters, prob_info::ProblemInfo)
    # create ODE problem(s)
    PROB = createPROB(par, prob_info)
    #set optimization problem
    OPT, searchPoints   = setOptProblem(par,prob_info,PROB)
    @show OPT
    # Perform search
    sol = map((x) -> myOptimizer(x,OPT),searchPoints)
    return sol
end
"""
Description createCallback
"""
function createCallback(par::ProblemParameters,prob_info::ProblemInfo)
    lX0     = length(par.x0)
    t2      = par.tspan
    detMin  = par.detMin
    stop    = par.stop
    col     = par.col
    mu      = prob_info.mu

    cRet   = ContinuousCallback((y,t,integrator) -> eConditionRet(y,t,integrator,mu),terminate!, nothing, save_positions = (false,true), rootfind=true) # Fee Return stoping

    if lX0 == 3

        if t2 < 0
            cDist   = ContinuousCallback((y,t,integrator) -> eCondition(y,t,integrator,mu),nothing, terminate!, save_positions = (false,true), rootfind=true)
            mDist   = ContinuousCallback((y,t,integrator) -> mCondition(y,t,integrator,mu),nothing, saveEvent!, save_positions = (false,true), rootfind=true)
        else
            cDist   = ContinuousCallback((y,t,integrator) -> eCondition(y,t,integrator,mu),terminate!,nothing, save_positions = (false,true), rootfind=true)
            mDist   = ContinuousCallback((y,t,integrator) -> mCondition(y,t,integrator,mu),saveEvent!, nothing, save_positions = (false,true), rootfind=true)
        end

        if detMin == 1
            call = CallbackSet(cDist,mDist)
        else
            call = CallbackSet(cDist)
        end

    else

        if t2 < 0
            if stop == 1
                cDist   = ContinuousCallback((y,t,integrator) -> gCondition(y,t,integrator,mu),nothing, (integrator) -> impulse!(integrator,col),   save_positions = (false,true), rootfind=true)
            else
                mDist   = ContinuousCallback((y,t,integrator) -> mCondition(y,t,integrator,mu),nothing, saveEvent!, save_positions = (false,true), rootfind=true)
                eDist   = ContinuousCallback((y,t,integrator) -> eCondition(y,t,integrator,mu),nothing, terminate!, save_positions = (false,true), rootfind=true)
            end
        else
            if stop == 1
                cDist   = ContinuousCallback((y,t,integrator) -> gCondition(y,t,integrator,mu), (integrator) -> impulse!(integrator,col),   nothing, save_positions = (false,true), rootfind=true)
            else
                mDist   = ContinuousCallback((y,t,integrator) -> mCondition(y,t,integrator,mu),saveEvent!, nothing, save_positions = (false,true), rootfind=true)
                eDist   = ContinuousCallback((y,t,integrator) -> eCondition(y,t,integrator,mu),terminate!, nothing, save_positions = (false,true), rootfind=true)
            end
        end

        if stop == 1
            call = CallbackSet(cDist)
        else
            iDist   = DiscreteCallback(iCondition, (integrator) -> impulseInt!(integrator,col),save_positions = (false,true) )
            if detMin == 0
                call = CallbackSet(eDist,iDist)
            else
                call = CallbackSet(eDist,iDist,mDist)
            end
        end
    end
    return call, cRet
end
"""
Description conditions and affect!'s
"""
function gCondition(y,t,integrator,mu::Float64)
    xp, yp, zp, vx, vy, vz = y
    #mu      = integrator.p[1]
    # the last parameter indicates where the minimum is to be looked for
    if integrator.p[end] == 1.0  # the minimum is with the secondary body when p[end] is "true"
        return (xp - 1.0 + mu)*vx +  yp*vy + zp*vz
    else
        return (xp + mu)*vx +  yp*vy + zp*vz
    end
end
function iCondition(y,t,integrator)
    integrator.p[4] == t
end
function eCondition(y::Array{Float64,1},t::Float64,integrator,mu::Float64)  #Float version
    xp, yp, zp, vx, vy, vz = y
    #mu      = integrator.p[1]
    tint    = integrator.p[4]
    if xp > 0.6 || t > tint
        return 0.0
    else
        return (xp + mu)*vx +  yp*vy + zp*vz
    end
end
function eCondition(y,t,integrator,mu::Float64)
    xp, yp, zp, vx, vy, vz = y
    #mu      = integrator.p[1]
    tint    = integrator.p[4]
    if xp > 0.6 || t > tint
        return convert(eltype(y),0.0)
    else
        return (xp + mu)*vx +  yp*vy + zp*vz
    end
end
function eConditionRet(y::Array{Float64,1},t::Float64,integrator,mu::Float64) # Float version
    xp, yp, zp, vx, vy, vz = y
    #mu = integrator.p[1]
    if xp > 0.6
        return 0.0
    else
        return (xp + mu)*vx +  yp*vy + zp*vz
    end
end
function eConditionRet(y,t,integrator,mu::Float64)
    xp, yp, zp, vx, vy, vz = y
    #mu = integrator.p[1]
    if xp > 0.6
        return convert(eltype(y),0.0)
    else
        return (xp + mu)*vx +  yp*vy + zp*vz
    end
end
function mCondition(y,t,integrator,mu::Float64)
    xp, yp, zp, vx, vy, vz = y
    #mu = integrator.p[1]
    cond = (xp - 1.0 + mu)*vx +  yp*vy + zp*vz
    return cond
end
function impulse!(integrator,col::Bool)
    imp = integrator.p[1:3]
    if integrator.p[end] == 1.0
        if col
            integrator.u[4:6] = integrator.u[4:6]*(1.0-imp[1]/norm(integrator.u[4:6]))
        else
            integrator.u[4:6] = integrator.u[4:6] - imp
        end
        integrator.p[end] = convert(eltype(integrator.u),0.0) # change the minimum search to the other body
    else
        terminate!(integrator)
    end
end
function saveEvent!(integrator)
    nothing
    #println("Moon encounter at t: ",integrator.t)
end
function impulseInt!(integrator,col::Bool)
    imp = integrator.p[1:3]
    if col
        integrator.u[4:6] = integrator.u[4:6]*(1.0-imp[1]/norm(integrator.u[4:6]))
    else
        integrator.u[4:6] = integrator.u[4:6] - imp
    end
end
"""
Description myOptimizer
"""
function myOptimizer(x::Array{Float64,1}, OPT)
    f, xf, res = optimize(OPT,x)
    println("The optimizer reached the value: ",f)
    println(res)
    return f, xf, res, x
end
"""
Description myOptimizer
"""
function filterGlobalSol(Gsol, par::ProblemParameters,prob_info::ProblemInfo)
    # create callbacks
    PROB = createPROB(par, prob_info)
    #
    #Eresults = DiffResults.GradientResult(par.x0)
    #for i in 1:length(Gsol)
    #    ForwardDiff.gradient!(Eresults, (x0::Vector) -> fAltitude(x0, par,prob_info, PROB) , Gsol[i][2])
    #    @show abs.( [DiffResults.value(Eresults)] ), Gsol[i][1]

        #@show abs.( fConstraint( Gsol[i][2], par,prob_info, PROB)  ), Gsol[i][1]
    #end
    filterSol = [ Gsol[i] for i in 1:length(Gsol) if abs.( fConstraint( Gsol[i][2], par,prob_info, PROB) ) < par.cTol*ones(par.rest[2]) ]

    sfSol = [filterSol[1]] # filtrar por diferencias
    for i in 2:length(filterSol)
        #println(filterSol[i][1])
        dif = 1
        for j in 1:length(sfSol)
            difF = abs(sfSol[j][1] - filterSol[i][1])/maximum([1,sfSol[j][1]])
            difX = norm(sfSol[j][2] - filterSol[i][2])/maximum([1,norm(sfSol[j][2])])
            if (difF < 0.1) && (difX < 0.1)
                dif = 0
            end
        end
        if dif == 1
            #println("new sol")
            sfSol = vcat(sfSol, filterSol[i])
        end
    end
    println("Number of solutions: ",length(sfSol))
    sfSol = sort!(sfSol, by = x -> x[1])
    return filterSol, sfSol
end
"""
Description plotSols
"""
function plotSols(sfSol, par::ProblemParameters, prob_info::ProblemInfo, orbit::Matrix{<:Real},guiP::Bool)
    # create ODE problem
    PROB = createPROB(par, prob_info)

    pyplot()
    ploty = plot(orbit[:,1],orbit[:,2],orbit[:,3], yflip=true, leg=true, lab="Halo Orbit", aspect_ratio=:equal)
    scatter!([-prob_info.mu],[0],[0],c=:blue,label="Earth",markersize=2)
    scatter!([1-prob_info.mu],[0],[0],c=:grey,label="Moon",markersize=1)

    y0 = par.y
    prob = PROB[1]
    for i in 1:length(sfSol)

        x = sfSol[i][2]

        ds, dp, ts, tp, eRet, ImpF, ImpS, ImpI, tImpF, tImpS, tImpI, fimp = calc_data(x, par, prob_info, PROB)

        if par.stop
            if length(x) == 3
                intProb           = remake(prob, u0 = y0 - vcat(zeros(3),x[1:3]) )
            elseif length(x) == 4
                intProb           = remake(prob, u0 = y0 - vcat(zeros(3),x[1:3]) , p = [x[4], 0.0, 0.0, 0.0, par.detMin] )
            else
                intProb           = remake(prob, u0 = y0 - vcat(zeros(3),x[1:3]) , p = [x[4], x[5], x[6], 0.0, par.detMin] )
            end
            tstop = [0.0]
        else
            if length(x) == 5
                intProb           = remake(prob, u0 = y0 - vcat(zeros(3),x[1:3]) , p = [x[4], 0.0, 0.0, x[5], par.detMin] )
            elseif length(x) == 7
                intProb           = remake(prob, u0 = y0 - vcat(zeros(3),x[1:3]) , p = [x[4], x[5], x[6], x[7], par.detMin] )
            end
            tstop = [x[end]]
        end
        if par.odeAlg
            sol         = solve(intProb, DP5(), abstol= par.atol, reltol= par.rtol, tstops =  tstop , dense=true, save_everystep=true)
        else
            sol         = solve(intProb, Vern7(), abstol= par.atol, reltol= par.rtol, tstops = tstop , dense=true, save_everystep=true)
        end
        if length(x) == 3
            deltav = @sprintf("%.2f (m/s)",(fimp + norm(x))*prob_info.r12*prob_info.w)
        elseif par.stop
            deltav = @sprintf("%.2f (m/s)",(fimp + norm(x[1:3]) + norm(x[4:end]))*prob_info.r12*prob_info.w)
        else
            deltav = @sprintf("%.2f (m/s)",(fimp + norm(x[1:3]) + norm(x[4:end-1]))*prob_info.r12*prob_info.w)
        end
        plot!(sol, vars=(1,2,3), leg=true, label=string(deltav), xlim=:auto)
    end
    #xlims!(-0.05,1.1)
    if guiP
        gui()
        return ploty
    end
end


text:
using AutoHashEquals

"""
    Compound <: RDate

Apply a list of rdates in order
"""
@auto_hash_equals struct Compound <: RDate
    parts::Vector{RDate}
end

apply(rdate::Compound, date::Dates.Date, cal_mgr::CalendarManager) = Base.foldl((x,y) -> apply(y, x, cal_mgr), rdate.parts, init=date)
multiply(rdate::Compound, count::Integer) = Compound(map(x -> multiply(x, count), rdate.parts))
combine(left::RDate, right::RDate) = Compound([left,right])
Base.:+(left::RDate, right::RDate) = combine(left, right)
Base.:-(left::RDate, right::RDate) = combine(left, -right)
Base.:-(rdate::Compound) = Compound(map(-, rdate.parts))
Base.:+(x::Compound, y::RDate) = Compound(vcat(x.parts, y))
Base.:+(x::RDate, y::Compound) = Compound(vcat(x, y.parts))

function Base.show(io::IO, rdate::Compound)
    for (i,part) in enumerate(rdate.parts)
        if i > 1 print(io, "+") end
        show(io, part)
    end
end

"""
    Repeat <: RDate

Repeat the application of an rdate n times.
"""
@auto_hash_equals struct Repeat <: RDate
    count::Int64
    part::RDate

    function Repeat(count::Int64, part::RDate)
        count > 0 || error("Repeat must use a positive count, not $count")
        new(count, part)
    end

    Repeat(part::RDate) = new(1, part)
end

function apply(rdate::Repeat, date::Dates.Date, cal_mgr::CalendarManager)
    for _ in 1:rdate.count
        date = apply(rdate.part, date, cal_mgr)
    end
    date
end

multiply(rdate::Repeat, count::Integer) = Repeat(count*rdate.count, rdate.part)
Base.:-(rdate::Repeat) = Repeat(rdate.count, -rdate.part)
Base.show(io::IO, rdate::Repeat) = (print(io, "$(rdate.count)*Repeat("), show(io, rdate.part), print(io,")"))

"""
    CalendarAdj(calendars, rdate::RDate, rounding::HolidayRoundingConvention)

Apply a calendar adjustment to an underlying rdate, applying an appropriate convention if
our final date falls on a holiday.

The calendars can only use alphanumeric characters, plus `/`, `-` and ` `.

In the string-form, you can apply a calendar adjustment using the `@` character and
provide `|` separated calendar names to apply it on. The convention is finally provided in
square brackets using its string-form name.

### Examples
```julia-repl
julia> cal_mgr = SimpleCalendarManager()
julia> setcalendar!(cal_mgr, "WEEKEND", WeekendCalendar())
julia> apply(rd"1d", Date(2019,9,27), cal_mgr)
2019-09-28
julia> apply(rd"1d@WEEKEND[NBD]", Date(2019,9,27), cal_mgr)
2019-09-30
julia> apply(rd"2m - 1d", Date(2019,7,23), cal_mgr)
2019-09-22
julia> apply(rd"(2m - 1d)@WEEKEND[PBD]", Date(2019,7,23), cal_mgr)
2019-09-20
```
"""
@auto_hash_equals struct CalendarAdj{R <: RDate, S <: HolidayRoundingConvention} <: RDate
    calendar_names::Vector{String}
    part::R
    rounding::S

    CalendarAdj(calendar_names, part::R, rounding::S) where {R <: RDate, S <: HolidayRoundingConvention} = new{R, S}(calendar_names, part, rounding)
    CalendarAdj(calendar_names::String, part::R, rounding::S) where {R<:RDate, S <: HolidayRoundingConvention} = new{R,S}(split(calendar_names, "|"), part, rounding)
end

function apply(rdate::CalendarAdj, date::Dates.Date, cal_mgr::CalendarManager)
    base_date = apply(rdate.part, date, cal_mgr)
    cal = calendar(cal_mgr, rdate.calendar_names)
    apply(rdate.rounding, base_date, cal)
end
multiply(rdate::CalendarAdj, count::Integer) = CalendarAdj(rdate.calendar_names, multiply(rdate.part, count), rdate.rounding)

Base.:-(x::CalendarAdj) = CalendarAdj(x.calendar_names, -x.part, x.rounding)
Base.show(io::IO, rdate::CalendarAdj) = print(io, "($(rdate.part))@$(join(rdate.calendar_names, "|"))[$(rdate.rounding)]")

"""
    Next(parts, inclusive::Bool = false)

Next is a mechanism through which we can find the next closest date in the future, given a list of
rdates to apply. We can choose whether today is also deemed a valid date.

This is commonly used in conjunction with rdates which don't necessarily always give a date in the
future, such as asking for the next Easter from today.

### Examples
```julia-repl
julia> RDates.Next([RDates.Easter(0), RDates.Easter(1)]) + Date(2019,1,1)
2019-04-21
julia> rd"Next(0E,1E)" + Date(2019,1,1)
2019-04-21
julia> RDates.Next([RDates.Easter(0), RDates.Easter(1)]) + Date(2019,4,21)
2020-04-12
julia> RDates.Next([RDates.Easter(0), RDates.Easter(1)], true) + Date(2019,4,21)
2019-04-21
julia> rd"Next!(0E,1E)" + Date(2019,4,21)
2019-04-21
```

!!! note
    The negation of `Next` will actually produce a `Previous`.

    ```julia-repl
    julia> -rd"Next(1d,2d)"
    Previous(-1d, -2d)
    julia> -3 * rd"Next(1d,2d)"
    Previous(-3d, -6d)
    ```

!!! warning
    While `Next` is a powerful operator, it does require application of every rdate
    every time, so can be expensive.

    When combining with ranges, it can often be useful to use an appropriate pivot
    point to start from instead.
"""
@auto_hash_equals struct Next <: RDate
    parts::Vector{RDate}
    inclusive::Bool

    Next(parts, inclusive::Bool = false) = new(parts, inclusive)
end

"""
    Previous(parts, inclusive::Bool = false)

Previous is a mechanism through which we can find the next closest date in the past, given a list of
rdates to apply. We can choose whether today is also deemed a valid date.

This is commonly used in conjunction with rdates which don't necessarily always give a date in the
future, such as asking for the previous Easter from today.

### Examples
```julia-repl
julia> RDates.Previous([RDates.Easter(0), RDates.Easter(-1)]) + Date(2019,12,31)
2019-04-21
julia> rd"Previous(0E,-1E)" + Date(2019,12,31)
2019-04-21
julia> RDates.Previous([RDates.Easter(0), RDates.Easter(-1)]) + Date(2019,4,21)
2018-04-01
julia> RDates.Previous([RDates.Easter(0), RDates.Easter(-1)], true) + Date(2019,4,21)
2019-04-21
julia> rd"Previous!(0E,-1E)" + Date(2019,4,21)
2019-04-21
```

!!! note
    The negation of `Previous` will actually produce a `Next`.

    ```julia-repl
    julia> -rd"Previous(-1d,-2d)"
    Next(1d, 2d)
    julia> -3 * rd"Previous(-1d,-2d)"
    Next(3d, 6d)
    ```

!!! warning
    While `Previous` is a powerful operator, it does require application of every rdate
    every time, so can be expensive.

    When combining with ranges, it can often be useful to use an appropriate pivot
    point to start from instead.
"""
@auto_hash_equals struct Previous <: RDate
    parts::Vector{RDate}
    inclusive::Bool

    Previous(parts, inclusive::Bool = false) = new(parts, inclusive)
end

function apply(rdate::Next, date::Dates.Date, cal_mgr::CalendarManager)
    dates = map(part -> apply(part, date, cal_mgr), rdate.parts)
    op = rdate.inclusive ? Base.:>= : Base.:>
    dates = [d for d in dates if op(d, date)]
    length(dates) > 0 || error("$rdate failed to apply to $date")
    min(dates...)
end

function apply(rdate::Previous, date::Dates.Date, cal_mgr::CalendarManager)
    dates = map(part -> apply(part, date, cal_mgr), rdate.parts)
    op = rdate.inclusive ? Base.:<= : Base.:<
    dates = [d for d in dates if op(d, date)]
    length(dates) > 0 || error("$rdate failed to apply to $date")
    max(dates...)
end

function multiply(rdate::Next, count::Integer)
    method = count >= 0 ? Next : Previous
    method(map(x -> multiply(x, count), rdate.parts), rdate.inclusive)
end

function multiply(rdate::Previous, count::Integer)
    method = count >= 0 ? Previous : Next
    method(map(x -> multiply(x, count), rdate.parts), rdate.inclusive)
end

Base.:-(rdate::Next) = Previous(map(-, rdate.parts), rdate.inclusive)
Base.:-(rdate::Previous) = Next(map(-, rdate.parts), rdate.inclusive)

function Base.show(io::IO, rdate::Next)
    incsign = rdate.inclusive ? "!" : ""
    print(io, "Next$incsign(")
    for (i,part) in enumerate(rdate.parts)
        show(io, part)
        if i < length(rdate.parts)
            print(io, ", ")
        end
    end
    print(io, ")")
end

function Base.show(io::IO, rdate::Previous)
    incsign = rdate.inclusive ? "!" : ""
    print(io, "Previous$incsign(")
    for (i,part) in enumerate(rdate.parts)
        show(io, part)
        if i < length(rdate.parts)
            print(io, ", ")
        end
    end
    print(io, ")")
end

ext: jl


text:
import Laplacians.condNumber
import Laplacians.sddmWrapLap

using StaticArrays

#=

A variations of the approxChol Laplacian solver
by Rasmus Kyng and Daniel A. Spielman.
It solves systems in non-singular BDD (Block Diagonally Dominant) Matrices.
This algorithm is an implementation of an approximate edge-by-edge elimination
algorithm inspired by the Approximate Gaussian Elimination algorithm of
Kyng and Sachdeva.

TODO:

- get rid of mat to analyze, and just build it in to llmatp
- Check on op norm, where faster and how call fastest.

 - Figure out what to do with singular Matrices
 - Make work with full matrices.

 - Test on approximate potential orths
 - Remove logging

=#

"""
    params = AC_Params(queue_forward, stag_test)

queue_forward=true if we use the forward queue (the original)
stag_test is a parameter for PCG, default 5.
"""
mutable struct AC_Params
    queue_forward::Bool
    stag_test::Integer
    useop::Bool
    perm::Bool
end

AC_Params() = AC_Params(true, 5, true, false)
AC_Params(b::Bool) = AC_Params(b, 5, true, false)
AC_Params(qf::Bool, uop::Bool) = AC_Params(qf, 5, uop, false)

#=
  For small matrices, the following is faster than "opnorm".
=#
_ournorm(x::StaticArrays.SArray) = sqrt(maximum(eigvals(x*x')))
ournorm(x::Real) = abs(x)
ournorm(mat::SMatrix{1,1}) = abs(mat[1])
ournorm(mat::SMatrix{2,2}) = _ournorm(mat)
ournorm(mat::SMatrix{3,3}) = _ournorm(mat)
ournorm(mat::SMatrix{4,4}) = _ournorm(mat)
ournorm(mat::SMatrix{5,5}) = _ournorm(mat)
ournorm(mat::SMatrix{6,6}) = _ournorm(mat)
ournorm(mat) = opnorm(mat)

function LLmata(a::SparseMatrixCSC{Tval,Tind}) where {Tind,Tval}
    n = size(a,1)
    m = nnz(a)

    degs = zeros(Tind,n)
    ptrs = copy(a.colptr)

    cols = a.colptr[1:n]
    llelems = Vector{LLa{Tind,Tval}}(undef, m)

    for i in 1:n
        degs[i] = a.colptr[i+1] - a.colptr[i]

        for ind in a.colptr[i]:(a.colptr[i+1]-1)
            j = a.rowval[ind]

            if j < i
                v = a.nzval[ind]

                if ptrs[i] < a.colptr[i+1]-1
                    llelems[ptrs[i]] = LLa(j, v, ptrs[i]+1, ptrs[j])
                else
                    llelems[ptrs[i]] = LLa(j, v, zero(Tind), ptrs[j])
                end

                if ptrs[j] < a.colptr[j+1]-1
                    llelems[ptrs[j]] = LLa(i, v, ptrs[j]+1, ptrs[i])
                else
                    llelems[ptrs[j]] = LLa(i, v, zero(Tind), ptrs[i])
                end

                ptrs[j] += 1
                ptrs[i] += 1
            end
        end

    end

    return LLmata{Tind,Tval}(n, degs, cols, llelems)
end

function print_ll_col(llmat::LLmata, i::Int)
    next = llmat.cols[i]
    while next != 0
      ll = llmat.lles[next]
      println("col $i, row $(ll.row) : $(ll.val)")

      next = ll.next
    end
end

#=============================================================

The approximate factorization

=============================================================#

function get_ll_col(llmat::LLmata{Tind,Tval},
    i,
    colspace::Vector{LLcol{Tind,Tval}}) where {Tind,Tval}

      ptr = llmat.cols[i]
      len = 0
      while ptr != 0

          if llmat.lles[ptr].val > zero(Tval)
              len = len+1

              item = LLcol(llmat.lles[ptr].row, ptr, llmat.lles[ptr].val)
              if (len > length(colspace))
                  push!(colspace,item)
              else
                  colspace[len] = item
              end
          end

          ptr = llmat.lles[ptr].next
      end

      return len
  end

#=
Add together entries of the column that are in the same row.
Decrease the degrees of vertices when their entries are added.
Returns the number of nonzero entries in the column after the compression, `ptr`.
=#

function compress_col!(llmat::LLmata{Tind,Tval},
    colspace::Vector{LLcol{Tind,Tval}},
    len::Tind,
    pq::acPQ{Tind},
    col::Int,
    compressions) where {Tind,Tval}

    o1 = Base.Order.ord(isless, x->x.row, false, Base.Order.Forward)

    sort!(colspace, one(len), len, QuickSort, o1)

    c = colspace

    ptr = 0
    currow = c[1].row
    curval = c[1].val
    curptr = c[1].ptr

    for i in 2:len

        if c[i].row != currow

            ptr = ptr+1
            c[ptr] = LLcol(currow, curptr, curval)

            currow = c[i].row
            curval = c[i].val
            curptr = c[i].ptr

        else

            curval = curval + c[i].val
            push!(compressions, Compression(col ,curptr, c[i].ptr))

            zeroval(llmat, llmat.lles[c[i].ptr].reverse)
            pqDec!(pq, currow)
        end

    end

    # emit the last row

    ptr = ptr+1
    c[ptr] = LLcol(currow, curptr, curval)

    o2 = Base.Order.ord(isless, x->x.val, false, Base.Order.Forward)
    sort!(colspace, one(ptr), ptr, QuickSort, o2)

    return ptr
end

const order_forward = Base.Order.ord(isless, identity, false, Base.Order.Forward)

"""
    a = mat_to_analyze(B::BlockCSC)

Produces an adjacency matrix of a graph whose edge weights are the norms
of matrices in B.  This could probably be accelerated.
"""
function mat_to_analyze(B::BlockCSC)

    mat = B.mat

    a = SparseMatrixCSC(mat.m, mat.n, copy(mat.colptr), copy(mat.rowval), ournorm.(mat.nzval))
    for i in 1:a.n
        a[i,i] = 0.0
    end
    dropzeros!(a)

    return a
end

function zeroval(llmat::LLmata{Tind,Tval}, revj::Tind) where {Tind, Tval}
    llmat.lles[revj] = LLa(llmat.lles[revj].row,
            zero(Tval),
            llmat.lles[revj].next ,
            llmat.lles[revj].reverse )
end

#=

Create edge (k,j), using the memory for (i,j)
=#
function rotate_edges(llmat::LLmata,k,j,ptr,revj,newEdgeVal)
    llmat.lles[revj] = LLa(k, newEdgeVal, llmat.lles[revj].next, ptr)

    # fix row j in col k
    khead = llmat.cols[k]
    llmat.cols[k] = ptr

    llmat.lles[ptr] = LLa(j, newEdgeVal, khead, revj)

end

function compute_elim(a::LLmata{Tind,Tval}; params = AC_Params()) where {Tind,Tval}

    #logstart(fn)

    n = a.n

    ijks = Vector{IJKop2{Tind,Tval}}(undef,0)

    #if params.queue_forward
        pq = acPQ(a.degs)
    #else
    #    pq = acPQ_rev(a.degs)
    #end

    colspace = Vector{LLcol{Tind,Tval}}(undef,n)
    cumspace = Vector{Tval}(undef,n)

    # elim_rank[i] is the iter in which vertex i is eliminated
    elim_rank = Vector{Tind}(undef,n)
    order_tmp = Vector{Tind}(undef, n)

    it = 1

    # DAS: might be computable from other terms, in which case can remove.
    diag_tmp = Vector{Tval}(undef,n)
    diag_adjusts = zeros(Tval, n)

    compressions = Vector{Compression}()

#    @inbounds while it < n
    while it < n

        i = pqPop!(pq)
        elim_rank[i] = it

        it = it + 1

        elim_vertex(a, i, ijks, colspace, cumspace,
        diag_tmp, diag_adjusts, order_tmp, pq, compressions, params)

    end

    last_col = pqPop!(pq)
    elim_rank[last_col] = n

    return ijks, diag_adjusts, elim_rank, last_col, compressions
end

function elim_vertex(llmat::LLmata{Tind,Tval},
    i::Int,
    ijks::Vector{IJKop2{Tind,Tval}},
    colspace::Vector{LLcol{Tind,Tval}},
    cumspace::Vector{Tval},
    diag_tmp::Vector{Tval},
    diag_adjusts::Vector{Tval},
    order_tmp::Vector{Tind},
    pq,
    compressions::Vector{Compression},
    params::AC_Params) where {Tind,Tval}

    len = get_ll_col(llmat, i, colspace)
    if len == 0
        return
    end

    len = compress_col!(llmat, colspace, len, pq, i, compressions)
    if len == 0
        return
    end

    #println("$(i): $(len)")

    # csum is the weighted degree of a vertex.
    # cumspace holds the cumulative sums of weights of edges, so that we can sample
    csum = zero(Tval)
    for ii in 1:len
        csum = csum + colspace[ii].val
        cumspace[ii] = csum
    end

    # we store adjustments to diagonals locally in diag_tmp, and then write them to
    # diagonal adjusts later.  This line sets the default adjusts, for when there are no incoming
    # edges to the vertex.
    for hoff in 2:len
        diag_tmp[hoff] = -(cumspace[hoff-1]) * colspace[hoff].val / csum
    end

    # Here we choose all the edges to be sampled.
    # This could be done in bulk, maybe faster, or in line with the next loop.
    for joffset in 1:(len-1)
        r = rand() * (csum - cumspace[joffset]) + cumspace[joffset]
        order_tmp[joffset] = searchsortedfirst(cumspace,r,joffset+1,len,order_forward)
    end

    # loop over the edges attached to vertex i, doing eliminations.
    for joffset in 1:(len-1)

        # are eliminating edge (i,j), and creating edge (j,k)

        ll = colspace[joffset]

        j = ll.row
        revj = llmat.lles[ll.ptr].reverse

        # (j,k) is the edge added.
        koff = order_tmp[joffset]
        k = colspace[koff].row

        pqInc!(pq, k)

        # this is the factor by which we multiply row i to elim row j
        op = - ll.val / csum

        # put in the new edge (j, k)
        oldEdgeVal = colspace[koff].val

        # this is what we really need to know for the factorization:
        # the probability that we chose this edge (j, k)
        probadd = oldEdgeVal / (csum-cumspace[joffset])

        # the following, commented out, form is equivalent
        # newEdgeVal = op * oldEdgeVal  / probadd
        newEdgeVal = (csum-cumspace[joffset]) * ll.val / csum

        # Compensate diagonal k for having sampled this edge
        diag_tmp[koff] += newEdgeVal

        rotate_edges(llmat,k,j,ll.ptr,revj,newEdgeVal)

        push!(ijks, IJKop2(i,j,k, probadd, ll.ptr, colspace[koff].ptr, revj))

    end # for joffset

    # adjust diagonals from sampling
    # takes a lot of time: all cache misses
    for hoff in 2:len
        h = colspace[hoff].row
        diag_adjusts[h] += diag_tmp[hoff]
    end

    # For the degre 1 case - last edge - just adjust diagonal

    ll = colspace[len]
    j = ll.row
    revj = llmat.lles[ll.ptr].reverse

    #= should firgure out why wanted it < n
    if it < n
        pqDec!(pq, j)
    end
    =#

    pqDec!(pq, j)

    # zero out that edge for the other side.
    zeroval(llmat, revj)

    # zero in k position indicates no edge
    push!(ijks, IJKop2(i,j,0, 1.0, ll.ptr, revj, 0))

end

function bdd_setup(B::BlockCSC{D,Tv,L}, elim_rank) where {D, Tv, L}
    mat = B.mat

    n = size(mat,1)
    cnt = count_nonzero_diags(mat)

    m = nnz(mat) - cnt

    matstore = Vector{StaticArrays.SArray{Tuple{D,D},Tv,2,L}}(undef, m)
    matstore = zeros(StaticArrays.SArray{Tuple{D,D},Tv,2,L}, m)
    dgs = Vector{StaticArrays.SArray{Tuple{D,D},Tv,2,L}}(undef, n)

    diag_free_ind = 0

    for i in 1:n

        for ind in mat.colptr[i]:(mat.colptr[i+1]-1)
            j = mat.rowval[ind]

            if i == j
                dgs[i] = mat.nzval[ind]
            else
                diag_free_ind += 1

                # SOMETHING IS WRONG IF NEED TO COMMENT THIS OUT!
                #if elim_rank[i] < elim_rank[j]
                    matstore[diag_free_ind] = mat.nzval[ind]
                #end
            end
        end
    end

    return matstore, dgs
end

"""
This makes a BlockCSC with structural zeros and the elements of bdd that we will use
in the elimination.  It allows us to do the factorizations just once.
"""
function prepare_triangular_factor(bdd::BlockCSC{D,Tv,L}, ijks, elim_rank) where {D,Tv,L}
    n = size(bdd,1)
    m = length(ijks)

    ti = Vector{Int}(undef, m+n)
    tj = Vector{Int}(undef, m+n)
    tv = zeros(typeof(bdd[1,1]), m+n)
    ptr = 1
    i = 0
    for it in 1:m
        ijk = ijks[it]
        if ijk.i != i
            i = ijk.i
            ti[ptr] = i
            tj[ptr] = i
            #tv[ptr] = la[i,i]
            ptr += 1
        end
        ti[ptr] = ijk.i
        tj[ptr] = ijk.j
        ptr += 1
    end
    lst = argmax(elim_rank)
    ti[ptr] = lst
    tj[ptr] = lst
    bdd_tri = BlockCSC(sparse(tj,ti,tv,n,n))

    aj, ai, av = findnz(bdd.mat)
    for it in 1:length(ai)
        if elim_rank[ai[it]] <= elim_rank[aj[it]]
            bdd_tri[aj[it], ai[it]] = av[it]
        #else
        #    bdd_tri[ai[it], aj[it]] = av[it]
        end
    end

    return bdd_tri
end

function factor_bdd_ijks_plt(ijks, diag_adjusts,
    elim_rank,
    matstore::Vector{StaticArrays.SArray{Tuple{D,D},Tv,2,L}}, dgs, comps) where {D,Tv,L}

    Tblock = StaticArrays.SArray{Tuple{D,D},Tv,2,L}

    n = length(dgs)
    colptr = Vector{Int}(undef, n + 1)
    rowval = Vector{Int}(undef, length(ijks) + n)
    nzval = Vector{Tblock}(undef, length(ijks) + n)
    colptr_ind = 1
    rowval_ind = 1

    E = one(Tblock)

    diaginv = E

    # true when are first entry in col.  Would be better to have a cscish ijks, and skip this.
    flag = true

    comps_ind = 1

    for it in 1:length(ijks)
        ijk = ijks[it]

        if flag

            while (comps_ind <= length(comps)) && (comps[comps_ind].col == ijk.i)
                matstore[comps[comps_ind].into] += matstore[comps[comps_ind].from]
                comps_ind += 1
            end

            dgs[ijk.i] += E*diag_adjusts[ijk.i]

            # place the diagonal entry
            colptr[colptr_ind] = rowval_ind
            colptr_ind += 1
            rowval[rowval_ind] = ijk.i
            nzval[rowval_ind] = dgs[ijk.i]
            rowval_ind += 1

            diaginv = inv(dgs[ijk.i])

            flag = false
        end

        op = - matstore[ijk.ptr_ji] * diaginv

        dgs[ijk.j] += op * (matstore[ijk.ptr_ji])'

       # push out an off-diagonal entry
       rowval[rowval_ind] = ijk.j

       nzval[rowval_ind] = matstore[ijk.ptr_ji]
       rowval_ind += 1

        if ijk.k != 0 # if not a degree 1 elim: again better to use cscish ijks.

            if elim_rank[ijk.j] < elim_rank[ijk.k]
                matstore[ijk.ptr_kj] = (op * ( matstore[ijk.ptr_ki])')' / ijk.prob
            else
                matstore[ijk.ptr_ji] = op * ( matstore[ijk.ptr_ki])'/ ijk.prob
            end

        else

            flag = true
        end
    end

    lst = argmax(elim_rank)

    colptr[colptr_ind] = rowval_ind
    rowval[rowval_ind] = lst

    nzval[rowval_ind] = dgs[lst] + E*diag_adjusts[lst]

    colptr[colptr_ind+1] = rowval_ind+1

    return PermLTCSC(colptr, rowval, nzval)

end

# this will be removed from the main code line
function AC_BDD_precon(bdd; verbose=false, params = AC_Params())

    t1 = time()
    t0 = t1

    a = mat_to_analyze(bdd)
    if verbose
        println("Time to construct mat to analyze ", time()-t1)
        t1 = time()
    end

    ijks, diag_adjusts, elim_rank, last_col, comps = compute_elim(LLmata(a))
    if verbose
        println("Time to compute elimination structure ", time()-t1)
        t1 = time()
    end

    matstore, dgs = bdd_setup(bdd, elim_rank)
    if verbose
        println("Time to prepare for elimination ", time()-t1)
        t1 = time()
    end

    ltcsc = factor_bdd_ijks_plt(ijks, diag_adjusts, elim_rank, matstore, dgs, comps)
    if verbose
        println("Time to eliminate ", time()-t1)
        t1 = time()
    end

    F(b) = solver(ltcsc, b)

    if verbose
        println("Total solver build time: ", time()-t0)
        println("Ratio of operator edges to original edges: ", 2 * length(ltcsc.nzval) / nnz(bdd))
    end

    return F
end

function solver(B, ltcsc;
    pcg_params = PCG_params())

    F!(b) = solver!(ltcsc, b)

    #=
    f1(b;tol=tol_,maxits=maxits_, maxtime=maxtime_,
        verbose=verbose_,
        pcgIts=pcgIts_) =
    pcg(Bp, b[elim_order], F, tol=tol, maxits=maxits, maxtime=maxtime,
        pcgIts=pcgIts,
        verbose=verbose,
        stag_test = params.stag_test)[elim_rank]
=#

    f1(b;tol=pcg_params.tol,
        maxits=pcg_params.maxits,
        maxtime=pcg_params.maxtime,
        verbose=pcg_params.verbose,
        pcgIts=pcg_params.pcgIts) =
    pcg_inplace(B, b, F!,
        tol=tol,
        maxits=maxits,
        maxtime=maxtime,
        pcgIts=pcgIts,
        verbose=verbose,
        stag_test = pcg_params.stag_test)

    return f1
end

function permed_solver(B, ltcsc, elim_rank::Vector{Int};
    pcg_params = PCG_params())

    elim_order = invperm(elim_rank)
    Bp = BlockCSC(B.mat[elim_order, elim_order])

    # note: the following does not copy, but directs pointer
    ltcsc_p = LTCSC(ltcsc.colptr, elim_rank[ltcsc.rowval], ltcsc.nzval)

    F(b) = solver(ltcsc_p, b)

    #=
    f1(b;tol=tol_,maxits=maxits_, maxtime=maxtime_,
        verbose=verbose_,
        pcgIts=pcgIts_) =
    pcg(Bp, b[elim_order], F, tol=tol, maxits=maxits, maxtime=maxtime,
        pcgIts=pcgIts,
        verbose=verbose,
        stag_test = params.stag_test)[elim_rank]
=#

    f1(b; tol=pcg_params.tol,
        maxits=pcg_params.maxits,
        maxtime=pcg_params.maxtime,
        verbose=pcg_params.verbose,
        pcgIts=pcg_params.pcgIts) =
    pcg_perm(elim_rank, Bp, copy(b), F,
        PCG_params(verbose, tol, maxits, maxtime, pcgIts, pcg_params.stag_test )
        )

    return f1
end

function pcg_perm(perm,
    mat::BlockCSC{D,Tv,L},
    b::Vector{SArray{Tuple{D},Tv,1,D}},
    pre::Function,
    pcg_params::PCG_params)  where {D,Tv,L}

    bp = b[perm]
    x = similar(b)

    x[perm] = pcg(mat, bp, pre, tol=pcg_params.tol,
        maxits=pcg_params.maxits,
        maxtime=pcg_params.maxtime,
        pcgIts=pcg_params.pcgIts,
        verbose=pcg_params.verbose,
        stag_test = pcg_params.stag_test)

    return x

end

"""
    solver = approxChoBDD(bdd); x = solver(b);
    solver = approxCholBDD(bdd; tol::Real=1e-6, maxits=1000, maxtime=Inf, verbose=false, pcgIts=Int[])

A heuristic by Daniel Spielman inspired by the linear system solver in https://arxiv.org/abs/1605.02353 by Rasmus Kyng and Sushant Sachdeva.  Whereas that paper eliminates vertices one at a time, this eliminates edges one at a time.  It is probably possible to analyze it.
It solves systems of equations in nonsingular Block Diagonally Dominant Matrices.
ADD NOTES ON HOW TO CONSTRUCT THESE.

For more info, see http://danspielman.github.io/Laplacians.jl/latest/usingSolvers/index.html
"""
function approxCholBDD(bdd;
    tol=1e-6,
    maxits=1000,
    maxtime=Inf,
    verbose=false,
    pcgIts=Int[],
    params = AC_Params())

    pcg_params = PCG_params(verbose,
        tol, maxits, maxtime, pcgIts,
        params.stag_test)

    tol_ =tol
    maxits_ =maxits
    maxtime_ =maxtime
    verbose_ =verbose
    pcgIts_ =pcgIts

    t1 = time()
    t0 = t1

    #M = sparse(bdd)

    a = mat_to_analyze(bdd)
    if verbose
        println("Time to construct mat to analyze ", time()-t1)
        t1 = time()
    end

    ijks, diag_adjusts, elim_rank, last_col, comps = compute_elim(LLmata(a))
    if verbose
        println("Time to compute elimination structure ", time()-t1)
        t1 = time()
    end

    matstore, dgs = bdd_setup(bdd, elim_rank)
    if verbose
        println("Time to prepare for elimination ", time()-t1)
        t1 = time()
    end

    ltcsc = factor_bdd_ijks_plt(ijks, diag_adjusts, elim_rank, matstore, dgs, comps)
    if verbose
        println("Time to eliminate ", time()-t1)
        t1 = time()
    end

    permute_outside = params.perm

    if permute_outside

        f1 = permed_solver(bdd, ltcsc, elim_rank,
            pcg_params = pcg_params)

        if verbose
            println("Time to permute ", time()-t1)
            t1 = time()
        end

    else

        f1 = solver(bdd,  ltcsc,
            pcg_params = pcg_params)

    end

    if verbose
        println("Total solver build time: ", time()-t0)
        println("Ratio of operator edges to original edges: ", 2 * length(ltcsc.nzval) / nnz(bdd))
    end

    return f1
end


text:
import Base: setindex!, getindex, haskey, pairs, keys, values

@enum DType begin
  FLOAT
  DOUBLE
  HALF
  UINT8
  INT8
  SHORT
  INT
  LONG
  BOOL
end

jltype2dtype(::Type{Float32}) = FLOAT
jltype2dtype(::Type{Float64}) = DOUBLE
jltype2dtype(::Type{Float16}) = HALF
jltype2dtype(::Type{UInt8}) = UINT8
jltype2dtype(::Type{Int8}) = INT8
jltype2dtype(::Type{Int16}) = SHORT
jltype2dtype(::Type{Int32}) = INT
jltype2dtype(::Type{Int64}) = LONG
jltype2dtype(::Type{Bool}) = BOOL

function dtype2jltype(t::DType)
  if t == FLOAT
    return Float32
  elseif t == DOUBLE
    return Float64
  elseif t == HALF
    return Float16
  elseif t == UINT8
    return UInt8
  elseif t == INT8
    return Int8
  elseif t == SHORT
    return Int16
  elseif t == INT
    return Int32
  elseif t == LONG
    return Int64
  elseif t == BOOL
    return Bool
  end
end

function bytewidth(t::DType)
  if t == FLOAT
    return 4
  elseif t == DOUBLE
    return 8
  elseif t == HALF
    return 2
  elseif t == UINT8
    return 1
  elseif t == INT8
    return 1
  elseif t == SHORT
    return 2
  elseif t == INT
    return 4
  elseif t == LONG
    return 8
  elseif t == BOOL
    return 1
  end
end

function thtype2dtype(defer)
  type = defer.head
  if type == Symbol("torch.FloatStorage")
    return FLOAT
  elseif type == Symbol("torch.DoubleStorage")
    return DOUBLE
  elseif type == Symbol("torch.HalfStorage")
    return HALF
  elseif type == Symbol("torch.ByteStorage")
    return UINT8
  elseif type == Symbol("torch.CharStorage")
    return INT8
  elseif type == Symbol("torch.ShortStorage")
    return SHORT
  elseif type == Symbol("torch.IntStorage")
    return INT
  elseif type == Symbol("torch.LongStorage")
    return LONG
  elseif type == Symbol("torch.BoolStorage")
    return BOOL
  else
    error("unknown type: $type")
  end
end

struct StorageManager
  data::Dict{String, Tuple{DType, Int, String, Array}}
end

StorageManager() = StorageManager(Dict())
setindex!(sm::StorageManager, value, key) = setindex!(sm.data, value, key)
getindex(sm::StorageManager, key) = getindex(sm.data, key)
haskey(sm, key) = haskey(sm.data, key)
pairs(sm::StorageManager) = pairs(sm.data)
keys(sm::StorageManager) = keys(sm.data)
values(sm::StorageManager) = values(sm.data)


text:
# auxiliary/stridedarray.jl
#
# Simple auxiliary methods to interface with StridedArray from Julia Base.

"""`numind(A)`

Returns the number of indices of a tensor-like object `A`, i.e. for a multidimensional array (`<:AbstractArray`) we have `numind(A) = ndims(A)`. Also works in type domain.
"""
numind(A::AbstractArray) = ndims(A)
numind{T<:AbstractArray}(::Type{T}) = ndims(T)

"""`similar_from_indices(T, indices, A, conjA=Val{:N})`

Returns an object similar to `A` which has an `eltype` given by `T` and dimensions/sizes corresponding to a selection of those of `op(A)`, where the selection is specified by `indices` (which contains integer between `1` and `numind(A)`) and `op` is `conj` if `conjA=Val{:C}` or does nothing if `conjA=Val{:N}` (default).
"""
function similar_from_indices{T,CA}(::Type{T}, indices, A::StridedArray, ::Type{Val{CA}}=Val{:N})
    dims = size(A)
    return similar(A,T,dims[indices])
end

"""`similar_from_indices(T, indices, A, B, conjA=Val{:N}, conjB={:N})`

Returns an object similar to `A` which has an `eltype` given by `T` and dimensions/sizes corresponding to a selection of those of `op(A)` and `op(B)` concatenated, where the selection is specified by `indices` (which contains integers between `1` and `numind(A)+numind(B)` and `op` is `conj` if `conjA` or `conjB` equal `Val{:C}` or does nothing if `conjA` or `conjB` equal `Val{:N}` (default).
"""
function similar_from_indices{T,CA,CB}(::Type{T}, indices, A::StridedArray, B::StridedArray, ::Type{Val{CA}}=Val{:N}, ::Type{Val{CB}}=Val{:N})
    dims = tuple(size(A)...,size(B)...)
    return similar(A,T,dims[indices])
end

"""`scalar(C)`

Returns the single element of a tensor-like object with zero dimensions, i.e. if `numind(C)==0`.
"""
scalar(C::StridedArray) = numind(C)==0 ? C[1] : throw(DimensionMismatch())


text:
module Bridges

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

const CI = MOI.ConstraintIndex

include("bridge.jl")
include("bridge_optimizer.jl")

# Variable bridges
include("Variable/Variable.jl")
# Constraint bridges
include("Constraint/Constraint.jl")
# Objective bridges
include("Objective/Objective.jl")

include("lazy_bridge_optimizer.jl")
include("debug.jl")

"""
    full_bridge_optimizer(model::MOI.ModelLike, ::Type{T}) where T

Returns a `LazyBridgeOptimizer` bridging `model` for every bridge defined in
this package and for the coefficient type `T`.
"""
function full_bridge_optimizer(model::MOI.ModelLike, T::Type)
    bridged_model = LazyBridgeOptimizer(model)
    Variable.add_all_bridges(bridged_model, T)
    Constraint.add_all_bridges(bridged_model, T)
    Objective.add_all_bridges(bridged_model, T)
    return bridged_model
end

print_num_bridges(io::IO, ::Variable.EmptyMap) = nothing
print_num_bridges(io::IO, ::Constraint.EmptyMap) = nothing
print_num_bridges(io::IO, ::Objective.EmptyMap) = nothing
function print_num_bridges(io::IO, B::Variable.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    return print(io, "\n$(indent)with $(n) variable bridge$(s(n))")
end
function print_num_bridges(io::IO, B::Constraint.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    return print(io, "\n$(indent)with $(n) constraint bridge$(s(n))")
end
function print_num_bridges(io::IO, B::Objective.Map)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    n = length(B)
    return print(io, "\n$(indent)with $(n) objective bridge$(s(n))")
end

function Base.show(io::IO, B::AbstractBridgeOptimizer)
    MOIU.print_with_acronym(io, summary(B))
    print_num_bridges(io, Variable.bridges(B))
    print_num_bridges(io, Constraint.bridges(B))
    print_num_bridges(io, Objective.bridges(B))
    if :model in propertynames(B)
        indent = " "^get(io, :indent, 0)
        print(io, "\n$(indent)with inner model ")
        show(IOContext(io, :indent => get(io, :indent, 0) + 2), B.model)
    end
end

end # module


text:
module CanopyOptics

###### Julia packages to import  ############
using Distributions            # Distributions for angular distributions of canopy elements
using FastGaussQuadrature      # Quadrature points for numerical integration
using Unitful                  # Units attached to some variable
using UnitfulEquivalences      # Spectral conversions
using SpecialFunctions:expint  # expint in Prospect
using DelimitedFiles           # File IO (Prospect csv File)
using DocStringExtensions      # Documentation
using UnPack                   # @unpack and stuff
using LazyArtifacts            # Artifacts
using LinearAlgebra            # Well, guess...
using Polynomials              # Polynomials for some empirical functions

# Filename for ProspectPro optical properties
const OPTI_2021 = artifact"Prospect" * "/dataSpec_PRO.csv";

# Needs to be more modular later:
FT = Float64
###### Own files to include #################
include("initialization/constants.jl")

include("types/canopy_types.jl")
include("types/angle_types.jl")
include("types/material_types.jl")

include("utils/quadrature.jl")
include("utils/canopy_angles.jl")
include("utils/fresnel.jl")

include("core/prospect.jl")
include("core/leafAngleRoutines.jl")

include("initialization/loadProspect.jl")
include("initialization/default_constructors.jl")

include("utils/dielectric.jl")

export prospect
export createLeafOpticalStruct, LeafProspectProProperties, LeafOpticalProperties, dirVector, dirVector_μ
export AbstractCanopyScatteringType, BiLambertianCanopyScattering, SpecularCanopyScattering
export PureIce, LiquidPureWater, LiquidSaltWater
export LeafProspectProProperties, LeafOpticalProperties, dielectric
# Functions:
export compute_Z_matrices, prospect, compute_specular_reflection, G

end # module


text:
function agent_step!(agent::Agents.AbstractAgent, model::Agents.ABM)
    if !agent.quarantined
        if model.clock % 10 != 0
            move_agent_random!(agent, model)
        else
            move_agent_home!(agent, model)
        end
    end
    update_fear!(agent, model)
    if model.behavior_activated
        update_socialnorm_memory!(agent, model)
        if model.clock % 10 == 0
            update_socialnorm!(agent, model)
            choose_strategy!(agent, model)
        end
        if model.behavior_policy
            adapt_to_policy!(agent, model)
        end
    end
    return agent
end

function move_agent_random!(agent::Agents.AbstractAgent, model::Agents.ABM)
    adjacent_places = setdiff(node_neighbors(agent.pos, model), model.homes)
    new_position = StatsBase.sample(adjacent_places)
    Agents.move_agent!(agent, new_position, model)
    return agent
end

function move_agent_home!(agent::Agents.AbstractAgent, model::Agents.ABM)
    Agents.move_agent!(agent, agent.home, model)
    return agent
end

function update_fear!(agent::Agents.AbstractAgent, model::Agents.ABM)
    incidence_influence = model.Id_count / model.agentcount
    for i in 1:10
        if Random.rand() < agent.neuroticism
            incidence_influence = pump(incidence_influence)
        end
    end
    agent.fear = StatsBase.mean([agent.fear, incidence_influence])
    return agent
end

function pump(x)
    return tanh(x) / tanh(1)
end

function update_socialnorm_memory!(agent::Agents.AbstractAgent, model::Agents.ABM)
    neighbor_behavior = [
        i.behavior
        for i in Agents.get_node_agents(agent.pos, model)
        if !(i === agent)
    ]
    if length(neighbor_behavior) > 0
        append!(
            agent.socialnorm_memory,
            neighbor_behavior
        )
    end
    return agent
end

function update_socialnorm!(agent::Agents.AbstractAgent, model::Agents.ABM)
    if length(agent.socialnorm_memory) > 0
        agent.socialnorm = StatsBase.mean(agent.socialnorm_memory)
    end
    agent.socialnorm_memory = Any[]
    return agent
end

function choose_strategy!(agent::Agents.AbstractAgent, model::Agents.ABM)
    weight = compute_weight(agent, model)
    costmatrix = compute_costmatrix(weight, model)
    if (b_is_dominant(costmatrix))
        agent.behavior = true
        agent.behavior_buffer = true
    elseif (notb_is_dominant(costmatrix))
        agent.behavior = false
        agent.behavior_buffer = false
    else
        agent.behavior = biased_coin(weight)
        agent.behavior_buffer = copy(agent.behavior)
    end
    return agent
end

function compute_weight(agent::Agents.AbstractAgent, model::Agents.ABM)
    return StatsBase.mean([
        agent.fear,
        agent.socialnorm,
        agent.trust_authorities * model.authority_recommendation
    ])
end

function compute_costmatrix(weight::Float64, model::Agents.ABM)
    costmatrix = (
        [(1 - weight) (1 - weight); weight weight]
        .* model.cost_matrix.player1
    )
    return costmatrix
end

function b_is_dominant(costmatrix::AbstractArray)
    is_dominant = (
        (costmatrix[1, 1] < costmatrix[2, 1])
        & (costmatrix[1, 2] < costmatrix[2, 2])
    )
    return is_dominant
end

function notb_is_dominant(costmatrix::AbstractArray)
    is_dominant = (
        (costmatrix[1, 1] > costmatrix[2, 1])
        & (costmatrix[1, 2] > costmatrix[2, 2])
    )
    return is_dominant
end

function biased_coin(weight::Float64)
    return Random.rand() < weight ? true : false
end

function adapt_to_policy!(agent::Agents.AbstractAgent, model::Agents.ABM)
    if agent.pos in model.public_spaces
        agent.behavior = true
    else
        agent.behavior = copy(agent.behavior_buffer)
    end
    return agent
end


text:
using Statistics

"""
    bgcf()
---
Create scenarios for Theo's `bgcf` program.
"""
function bgcf()
    subtitle("A setup for Theo's bgcf")
    item("Parameters")
    isfile(joinpath(abg_rst, "toy.G")) || simu_qp()
    nlc, nid = countlines.(joinpath.(abg_dat, ["toy.bim", "toy.fam"]))
    open(joinpath(abg_rst, "bgcf.ph"), "w") do io
        for ph in eachline(joinpath(abg_rst, "toy.ph"))
            println(io, "$ph 1 1")
            done
        end
    end
    seed = rand(1:2^31)
    write(joinpath(abg_rst, "dseed.rec"), "$seed\n")
    giv = reshape(reinterpret(Float64, read(joinpath(abg_rst, "toy.giv"))), (nid, nid))
    open(joinpath(abg_rst, "g.giv"), "w") do io
        for i in 1:nid
            for j in 1:nid
                ele = giv[i, j]
                println(io, "$i $j $ele")
            end
        end
    end
    cp(joinpath(abg_dat, "toy.bim"), joinpath(abg_rst, "toy.bim"), force=true)
    cp(joinpath(abg_dat, "toy.bed"), joinpath(abg_rst, "toy.bed"), force=true)
    nchr = 1
    _π = 0.1
    ve = 0.2
    vsnp = 0.8 * 10 / nlc
    vu = 0.8
    cycle = 1000
    burnin = 100
    nworker = 1
    ithin = 0
    ispeed = 0

    inp_str=["&inputvalues",
               "Nanim=$nid",
               "Nchr=$nchr",
               "Nsnptotal=$nlc",
               "bimfil=toy.bim",
               "bedfil=toy.bed",
               "phenfil=bgcf.ph", # line: y weight class; weight=1 → equal weight
               "PI_value=$_π", # for BayesC analysis. π ∈ [0, 1].  π < 0 → estimate π
               "Ve=$ve",
               "Vsnp=$vsnp",  # Vsnp matrix; SNP ∼ N(0, Vsnp) 1% of Vg
               "Vu=$vu", # G=polygenic covariance matrix, 70% of Vg.  All 0 → no polygenic variance.
               "Ncycle=$cycle",   # Gibbs cycles
               "Nburnin=$burnin", # burnin cycles
               "Nwork=$nworker",  # number of Gibbs sampling threads
               "Ithin=$ithin", # thinning output (prt each ithin sample).  ithin = 0 → no sample output, only smry stats
               "Ispeedup=$ispeed", # reduce computations by factor ispeedup.  ispeedup = 0 → no reduction.
               "/"]
    write(joinpath(abg_rst, "bgcf.inp"), join(inp_str, '\n'), '\n')
    done()
end


text:
module adapters_odbc_mariadb_repo_test

using Test # @test
using Octo.Adapters.ODBC # Repo Schema
using MariaDB_Connector_C_jll

Repo.debug_sql()

PLUGIN_DIR = joinpath(MariaDB_Connector_C_jll.artifact_dir, "lib", "mariadb", "plugin")

Repo.connect(
    adapter  = Octo.Adapters.ODBC,
    dsn      = "Driver={ODBC_Test_MariaDB};SERVER=127.0.0.1;PLUGIN_DIR=$PLUGIN_DIR;Option=67108864;CHARSET=utf8mb4;USER=root;PWD="
)

# Repo.execute([DROP DATABASE IF EXISTS :mysqltest])
# Repo.execute([CREATE DATABASE :mysqltest])
Repo.execute([USE :mysqltest])

# https://discourse.julialang.org/t/odbc-example-for-prepared-statement-doesnt-work
struct CoolTable
end
Schema.model(CoolTable, table_name="cool_table", primary_key="ID")

Repo.execute([DROP TABLE IF EXISTS CoolTable])
Repo.execute(Raw("""
    CREATE TABLE cool_table (
        ID SERIAL,
        col1 INTEGER,
        col2 DOUBLE,
        col3 VARCHAR(255),
        PRIMARY KEY (ID)
    )"""))

Repo.insert!(CoolTable, [(col1=1, col2=4.0, col3="hey"),
                         (col1=2, col2=5.0, col3="there"),
                         (col1=3, col2=6.0, col3="sailor")])

df = Repo.query(CoolTable)
@info :df df
@test size(df) == (3,)
@test df[1].col3 == "hey"

struct Employee
end
Schema.model(Employee, table_name="Employee", primary_key="ID")

Repo.execute([DROP TABLE IF EXISTS Employee])
Repo.execute(Raw("""
    CREATE TABLE Employee (
        ID SERIAL,
        Name VARCHAR(255),
        Salary DOUBLE,
        PRIMARY KEY (ID)
    )"""))

changes = (Name="John", Salary=10000.50)
Repo.insert!(Employee, changes)
changes = (Name="Cloris", Salary=20000.50)
Repo.insert!(Employee, changes)
df = Repo.query(Employee)
@info :df df
@test size(df) == (2,)
@test df[1].Name == "John"

Repo.disconnect()

end # module adapters_odbc_mariadb_repo_test


text:
for T in (:Int8, :Int16, :Int32, :Int64, :Int128, :Float16, :Float32, :Float64)
    @eval begin
        promote_type(::Type{ArbComplex}, ::Type{$T}) = ArbComplex
        promote_type(::Type{ArbReal}, ::Type{$T}) = ArbReal
        promote_type(::Type{ArbFloat}, ::Type{$T}) = ArbFloat

        promote_type(::Type{ArbComplex{P}}, ::Type{$T}) where {P} = ArbComplex{P}
        promote_type(::Type{ArbReal{P}}, ::Type{$T}) where {P} = ArbReal{P}
        promote_type(::Type{ArbFloat{P}}, ::Type{$T}) where {P} = ArbFloat{P}

        promote_rule(::Type{ArbComplex}, ::Type{$T}) = ArbComplex
        promote_rule(::Type{ArbReal}, ::Type{$T}) = ArbReal
        promote_rule(::Type{ArbFloat}, ::Type{$T}) = ArbFloat

        promote_rule(::Type{ArbComplex{P}}, ::Type{$T}) where {P} = ArbComplex{P}
        promote_rule(::Type{ArbReal{P}}, ::Type{$T}) where {P} = ArbReal{P}
        promote_rule(::Type{ArbFloat{P}}, ::Type{$T}) where {P} = ArbFloat{P}
    end
end

promote_type(::Type{ArbComplex}, ::Type{ArbReal}) = ArbComplex
promote_type(::Type{ArbComplex}, ::Type{ArbFloat}) = ArbComplex
promote_type(::Type{ArbReal}, ::Type{ArbFloat}) = ArbReal

promote_rule(::Type{ArbComplex}, ::Type{ArbReal}) = ArbComplex
promote_rule(::Type{ArbComplex}, ::Type{ArbFloat}) = ArbComplex
promote_rule(::Type{ArbReal}, ::Type{ArbFloat}) = ArbReal

promote_type(::Type{ArbComplex{P}}, ::Type{ArbReal{P}}) where {P} = ArbComplex{P}
promote_type(::Type{ArbComplex{P}}, ::Type{ArbFloat{P}}) where {P} = ArbComplex{P}
promote_type(::Type{ArbReal{P}}, ::Type{ArbFloat{P}}) where {P} = ArbReal{P}

promote_rule(::Type{ArbComplex{P}}, ::Type{ArbReal{P}}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{ArbFloat{P}}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbReal{P}}, ::Type{ArbFloat{P}}) where {P} = ArbReal{P}

promote_rule(::Type{ArbComplex{P}}, ::Type{Mag}) where {P} = ArbComplex{P}
promote_rule(::Type{ArbReal{P}}, ::Type{Mag}) where {P} = ArbReal{P}
promote_rule(::Type{ArbFloat{P}}, ::Type{Mag}) where {P} = ArbFloat{P}

promote_type(::Type{ArbFloat{P}}, ::Type{BigInt}) where {P} = ArbFloat{P}
promote_type(::Type{ArbReal{P}}, ::Type{BigInt}) where {P} = ArbReal{P}
promote_type(::Type{ArbComplex{P}}, ::Type{BigInt}) where {P} = ArbComplex{P}

promote_type(::Type{ArbFloat{P}}, ::Type{Rational{I}}) where {P,I} = ArbFloat{P}
promote_type(::Type{ArbReal{P}}, ::Type{Rational{I}}) where {P,I} = ArbReal{P}
promote_type(::Type{ArbComplex{P}}, ::Type{Rational{I}}) where {P,I} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{BigInt}) where {P} = promote_type(ArbFloat{P}, BigInt)
promote_rule(::Type{ArbReal{P}}, ::Type{BigInt}) where {P} = promote_type(ArbReal{P}, BigInt)
promote_rule(::Type{ArbComplex{P}}, ::Type{BigInt}) where {P} = promote_type(ArbComplex{P}, BigInt)

promote_rule(::Type{ArbFloat{P}}, ::Type{Rational{I}}) where {P,I} = promote_type(ArbFloat{P}, Rational{I})
promote_rule(::Type{ArbReal{P}}, ::Type{Rational{I}}) where {P,I} = promote_type(ArbReal{P}, Rational{I})
promote_rule(::Type{ArbComplex{P}}, ::Type{Rational{I}}) where {P,I} = promote_type(ArbComplex{P}, Rational{I})

promote_type(::Type{BigFloat},::Type{ArbFloat{P}}) where {P} = ArbFloat{P}
promote_type(::Type{BigFloat},::Type{ArbReal{P}}) where {P} = ArbReal{P}
promote_type(::Type{BigFloat},::Type{ArbComplex{P}}) where {P} = ArbComplex{P}
promote_type(::Type{ArbFloat{P}}, ::Type{BigFloat}) where {P} = ArbFloat{P}
promote_type(::Type{ArbReal{P}}, ::Type{BigFloat}) where {P} = ArbReal{P}
promote_type(::Type{ArbComplex{P}}, ::Type{BigFloat}) where {P} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbFloat{P}, BigFloat)
promote_rule(::Type{ArbReal{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbReal{P}, BigFloat)
promote_rule(::Type{ArbComplex{P}}, ::Type{BigFloat}) where {P} = promote_type(ArbComplex{P}, BigFloat)

promote_rule(::Type{ArbFloat{P}}, ::Type{I}) where {P, I<:Integer} = ArbFloat{P}
promote_rule(::Type{ArbReal{P}}, ::Type{I}) where {P, I<:Integer} = ArbReal{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{I}) where {P, I<:Integer} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbFloat{P}
promote_rule(::Type{ArbReal{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbReal{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{F}) where {P, F<:AbstractFloat} = ArbComplex{P}
promote_rule(::Type{ArbComplex{P}}, ::Type{C}) where {P, C<:Complex} = ArbComplex{P}

promote_rule(::Type{ArbFloat{P}}, x::ArbFloat{Q}) where {P,Q} =
    P>Q ? ArbFloat{P} : ArbFloat{Q}
promote_rule(::Type{ArbReal{P}}, x::ArbReal{Q}) where {P,Q} =
    P>Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbComplex{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}
promote_rule(::Type{ArbFloat{P}}, x::ArbReal{Q}) where {P,Q} =
    P>Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbFloat{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}
promote_rule(::Type{ArbReal{P}}, x::ArbComplex{Q}) where {P,Q} =
    P>Q ? ArbComplex{P} : ArbComplex{Q}

promote_rule(::Type{ArbFloat{P}},::Type{ArbFloat{Q}}) where {P,Q} = P<Q ? ArbFloat{P} : ArbFloat{Q}
promote_rule(::Type{ArbReal{P}},::Type{ArbReal{Q}}) where {P,Q} = P<Q ? ArbReal{P} : ArbReal{Q}
promote_rule(::Type{ArbComplex{P}},::Type{ArbComplex{Q}}) where {P,Q} = P<Q ? ArbComplex{P} : ArbComplex{Q}

promote_rule(::Type{<:AbstractIrrational}, ::Type{ArbFloat{P}}) where {P} = ArbFloat{P}
promote_rule(::Type{<:AbstractIrrational}, ::Type{ArbReal{P}}) where {P} = ArbReal{P}
promote_rule(::Type{<:AbstractIrrational}, ::Type{ArbComplex{P}}) where {P} = ArbComplex{P}

convert(::Type{ArbFloat}, x::ArbFloat{P}) where {P} = x
convert(::Type{ArbReal}, x::ArbReal{P}) where {P} = x
convert(::Type{ArbComplex}, x::ArbComplex{P}) where {P} = x
convert(::Type{ArbFloat}, x::ArbReal{P}) where {P} = ArbFloat{P}(x)
convert(::Type{ArbFloat}, x::ArbComplex{P}) where {P} = ArbFloat{P}(real(x))
convert(::Type{ArbReal}, x::ArbComplex{P}) where {P} = real(x)
convert(::Type{ArbReal}, x::ArbFloat{P}) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex}, x::ArbFloat{P}) where {P} = ArbComplex{P}(x)

convert(::Type{ArbReal{P}}, x::ArbComplex{P}) where {P} = real(x)
convert(::Type{ArbFloat{P}}, x::ArbComplex{P}) where {P} = midpoint(real(x), ArbFloat{P})
convert(::Type{Mag}, x::ArbComplex{P}) where {P} = radius(real(x), Mag)

convert(::Type{ArbComplex{P}}, x::ArbReal{P}) where {P} = ArbComplex{P}(x)
convert(::Type{ArbComplex{P}}, x::ArbFloat{P}) where {P} = ArbComplex{P}(ArbReal{P}(x))
convert(::Type{ArbComplex{P}}, x::Mag) where {P} = ArbComplex{P}(ArbFloat{P}(x))

convert(::Type{ArbFloat{P}}, x::ArbReal{P}) where {P} = midpoint(x, ArbFloat{P})
convert(::Type{Mag}, x::ArbReal{P}) where {P} = radius(x, Mag)

convert(::Type{ArbReal{P}}, x::ArbFloat{P}) where {P} = ArbReal{P}(x)
convert(::Type{ArbReal{P}}, x::Mag) where {P} = ArbReal{P}(ArbFloat{P}(x))

convert(::Type{Mag}, x::ArbFloat{P}) where {P} = Mag(x)
convert(::Type{ArbFloat{P}}, x::Mag) where {P} = ArbFloat{P}(x)

convert(::Type{ArbFloat{P}}, x::BigFloat) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::BigFloat) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::BigFloat) where {P} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::BigInt) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::BigInt) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::BigInt) where {P} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::I) where {P, I<:Integer} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::I) where {P, I<:Integer} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::I) where {P, I<:Integer} = ArbComplex{P}(x)

convert(::Type{ArbFloat{P}}, x::Rational{I}) where {P, I} = ArbFloat{P}(x.num) / ArbFloat{P}(x.den)
convert(::Type{ArbReal{P}}, x::Rational{I}) where {P, I} = ArbReal{P}(x.num) / ArbReal{P}(x.den)
convert(::Type{ArbComplex{P}}, x::Rational{I}) where {P, I} = ArbComplex{P}(ArbReal{P}(x))

function rationalize(::Type{I}, x::ArbFloat{P}) where {P,I<:Integer}
    oldprec = precision(BigFloat)
    setprecision(BigFloat, P)
    y = BigFloat(x)
    z = rationalize(I, y)
    setprecision(BigFloat, oldprec)
    return z
end

rationalize(::Type{I}, x::ArbReal{P}) where {P, I<:Integer} = rationalize(I, ArbFloat{P}(x))
rationalize(::Type{I}, x::ArbComplex{P}) where {P, I<:Integer} = Complex{I}(rationalize(I, ArbFloat{P}(real(x))), rationalize(I, ArbFloat{P}(imag(x))))

rationalize(x::ArbFloat{P}) where {P} = rationalize(Int64, x)
rationalize(x::ArbReal{P}) where {P} = rationalize(Int64, ArbFloat{P}(x))
rationalize(x::ArbComplex{P}) where {P} = Complex{Int64}(rationalize(Int64, ArbFloat{P}(real(x))), rationalize(Int64, ArbFloat{P}(imag(x))))

convert(::Type{Rational{I}}, x::ArbFloat{P}) where {P, I} = rationalize(I, x)
convert(::Type{Rational{I}}, x::ArbReal{P}) where {P, I} = rationalize(I, x)
convert(::Type{Rational{I}}, x::ArbComplex{P}) where {P, I} = rationalize(I, x)

convert(::Type{ArbFloat{P}}, x::AbstractFloat) where {P} = ArbFloat{P}(x)
convert(::Type{ArbReal{P}}, x::AbstractFloat) where {P} = ArbReal{P}(x)
convert(::Type{ArbComplex{P}}, x::AbstractFloat) where {P} = ArbComplex{P}(x)
convert(::Type{ArbComplex{P}}, x::Complex) where {P} = ArbComplex{P}(x)


text:
################################################################################
#                        BGV Scheme definition
################################################################################

struct BGVParams <: SHEShemeParams
    # The Cypertext ring over which operations are performed
    ℛ
    # The plaintext ring.
    ℛplain
    σ
end
scheme_name(p::Type{BGVParams}) = "BGV"
plaintext_modulus(p::BGVParams) = modulus(base_ring(p.ℛplain))

BGVParams(ring, p::Integer, σ) =
    BGVParams(ring, plaintext_space(ring, p), σ)

ℛ_plain(p::BGVParams) = p.ℛplain
ℛ_cipher(p::BGVParams) = p.ℛ

π⁻¹(params::BGVParams, plaintext) = convert(params.ℛ, params.ℛplain(plaintext))
function π(params::BGVParams, b)
    @fields_as_locals params::BGVParams
    ℛplain(map(x->coefftype(ℛplain)(convert(Integer, mod(SignedMod(x), plaintext_modulus(params)))), NTT.coeffs_primal(b)))
end

struct ShiftedDiscreteNormal
    p::Int
    dn::DiscreteNormal
end
Base.rand(rng::AbstractRNG, d::ShiftedDiscreteNormal) = d.p*rand(rng, d.dn)

𝒩(params::BGVParams) = RingSampler(params.ℛ, ShiftedDiscreteNormal(plaintext_modulus(params), DiscreteNormal(0, params.σ)))
𝒢(params::BGVParams) = RingSampler(params.ℛ, DiscreteNormal(0, params.σ))


text:
# Tests for Label and continuous output combiners
function t_combiners()

# test labels
L = [ 1 1 2 0;
      0 0 1 1;
      1 0 1 2 ]

# train labels
TL =( [1 2 0 0; # estimate labels from 3 ensemble members
      1 0 1 1;
      1 1 2 1],
     [1,2,0,0] ) # true labels

# Test for vote combiner
Test.@test all((L |> j4pr.votecombiner(size(L,1)))[:,1:3] .== [1 0 1])

# Test for weighted vote combiner
wcmb=TL|> j4pr.wvotecombiner(3) # 1'st line closest, Inf weight to 1'st ensemble memeber
Test.@test all((L |> wcmb) .== L[1:1,:])# the result is identical to the first line of L

# Test for Naive Bayes combiner
nbcmb=TL |> j4pr.naivebayescombiner(3)
Test.@test all((L |> nbcmb) .== [1 1 1 1])

# Test for Naive Bayes combiner - add 2 more label output combinations
#to L that were seen in training
bkscmb=TL |> j4pr.bkscombiner(3)
Test.@test all(([L [1,1,1] [2,0,1]] |> bkscmb) .== [0 0 0 0 1 2])

# Test for the mean combiner (3 outputs, 1 continuous output)
Test.@test L |> j4pr.meancombiner(3,1) ≈ [2/3 1/3 4/3 1.0]

# test for the weighted mean combiner
Test.@test L |> j4pr.wmeancombiner(3,1,[1.0,2,3]) ≈ [(1*1.0+3*1)/3 (1*1.0)/3 (2*1+1*2+1*3)/3 (1*2+2*3)/3]

# Test for the product combiner
Test.@test all((L |> j4pr.productcombiner(3,1)) .== [0 0 2 0])

# Test for the median combiner
Test.@test all((L |> j4pr.mediancombiner(3,1)) .== [1 0 1 1])

# Test for the DT combiner
TL=([ 0 0 1;
      1 0 0;
      0.3 0.2 1;
      0.7 0.8 0], # output for 3 samples, 2 classes and 2 ensemble members
    [0,0,1])	  # labels
L =j4pr.mat([0,1,0.2,0.8])
dtcomb=j4pr.dtcombiner(TL,2,2)
Test.@test all(vec(L |> dtcomb) .≈ [1-1/4*Distances.euclidean(dtcomb.x.data.DT[:,:,1],[0 0.2;1 0.8]),
					1-1/4*Distances.euclidean(dtcomb.x.data.DT[:,:,2],[0 0.2;1 0.8])])

end


text:
include("Transform/dct.jl")
"""
Some runtime parameters. The default value is after the `=`

## step 1:
```julia
bm3d_config.thr_patchSize = CartesianIndex(8, 8)	size of patches;
bm3d_config.thr_searchStride = CartesianIndex(3, 3)	In order to speed up the processing,the loop over the pixels of the imageis done with a stepp(integer) in row and column;
bm3d_config.thr_searchWin = CartesianIndex(19, 19)	search window size;
bm3d_config.thr_nMatch = 31	maximum number of similar patches kept;
bm3d_config.thr_thresh3D = 2.7	3D group hreshold;
bm3d_config.thr_distFunction = Euclidean2	A function that measures the distance between blocks;
bm3d_config.thr_group_transform! = dct!	A function that 3D transform in group
bm3d_config.thr_group_itransform! = idct!	A function that 3D inverse transform in group;
```

## step2:
```julia
bm3d_config.wie_patchSize = CartesianIndex(8, 8)	size of patches;
bm3d_config.wie_searchStride = CartesianIndex(3, 3)	In order to speed up the processing,the loop over the pixels of the imageis done with a stepp(integer) in row and column;
bm3d_config.wie_searchWin = CartesianIndex(11, 11)	search window size;
bm3d_config.wie_nMatch = 15	maximum number of similar patches kept;
bm3d_config.wie_distFunction = Euclidean2	A function that measures the distance between blocks;
bm3d_config.wie_group_transform! = dct!	A function that 3D transform in group;
bm3d_config.wie_group_itransform! = idct!	A function that 3D inverse transform in group;
```

## others
bm3d_config.show_info = false	show some running message

note:

The 3D group data struct:
```
G3D[k, :, :]	# k-th matching block.
size(G3D)	# (nMatch, patchSize[1], patchSize[2])
```
"""
mutable struct bm3d_config
	show_info::Bool

	thr_patchSize::CartesianIndex{2}
	thr_searchStride::CartesianIndex{2}
	thr_searchWin::CartesianIndex{2}
	thr_nMatch::Int
	thr_thresh3D::Float64
	thr_distFunction::Function
	thr_group_transform!::Function
	thr_group_itransform!::Function

	wie_patchSize::CartesianIndex{2}
	wie_searchStride::CartesianIndex{2}
	wie_searchWin::CartesianIndex{2}
	wie_nMatch::Int
	wie_distFunction::Function
	wie_group_transform!::Function
	wie_group_itransform!::Function

	bm3d_config() = new(
		false,
		CartesianIndex(8, 8),
		CartesianIndex(2, 2),
		CartesianIndex(24, 24),
		31,
		2.6,
		Euclidean2,
		dct!,
		idct!,

		CartesianIndex(8, 8),
		CartesianIndex(2, 2),
		CartesianIndex(19, 19),
		15,
		Euclidean2,
		dct!,
		idct!,
	)
end

"""
Euclidean2(a, b)

(Euclidean distance)^2
"""
function Euclidean2(a::AbstractArray{T, N}, b::AbstractArray{T, N}) where {T, N}
	norm(a .- b)^2
end

function print(config::bm3d_config)
	@info "========================================"
	@info "|       bm3d config                    |"
	@info "========================================"
	@info "| thr_patchSize: $(config.thr_patchSize)"
	@info "| thr_searchStride: $(config.thr_searchStride)"
	@info "| thr_searchWin: $(config.thr_searchWin)"
	@info "| thr_nMatch: $(config.thr_nMatch)"
	@info "| thr_thresh3D: $(config.thr_thresh3D)"
	@info "| thr_distFunction: $(config.thr_distFunction)"
	@info "| thr_group_transform!: $(config.thr_group_transform!)"
	@info "| thr_group_itransform!: $(config.thr_group_itransform!)"
	@info "| wie_patchSize: $(config.wie_patchSize)"
	@info "| wie_searchStride: $(config.wie_searchStride)"
	@info "| wie_searchWin: $(config.wie_searchWin)"
	@info "| wie_nMatch: $(config.wie_nMatch)"
	@info "| wie_distFunction: $(config.wie_distFunction)"
	@info "| wie_group_transform!: $(config.wie_group_transform!)"
	@info "| wie_group_itransform!: $(config.wie_group_itransform!)"
	@info "========================================"
end


text:
using SIIPExamples

using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
using D3TypeTrees

using Dates
using DataFrames

using Cbc #solver

rts_dir = SIIPExamples.download("https://github.com/GridMod/RTS-GMLC")
rts_src_dir = joinpath(rts_dir, "RTS_Data", "SourceData")
rts_siip_dir = joinpath(rts_dir, "RTS_Data", "FormattedData", "SIIP");

rawsys = PowerSystems.PowerSystemTableData(
    rts_src_dir,
    100.0,
    joinpath(rts_siip_dir, "user_descriptors.yaml"),
    timeseries_metadata_file = joinpath(rts_siip_dir, "timeseries_pointers.json"),
    generator_mapping_file = joinpath(rts_siip_dir, "generator_mapping.yaml"),
);
sys = System(rawsys; time_series_resolution = Dates.Hour(1));

branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticTransformer),
    :TT => DeviceModel(TapTransformer, StaticTransformer),
)

devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroROR => DeviceModel(HydroDispatch, FixedOutput),
    :Hydro => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
    :RenFx => DeviceModel(RenewableFix, FixedOutput),
    :ILoads => DeviceModel(InterruptibleLoad, InterruptiblePowerLoad),
)

services = Dict(
    :ReserveUp => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
    :ReserveDown => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
)

template_uc = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services);

solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)

horizon = 24;
interval = Dates.Hour(24);
transform_single_time_series!(sys, horizon, interval)

op_problem = OperationsProblem(
    GenericOpProblem,
    template_uc,
    sys;
    optimizer = solver,
    horizon = horizon,
)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl