text:
export Day4, Day4_1, Day4_2
abstract type Day4 <: Day end
abstract type Day4_1 <: Day4 end
abstract type Day4_2 <: Day4 end

solve(::Type{T}, lines::Vector) where T <: Day4 = begin
    (lines |> topassports .|> p -> validpassport(T, p)) |> sum
end

topassports(lines) = begin
    passports = Dict{String,String}[]
    parts = split_on_empty_line(lines)
    for part in parts
        p = (
                join(part, " ") |>
                split .|>
                l -> split(l, ":")
                ) |> Dict{String,String}
            push!(passports, p)
    end
    passports
end

hgt_test(hgt) = begin
    m=match(r"^(\d+)(.+)", hgt)
    if m == Nothing
        return false
    end
    n = parse(Int, m[1])
    t = m[2]
    if "cm" == t
        low, high = 150, 193
    elseif "in" == t
        low, high = 59, 76
    else
        return false
    end
    low <= n <= high
end

rules(::Type{Day4_1}, passport) = begin
    keyslength = keys(passport) |> length
    Dict([
        "length"=> (keyslength == 8) | ((keyslength == 7) & ("cid" ∉ keys(passport))),
    ])
end

rules(::Type{Day4_2}, passport) = begin
    r1 = rules(Day4_1, passport)
    r2 = Dict([
        "byr"=>1920 <= parse(Int, get(passport, "byr", "0")) <= 2002,
        "iyr"=>2010 <= parse(Int, get(passport, "iyr", "0")) <= 2020,
        "eyr"=>2020 <= parse(Int, get(passport, "eyr", "0")) <= 2030,
        "hcl"=>occursin(r"^#[a-f0-9]{6}$", get(passport, "hcl", "")),
        "ecl"=>get(passport, "ecl", "") ∈ ("amb", "blu", "brn", "gry", "grn", "hzl", "oth"),
        "pid"=>occursin(r"^[0-9]{9}$", get(passport, "pid", "")),
        "hgt"=>hgt_test(get(passport, "hgt", "0cm"))
    ])
    merge(r1,r2)
end

validpassport(::Type{T}, passport::Dict) where T <: Day4 = begin
    r = rules(T, passport)
    (false in values(r)) ? false : true
end


text:
# */ --------------------------------------------------------------------------------------------- /* #
# */ --- help functions -------------------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function circ(h::Int64,H::UnitRange{Int64})::Int64

    (h==minimum(H)) ? idx = maximum(H) : idx = h-1

    return idx
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- optimization models --------------------------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function MP!(m::JuMP.Model,ms::ms_type,mp::mp_type,ps::ps_type,pp::pp_type,unc::u_type)::JuMP.Model

    # problem structure
    # min   f + ∑ᵢ πᵢ g(xᵢ,cᵢ)

    xr,cr=1:unc.nx,1:unc.nc
    # */ -- master variables --------------------------------------- /* #
    @variable(m, f, container=Array) # investment-only cost
    @variable(m, x0[ms.P,ms.I0] >= .0, container=Array) # newly installed capacity
    @variable(m, x[xr,ms.I], container=Array) # rhs paramters for operational subproblem (eg, accumulated capacity)
    @variable(m, β[ms.I] >= .0, container=Array) # operational cost of node i
    # */ -- subproblems variables ---------------------------------- /* #
    @variable(m, ϕ[ms.I,cr], container=Array) # operational cost dependent on uncertain cost cᵢ
    @variable(m, c0[ms.I], container=Array) # operational cost independent of uncertain costs c
    @variable(m, yG[ms.I,ps.G,ps.S,ps.H] >= 0, container=Array) # generation level of conventional generation
    @variable(m, yI[ms.I,ps.B,ps.S,ps.H] >= 0, container=Array) # charging power of storage generation
    @variable(m, yO[ms.I,ps.B,ps.S,ps.H] >= 0, container=Array) # discharging power of storage generation
    @variable(m, yL[ms.I,ps.B,ps.S,ps.H] >= 0, container=Array) # energy level of storage generation
    @variable(m, yS[ms.I,ps.S,ps.H] >= 0, container=Array) # load shedding
    @variable(m, x1[ms.I,xr], container=Array) # values of rhs parameters in the subproblem (eg, capacity)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, f + mp.κ*sum(mp.π[i]*β[i] for i in ms.I) ) # investment plus operational cost (10⁶£)
    # */ ----------------------------------------------------------- /* #

    # */ -- master constraints ------------------------------------- /* #
    @constraint(m, cm01, f >= exp10(-6)*sum(mp.π0[i0]*sum(mp.ci[p,i0]*x0[p,i0] for p in ms.P) for i0 in ms.I0) + exp10(-6)*mp.κ*sum(mp.π[i]*sum(mp.cf[p]*x[p,i] for p in ms.P) for i in ms.I) ) # compute investment-only cost
    @constraint(m, cm02[p=ms.P,i=ms.I], x[p,i] == mp.xh[p,i] + sum(x0[p,i0] for i0 in mp.map[i]) ) # compute accumulated capacity at node i
    @constraint(m, cm03[p=ms.P,i=ms.I], x[p,i] >= .0 ) # minimum accumulated capacity
    @constraint(m, cm04[p=ms.P,i=ms.I], x[p,i] <= mp.xm[p] ) # maximum accumulated capacity
    for i in ms.I for j in 1:unc.nh
        fix(x[unc.nx0+j,i],unc.h[i,j];force=true) # set rhs paramters to dependent on subproblem i (eg, CO₂ budget parameter)
    end end
    # */ -- subproblems constraints -------------------------------- /* #
    @constraint(m, cs00[i=ms.I], β[i] >= c0[i] + unc.c[i,1]*ϕ[i,1] + unc.c[i,2]*ϕ[i,2]) # compute generation cost (10⁶£) of subproblem i
    @constraint(m, cs01[i=ms.I], c0[i] == exp10(-6)*sum(sum(sum((pp.cvg[g]+pp.cfg[g]/pp.ηg[g])*yG[i,g,s,h] for g in ps.G0) + pp.cvg[ps.gn]*yG[i,ps.gn,s,h] + pp.cs*yS[i,s,h] for h in ps.H)*pp.α for s in ps.S)) # compute operational cost independent of uncertain costs c
    @constraint(m, cs02[i=ms.I], ϕ[i,1] == exp10(-6)*sum(sum(sum(pp.eg[g]/pp.ηg[g]*yG[i,g,s,h] for g in ps.G) for h in ps.H)*pp.α for s in ps.S)) # compute cost dependent on CO₂ emission cost
    @constraint(m, cs03[i=ms.I], ϕ[i,2] == exp10(-6)*sum(sum(1/pp.ηg[ps.gn]*yG[i,ps.gn,s,h] for h in ps.H)*pp.α for s in ps.S)) # compute cost dependent on nuclear fuel cost
    @constraint(m, cs04[i=ms.I,b=ps.B,s=ps.S,h=ps.H], yL[i,b,s,h] - yL[i,b,s,circ(h,ps.H)] == pp.ηb[b]*yI[i,b,s,h] - yO[i,b,s,h]) # storage energy balance
    @constraint(m, cs05[i=ms.I,g=ps.G,s=ps.S,h=ps.H], yG[i,g,s,h] - yG[i,g,s,circ(h,ps.H)] <= pp.rg[g]*x1[i,g]) # upward ramping limitation on conventional generation
    @constraint(m, cs06[i=ms.I,g=ps.G,s=ps.S,h=ps.H], yG[i,g,s,circ(h,ps.H)] - yG[i,g,s,h] <= pp.rg[g]*x1[i,g]) # downward ramping limitation on conventional generation
    @constraint(m, cs07[i=ms.I,g=ps.G,s=ps.S,h=ps.H], yG[i,g,s,h] <= x1[i,g]) # maximum capacity limitation on conventional generation
    @constraint(m, cs08[i=ms.I,b=ps.B,s=ps.S,h=ps.H], yI[i,b,s,h] <= pp.pb[b]*x1[i,b+ps.b0]) # maximum capacity limitation on storage charging power
    @constraint(m, cs09[i=ms.I,b=ps.B,s=ps.S,h=ps.H], yO[i,b,s,h] <= pp.pb[b]*x1[i,b+ps.b0]) # maximum capacity limitation on storage discharging power
    @constraint(m, cs10[i=ms.I,b=ps.B,s=ps.S,h=ps.H], yL[i,b,s,h] <= x1[i,b+ps.b0]) # maximum capacity limitation on storage energy level
    @constraint(m, cs11[i=ms.I], sum(yG[i,g,s,h]*pp.eg[g]/pp.ηg[g] for g in ps.G for s in ps.S for h in ps.H) <= pp.lco*x1[i,unc.nx0+1]) # CO₂ emission budget
    @constraint(m, cs12[i=ms.I,s=ps.S,h=ps.H], sum(yG[i,g,s,h] for g in ps.G) + sum(yO[i,b,s,h]-yI[i,b,s,h] for b in ps.B) + yS[i,s,h] >= -x1[i,unc.nx0+2]*pp.pd[s,h] - sum(pp.pr[r,s,h]*x1[i,r+ps.r0] for r in ps.R)) # energy demand
    @constraint(m, cs13[i=ms.I,n=1:unc.nx0], x1[i,n]*exp10(-3) == x[n,i]) # generation to level of decisions {xᵢ, ∀i} set by the master problem
    @constraint(m, cs14[i=ms.I,n=1:unc.nh ], x1[i,unc.nx0+n]   == x[unc.nx0+n,i]) # rhs uncertain parameters (energy level and CO₂ budget) specific to problem i
    # */ ----------------------------------------------------------- /* #

    return m
end

function RMP!(m::JuMP.Model,ms::ms_type,mp::mp_type,unc::u_type)::JuMP.Model

    # problem structure
    # min   f + ∑ᵢ πᵢ βᵢ)
    # s.t.  f = q₀ᵀx₀ + ∑ᵢ qᵢᵀxᵢ
    #      (x₀,x₁,..,xᵢ,..,xₙ) ∈ 𝕏
    #       βᵢ ∈ Θᵢ,  ∀i

    xr=1:unc.nx
    # */ -- variables ---------------------------------------------- /* #
    @variable(m, f, container=Array) # investment-only cost
    @variable(m, x0[ms.P,ms.I0] >= .0, container=Array) # newly installed capacity
    @variable(m, x[xr,ms.I], container=Array) # rhs paramters for operational subproblem (eg, accumulated capacity)
    @variable(m, β[ms.I] >= .0, container=Array) # operational cost of node i
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, f + mp.κ*sum(mp.π[i]*β[i] for i in ms.I) ) # investment plus operational cost (10⁶£)
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    @constraint(m, c01, f >= exp10(-6)*sum(mp.π0[i0]*sum(mp.ci[p,i0]*x0[p,i0] for p in ms.P) for i0 in ms.I0) + exp10(-6)*mp.κ*sum(mp.π[i]*sum(mp.cf[p]*x[p,i] for p in ms.P) for i in ms.I) ) # compute investment-only cost
    @constraint(m, c02[p=ms.P,i=ms.I], x[p,i] == mp.xh[p,i] + sum(x0[p,i0] for i0 in mp.map[i]) ) # compute accumulated capacity at node i
    @constraint(m, c03[p=ms.P,i=ms.I], x[p,i] >= .0 ) # minimum accumulated capacity
    @constraint(m, c04[p=ms.P,i=ms.I], x[p,i] <= mp.xm[p] ) # maximum accumulated capacity
    for i in ms.I for j in 1:unc.nh
        fix(x[unc.nx0+j,i],unc.h[i,j];force=true) # set rhs paramters to dependent on subproblem i (eg, CO₂ budget parameter)
    end end
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

function SP!(m::JuMP.Model,ps::ps_type,pp::pp_type,unc::u_type)::JuMP.Model

    # problem structure
    # min   c₀ + cᵀϕ
    # s.t.  c₀  = C₀y
    #       ϕ   = C y
    #       y ∈ 𝕐
    #       A y ≦ B x
    #       x = xf :(λ)

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ[1:unc.nc], container=Array) # operational cost dependent on uncertain cost cᵢ
    @variable(m, c0, container=Array) # operational cost independent of uncertain costs c
    @variable(m, yG[ps.G,ps.S,ps.H] >= 0, container=Array) # generation level of conventional generation
    @variable(m, yI[ps.B,ps.S,ps.H] >= 0, container=Array) # charging power of storage generation
    @variable(m, yO[ps.B,ps.S,ps.H] >= 0, container=Array) # discharging power of storage generation
    @variable(m, yL[ps.B,ps.S,ps.H] >= 0, container=Array) # energy level of storage generation
    @variable(m, yS[ps.S,ps.H] >= 0, container=Array) # load shedding
    @variable(m, x0[1:unc.nx], container=Array) # values of rhs parameters in the subproblem (eg, capacity)
    @variable(m, x[1:unc.nx], container=Array) # values of decisions {xᵢ, ∀i} set by the master problem
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, c0 ) # generation cost (10⁶£)
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    @constraint(m, c01, c0 == exp10(-6)*sum(sum(sum((pp.cvg[g]+pp.cfg[g]/pp.ηg[g])*yG[g,s,h] for g in ps.G0) + pp.cvg[ps.gn]*yG[ps.gn,s,h] + pp.cs*yS[s,h] for h in ps.H)*pp.α for s in ps.S)) # compute operational cost independent of uncertain costs c
    @constraint(m, c02, ϕ[1] == exp10(-6)*sum(sum(sum(pp.eg[g]/pp.ηg[g]*yG[g,s,h] for g in ps.G) for h in ps.H)*pp.α for s in ps.S)) # compute cost dependent on CO₂ emission cost
    @constraint(m, c03, ϕ[2] == exp10(-6)*sum(sum(1/pp.ηg[ps.gn]*yG[ps.gn,s,h] for h in ps.H)*pp.α for s in ps.S)) # compute cost dependent on nuclear fuel cost
    @constraint(m, c04[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] - yL[b,s,circ(h,ps.H)] == pp.ηb[b]*yI[b,s,h] - yO[b,s,h]) # storage energy balance
    @constraint(m, c05[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] - yG[g,s,circ(h,ps.H)] <= pp.rg[g]*x0[g]) # upward ramping limitation on conventional generation
    @constraint(m, c06[g=ps.G,s=ps.S,h=ps.H], yG[g,s,circ(h,ps.H)] - yG[g,s,h] <= pp.rg[g]*x0[g]) # downward ramping limitation on conventional generation
    @constraint(m, c07[g=ps.G,s=ps.S,h=ps.H], yG[g,s,h] <= x0[g]) # maximum capacity limitation on conventional generation
    @constraint(m, c08[b=ps.B,s=ps.S,h=ps.H], yI[b,s,h] <= pp.pb[b]*x0[b+ps.b0]) # maximum capacity limitation on storage charging power
    @constraint(m, c09[b=ps.B,s=ps.S,h=ps.H], yO[b,s,h] <= pp.pb[b]*x0[b+ps.b0]) # maximum capacity limitation on storage discharging power
    @constraint(m, c10[b=ps.B,s=ps.S,h=ps.H], yL[b,s,h] <= x0[b+ps.b0]) # maximum capacity limitation on storage energy level
    @constraint(m, c11, sum(yG[g,s,h]*pp.eg[g]/pp.ηg[g] for g in ps.G for s in ps.S for h in ps.H) <= pp.lco*x0[unc.nx0+1]) # CO₂ emission budget
    @constraint(m, c12[s=ps.S,h=ps.H], sum(yG[g,s,h] for g in ps.G) + sum(yO[b,s,h]-yI[b,s,h] for b in ps.B) + yS[s,h] >= -x0[unc.nx0+2]*pp.pd[s,h] - sum(pp.pr[r,s,h]*x0[r+ps.r0] for r in ps.R)) # energy demand
    @constraint(m, c13[n=1:unc.nx0], x0[n]*exp10(-3) == x[n]) # generation to level of decisions {xᵢ, ∀i} set by the master problem
    @constraint(m, c14[n=1:unc.nh ], x0[unc.nx0+n]   == x[unc.nx0+n]) # rhs uncertain parameters (energy level and CO₂ budget) specific to problem i
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

function LBoracle!(m::JuMP.Model,unc::u_type,J::Int64)::JuMP.Model

    # problem structure
    # min   ϕ + γᵀcᵢ
    # s.t.  ϕ + γᵀcₛ ≧ θₛ + λᵀₛ(x-xₛ),   ∀s
    #       x = xᵢ

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ, container=Array)
    @variable(m, γ[1:unc.nc] >= .0, container=Array)
    @variable(m, x[1:unc.nx], container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Min, 0. )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

function UBoracle!(m::JuMP.Model,unc::u_type,J::Int64)::JuMP.Model

    # problem structure
    # max   ϕ + γᵀxᵢ
    # s.t.  ϕ + γᵀxₛ ≦ θₛ + ϕᵀₛ(c-cₛ),   ∀s
    #       c = cᵢ

    # */ -- variables ---------------------------------------------- /* #
    @variable(m, ϕ, container=Array)
    @variable(m, γ[1:unc.nx] <= .0, container=Array)
    @variable(m, c[1:unc.nc], container=Array)
    # */ ----------------------------------------------------------- /* #

    # */ -- obj function ------------------------------------------- /* #
    @objective(m, Max, 0. )
    # */ ----------------------------------------------------------- /* #

    # */ -- constraints -------------------------------------------- /* #
    # */ ----------------------------------------------------------- /* #
    optimize!(m)

    return m
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : master -------------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function gen_kiidx(u::u_type,w::Int64)::Tuple{Array{BitArray{1},1},Array{Array{Int64,1},1}}

    # function to generate (kidx,iidx) for given number of subproblems w solved at each iteration
    # iidx -> array of subsets {ℰₖ, 𝑘=1,..,w}
    # kidx -> array of 0/1 association with subset ℰₖ

    if w == 1 # (trivial case with one subproblem) ℰₖ := ℰ
        kidx = Array{BitArray{1},1}(undef,0)
        push!(kidx,zeros(Int64,u.ni).==0)
        iidx = Array{Array{Int64,1},1}(undef,0)
        push!(iidx,[i for i in 1:u.ni])
    else
        x = zeros(u.nh+u.nc,u.ni) # generate empty matrix x for (H,C) values of each subproblem
        for j in 1:u.nh
            uh = u.h[:,j] # collect values of H (rhs uncertain params)
            hmin,hmax = minimum(uh),maximum(uh) # evaluate min & max values
            (hmax > hmin) ? x[j,:] .= (uh.-hmin)./(hmax-hmin) :  x[j,:] .= zeros(u.ni) # set normalized H values into x
        end
        for j in 1:u.nc
            uc = u.c[:,j] # collect values of C (uncertain cost coeff)
            cmin,cmax = minimum(uc),maximum(uc) # evaluate min & max values
            (cmax > cmin) ? x[j+u.nh,:] .= (uc.-cmin)./(cmax-cmin) : x[j+u.nh,:] .= zeros(u.ni) # set normalized C values into x
        end
        kass = kmeans(x,w).assignments # run kmeans clustering on x to obtain w subsets
        kidx = Array{BitArray{1},1}(undef,0) # generate empty kidx
        for j in 1:w
            push!(kidx,kass.==j) # for each subset j = 1,..,w generate a bitarry with 1 if subproblem i belongs to subset j (and 0 otherwise)
        end
        idx = [i for i in 1:u.ni] # generate array of indices of subproblems
        iidx = Array{Array{Int64,1},1}(undef,0) # generate empty iidx
        for j in 1:w
            push!(iidx,idx[kidx[j]]) # for each subset j = 1,..,w generate an array with all the subproblems that belongs to subset j
        end
    end

    return kidx,iidx
end

function gen_D__type(ms::ms_type,mp::mp_type,ps::ps_type,pp::pp_type,u::u_type,c::Int64,a::Int64,q::Float64,w::Int64,j::Int64,e::Float64)::D__type

    # function to generate D__type (structure of problem data) for given ms (master problem sets), mp (master problem params),
    # ps (subproblem sets), pp (subproblem params), u (uncertainty params), c (case), a (algorithm), q (tolerance damping),
    # w (number of exact subproblems), j (maximum number of iteration), and e (convergence tolerance)

    k,i = gen_kiidx(u,w) # generate (kidx,iidx) for given uncertain data u and given number of subproblems w solved at each iteration

    return D__type(ms,mp,ps,pp,u,c,a,q,w,k,i,j,e)
end

function gen_R__type(d::D__type)::R__type

    # function to generate R__type (structure of relax master problem) given data D__type (problem data)

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env))) # generate empty JuMP optimization model
    JuMP.set_optimizer_attribute(m,"OutputFlag",0) # set JuMP model silent
    RMP!(m,d.ms,d.mp,d.unc) # build relaxed master problem
    v = Array{JuMP.VariableRef,2}(undef,d.unc.nx+1,d.unc.ni) # generate empty array of vars
    for i in 1:d.unc.ni v[:,i] .= vcat([m[:β][i]],m[:x][:,i]) end # set reference to variablems {(βᵢ,xᵢ), ∀i}
    c = Array{JuMP.ConstraintRef,2}(undef,d.J,d.unc.ni) # generate empty array of cons

    return R__type(m,v,c)
end

function gen_T1_type(d::D__type)::T1_type

    # function to generate T1_type (structure of temporary data, algorithm 1 ) given data D__type (problem data)

    x = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    β = zeros(d.unc.ni) # generate empty array for variables value {βᵢ, ∀i}
    θ = zeros(d.unc.ni) # generate empty array for optimal objectives {θᵢ, ∀i}
    λ = zeros(d.unc.nx,d.unc.ni) # generate empty array for subgradient {λᵢ, ∀i} wrt xᵢ

    return T1_type(x,c,β,θ,λ,0.)
end

function gen_T3_type(d::D__type)::T3_type

    # function to generate T1_type (structure of temporary data, algorithm 3) given data D__type (problem data)

    x  = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c  = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    β  = zeros(d.unc.ni) # generate empty array for variables value {βᵢ, ∀i}
    θp = zeros(d.unc.ni) # generate empty array for valid upper bounds on objectives {θᵢ, ∀i}
    θd = zeros(d.unc.ni) # generate empty array for valid lower bounds on objectives {θᵢ, ∀i}
    λ  = zeros(d.unc.nx,d.unc.ni) # generate empty array for valid subgradient {λᵢ, ∀i} wrt xᵢ

    return T3_type(x,c,β,θp,θd,λ,0.)
end

function gen_T2_type(d::D__type)::T2_type

    # function to generate T1_type (structure of temporary data, algorithm 2) given data D__type (problem data)

    x  = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c  = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    β  = zeros(d.unc.ni) # generate empty array for variables value {βᵢ, ∀i}
    θl = zeros(d.unc.ni) # generate empty array for valid lower bounds on objectives {θᵢ, ∀i}
    λl = zeros(d.unc.nx,d.unc.ni) # generate empty array for valid subgradient {λᵢ, ∀i} wrt xᵢ
    θu = zeros(d.unc.ni) # generate empty array for valid upper bounds on objectives {θᵢ, ∀i}
    ϕu = zeros(d.unc.nc,d.unc.ni) # generate empty array for valid subgradient {ϕᵢ, ∀i} wrt cᵢ
    Ie = zeros(Int64,d.w) # generate empty array for indices of w subproblems

    return T2_type(x,c,β,θl,λl,θu,ϕu,Ie,0.)
end

function gen_H__type(d::D__type)::H__type

    # function to generate H_type (structure of stored data) given data D__type (problem data)

    L = zeros(d.J) # generate empty array for lower bound values
    U = zeros(d.J) # generate empty array for upper bound values
    t = zeros(d.J) # generate empty array for time spend to perform each iteration
    T = zeros(d.J,6) # generate empty array for time spend to perform each iteration, split by steps

    return H__type(L,U,t,T,0)
end

function gen_M__type(d::D__type)::M__type

    # function to generate M_type (structure of exact solutions stored) given data D__type (problem data)

    θ = zeros(nprocs()) # generate empty array of optimal objective
    λ = zeros(nprocs(),d.unc.nx) # generate empty array of subgradient wrt x
    ϕ = zeros(nprocs(),d.unc.nc) # generate empty array of subgradient wrt c
    x = zeros(nprocs(),d.unc.nx) # generate empty array of values of x
    c = zeros(nprocs(),d.unc.nc) # generate empty array of values of c

    return M__type(0,θ,λ,ϕ,x,c)
end

function gen_B0_type(c::Int64)::JuMP.Model

    # function to generate B0_type (Benders master structure of algorithm 0) given c (case)

    ms,mp,ps,pp,unc = load_data(c) # generate data structures ms,mp,ps,pp,unc for given value of c (case)
    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env))) # generate empty JuMP model
    println(" Attempting to build JuMP model of deterministic equivalent")
    println(" ")
    MP!(m,ms,mp,ps,pp,unc) # build master problem
    println(" JuMP Model build successfully!")
    println(" ")

    return m
end

function gen_B1_type(c::Int64,a::Int64,j::Int64,e::Float64)::B1_type

    # function to generate B1_type (Benders master structure of algorithm 1 ) given c (case), a (algorithm),
    # j (maximum number of iteration), and e (convergence tolerance).

    ms,mp,ps,pp,unc = load_data(c) # generate data structures ms,mp,ps,pp,unc for given value of c (case)
    d = gen_D__type(ms,mp,ps,pp,unc,c,a,.0,1,j,e) # generate D__type (structure of problem data)
    r = gen_R__type(d) # generate R__type (structure of relax master problem)
    t = gen_T1_type(d) # generate T1_type (structure of temporary data)
    h = gen_H__type(d) # generate H__type (structure of stored data)

    return B1_type(r,t,h,d)
end

function gen_B3_type(c::Int64,a::Int64,q::Float64,j::Int64,e::Float64)::B3_type

    # function to generate B3_type (Benders master structure of algorithm 3) given c (case), a (algorithm), q (tolerance damping),
    # j (maximum number of iteration), and e (convergence tolerance).

    ms,mp,ps,pp,unc = load_data(c) # generate data structures ms,mp,ps,pp,unc for given value of c (case)
    d = gen_D__type(ms,mp,ps,pp,unc,c,a,q,1,j,e) # generate D__type (structure of problem data)
    r = gen_R__type(d) # generate R__type (structure of relax master problem)
    t = gen_T3_type(d) # generate T3_type (structure of temporary data)
    h = gen_H__type(d) # generate H__type (structure of stored data)

    return B3_type(r,t,h,d)
end

function gen_B2_type(c::Int64,a::Int64,w::Int64,j::Int64,e::Float64)::B2_type

    # function to generate B2_type (Benders master structure of algorithm 2) given c (case), a (algorithm), w (number of exact subproblems),
    # j (maximum number of iteration), and e (convergence tolerance).

    ms,mp,ps,pp,unc = load_data(c) # generate data structures ms,mp,ps,pp,unc for given value of c (case)
    d = gen_D__type(ms,mp,ps,pp,unc,c,a,.0,w,j,e) # generate D__type (structure of problem data)
    r = gen_R__type(d) # generate R__type (structure of relax master problem)
    t = gen_T2_type(d) # generate T2_type (structure of temporary data)
    h = gen_H__type(d) # generate H__type (structure of stored data)
    m = gen_M__type(d) # generate M__type (structure of exact solutions stored)

    return B2_type(r,t,h,m,d)
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : subproblem ---------------------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function gen_d__type(d::D__type)::d__type

    # function to generate d__type (structure of subproblem data) for given D__type (structure of problem data)

    s = d.ps # set subproblem sets
    p = d.pp # set subproblem parameters
    u = d.unc # set uncertainty parameters
    a = d.algm # set algorithm
    j = Int64(round(d.J*d.w;digits=0)) # set maximum number of exact solutions to store

    return d__type(s,p,u,a,j)
end

function gen_Lb_type(d::d__type)::O__type

    # function to generate O__type (structure of oracle for lower bound) for given d__type (structure of subproblem data)

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env))) # generate empty JuMP model
    JuMP.set_optimizer_attribute(m,"OutputFlag",0) # set model silent
    LBoracle!(m,d.unc,d.J+1) # build lower bound oracle
    v = vcat(m[:ϕ],m[:γ]) # set reference to variables (ϕ,γ)
    c = Array{JuMP.ConstraintRef,1}(undef,d.J+1) # generate empty reference to new constraints
    o = v'*ones(length(v)) # generate reference for objective function
    e = v'*ones(length(v)) # generate reference for objective function
    h = zeros(d.J+1) # generate empty array for rhs values of old constraints (old exact solutions)

    return O__type(m,c,v,o,e,h,0)
end

function gen_Ub_type(d::d__type)::O__type

    # function to generate O__type (structure of oracle for upper bound) for given d__type (structure of subproblem data)

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env))) # generate empty JuMP model
    JuMP.set_optimizer_attribute(m,"OutputFlag",0) # set model silent
    UBoracle!(m,d.unc,d.J+1) # build upper bound oracle
    v = vcat(m[:ϕ],m[:γ]) # set reference to variables (ϕ,γ)
    c = Array{JuMP.ConstraintRef,1}(undef,d.J+1) # generate empty reference to new constraints
    o = v'*ones(length(v)) # generate reference for objective function
    e = v'*ones(length(v)) # generate reference for objective function
    h = zeros(d.J+1) # generate empty array for rhs values of old constraints (old exact solutions)

    return O__type(m,c,v,o,e,h,0)
end

function gen_E__type(d::d__type)::E__type

    # function to generate E__type (structure of exact subproblem) for given d__type (structure of subproblem data)

    m = JuMP.Model(JuMP.optimizer_with_attributes(()->Gurobi.Optimizer(gurobi_env))) # generate empty JuMP model
    JuMP.set_optimizer_attribute(m,"OutputFlag",0) # set model silent
    JuMP.set_optimizer_attribute(m,"Method",2) # set solution method (barrier)
    d.algm == 3 ? JuMP.set_optimizer_attribute(m,"Crossover",0) : nothing # deactivate crossover (if algorithm 3)
    SP!(m,d.ps,d.pp,d.unc) # build subproblem
    v = vcat(m[:c0],m[:ϕ][:]) # set reference to variables (c₀,ϕ)
    o = v'*ones(length(v)) # generate reference for objective function

    return E__type(m,v,o)
end

function gen_M__type(d::d__type)::M__type

    # function to generate M__type (structure of exact solutions stored) for given d__type (structure of subproblem data)

    θ = zeros(d.J+1) # generate empty array for optimal objective of exact solutions
    λ = zeros(d.J+1,d.unc.nx) # generate empty array for subgradient wrt x of exact solutions
    ϕ = zeros(d.J+1,d.unc.nc) # generate empty array for subgradient wrt c of exact solutions
    x = zeros(d.J+1,d.unc.nx) # generate empty array for values of x of exact solutions
    c = zeros(d.J+1,d.unc.nc) # generate empty array for values of c of exact solutions

    return M__type(0,θ,λ,ϕ,x,c)
end

function gen_Q__type(d::d__type)::Array{Q__type,1}

    # function to generate Q__type (structure of extended exact solutions stored) for given d__type (structure of subproblem data)

    return Array{Q__type,1}(undef,0)
end

function gen_t1_type(d::d__type)::t1_type

    # function to generate t1_type (structure of temporary data, algorithm 1 ) for given d__type (structure of subproblem data)

    x = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    θ = zeros(d.unc.ni) # generate empty array for optimal objectives {θᵢ, ∀i}
    λ = zeros(d.unc.nx,d.unc.ni) # generate empty array for subgradient {λᵢ, ∀i} wrt xᵢ

   return t1_type(x,c,θ,λ,0)
end

function gen_t3_type(d::d__type)::t3_type

    # function to generate t3_type (structure of temporary data, algorithm 3) for given d__type (structure of subproblem data)

    x  = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c  = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    θp = zeros(d.unc.ni) # generate empty array for valid upper bounds on objectives {θᵢ, ∀i}
    θd = zeros(d.unc.ni) # generate empty array for valid power bounds on objectives {θᵢ, ∀i}
    λ  = zeros(d.unc.nx,d.unc.ni) # generate empty array for valid subgradients {λᵢ, ∀i} wrt xᵢ

   return t3_type(x,c,θp,θd,λ,0)
end

function gen_t2_type(d::d__type)::t2_type

    # function to generate t2_type (structure of temporary data, algorithm 2) for given d__type (structure of subproblem data)

    x  = zeros(d.unc.nx,d.unc.ni) # generate empty array for decisions {xᵢ, ∀i}
    c  = convert(Array{Float64,2},d.unc.c') # set values of parameters {cᵢ, ∀i}
    θl = zeros(d.unc.ni) # generate empty array for valid power bounds on objectives {θᵢ, ∀i}
    θu = zeros(d.unc.ni) # generate empty array for valid upper bounds on objectives {θᵢ, ∀i}
    λl = zeros(d.unc.nx,d.unc.ni) # generate empty array for valid subgradients {λᵢ, ∀i} wrt xᵢ
    ϕu = zeros(d.unc.nc,d.unc.ni) # generate empty array for valid subgradients {ϕᵢ, ∀i} wrt cᵢ

   return t2_type(x,c,θl,θu,λl,ϕu,0)
end

function gen_S1_type(D::D__type)::S1_type

    # function to generate S1_type (structure of Benders subproblem, algorithm 1 ) for given d__type (structure of subproblem data)

    d = gen_d__type(D) # generate d__type (structure of subproblem data)
    e = gen_E__type(d) # generate E__type (structure of exact subproblem)
    t = gen_t1_type(d) # generate t1_type (structure of temporary data)

    return S1_type(e,d,t)
end

function gen_S3_type(D::D__type)::S3_type

    # function to generate S3_type (structure of Benders subproblem, algorithm 3) for given d__type (structure of subproblem data)

    d = gen_d__type(D) # generate d__type (structure of subproblem data)
    e = gen_E__type(d) # generate E__type (structure of exact subproblem)
    t = gen_t3_type(d) # generate t3_type (structure of temporary data)

    return S3_type(e,d,t)
end

function gen_S2_type(D::D__type)::S2_type

    # function to generate S2_type (structure of Benders subproblem, algorithm 2) for given d__type (structure of subproblem data)

    d = gen_d__type(D) # generate d__type (structure of subproblem data)
    e = gen_E__type(d) # generate E__type (structure of exact subproblem)
    l = gen_Lb_type(d) # generate O__type (structure of oracle for lower bound)
    u = gen_Ub_type(d) # generate O__type (structure of oracle for upper bound)
    m = gen_M__type(d) # generate M__type (structure of exact solutions stored)
    q = gen_Q__type(d) # generate Q__type (structure of extended exact solutions stored)
    t = gen_t2_type(d) # generate t2_type (structure of temporary data)

    return S2_type(e,l,u,m,q,d,t)
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- functions to generate structs : master and subproblem ----------------------------------- /* #
# */ --------------------------------------------------------------------------------------------- /* #

function generate_structures(c::Int64,a::Int64,q::Float64,w::Int64,j::Int64,e::Float64)

    # function to generate B__struct and S__struct for given values of c (case), a (algorithm), q (tolerance damping),
    # w (number of exact subproblems), j (maximum number of iteration), and e (convergence tolerance)

    if a == 0
        B = gen_B0_type(c) # generate B__struct (Benders master structure)
        S = nothing # empty subproblem structure
    end
    if a == 1 # (if Stand_Bend)
        B = gen_B1_type(c,a,j,e) # generate B__struct (Benders master structure)
        S = gen_S1_type(B.data) # generate S__struct (Benders subproblem structure)
    end
    if a == 2 # (if Adapt_Bend)
        B = gen_B2_type(c,a,w,j,e) # generate B__struct (Benders master structure)
        S = gen_S2_type(B.data) # generate S__struct (Benders subproblem structure)
    end
    if a == 3 # (if Zaker_Bend)
        B = gen_B3_type(c,a,q,j,e) # generate B__struct (Benders master structure)
        S = gen_S3_type(B.data) # generate S__struct (Benders subproblem structure)
    end

    return B,S
end

# */ --------------------------------------------------------------------------------------------- /* #
# */ --- lower level functions : subproblem ------------------------------------------------------ /* #
# */ --------------------------------------------------------------------------------------------- /* #

function set_Iex!(b::B2_type)::B2_type

    # function to update B2_type.temp.Ie (index of w subproblems solved exactly @ iter k)
    # for given B2_type (Benders master problem) -> algorithm 2 (Adapt_Bend)

    for k in 1:b.data.w # for each k = 1,..,w (for each subset ℰₖ)
        ie = findmax(b.data.mp.π[b.data.kidx[k]].*(b.temp.θu[b.data.kidx[k]].-b.temp.θl[b.data.kidx[k]]))[2] # select subproblem ie ∈ ℰₖ for which the gap is the largest
        b.temp.Ie[k] = b.data.iidx[k][ie] # store ie in B2_type.temp.Ie
    end

    return b
end

function solv_exact!(s::S1_type)::S1_type

    # function to update and solve S1_type.ex.m (exact subproblem) and store exact solution in S1_type.temp
    # for given S1_type (Benders subproblem) -> algorithm 1 (Stand_Bend)

    s.ex.objf =  s.ex.vars'*vcat([1.],s.temp.c[:,s.temp.i]) # compute new objective of subproblem i (c₀ + cᵢᵀϕ)
    @objective(s.ex.m, Min, s.ex.objf) # update objective in JuMP model (c₀ + cᵢᵀϕ)
    fix.(s.ex.m[:x],s.temp.x[:,s.temp.i];force=true) # fix values of variables xᵢ in the model to decisions xᵢ set by the master problem
    optimize!(s.ex.m) # solve the JuMP model to optimality
    s.temp.θ[s.temp.i] = objective_value(s.ex.m) # store optimal objective θᵢ
    s.temp.λ[:,s.temp.i] .= dual.(FixRef.(s.ex.m[:x])) # store subgradient λᵢ wrt xᵢ

    return s
end

function solv_exact!(s::S3_type)::S3_type

    # function to update and solve S3_type.ex.m (exact subproblem) and store exact solution in S3_type.temp
    # for given S3_type (Benders subproblem) -> algorithm 3 (Zaker_Bend)

    s.ex.objf =  s.ex.vars'*vcat([1.],s.temp.c[:,s.temp.i]) # compute new objective of subproblem i (c₀ + cᵢᵀϕ)
    @objective(s.ex.m, Min, s.ex.objf) # update objective in JuMP model (c₀ + cᵢᵀϕ)
    fix.(s.ex.m[:x],s.temp.x[:,s.temp.i];force=true) # fix values of variables xᵢ in the model to decisions xᵢ set by the master problem
    optimize!(s.ex.m) # solve the JuMP model up to optimality tolerance δ
    s.temp.θp[s.temp.i] = objective_value(s.ex.m) # save suboptimal primal objective θᵢ (upper bound)
    s.temp.λ[:,s.temp.i] .= dual.(FixRef.(s.ex.m[:x])) # store subgradient λᵢ wrt xᵢ
    s.temp.θd[s.temp.i] = s.temp.λ[:,s.temp.i]'*value.(s.ex.m[:x]) # save suboptimal dual objective θᵢ (lower bound)

    return s
end

function solv_exact!(s::S2_type)::S2_type

    # function to update and solve S2_type.ex.m (exact subproblem) and store exact solution in S2_type.temp
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.ex.objf =  s.ex.vars'*vcat([1.],s.temp.c[:,s.temp.i]) # compute new objective of subproblem i (c₀ + cᵢᵀϕ)
    @objective(s.ex.m, Min, s.ex.objf) # update objective in JuMP model (c₀ + cᵢᵀϕ)
    fix.(s.ex.m[:x],s.temp.x[:,s.temp.i];force=true) # fix values of variables xᵢ in the model to decisions xᵢ set by the master problem
    optimize!(s.ex.m) # solve the JuMP model to optimality
    s.temp.θl[s.temp.i] = objective_value(s.ex.m) # store optimal objective θᵢ
    s.temp.λl[:,s.temp.i] .= dual.(FixRef.(s.ex.m[:x])) # store subgradient λᵢ wrt xᵢ
    s.temp.ϕu[:,s.temp.i] .= value.(s.ex.m[:ϕ]) # store subgradient ϕᵢ wrt cᵢ

    return s
end

function extended_sol!(s::S2_type)::S2_type

    # function to store new extended exact solution to S2_type.q (extended exact solutions stored)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    ϕ  = value.(s.ex.m[:ϕ] ) # extract the operational cost dependent on uncertain cost cᵢ
    c0 = value.(s.ex.m[:c0]) # extract the operational cost independent of uncertain costs c
    yG = value.(s.ex.m[:yG]) # extract the generation level of conventional generation
    yI = value.(s.ex.m[:yI]) # extract the charging power of storage generation
    yO = value.(s.ex.m[:yI]) # extract the discharging power of storage generation
    yL = value.(s.ex.m[:yL]) # extract the energy level of storage generation
    yS = value.(s.ex.m[:yS]) # extract the load shedding
    x0 = value.(s.ex.m[:x0]) # extract the values of rhs parameters in the subproblem (eg, capacity)
    push!(s.q,Q__type(ϕ,c0,yG,yI,yO,yL,yS,x0)) # add new extended exact solution

    return s
end

function update_m_!(s::S2_type)::S2_type

    # function to add a new exact solution to S2_type.m (exact solutions stored)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.m.n += 1 # update number of solution stored (+1)
    s.m.θ[s.m.n] = s.temp.θl[s.temp.i] # store optimal solution θᵢ
    s.m.λ[s.m.n,:] .= s.temp.λl[:,s.temp.i] # store subgradient λᵢ wrt xᵢ
    s.m.ϕ[s.m.n,:] .= s.temp.ϕu[:,s.temp.i] # store subgradient ϕᵢ wrt cᵢ
    s.m.x[s.m.n,:] .= s.temp.x[:,s.temp.i] # store xᵢ
    s.m.c[s.m.n,:] .= s.temp.c[:,s.temp.i] # store cᵢ

    return s
end

function update_lb!(s::S2_type)::S2_type

    # function to add a new exact solution to S2_type.olb (oracle for lower bound)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.olb.cons[s.m.n] = @constraint(s.olb.m, s.olb.m[:ϕ] + s.olb.m[:γ]'*s.m.c[s.m.n,:] >= 0. ) # add constraint (ϕ + γᵀcₙ >= 0) to the oracle
    s.olb.help[1:s.m.n] .= s.m.θ[1:s.m.n] .- sum!(ones(s.m.n),s.m.λ[1:s.m.n,:].*s.m.x[1:s.m.n,:]) # compute vector h used when solving the lb oracle (hₛ = θₛ - λₛᵀxₛ, ∀s=1,..,n)
    s.olb.n += 1 # update number of solution stored (+1)

    return s
end

function update_ub!(s::S2_type)::S2_type

    # function to add a new exact solution in S2_type.olb (oracle for upper bound)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.oub.cons[s.m.n] = @constraint(s.oub.m, s.oub.m[:ϕ] + s.oub.m[:γ]'*s.m.x[s.m.n,:] <= 0. ) # add constraint (ϕ + γᵀxₙ <= 0) to the oracle
    s.oub.help[1:s.m.n] .= s.m.θ[1:s.m.n] .- sum!(ones(s.m.n),s.m.ϕ[1:s.m.n,:].*s.m.c[1:s.m.n,:]) # compute vector h used when solving the ub oracle (hₛ = θₛ - ϕₛᵀcₛ, ∀s)
    s.oub.n += 1 # update number of solution stored (+1)

    return s
end

function update_s!(s::S2_type)::S2_type

    # function to add a new exact solution to S2_type.m (exact solutions stored), S2_type.olb (oracle for lower bound), and S2_type.olb (oracle for upper bound)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    update_m_!(s) # store new exact solution in S2_type.m
    update_lb!(s) # store new exact solution in lower bound oracle
    update_ub!(s) # store new exact solution in upper bound oracle

    return s
end

function run_lb!(s::S2_type)::S2_type

    # function to solve S2_type.olb (oracle for lower bound) and store the solutions in S2_type.temp (temporary data)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.olb.objf = s.olb.vars'*vcat([1.],s.temp.c[:,s.temp.i]) # compute objective function for subproblem i (ϕ + γᵀcᵢ)
    @objective(s.olb.m, Min, s.olb.objf) # set objective function
    set_normalized_rhs.(s.olb.cons[1:s.olb.n], s.olb.help[1:s.olb.n].+ s.m.λ[1:s.olb.n,:]*s.temp.x[:,s.temp.i]) # update the rhs of the constraints (ϕ + γᵀcₛ >= hₛ + λₛᵀxᵢ, ∀s)
    optimize!(s.olb.m) # solve the JuMP model to optimality
    s.temp.θl[s.temp.i] = objective_value(s.olb.m) # store the valid lower bound θᵢ
    s.temp.λl[:,s.temp.i] .= s.m.λ[1:s.m.n,:]'*dual.(s.olb.cons[1:s.olb.n]) # store the valid subgradient λᵢ wrt xᵢ

    return s
end

function run_ub!(s::S2_type)::S2_type

    # function to solve S2_type.oub (oracle for upper bound) and store the solutions in S2_type.temp (temporary data)
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.oub.objf = s.oub.vars'*vcat([1.],s.temp.x[:,s.temp.i]) # compute objective function for subproblem i (ϕ + γᵀxᵢ)
    @objective(s.oub.m, Max, s.oub.objf) # set objective function
    set_normalized_rhs.(s.oub.cons[1:s.oub.n], s.oub.help[1:s.oub.n].+ s.m.ϕ[1:s.oub.n,:]*s.temp.c[:,s.temp.i])  # update the rhs of the constraints (ϕ + γᵀxₛ <= hₛ + ϕₛᵀcᵢ, ∀s)
    optimize!(s.oub.m) # solve the JuMP model to optimality
    s.temp.θu[s.temp.i] = objective_value(s.oub.m) # store the valid upper bound θᵢ
    s.temp.ϕu[:,s.temp.i] .= s.m.ϕ[1:s.m.n,:]'*dual.(s.oub.cons[1:s.oub.n]) # store the valid subgradient ϕᵢ wrt cᵢ

    return s
end

function run_oracles!(s::S2_type)::S2_type

    # function to call S2_type.olb (oracle for lower bound) and S2_type.oub (oracle for upper bound) to generate valid bounds on subproblem i
    # for given S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    run_lb!(s) # call the oracle for lower bound at (xᵢ,cᵢ)
    run_ub!(s) # call the oracle for upper bound at (xᵢ,cᵢ)

    return s
end

# */ --------------------------------------------------------------------------------------------- /* #

function step_a!(b::B1_type,s::S1_type)::Tuple{B1_type,S1_type}

    # function to perform 'step_a' of the Benders algorithm (solve the relaxed master problem and store decisions {xᵢ, ∀i})
    # for given B1_type (Benders master problem) and S1_type (Benders subproblem) -> algorithm 1 (Stand_Bend)

    optimize!(b.rmp.m) # solve the relaxed master problem to optimality
    b.temp.β .= value.(b.rmp.m[:β]) # store optimal values of {βᵢ, ∀i}
    b.temp.x .= value.(b.rmp.m[:x]) # store opitmal decisions {xᵢ, ∀i}
    s.temp.x .= b.temp.x # send decisions {xᵢ, ∀i} to the Benders subproblem

    return b,s
end

function step_a!(b::B3_type,s::S3_type)::Tuple{B3_type,S3_type}

    # function to perform 'step_a' of the Benders algorithm (solve the relaxed master problem and store decisions {xᵢ, ∀i})
    # for given B3_type (Benders master problem) and S3_type (Benders subproblem) -> algorithm 3 (Zaker_Bend)

    optimize!(b.rmp.m) # solve the relaxed master problem to optimality
    b.temp.β .= value.(b.rmp.m[:β]) # store optimal values of {βᵢ, ∀i}
    b.temp.x .= value.(b.rmp.m[:x]) # store opitmal decisions {xᵢ, ∀i}
    s.temp.x .= b.temp.x # send decisions {xᵢ, ∀i} to the Benders subproblem

    return b,s
end

function step_a!(b::B2_type,s::S2_type)::Tuple{B2_type,S2_type}

    # function to perform 'step_a' of the Benders algorithm (solve the relaxed master problem and store decisions {xᵢ, ∀i})
    # for given B2_type (Benders master problem) and S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    optimize!(b.rmp.m) # solve the relaxed master problem to optimality
    b.temp.β .= value.(b.rmp.m[:β]) # store optimal values of {βᵢ, ∀i}
    b.temp.x .= value.(b.rmp.m[:x]) # store opitmal decisions {xᵢ, ∀i}
    s.temp.x .= b.temp.x # send decisions {xᵢ, ∀i} to the Benders subproblem

    return b,s
end

function step_b!(b::B1_type)::B1_type

    # function to perform 'step_b' of the Benders algorithm (compute lower bound on optimal objective)
    # for given B1_type (Benders master problem) -> algorithm 1 (Stand_Bend)

    b.hist.L[b.hist.k] = exp10(6)*objective_value(b.rmp.m) # compute and store the lower bound Lₖ on the optimal objective @ iter k

    return b
end

function step_b!(b::B3_type)::B3_type

    # function to perform 'step_b' of the Benders algorithm (compute lower bound on optimal objective)
    # for given B3_type (Benders master problem) -> algorithm 3 (Zaker_Bend)

    b.hist.L[b.hist.k] = exp10(6)*objective_value(b.rmp.m) # compute and store the lower bound Lₖ on the optimal objective @ iter k

    return b
end

function step_b!(b::B2_type)::B2_type

    # function to perform 'step_b' of the Benders algorithm (compute lower bound on optimal objective)
    # for given B2_type (Benders master problem) -> algorithm 2 (Adapt_Bend)

    b.hist.L[b.hist.k] = exp10(6)*objective_value(b.rmp.m) # compute and store the lower bound Lₖ on the optimal objective @ iter k

    return b
end

function step_c!(b::B1_type,s::S1_type)::Tuple{B1_type,S1_type}

    # function to perform 'step_c' of the Benders algorithm (solve each subproblem exactly and store optimal solutions)
    # for given B1_type (Benders master problem) and S1_type (Benders subproblem) -> algorithm 1 (Stand_Bend)

    for s.temp.i in 1:b.data.unc.ni
        solv_exact!(s) # solve subproblem i exactly and store solutions
    end
    b.temp.θ .= s.temp.θ # send optimal objectives {θᵢ, ∀i} to the Benders master problem
    b.temp.λ .= s.temp.λ # send subgradients {λᵢ, ∀i} to the Benders master problem

    return b,s
end

function step_c!(b::B3_type,s::S3_type)::Tuple{B3_type,S3_type}

    # function to perform 'step_c' of the Benders algorithm (solve each subproblem to tolerance δ and store optimal solutions)
    # for given B3_type (Benders master problem) and S3_type (Benders subproblem) -> algorithm 3 (Zaker_Bend)

    δ = max(exp10(-8),b.data.q^(-b.hist.k+1.)) # compute optimality tolerance δ = q^(-k+1), where q is the damping factor and k is the current iter number
    JuMP.set_optimizer_attribute(s.ex.m,"BarConvTol",δ) # set optimality tolerance δ
    for s.temp.i in 1:b.data.unc.ni
        solv_exact!(s) # solve subproblem i up to tolerance δ and store solutions
    end
    b.temp.θp .= s.temp.θp # send primal objectives {θᵢ, ∀i} to the Benders master problem
    b.temp.θd .= s.temp.θd # send dual objectives {θᵢ, ∀i} to the Benders master problem
    b.temp.λ  .= s.temp.λ # send subgradients {λᵢ, ∀i} to the Benders master problem

    return b,s
end

function step_c!(b::B2_type,s::S2_type)::Tuple{B2_type,S2_type}

    # function to perform 'step_c' of the Benders algorithm (solve w subproblem exactly and store optimal solutions)
    # for given B2_type (Benders master problem) and S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    set_Iex!(b) # select the w subproblems to solve exactly
    for s.temp.i in b.temp.Ie
        solv_exact!(s) # solve subproblem i exactly and store solutions
        update_s!(s) # send the new exact solution to the adaptive oracles
        extended_sol!(s) # store the new extended exact solution
    end

    return b,s
end

function step_d!()

end

function step_d!(b::B2_type,s::S2_type)::Tuple{B2_type,S2_type}

    # function to perform 'step_d' of the Benders algorithm (solve the adaptive oracles to generate inexact but valid solutions)
    # for given B2_type (Benders master problem) and S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    for s.temp.i in 1:b.data.unc.ni
        run_oracles!(s) # solve the adaptive oracles for subproblem i and store the solutions
    end
    b.temp.θl .= s.temp.θl # send valid lower bound {θᵢ, ∀i} to the Benders master problem
    b.temp.θu .= s.temp.θu # send valid upper bound {θᵢ, ∀i} to the Benders master problem
    b.temp.λl .= s.temp.λl # send valid subgradients {λᵢ, ∀i} to the Benders master problem
    b.temp.ϕu .= s.temp.ϕu # send valid subgradients {ϕᵢ, ∀i} to the Benders master problem

    return b,s
end

function step_e!(b::B1_type)::B1_type

    # function to perform 'step_e' of the Benders algorithm (compute upper bound on optimal objective)
    # for given B1_type (Benders master problem) -> algorithm 1 (Stand_Bend)

    if (b.hist.k==1)
        b.hist.U[b.hist.k] = exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θ) # compute and store the upper bound Uₖ on the optimal objective @ iter k = 1
    else
        b.hist.U[b.hist.k] = min(b.hist.U[b.hist.k-1],exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θ)) # compute and store the upper bound Uₖ on the optimal objective @ iter k > 1
    end

    return b
end

function step_e!(b::B3_type)::B3_type

    # function to perform 'step_e' of the Benders algorithm (compute upper bound on optimal objective)
    # for given B3_type (Benders master problem) -> algorithm 3 (Zaker_Bend)

    if (b.hist.k==1)
        b.hist.U[b.hist.k] = exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θp) # compute and store the upper bound Uₖ on the optimal objective @ iter k = 1
    else
        b.hist.U[b.hist.k] = min(b.hist.U[b.hist.k-1],exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θp)) # compute and store the upper bound Uₖ on the optimal objective @ iter k > 1
    end

    return b
end

function step_e!(b::B2_type)::B2_type

    # function to perform 'step_e' of the Benders algorithm (compute upper bound on optimal objective)
    # for given B2_type (Benders master problem) -> algorithm 2 (Adapt_Bend)

    if (b.hist.k==1)
        b.hist.U[b.hist.k] = exp10(6)*(value(b.rmp.m[:f])+b.data.mp.κ*b.data.mp.π'*b.temp.θu) # compute and store the upper bound Uₖ on the optimal objective @ iter k = 1
    else
        b.hist.U[b.hist.k] = min(b.hist.U[b.hist.k-1], exp10(6)*(value(b.rmp.m[:f])+(b.data.mp.κ*b.data.mp.π'*b.temp.θu))) # compute and store the upper bound Uₖ on the optimal objective @ iter k > 1
    end

    return b
end

function step_f!(b::B1_type)::B1_type

    # function to perform 'step_f' of the Benders algorithm (compute optimality gap Δ and update the relaxed master problem)
    # for given B1_type (Benders master problem) -> algorithm 1 (Stand_Bend)

    b.temp.Δ = (b.hist.U[b.hist.k]-b.hist.L[b.hist.k])/b.hist.U[b.hist.k]*100 # compute optimality gap Δ = (Uₖ-Lₖ)/Uₖ @ iter k
    for i in 1:b.data.unc.ni
        @constraint(b.rmp.m,b.rmp.m[:β][i] >= b.temp.θ[i]+b.temp.λ[:,i]'*(b.rmp.m[:x][:,i].-b.temp.x[:,i])) # add constraints βᵢ >= θᵢ + λᵢᵀ(𝑥ᵢ-xᵢ) to the relaxed master problem
    end

    return b
end

function step_f!(b::B3_type)::B3_type

    # function to perform 'step_f' of the Benders algorithm (compute optimality gap Δ and update the relaxed master problem)
    # for given B3_type (Benders master problem) -> algorithm 3 (Zaker_Bend)

    b.temp.Δ = (b.hist.U[b.hist.k]-b.hist.L[b.hist.k])/b.hist.U[b.hist.k]*100 # compute optimality gap Δ = (Uₖ-Lₖ)/Uₖ @ iter k
    for i in 1:b.data.unc.ni
        @constraint(b.rmp.m,b.rmp.m[:β][i] >= b.temp.θd[i]+b.temp.λ[:,i]'*(b.rmp.m[:x][:,i].-b.temp.x[:,i])) # add constraints βᵢ >= θᵢ + λᵢᵀ(𝑥ᵢ-xᵢ) to the relaxed master problem
    end

    return b
end

function step_f!(b::B2_type)::B2_type

    # function to perform 'step_f' of the Benders algorithm (compute optimality gap Δ and update the relaxed master problem)
    # for given B2_type (Benders master problem) -> algorithm 2 (Adapt_Bend)

    b.temp.Δ = (b.hist.U[b.hist.k]-b.hist.L[b.hist.k])/b.hist.U[b.hist.k]*100 # compute optimality gap Δ = (Uₖ-Lₖ)/Uₖ @ iter k
    for i in 1:b.data.unc.ni
        @constraint(b.rmp.m,b.rmp.m[:β][i] >= b.temp.θl[i]+b.temp.λl[:,i]'*(b.rmp.m[:x][:,i].-b.temp.x[:,i])) # add constraints βᵢ >= θᵢ + λᵢᵀ(𝑥ᵢ-xᵢ) to the relaxed master problem
    end

    return b
end

# */ --------------------------------------------------------------------------------------------- /* #

function step_0!(b::B2_type,s::S2_type)::Tuple{B2_type,S2_type}

    # function to perform 'step_0' of the Benders algorithm (solve subproblem exactly at special point (x_min,c_min) to make the adaptive oracles always feasible)
    # for given B2_type (Benders master problem) and S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    s.temp.i = 1 # set i = 1
    s.temp.x[:,1] .= round.(vcat(minimum(b.data.mp.xh,dims=2)[:]*.99,minimum(b.data.unc.h,dims=1)[:].*[.99,1.01]); digits=4) # set x₁ = x_min
    s.temp.c[:,1] .= round.(minimum(b.data.unc.c,dims=1)[:]*.99; digits=4) # set c₁ = c_min
    solv_exact!(s) # solve subproblem 1 exactly (at point (x_min,c_min)) and store solutions
    update_s!(s) # send the new exact solution to the adaptive oracles
    s.temp.x .= vcat(b.data.mp.xh,b.data.unc.h') # set decisions {xᵢ, ∀i} on historial data
    s.temp.c .= b.data.unc.c' # set {cᵢ, ∀i}
    for s.temp.i in 1:b.data.unc.ni
        run_oracles!(s) # solve the adaptive oracles for subproblem i and store the solutions
    end
    b.temp.θl .= s.temp.θl # send valid lower bound {θᵢ, ∀i} to the Benders master problem
    b.temp.θu .= s.temp.θu # send valid upper bound {θᵢ, ∀i} to the Benders master problem
    b.temp.λl .= s.temp.λl # send valid subgradients {λᵢ, ∀i} to the Benders master problem
    b.temp.ϕu .= s.temp.ϕu # send valid subgradients {ϕᵢ, ∀i} to the Benders master problem

    return b,s
end

function step!(b::B1_type,s::S1_type)::Tuple{B1_type,S1_type}

    # function to perform a Benders 'step'
    # for given B1_type (Benders master problem) and S1_type (Benders subproblem) -> algorithm 1 (Stand_Bend)

    b.hist.k += 1                                # update iter k number
    b.hist.T[b.hist.k,1] = @elapsed step_a!(b,s) # solve the relaxed master problem and store decisions {xᵢ, ∀i}
    b.hist.T[b.hist.k,2] = @elapsed step_b!(b)   # compute lower bound on optimal objective
    b.hist.T[b.hist.k,3] = @elapsed step_c!(b,s) # solve each subproblem exactly and store optimal solutions
    b.hist.T[b.hist.k,4] = @elapsed step_d!( )   #
    b.hist.T[b.hist.k,5] = @elapsed step_e!(b)   # compute upper bound on optimal objective
    b.hist.T[b.hist.k,6] = @elapsed step_f!(b)   # compute optimality gap Δ and update the relaxed master problem

    return b,s
end

function step!(b::B3_type,s::S3_type)::Tuple{B3_type,S3_type}

    # function to perform a Benders 'step'
    # for given B3_type (Benders master problem) and S3_type (Benders subproblem) -> algorithm 3 (Zaker_Bend)

    b.hist.k += 1                                # update iter k number
    b.hist.T[b.hist.k,1] = @elapsed step_a!(b,s) # solve the relaxed master problem and store decisions {xᵢ, ∀i}
    b.hist.T[b.hist.k,2] = @elapsed step_b!(b)   # compute lower bound on optimal objective
    b.hist.T[b.hist.k,3] = @elapsed step_c!(b,s) # solve each subproblem to tolerance δ and store optimal solutions
    b.hist.T[b.hist.k,4] = @elapsed step_d!( )   #
    b.hist.T[b.hist.k,5] = @elapsed step_e!(b)   # compute upper bound on optimal objective
    b.hist.T[b.hist.k,6] = @elapsed step_f!(b)   # compute optimality gap Δ and update the relaxed master problem

    return b,s
end

function step!(b::B2_type,s::S2_type)::Tuple{B2_type,S2_type}

    # function to perform a Benders 'step'
    # for given B2_type (Benders master problem) and S2_type (Benders subproblem) -> algorithm 2 (Adapt_Bend)

    b.hist.k += 1                                # update iter k number
    b.hist.T[b.hist.k,1] = @elapsed step_a!(b,s) # solve the relaxed master problem and store decisions {xᵢ, ∀i}
    b.hist.T[b.hist.k,2] = @elapsed step_b!(b)   # compute lower bound on optimal objective
    b.hist.T[b.hist.k,3] = @elapsed step_c!(b,s) # solve w subproblem exactly and store optimal solutions
    b.hist.T[b.hist.k,4] = @elapsed step_d!(b,s) # solve the adaptive oracles to generate inexact but valid solutions
    b.hist.T[b.hist.k,5] = @elapsed step_e!(b)   # compute upper bound on optimal objective
    b.hist.T[b.hist.k,6] = @elapsed step_f!(b)   # compute optimality gap Δ and update the relaxed master problem

    return b,s
end

# */ --------------------------------------------------------------------------------------------- /* #

function solve_deterministic!(b::JuMP.Model,c::Int64;print_output::Int64=1,save_results::Int64=1)::JuMP.Model

    # function to solve the deterministic equivalent of the stochastic investment planning problem
    # for given Model (master problem) and c (case)

    print_output == 1 ? print0(b,c) : nothing
    println(" Sending JuMP model to the Optimizer")
    println(" ")
    print_output == 1 ? nothing : JuMP.set_optimizer_attribute(b,"OutputFlag",0)
    optimize!(b)
    print_output == 1 ? print_summary_a(b,c) : nothing

    return b
end

function solve_Benders!(b::B1_type,s::S1_type;print_output::Int64=1)::Tuple{B1_type,S1_type}

    # function to run Stand_Bend algorithm to solve the investment planning problem
    # for given B1_type (Benders master problem) and S1_type (Benders subproblem) -> algorithm 1 (Stand_Bend)

    print_output == 1 ? print0(b) : nothing
    while (b.hist.k < b.data.J) # check the iteration number k is lower than the limit J
        step!(b,s) # perform a Benders step
        print_output == 1 ? try print_info(b) catch; end : nothing
        (b.temp.Δ <= b.data.δ) ? break : nothing # check if the optimality gap Δₖ is lower than the targer tolerance δ
    end
    print_output == 1 ? print_summary(b) : nothing

    return b,s
end

function solve_Benders!(b::B3_type,s::S3_type;print_output::Int64=1)::Tuple{B3_type,S3_type}

    # function to run Zaker_Bend algorithm to solve the investment planning problem
    # for given B3_type (Benders master problem) and S3_type (Benders subproblem) -> algorithm 3 (Zaker_Bend)

    print_output == 1 ? print0(b) : nothing
    while (b.hist.k < b.data.J) # check the iteration number k is lower than the limit J
        step!(b,s) # perform a Benders step
        print_output == 1 ? try print_info(b) catch; end : nothing
        (b.temp.Δ <= b.data.δ) ? break : nothing # check if the optimality gap Δₖ is lower than the targer tolerance δ
    end
    print_output == 1 ? print_summary(b) : nothing

    return b,s
end

function solve_Benders!(b::B2_type,s::S2_type;print_output::Int64=1)::Tuple{B2_type,S2_type}

    print_output == 1 ? print0(b) : nothing
    step_0!(b,s) # perform step0 to solve the special point (x_min,c_min)
    while (b.hist.k < b.data.J) # check the iteration number k is lower than the limit J
        step!(b,s) # perform a Benders step
        print_output == 1 ? try print_info(b) catch; end : nothing
        (b.temp.Δ <= b.data.δ) ? break : nothing # check if the optimality gap Δₖ is lower than the targer tolerance δ
    end
    print_output == 1 ? print_summary(b) : nothing

    return b,s
end

# */ --------------------------------------------------------------------------------------------- /* #

ext: jl


text:
using CSPSolver

include("test_kenken.jl")
include("test_sudoku.jl")

max_stars_repo_name: gindij/CSPSolver

max_issues_repo_name: gindij/CSPSolver

max_forks_repo_name: gindij/CSPSolver


text:
# Routine for finding the roots of the boundary equation
include("roots.jl")
# Routine for computing the normalization constant of the quark wavefunction
include("norms.jl")
# Routine for evaluating the coefficients involved in the NJL mass gap equations
include("coefs.jl")
# Routine for solving the nonlinear system of differential equations, namely the NJL equations
include("solvenjl.jl")
# Routine for computing various observables
# include("observables.jl")

function njlradius(R)

    @info "Running for a fixed radius: " R

    xroot, lroot, ϵroot, κ = roots(R)
    norm = norms(R, xroot, lroot, ϵroot)
    coef, avgcoef, energ, j = coefs(R, xroot, lroot, ϵroot, κ, norm)
    mass = solvenjl(κ, coef, energ, j)
    # Mvavg, Msavg, Mr, Mavg, qq̅, ψ̅ψavg, ε = observables(R, xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass)

    # @info "Done!"

    return xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass

end

function saveradius(R, xroot, lroot, ϵroot, κ, norm, coef, energ, j, mass; foldername="results/")

    # @info "Saving results to file"

    # Saving the relevant results to file
    !isdir(foldername) && mkdir(foldername)
    fn = sprint(show, R, context=:compact=>true)
    filename = joinpath(foldername, "$fn.jld")
    save(filename, "xroot", xroot, "lroot", lroot, "ϵroot", ϵroot, "κ", κ, "norm", norm, "coef", coef, "energ", energ, "j", j, "mass", mass)

    # @info "Done!"

end

max_stars_repo_name: avramescudana/QuarkDroplets

max_issues_repo_name: avramescudana/QuarkDroplets

max_forks_repo_name: avramescudana/QuarkDroplets


text:
# Functions for computing indicators on solutions

# true is sol1 dominates sol2
function INDICATOR_is_dominating(sol1::Sol, sol2::Sol)

    debug && DEBUG_correct_solutions([sol1,sol2])

    return (sol1.z[1] >= sol2.z[1] && sol1.z[2] > sol2.z[2]) || (sol1.z[1] > sol2.z[1] && sol1.z[2] >= sol2.z[2])
end

# evaluates the rank of each solution in the given population
function INDICATOR_fast_non_dominated_sort!(population::Vector{Sol})

    debug && DEBUG_feasible_solutions(population)

    # true if no individual in population has rank idx_rank
    function rank_empty(population::Vector{Sol}, idx_rank::Int)
        for p in population
            if p.rank == idx_rank
                return false
            end
        end
        return true
    end

    # initialize the rank of each sol in the population
    for p in population
        p.rank = 0
    end

    n = length(population)
    S = [ [] for idx in 1:n] # S[i] list of indices of sols dominated by pop[i]
    domination_count = zeros(Int, n) # ...[i] number of sols that dominates pop[i]

    for idx_p in 1:n
        p = population[idx_p]
        for idx_q in idx_p+1:n
            q = population[idx_q]
            if INDICATOR_is_dominating(p, q) || p.z == q.z
                push!(S[idx_p], idx_q)
                domination_count[idx_q] += 1
            elseif INDICATOR_is_dominating(q, p)
                push!(S[idx_q], idx_p)
                domination_count[idx_p] += 1
            end
        end

        if domination_count[idx_p] == 0
            p.rank = 1
        end
    end

    idx_rank = 1

    while !rank_empty(population, idx_rank)
        for idx_p in 1:n
            if population[idx_p].rank == idx_rank
                for idx_q in S[idx_p]
                    domination_count[idx_q] -= 1
                    if domination_count[idx_q] == 0
                        population[idx_q].rank = idx_rank + 1
                    end
                end
            end
        end
        idx_rank += 1
    end

    # compute explicite fronts
    nb_ranks = idx_rank - 1
    F = Vector{Vector{Sol}}(undef, nb_ranks)
    for idx in 1:nb_ranks
        F[idx] = Vector{Sol}(undef, 0)
        for p in population
            if p.rank == idx
                push!(F[idx], p)
            end
        end
    end

    return F

end

# evaluates the crowding distance of each solution in the given front
# `Assertion` : front is a front, hence solutions are mutually non dominated
function INDICATOR_crowding_distance!(front::Vector{Sol})

    debug && DEBUG_feasible_solutions(front)
    debug && DEBUG_is_front(front)

    n = length(front)

    if n == 1

        front[1].crowding = UTIL_inf()

    else

        # initialize crowding distance for each solution
        for sol in front
            sol.crowding = 0
        end

        # permutation to sort front with regards to the first objective
        perm = sortperm(front, by = sol -> sol.z[1]) # O(nlog(n)))

        # set crowding of extreme sols to Infinity to always select them
        front[perm[1]].crowding = UTIL_inf()
        front[perm[n]].crowding = UTIL_inf()

        # first objective
        for idx in 2:n-1
            front[perm[idx]].crowding += front[perm[idx+1]].z[1] - front[perm[idx-1]].z[1]
        end

        # second objective
        for idx in 2:n-1
            front[perm[n+1-idx]].crowding += front[perm[n+1-(idx+1)]].z[2] - front[perm[n+1-(idx-1)]].z[2]
        end

    end

end

# evaluates the hypervolume contribution of each solution in the given front
# `Assertion` : front is a front, hence solutions are mutually non dominated
function INDICATOR_hypervolume_contribution!(front::Vector{Sol})

    debug && DEBUG_feasible_solutions(front)
    debug && DEBUG_is_front(front)

    N = length(front)

    # sort front wrt z1
    perm = sortperm(front, by = sol -> sol.z[1])

    # extreme solutions get an infinite hypervolume contribution
    front[perm[1]].hv_contrib = UTIL_inf()
    front[perm[N]].hv_contrib = UTIL_inf()

    for i in 2:N-1
        front[perm[i]].hv_contrib = (front[perm[i]].z[1] - front[perm[i-1]].z[1]) + (front[perm[i]].z[2] - front[perm[i+1]].z[2])
    end

end


text:
"""
    SoilWaterParameterizations

van Genuchten, Brooks and Corey, and Haverkamp parameters for and formulation of
  - hydraulic conductivity
  - matric potential

Hydraulic conductivity can be chosen to be dependent or independent of
impedance, viscosity and moisture.

Functions for hydraulic head, effective saturation, pressure head, matric
potential, and the relationship between augmented liquid fraction and liquid
fraction are also included.
"""
module SoilWaterParameterizations

using DocStringExtensions
using UnPack

export AbstractImpedanceFactor,
    NoImpedance,
    IceImpedance,
    impedance_factor,
    AbstractViscosityFactor,
    ConstantViscosity,
    TemperatureDependentViscosity,
    viscosity_factor,
    AbstractMoistureFactor,
    MoistureDependent,
    MoistureIndependent,
    moisture_factor,
    AbstractHydraulicsModel,
    vanGenuchten,
    BrooksCorey,
    Haverkamp,
    hydraulic_conductivity,
    effective_saturation,
    pressure_head,
    hydraulic_head,
    matric_potential,
    volumetric_liquid_fraction,
    inverse_matric_potential

"""
    AbstractImpedanceFactor{FT <: AbstractFloat}

"""
abstract type AbstractImpedanceFactor{FT <: AbstractFloat} end

"""
    AbstractViscosityFactor{FT <: AbstractFloat}
"""
abstract type AbstractViscosityFactor{FT <: AbstractFloat} end

"""
    AbstractMoistureFactor{FT <:AbstractFloat}
"""
abstract type AbstractMoistureFactor{FT <: AbstractFloat} end

"""
    AbstractsHydraulicsModel{FT <: AbstractFloat}

Hydraulics model is used in the moisture factor in hydraulic
conductivity and in the matric potential. The single hydraulics model
choice sets both of these.
"""
abstract type AbstractHydraulicsModel{FT <: AbstractFloat} end

"""
    vanGenuchten{FT} <: AbstractHydraulicsModel{FT}

The necessary parameters for the van Genuchten hydraulic model;
defaults are for Yolo light clay.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct vanGenuchten{FT} <: AbstractHydraulicsModel{FT}
    "Exponent parameter - used in matric potential"
    n::FT
    "used in matric potential. The inverse of this carries units in
     the expression for matric potential (specify in inverse meters)."
    α::FT
    "Exponent parameter - determined by n, used in hydraulic conductivity"
    m::FT
    function vanGenuchten{FT}(; n::FT = FT(1.43), α::FT = FT(2.6)) where {FT}
        new(n, α, FT(1) - FT(1) / FT(n))
    end
end

"""
    BrooksCorey{FT} <: AbstractHydraulicsModel{FT}

The necessary parameters for the Brooks and Corey hydraulic model.

Defaults are chosen to somewhat mirror the Havercamp/vG Yolo light
clay hydraulic conductivity/matric potential.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct BrooksCorey{FT} <: AbstractHydraulicsModel{FT}
    "ψ_b - used in matric potential. Units of meters."
    ψb::FT = FT(0.1656)
    "Exponent used in matric potential and hydraulic conductivity."
    m::FT = FT(0.5)
end

"""
    Haverkamp{FT} <: AbstractHydraulicsModel{FT}

The necessary parameters for the Haverkamp hydraulic model for Yolo light
 clay.

Note that this only is used in creating a hydraulic conductivity function,
 and another formulation for matric potential must be used.

# Fields
$(DocStringExtensions.FIELDS)
"""
struct Haverkamp{FT} <: AbstractHydraulicsModel{FT}
    "exponent in conductivity"
    k::FT
    "constant A (units of cm^k) using in conductivity. Our sim is in meters"
    A::FT
    "Exponent parameter - using in matric potential"
    n::FT
    "used in matric potential. The inverse of this carries units in the
     expression for matric potential (specify in inverse meters)."
    α::FT
    "Exponent parameter - determined by n, used in hydraulic conductivity"
    m::FT
    function Haverkamp{FT}(;
        k::FT = FT(1.77),
        A::FT = FT(124.6 / 100.0^1.77),
        n::FT = FT(1.43),
        α::FT = FT(2.6),
    ) where {FT}
        new(k, A, n, α, FT(1) - FT(1) / FT(n))
    end
end

"""
    MoistureIndependent{FT} <: AbstractMoistureFactor{FT} end

Moisture independent moisture factor.
"""
struct MoistureIndependent{FT} <: AbstractMoistureFactor{FT} end

"""
    MoistureDependent{FT} <: AbstractMoistureFactor{FT} end

Moisture dependent moisture factor.
"""
struct MoistureDependent{FT} <: AbstractMoistureFactor{FT} end

"""
    moisture_factor(
        mm::MoistureDependent{FT},
        hm::vanGenuchten{FT},
        S_l::FT,
    ) where {FT}

Returns the moisture factor of the hydraulic conductivy assuming a
MoistureDependent and van Genuchten hydraulic model.
"""
function moisture_factor(
    mm::MoistureDependent{FT},
    hm::vanGenuchten{FT},
    S_l::FT,
) where {FT}
    n = hm.n
    m = hm.m
    if S_l < FT(1)
        K = sqrt(S_l) * (FT(1) - (FT(1) - S_l^(FT(1) / m))^m)^FT(2)
    else
        K = FT(1)
    end
    return K
end

"""
    moisture_factor(
        mm::MoistureDependent{FT},
        hm::BrooksCorey{FT},
        S_l::FT,
    ) where {FT}

Returns the moisture factor of the hydraulic conductivy assuming a
MoistureDependent and Brooks/Corey hydraulic model.
"""
function moisture_factor(
    mm::MoistureDependent{FT},
    hm::BrooksCorey{FT},
    S_l::FT,
) where {FT}
    ψb = hm.ψb
    m = hm.m

    if S_l < 1
        K = S_l^(FT(2) * m + FT(3))
    else
        K = FT(1)
    end
    return K
end

"""
    moisture_factor(
        mm::MoistureDependent{FT},
        hm::Haverkamp{FT},
        S_l::FT,
    ) where {FT}

Returns the moisture factor of the hydraulic conductivy assuming a
MoistureDependent and Haverkamp hydraulic model.
"""
function moisture_factor(
    mm::MoistureDependent{FT},
    hm::Haverkamp{FT},
    S_l::FT,
) where {FT}
    k = hm.k
    A = hm.A
    if S_l < 1
        ψ = matric_potential(hm, S_l)
        K = A / (A + abs(ψ)^k)
    else
        K = FT(1)
    end
    return K
end

"""
    moisture_factor(mm::MoistureIndependent{FT},
                    hm::AbstractHydraulicsModel{FT},
                    S_l::FT,
    ) where {FT}
Returns the moisture factor in hydraulic conductivity when a
MoistureIndependent model is chosen. Returns 1.

Note that the hydraulics model and S_l are not used, but are included
as arguments to unify the function call.
"""
function moisture_factor(
    mm::MoistureIndependent{FT},
    hm::AbstractHydraulicsModel{FT},
    S_l::FT,
) where {FT}
    Factor = FT(1.0)
    return Factor
end

"""
    ConstantViscosity{FT} <: AbstractViscosityFactor{FT}

A model to indicate a constant viscosity - independent of temperature -
factor in hydraulic conductivity.
"""
struct ConstantViscosity{FT} <: AbstractViscosityFactor{FT} end

"""
    TemperatureDependentViscosity{FT} <: AbstractViscosityFactor{FT}

The necessary parameters for the temperature dependent portion of hydraulic
conductivity.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct TemperatureDependentViscosity{FT} <:
                   AbstractViscosityFactor{FT}
    "Empirical coefficient"
    γ::FT = FT(2.64e-2)
    "Reference temperature"
    T_ref::FT = FT(288.0)
end

"""
    viscosity_factor(
        vm::ConstantViscosity{FT},
        T::FT,
    ) where {FT}

Returns the viscosity factor when we choose no temperature dependence, i.e.
a constant viscosity. Returns 1.

T is included as an argument to unify the function call.
"""
function viscosity_factor(vm::ConstantViscosity{FT}, T::FT) where {FT}
    Theta = FT(1.0)
    return Theta
end

"""
    viscosity_factor(
        vm::TemperatureDependentViscosity{FT},
        T::FT,
    ) where {FT}

Returns the viscosity factor when we choose a TemperatureDependentViscosity.
"""
function viscosity_factor(
    vm::TemperatureDependentViscosity{FT},
    T::FT,
) where {FT}
    γ = vm.γ
    T_ref = vm.T_ref
    factor = FT(γ * (T - T_ref))
    Theta = FT(exp(factor))
    return Theta
end

"""
    NoImpedance{FT} <: AbstractImpedanceFactor{FT}

A model to indicate to dependence on ice for the hydraulic conductivity.
"""
struct NoImpedance{FT} <: AbstractImpedanceFactor{FT} end

"""
    IceImpedance{FT} <: AbstractImpedanceFactor{FT}

The necessary parameters for the empirical impedance factor due to ice.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef struct IceImpedance{FT} <: AbstractImpedanceFactor{FT}
    "Empirical coefficient from Hansson 2014. "
    Ω::FT = FT(7)
end

"""
    impedance_factor(
        imp::NoImpedance{FT},
        θ_i::FT,
        θ_l::FT,
    ) where {FT}

Returns the impedance factor when no effect due to ice is desired.
Returns 1.

The other arguments are included to unify the function call.
"""
function impedance_factor(imp::NoImpedance{FT}, θ_i::FT, θ_l::FT) where {FT}
    gamma = FT(1.0)
    return gamma
end

"""
    impedance_factor(
        imp::IceImpedance{FT},
        θ_i::FT,
        θ_l::FT,
    ) where {FT}

Returns the impedance factor when an effect due to the fraction of
ice is desired.
"""
function impedance_factor(imp::IceImpedance{FT}, θ_i::FT, θ_l::FT) where {FT}
    Ω = imp.Ω
    f_ice = θ_i / (θ_i + θ_l)
    gamma = FT(10.0^(-Ω * f_ice))
    return gamma
end

"""
    hydraulic_conductivity(
        impedance::AbstractImpedanceFactor{FT},
        viscosity::AbstractViscosityFactor{FT},
        moisture::AbstractMoistureFactor{FT},
        hydraulics::AbstractHydraulicsModel{FT},
        θ_i::FT,
        porosity::FT,
        T::FT,
        S_l::FT,
    ) where {FT}

Returns the hydraulic conductivity.
"""
function hydraulic_conductivity(
    impedance::AbstractImpedanceFactor{FT},
    viscosity::AbstractViscosityFactor{FT},
    moisture::AbstractMoistureFactor{FT},
    hydraulics::AbstractHydraulicsModel{FT},
    θ_i::FT,
    porosity::FT,
    T::FT,
    S_l::FT,
) where {FT}
    K = FT(
        viscosity_factor(viscosity, T) *
        impedance_factor(impedance, θ_i, porosity * S_l) *
        moisture_factor(moisture, hydraulics, S_l),
    )
    return K
end

"""
    hydraulic_head(z,ψ)

Return the hydraulic head.

The hydraulic head is defined as the sum of vertical height z and
pressure head ψ; meters.
"""
hydraulic_head(z, ψ) = z + ψ

"""
    volumetric_liquid_fraction(
        ϑ_l::FT,
        eff_porosity::FT,
    ) where {FT}

Compute the volumetric liquid fraction from the effective porosity and the augmented liquid
fraction.
"""
function volumetric_liquid_fraction(ϑ_l::FT, eff_porosity::FT) where {FT}
    if ϑ_l < eff_porosity
        θ_l = ϑ_l
    else
        θ_l = eff_porosity
    end
    return θ_l
end

"""
    effective_saturation(
        porosity::FT,
        ϑ_l::FT,
        θ_r::FT,
    ) where {FT}

Compute the effective saturation of soil.

`ϑ_l` is defined to be larger than `θ_r`. If `ϑ_l-θ_r` is negative,
hydraulic functions that take it as an argument will return
imaginary numbers, resulting in domain errors. Exit in this
case with an error.
"""
function effective_saturation(porosity::FT, ϑ_l::FT, θ_r::FT) where {FT}
    ϑ_l < θ_r && error("Effective saturation is negative")
    S_l = (ϑ_l - θ_r) / (porosity - θ_r)
    return S_l
end

"""
    pressure_head(
        model::AbstractHydraulicsModel{FT},
        param_functions::PS,
        S_s::FT,
        ϑ_l::FT,
        θ_i::FT,
    ) where {FT,PS}

Determine the pressure head in both saturated and unsaturated soil.

If ice is present, it reduces the volume available for liquid water.
The augmented liquid fraction changes behavior depending on if this
volume is full of liquid water vs not. Therefore, the region of saturated
vs unsaturated soil depends on porosity - θ_i, not just on porosity.
If the liquid water is unsaturated, the usual matric potential expression
is treated as unaffected by the presence of ice.
"""
function pressure_head(
    model::AbstractHydraulicsModel{FT},
    param_functions::PS,
    ϑ_l::FT,
    θ_i::FT,
) where {FT, PS}
    porosity = param_functions.porosity
    S_s = param_functions.S_s
    θ_r = param_functions.θ_r
    eff_porosity = porosity - θ_i
    if ϑ_l < eff_porosity
        S_l = effective_saturation(porosity, ϑ_l, θ_r)
        ψ = matric_potential(model, S_l)
    else
        ψ = (ϑ_l - eff_porosity) / S_s
    end
    return ψ
end

"""
    matric_potential(
            model::vanGenuchten{FT},
            S_l::FT
    ) where {FT}

Compute the van Genuchten function for matric potential.
"""
function matric_potential(model::vanGenuchten{FT}, S_l::FT) where {FT}
    @unpack n, m, α = model
    ψ_m = -((S_l^(-FT(1) / m) - FT(1)) * α^(-n))^(FT(1) / n)
    return ψ_m
end

"""
    matric_potential(
            model::Haverkamp{FT},
            S_l::FT
    ) where {FT}

Compute the van Genuchten function as a proxy for the Haverkamp model
matric potential (for testing purposes).
"""
function matric_potential(model::Haverkamp{FT}, S_l::FT) where {FT}
    @unpack n, m, α = model

    ψ_m = -((S_l^(-FT(1) / m) - FT(1)) * α^(-n))^(FT(1) / n)
    return ψ_m
end

"""
    matric_potential(
            model::BrooksCorey{FT},
            S_l::FT
    ) where {FT}

Compute the Brooks and Corey function for matric potential.
"""
function matric_potential(model::BrooksCorey{FT}, S_l::FT) where {FT}
    @unpack ψb, m = model
    ψ_m = -ψb * S_l^(-FT(1) / m)
    return ψ_m
end

"""
    inverse_matric_potential(
        model::vanGenuchten{FT},
        ψ::FT
    ) where {FT}

Compute the effective saturation given the matric potential, using
the van Genuchten formulation.
"""
function inverse_matric_potential(model::vanGenuchten{FT}, ψ::FT) where {FT}

    ψ > 0 && error("Matric potential is positive")

    @unpack n, m, α = model
    S = (FT(1) + (α * abs(ψ))^n)^(-m)
    return S
end

"""
    inverse_matric_potential(
        model::Haverkamp{FT}
        ψ::FT
    ) where {FT}

Compute the effective saturation given the matric potential using the
Haverkamp hydraulics model. This model uses the van Genuchten
formulation for matric potential.
"""
function inverse_matric_potential(model::Haverkamp{FT}, ψ::FT) where {FT}
    ψ > 0 && error("Matric potential is positive")
    @unpack n, m, α = model
    S = (FT(1) + (α * abs(ψ))^n)^(-m)
    return S
end

"""
    inverse_matric_potential(
        model::BrooksCorey{FT}
        ψ::FT
    ) where {FT}

Compute the effective saturation given the matric potential using the
Brooks and Corey formulation.
"""
function inverse_matric_potential(model::BrooksCorey{FT}, ψ::FT) where {FT}
    ψ > 0 && error("Matric potential is positive")

    @unpack ψb, m = model
    S = (-ψ / ψb)^(-m)
    return S
end

end #Module


text:
"""
    IVLogitShare <: GMMModel

An `IVLogitShare` model consists of outcomes, `y` ∈ (0,1),  regressors
`x` and instruments `z`.  The moment condition is

``E[ (\\log(y/(1-y)) - xβ)z ] = 0``

The dimensions of `x`, `y`, and `z` must be such that
`length(y) == size(x,1) == size(z,1)`
and
`size(x,2) ≤ size(z,2)`.
"""
struct IVLogitShare <: GMMModel
  x::Matrix{Float64}
  y::Vector{Float64}
  z::Matrix{Float64}
end

"""
    IVLogitShare(n::Integer, β::AbstractVector,
                      π::AbstractMatrix, ρ)

Simulate an IVLogitShare model.

# Arguments

- `n` number of observations
- `β` coefficients on `x`
- `π` first stage coefficients `x = z*π + v`
- `ρ` correlation between x[:,1] and structural error.

Returns an IVLogitShare GMMModel.
"""
function IVLogitShare(n::Integer, β::AbstractVector,
                      π::AbstractMatrix, ρ)
  z = randn(n, size(π)[1])
  endo = randn(n, length(β))
  x = z*π .+ endo
  ξ = rand(Normal(0,sqrt((1.0-ρ^2))),n).+endo[:,1]*ρ
  y = cdf.(Logistic(), x*β .+ ξ)
  return(IVLogitShare(x,y,z))
end

number_parameters(model::IVLogitShare) = size(model.x,2)
number_observations(model::IVLogitShare) = length(model.y)
number_moments(model::IVLogitShare) = size(model.z,2)

function get_gi(model::IVLogitShare)
  function(β)
    ξ = quantile.(Logistic(), model.y) .- model.x*β
    ξ.*model.z
  end
end

function gel_jump_problem(model::IVLogitShare)
  n = number_observations(model)
  d = number_parameters(model)
  k = number_moments(model)
  Ty = quantile.(Logistic(),model.y)
  m = Model()
  @variable(m, 0.0 <= p[1:n] <= 1.0)
  @variable(m, θ[1:d])
  @constraint(m, prob,sum(p)==1.0)
  @constraint(m, momentcon[i=1:k], dot((Ty - model.x*θ).*model.z[:,i],p)==0.0)
  @NLobjective(m,Max, sum(log(p[i]) for i in 1:n))
  return(m)
end


text:
function Base.size( ba::BOSSArray{T,N}) where {T,N}
    ([typemax(Int64) for i in 1:N]...)
end

function Base.eltype( ba::BOSSArray{T,N} ) where {T,N}
    T
end

function Base.ndims( ba::BOSSArray{T,N} ) where {T,N}
    N
end

function Base.reshape( ba::BOSSArray{T,N}, newShape ) where {T,N}
    warn("BOSSArray do not support reshaping")
end

"""
transform int to a UnitRange
"""
function idx2unitrange( idx::Union{UnitRange, Int} )
    if isa(idx, Int)
        return idx:idx
    else
        return idx
    end
end

function Base.getindex( ba::BOSSArray{T,3}, idxes::Union{UnitRange, Int} ... ) where T
    idxes = map(idx2unitrange, idxes)
    # construct the url
    # note that the start should -1 to match the coordinate system of numpy
    urlPath = "$(ba.urlPrefix)/cutout/$(ba.collectionName)"*
                    "/$(ba.experimentName)"*
                    "/$(ba.channelName)/$(ba.resolutionLevel)"*
                    "/$(idxes[1].start-1):$(idxes[1].stop)"*
                    "/$(idxes[2].start-1):$(idxes[2].stop)"*
                    "/$(idxes[3].start-1):$(idxes[3].stop)"
    resp = Requests.get(URI(urlPath); headers = ba.headers)
    data = Blosc.decompress(T, resp.data)
    data = reshape(data, map(length, idxes))
    return data
end

function Base.getindex( ba::BOSSArray{T,4}, idxes::Union{UnitRange, Int} ... ) where T
    idxes = map(idx2unitrange, idxes)
    # construct the url
    # note that the start should -1 to match the coordinate system of numpy
    urlPath = "$(ba.urlPrefix)/cutout/$(ba.collectionName)"*
                    "/$(ba.experimentName)"*
                    "/$(ba.channelName)/$(ba.resolutionLevel)"*
                    "/$(idxes[1].start-1):$(idxes[1].stop)"*
                    "/$(idxes[2].start-1):$(idxes[2].stop)"*
                    "/$(idxes[3].start-1):$(idxes[3].stop)"*
                    "/$(idxes[4].start-1):$(idxes[4].stop)"
    resp = Requests.get(URI(urlPath); headers = ba.headers)
    data = Blosc.decompress(T, resp.data)
    data = reshape(data, map(length, idxes))
    return data
end

function Base.setindex!(ba::BOSSArray{T,3}, buffer::AbstractArray,
                         idxes::Union{UnitRange, Int}...) where T
    @assert ndims(buffer)==3
    buffer = convert(Array{T,3}, buffer)
    idxes = map(idx2unitrange, idxes)
    # construct the url
    # note that the start should -1 to match the coordinate system of numpy
    urlPath = "$(ba.urlPrefix)/cutout/$(ba.collectionName)"*
                    "/$(ba.experimentName)"*
                    "/$(ba.channelName)/$(ba.resolutionLevel)"*
                    "/$(idxes[1].start-1):$(idxes[1].stop)"*
                    "/$(idxes[2].start-1):$(idxes[2].stop)"*
                    "/$(idxes[3].start-1):$(idxes[3].stop)"
    data = Blosc.compress( buffer )
    local resp
    for i in 1:3
        resp = Requests.post(URI(urlPath); data = data, headers = ba.headers)
        if statuscode(resp) == 201
            break
        else
            sleep(3)
        end
    end

    if statuscode(resp) != 201
        error("writting to boss have error: $resp")
    end
    return resp
end


text:
#=
2520 is the smallest number that can be divided by each of the numbers from 1 to 10 without any remainder.

What is the smallest positive number that is evenly divisible by all of the numbers from 1 to 20?
=#

# Ugly but fast.
function calc()
  current = 20
  maxdivisor = 20
  while true
    maxforcurrent = 1
    for i in 1:maxdivisor
      if rem(current, i) == 0
        maxforcurrent = i
      else
        break
      end
    end
    if maxforcurrent == maxdivisor
      break
    end
    current += 20
  end
  current
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler


text:
c = '\u25b3'

if ismatch(r"UTF", get(ENV, "LANG", ""))
    println("This output device supports Unicode: ", c)
else
    println("This output device does not support Unicode.")
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export iauTdbtt
"""
Time scale transformation:  Barycentric Dynamical Time, TDB, to
Terrestrial Time, TT.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  canonical.

Given:
   tdb1,tdb2  double    TDB as a 2-part Julian Date
   dtr        double    TDB-TT in seconds

Returned:
   tt1,tt2    double    TT as a 2-part Julian Date

Returned (function value):
               int       status:  0 = OK

Notes:

   1. tdb1+tdb2 is Julian Date, apportioned in any convenient way
   between the two arguments, for example where tdb1 is the Julian
   Day Number and tdb2 is the fraction of a day.  The returned
   tt1,tt2 follow suit.

   2. The argument dtr represents the quasi-periodic component of the
   GR transformation between TT and TCB.  It is dependent upon the
   adopted solar-system ephemeris, and can be obtained by numerical
   integration, by interrogating a precomputed time ephemeris or by
   evaluating a model such as that implemented in the SOFA function
   iauDtdb.   The quantity is dominated by an annual term of 1.7 ms
   amplitude.

   3. TDB is essentially the same as Teph, the time argument for the
   JPL solar system ephemerides.

References:

   McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
   IERS Technical Note No. 32, BKG (2004)

   IAU 2006 Resolution 3

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# int iauTdbtt(double tdb1, double tdb2, double dtr,
#              double *tt1, double *tt2 )

function iauTdbtt(tdb1::Real, tdb2::Real, dtr::Real)
   ref_tt1 = Ref{Float64}(0.0)
   ref_tt2 = Ref{Float64}(0.0)

   status = ccall((:iauTdbtt, libsofa_c), Cint,
                   (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
                   convert(Float64, tdb1), convert(Float64, tdb2), convert(Float64, dtr),
                   ref_tt1, ref_tt2)

   return status, ref_tt1[], ref_tt2[]
end


text:
__precompile__()

module Baseline
using Flux

using Common: minibatch, vecvectomat!, vecvectomat, copymodel, copymodel!

#using BSON: @save
using JLD2
using LinearAlgebra
using Distributed
using Random

abstract type AbstractBaseline end

struct Linear{DT<:AbstractFloat} <: AbstractBaseline
    nfeat::Int
    T::Int
    K::Int
    reg_coeff::DT
    coeff::AbstractVector{DT} # states x N

    # workspace
    featnd::Array{DT,3}
    A::AbstractMatrix{DT}
    bline::AbstractVector{DT}
    function Linear{DT}(reg::AbstractFloat, ns::Integer, T::Integer, K::Integer) where DT<:AbstractFloat
        nfeat = 2*ns+3+1
        featnd=zeros(DT, nfeat, T, K)
        A = zeros(DT, nfeat, nfeat)
        bline = zeros(DT, T*K)
        return new(nfeat, T, K, reg, zeros(DT, nfeat),
                   featnd, A, bline)
    end
end

struct Quadratic{DT<:AbstractFloat} <: AbstractBaseline
    nfeat::Int
    T::Int
    K::Int
    reg_coeff::DT
    coeff::AbstractVector{DT} # states x N

    # workspace
    featnd::Array{DT,3}
    A::AbstractMatrix{DT}
    bline::AbstractVector{DT}
    function Quadratic{DT}(reg::AbstractFloat, ns::Integer, T::Integer, K::Integer) where DT<:AbstractFloat
        nfeat = ns+convert(Int,ns*(ns+1)/2)+ns+ns+4+1
        featnd=zeros(DT, nfeat, T, K)
        gettimefeatures!(featnd)

        A = zeros(DT, nfeat, nfeat)
        bline = zeros(DT, T*K)
        return new(nfeat, T, K, reg, zeros(DT, nfeat),
                   featnd, A, bline)
    end
end

struct NN{DT<:AbstractFloat} <: AbstractBaseline
    nhidden::Int
    nfeat::Int
    T::Int
    K::Int

    features::Function
    loss::Function
    opt::Function
    epochs::Int

    #workspace
    featnd::Array{DT,2}
    m::Flux.Chain
    bline::AbstractVector{DT}
    function NN{DT}(ns::Integer, T::Integer, K::Integer;
                    nhidden=64,
                    epochs=3,
                    step=0.001,
                    bias=0.0, # l2 regularization is broken when bias vec is 0...
                    priorreg=false,
                    featurefunc=getNNfeaturesNoTime
                   ) where DT<:AbstractFloat

        nfeat = size(featurefunc(rand(ns,2)), 1)
        featnd = zeros(DT, nfeat, T*K)
        m = Chain(Dense(nfeat, nhidden, tanh),
                  #Dense(nhidden, nhidden, tanh),
        Dense(nhidden, nhidden, tanh),
        Dense(nhidden, 1, initb=(x)->fill(bias,x)) )
        #Dense(nhidden, 1, (x)->bias+x))

        #m[end].b.data .= bias
        m0 = deepcopy(params(m))

        bline = zeros(DT, T*K)
        #opt = ADAMW(params(m), step)

        if priorreg
            println("doing Prior function")
            my_vecnorm(m) = vecnorm(m.+0.00001)
            opt = ADAM(params(m), step)
            return new(nhidden, nfeat, T, K,
                       featurefunc,
                       (a,b)->Flux.mse(m(a), b)+sum(vecnorm, params(m)-m0.+0.00001)/length(b),
                       #(a,b)->sum(abs.(m(a)-b))/length(b)+sum(my_vecnorm, params(m)-m0)/length(b),
            opt,
            epochs,
            featnd, m, bline)
        else
            println("doing MSE Loss function")
            opt = ADAM(params(m), step)
            return new(nhidden, nfeat, T, K,
                       featurefunc,
                       (a,b)->Flux.mse(m(a), b), #+0.1*(sum(vecnorm, params(m).+0.00001)/length(b)),
                       #(a,b)->sum(abs.(m(a)-b))/length(b),
            opt,
            epochs,
            featnd, m, bline)
        end
    end
end

######################################################################## common

function gettimefeatures!(feat::Array{DT,3}) where DT<:AbstractFloat
    T    = size(feat, 2)
    K = size(feat, 3)
    al = convert(Array{DT}, linspace(1, T, T)') ./ DT(T)
    for t=1:K
        feat[end-4,:,t] = al
        feat[end-3,:,t] = al.^2
        feat[end-2,:,t] = al.^3
        feat[end-1,:,t] = al.^4
        feat[end,:,t]   .= DT(1.0)
    end
end

function getlinearfeatures(obs::AbstractMatrix{Float64})
    T = size(obs,2)
    # TODO add clipping
    o = max.(obs, -10.0)
    o = min.(o, 10.0)
    al = linspace(1, T, T)' / T
    feat = vcat(o, o.^2, al, al.^2, al.^3, ones(1, T))
end

function getquadfeatures!(feat::Array{DT,3},
                          obs::AbstractMatrix{DT},
                          t::Integer) where DT<:AbstractFloat
    n = size(obs, 1)
    T = size(feat, 2)

    feat[1:n,:,t] = obs[:,(t-1)*T+1:t*T] #implicite conversion between 32, 64
    feat[1:n,:,t] = max.(feat[1:n,:,t], DT(-10.0))
    feat[1:n,:,t] = min.(feat[1:n,:,t], DT(10.0))
    feat[1:n,:,t] /= DT(10.0)

    qfsize = convert(Int, n*(n+1)/2)
    qf = view(feat, (n+1):(n+qfsize) , :, t) #Matrix{Float64}(qfsize, T)
    k = 1
    @inbounds for i=1:n
        for j=i:n
            for x=1:T
                qf[k,x] = feat[i,x,t].*feat[j,x,t]
            end
            k += 1
        end
    end

    i0 = n+qfsize
    i1 = i0+n
    @inbounds for j = 1:n
        for i=1:T
            feat[i0+j,i,t] = sin(pi*feat[j,i,t])
            feat[i1+j,i,t] = cos(pi*feat[j,i,t])
        end
    end
end

function getNNfeatures(obs::AbstractMatrix{Float64})
    T = size(obs,2)
    o = max.(obs, -10.0)
    o = min.(o, 10.0)
    al = linspace(1, T, T)' / 1000.0
    feat = vcat(o, o.^2, al, al.^2, al.^3, ones(1, T))
    #feat = vcat(o, al, al.^2, al.^3, ones(1, T))
end

function getNNfeaturesNoTime(obs::AbstractMatrix{Float64})
    T = size(obs,2)
    o = max.(obs, -10.0)
    o = min.(o, 10.0)
    feat = vcat(o, o.^2, sin.(o), cos.(o), ones(1, T))
end

function prefit!(b::AbstractBaseline, returns::AbstractVector{DT}) where DT<:AbstractFloat
    featmat = reshape(b.featnd, b.nfeat, b.K*b.T)
    #b.A[:] = (featmat*featmat')
    mul!(b.A, featmat, transpose(featmat))

    return featmat*returns
end

# returns: T*K
function fit!(b::AbstractBaseline, returns::AbstractVector{DT}) where DT<:AbstractFloat

    target = prefit!(b, returns) # sets A matrix, returns target

    for i=1:b.nfeat # set diag
        b.A[i,i] += b.reg_coeff
    end
    #b.coeff[:] = b.A\(target)
    if isposdef(b.A)
        A_ldiv_B!(b.coeff, cholfact!(b.A), target) # needs A to be factorized
    else
        A_ldiv_B!(b.coeff, lufact!(convert(Array{Float64}, b.A)), target) # needs A to be factorized
    end
end

function fit!(b::NN, returns::AbstractVector{DT}) where DT<:AbstractFloat
    nsamples = length(returns)
    nbatch = 50

    X = b.featnd |> gpu
    Y = returns |> gpu

    rndidx = Flux.chunk(shuffle(1:nsamples), nbatch)

    datas = [ (X[:,i], Y[i]') for i in rndidx ]

    #opt = ADAM(params(b.m))
    #evalcb() = @show( norm(b.loss(b.featnd, returns')) / norm(returns) )
    evalcb() = @show( norm(b.loss(datas[1]...)) / norm(datas[1][2]) )
    evalcb(x, y) = norm(b.loss(x, y')) / norm(y)

    for i=1:b.epochs
        Flux.train!(b.loss, datas, b.opt)#, cb = Flux.throttle(evalcb, 1))
        #if evalcb(b.featnd, returns) < 0.5
        #if evalcb() < 0.5
        #   break
        #end
    end
end

function bellmanfit!(b::NN,
                     features::AbstractMatrix{DT},
                     nstep_features::AbstractMatrix{DT},
                     returns::AbstractVector{DT},
                     discount::Float64;
                     batch=16, epochs=20,
                     nbatch=ceil(Int, N/batch),
                     rndidx=minibatch(length(returns), batch)) where DT<:AbstractFloat
    N = length(returns)

    X = b.features(features) #|> gpu
    Y = returns #|> gpu
    Z = nstep_features #|> gpu

    #evalcb() = @show( norm(b.loss(datas[1]...)) / norm(datas[1][2]) )
    #evalcb(x, y) = norm(b.loss(x, y')) / norm(y)

    idx = 1:batch
    values = zeros(DT, batch)
    #for i=1:epochs
    for n=1:nbatch
        #idx = rndidx[n]
        values[:] = discount*Baseline.predict(b, Z[:,idx])
        data = [( @view( X[:,idx] ), (Y[idx]+values)' )]

        Flux.train!(b.loss, data, b.opt)

        idx += batch
    end
    #   rndidx = minibatch(N, batch)
    #end
end

######################################################################## linear

function predict!(b::Linear, obs::AbstractMatrix{T}) where T<:AbstractFloat
    nthread = min.(b.K, Threads.nthreads()) # cant have more threads than data
    Threads.@threads for tid=1:nthread
        #for tid=1:nthread
        thread_range = Distributed.splitrange(b.K, nthread)[tid]

        for t=thread_range
            b.featnd[:,:,t] = getlinearfeatures(obs[:,(t-1)*b.T+1:t*b.T])
        end
    end
    b.bline[:] = reshape(b.featnd, b.nfeat, b.K*b.T)'*b.coeff
end

##################################################################### quadratic
# TODO this is hacky and this whole module is a disgrace
function predict(b::Quadratic, obs::AbstractMatrix{DT}) where DT<:AbstractFloat
    nsamples = size(obs, 2)
    featnd=zeros(DT, b.nfeat, 1, nsamples)
    gettimefeatures!(featnd)
    for t=1:nsamples
        getquadfeatures!(featnd, obs, t)
    end
    values = reshape(featnd, b.nfeat, nsamples)'*b.coeff
    return values
end

function predict!(b::Quadratic, obs::AbstractMatrix{DT}) where DT<:AbstractFloat
    nthread = min.(b.K, Threads.nthreads()) # cant have more threads than data
    Threads.@threads for tid=1:nthread
        #for tid=1:nthread
        thread_range = Distributed.splitrange(b.K, nthread)[tid]

        for t=thread_range
            getquadfeatures!(b.featnd, obs, t)
        end
    end
    b.bline[:] = reshape(b.featnd, b.nfeat, b.K*b.T)'*b.coeff
end

##################################################################### Neural Network

function predict(b::NN, obs::AbstractMatrix{DT}) where DT<:AbstractFloat
    return b.m( b.features(obs) ).data[1,:]
end

function predict!(b::NN, obs::AbstractMatrix{DT}) where DT<:AbstractFloat
    b.featnd[:,:] = b.features(obs)
    b.bline .= b.m(b.featnd).data[1,:]
end

end # module


text:
using FiniteDifference
using Test

@testset "FiniteDifference.jl" begin
    # Write your tests here.
end