text:
export ExperienceBasedSampleModel, sample

import StatsBase: sample

"""
    ExperienceBasedSampleModel() -> ExperienceBasedSampleModel

Generate a transition based on previous experiences.
"""
mutable struct ExperienceBasedSampleModel <: AbstractEnvironmentModel
    experiences::Dict{
        Any,
        Dict{Any,NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}}},
    }
    sample_count::Int
    ExperienceBasedSampleModel() =
        new(
            Dict{
                Any,
                Dict{
                    Any,
                    NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}},
                },
            }(),
            0,
        )
end

function extract_experience(t::AbstractTrajectory, m::ExperienceBasedSampleModel)
    if length(t) > 0
        get_trace(t, :state)[end],
        get_trace(t, :action)[end],
        get_trace(t, :reward)[end],
        get_trace(t, :terminal)[end],
        get_trace(t, :next_state)[end]
    else
        nothing
    end
end

RLBase.update!(m::ExperienceBasedSampleModel, ::Nothing) = nothing

function RLBase.update!(m::ExperienceBasedSampleModel, transition::Tuple)
    s, a, r, d, s′ = transition
    if haskey(m.experiences, s)
        m.experiences[s][a] = (reward = r, terminal = d, nextstate = s′)
    else
        m.experiences[s] = Dict{
            Any,
            NamedTuple{(:reward, :terminal, :nextstate),Tuple{Float64,Bool,Any}},
        }(a => (reward = r, terminal = d, nextstate = s′))
    end
end

function sample(model::ExperienceBasedSampleModel)
    s = rand(keys(model.experiences))
    a = rand(keys(model.experiences[s]))
    model.sample_count += 1
    s, a, model.experiences[s][a]...
end

ext: jl


text:
# Sudoku Solver using backtracking

function ispossible(m::Array{Int, 4}, i::CartesianIndex{4}, n)
    return !(
        n ∈ m[  : ,   : , i[3], i[4]] || # column
        n ∈ m[i[1], i[2],   : ,   : ] || # line
        n ∈ m[  : , i[2],   : , i[4]]    # cell
    )
end

function sudokusolve!(::Array{Int, 4}, ::Nothing)
    # after `CartesianIndex(3,3,3,3)`, `iterate` returns `nothing`
    return true
end

function sudokusolve!(m::Array{Int, 4}, t::Tuple)
    return sudokusolve!(m, t[1])
end

function sudokusolve!(m::Array{Int, 4}, index::CartesianIndex{4})
    next_index = iterate(CartesianIndices(m), index)

    if m[index] != 0
        return sudokusolve!(m, next_index)
    end

    for n in 1:9
        if ispossible(m, index, n)
            m[index] = n
            if sudokusolve!(m, next_index)
                return true
            end
        end
    end

    m[index] = 0
    return false
end

function sudokusolve!(m::Array{Int, 2})
    sudokusolve!(reshape(m, 3, 3, 3, 3), CartesianIndex(1,1,1,1)) # reshape for easier manipulation
end

function string2matrix(s::String)
    return permutedims(reshape(parse.(Int, split(s, "")), 9, 9))
end

m = string2matrix("004300209005009001070060043006002087190007400050083000600000105003508690042910300")

println("\nOriginal :")
show(stdout, MIME("text/plain"), m)

sudokusolve!(m)

println("\n\nSolved :")
show(stdout, MIME("text/plain"), m)


text:
using NamedArrays, StatsPlots, JLD2, Measures

sumdimdrop(x::AbstractArray; dims) = dropdims(sum(x, dims=dims), dims=dims)

# Convert a JuMPArray to an AxisArray. Uses current JuMPArray internals, will need a rewrite in the next JuMP version.
AxisArrays.AxisArray(ja::JuMP.JuMPArray) = AxisArray(ja.innerArray, ja.indexsets...)

# Convert a JuMPDict to an AxisArray. Uses current JuMPDict internals, will need a rewrite in the next JuMP version.
# function AxisArrays.AxisArray(jd::JuMP.JuMPDict)
#   d = jd.tupledict
#   k = keys(d)
#   ndims = length(iterate(k)[1])
#   sets = [unique(getindex.(k,dim)) for dim=1:ndims]
#   a = AxisArray(zeros(length.(sets)...), sets...)
#   for (key, value) in d
#       a[key...] = value
#   end
#   return a
# end

# Convert a JuMPDict to a Dict. Uses current JuMPDict internals, will need a rewrite in the next JuMP version.
getdict(jd::JuMP.JuMPDict) = jd.tupledict

function readresults(model::ModelInfo, status::Symbol)
    @unpack REGION, TECH, CLASS, HOUR, techtype, STORAGECLASS = model.sets
    @unpack Systemcost, CO2emissions, FuelUse, Electricity, Charging, StorageLevel, Transmission, TransmissionCapacity, Capacity = model.vars
    @unpack demand, classlimits, hydrocapacity = model.params
    storagetechs = [k for k in TECH if techtype[k] == :storage]

    params = Dict(:demand => demand, :classlimits => classlimits, :hydrocapacity => hydrocapacity)

    cost = AxisArray(getvalue(Systemcost))
    emis = AxisArray(getvalue(CO2emissions))
    fuel = AxisArray(getvalue(FuelUse))
    # getting Electricity for all set combos is slowest, so let's optimize storage format and use faster internal function call
    elec = Dict((k,c) => [JuMP._getValue(Electricity[r,k,c,h]) for h in HOUR, r in REGION] for k in TECH for c in CLASS[k]);
    charge = AxisArray(getvalue(Charging))
    # oops, StorageLevel is also slow
    storage = Dict((k,c) => [JuMP._getValue(StorageLevel[r,k,c,h]) for h in HOUR, r in REGION] for k in storagetechs for c in STORAGECLASS[k]);
    transmission = AxisArray(getvalue(Transmission))
    transmissioncapac = AxisArray(getvalue(TransmissionCapacity))
    capac = getdict(getvalue(Capacity))

    return Results(status, model.options, model.hourinfo, model.sets, params, cost, emis, fuel, elec, charge, storage, transmission, transmissioncapac, capac)
end

function saveresults(results::Results, runname; resultsfile="", group="", compress=true)
    isempty(resultsfile) && return nothing
    if !isempty(group) && group[end] != '/'
        group *= "/"
    end
    runname = "$group$runname"
    JLD2.jldopen(resultsfile, "a+", compress=compress) do file
        if haskey(file, runname)
            @warn "The run $runname already exists in $resultsfile (new run not saved to disk). "
        else
            file[runname] = results
        end
    end
    return nothing
end

function listresults(; resultsfile="results.jld2", group="")
    if !isempty(group) && group[end] != '/'
        group *= "/"
    end
    JLD2.jldopen(resultsfile, "r") do file
        return keys(file)
    end
end

function loadresults(; resultsfile="results.jld2", group="", loadoptions...)
    options = merge(defaultoptions(), loadoptions)
    loadresults(autorunname(options); resultsfile=resultsfile, group=group)
end

function loadresults(runname::String; resultsfile="results.jld2", group="")
    if !isempty(group) && group[end] != '/'
        group *= "/"
    end
    runname = "$group$runname"
    results = nothing
    JLD2.jldopen(resultsfile, "r") do file
        if haskey(file, runname)
            results = file[runname]
        else
            println("\nThe run $runname does not exist in $resultsfile.")
        end
    end
    return results
end

const CHARTTECHS = Dict(
    :palette => [
        #RGB([68,131,208]/255...),  #hydroRoR
        RGB([216,137,255]/255...),  #nuclear            or RGBA{Float64}(0.76444,0.444112,0.824298,1.0)
        RGB([119,112,71]/255...),   #coal
        #RGB([164,155,104]/255...), #coal CCS
        RGB([199,218,241]/255...),  #wind
        RGB([149,179,215]/255...),  #wind offshore
        RGB([214,64,64]/255...),    #solarCSP
        RGB([255,255,64]/255...),   #solarPV
        RGB([240,224,0]/255...),    #PVrooftop
        RGB([255,192,0]/255...),    #CCGT
        RGB([99,172,70]/255...),    #bioCCGT
        RGB([100,136,209]/255...),  #hydro
        RGB([144,213,93]/255...),   #bioGT
        RGB([148,138,84]/255...),   #gasGT
        RGB([157,87,205]/255...)    #battery
    ],

    :labellookup => Dict(:nuclear => "nuclear", :coal => "coal", :wind => "onshore wind", :offwind => "offshore wind",
        :csp => "CSP", :pv => "PV plant", :pvroof => "PV rooftop", :gasCCGT => "gas CCGT", :bioCCGT => "bio CCGT",
        :hydro => "hydro", :bioGT => "bio GT", :gasGT => "gas GT", :battery => "battery"),

    :displaytechs => [:nuclear, :coal, :wind, :offwind, :csp, :pv, :pvroof, :gasCCGT, :bioCCGT, :hydro, :bioGT, :gasGT, :battery]
)

function analyzeresults(results::Results)
    @unpack REGION, FUEL, TECH, CLASS, HOUR, techtype, STORAGECLASS = results.sets
    @unpack demand, classlimits, hydrocapacity = results.params
    @unpack CO2emissions, FuelUse, Electricity, Transmission, Capacity, TransmissionCapacity, Charging, StorageLevel, Systemcost = results

    hoursperperiod = results.hourinfo.hoursperperiod

    capacmat = [sum(Capacity[r,k,c] for c in CLASS[k]) for k in TECH, r in REGION]
    capac = NamedArray([capacmat sum(capacmat, dims=2)], (TECH, [REGION; :TOTAL]))
    elec = zeros(length(HOUR), length(TECH), length(REGION))
    for (i,k) in enumerate(TECH)
        elec[:,i,:] = sum(Electricity[k,c] for c in CLASS[k])
    end
    annualelec = NamedArray([sumdimdrop(elec,dims=1) sumdimdrop(elec, dims=(1,3))], (TECH, [REGION; :TOTAL]), (:TECH, :REGION))
    charge = [Charging[r,:battery,h] for h in HOUR, r in REGION]
    storagetechs = [k for k in TECH if techtype[k] == :storage]
    storage = zeros(length(HOUR), length(storagetechs), length(REGION))
    for (i,k) in enumerate(storagetechs)
        storage[:,i,:] = sum(StorageLevel[k,c] for c in STORAGECLASS[k])
    end
    existingstoragelevel = NamedArray(storage, (collect(HOUR), storagetechs, REGION), (:HOUR, :TECH, :REGION))
    tcapac = NamedArray([TransmissionCapacity[r1,r2] for r1 in REGION, r2 in REGION], (REGION,REGION))

    # @unpack ElecDemand = model.constraints
    # prices = [getdual(ElecDemand[r,h]) for r in REGION, h in HOUR]
    # prices = NamedArray([getdual(ElecDemand[r,h]) for r in REGION, h in HOUR], (REGION, collect(HOUR)))       # €/kWh

    plotly()

    @unpack palette, labellookup, displaytechs = CHARTTECHS
    techlabels = [labellookup[k] for r=1:1, k in displaytechs]
    displayorder = [i for (i,k) in enumerate(TECH), d in displaytechs if d == k]

    function chart(country::Symbol; plotstoragetech=:none, plotbatterycharge=false, plotbatterylevel=false, optionlist...)
        if country == :BARS
            existinghydro = vec(sum(Electricity[:hydro,:x0], dims=1))
            regcost = Systemcost ./ (vec(sum(annualelec, dims=1)[1:end-1]) - existinghydro) * 1000
            totcost = sum(Systemcost) / (sum(annualelec[:,:TOTAL]) - sum(existinghydro)) * 1000
            lcoe = NamedArray(collect([regcost; totcost]'), (["system cost (€/MWh)"], [REGION; :TOTAL]))
            println("Regional system cost per MWh generated (€/MWh):")
            display(round.(lcoe, digits=2))

            lr = length(REGION)
            stackedbar(String.(REGION), collect(annualelec[displayorder,1:end-1]'/1000); labels=techlabels, size=(340+70*lr,550), left_margin=25px,
                legend=:outertopright, line=0, tickfont=10, legendfont=14, guidefont=14, color_palette=palette, ylabel="TWh/year", yformatter=:plain, optionlist...)

            xpos = (1:lr)' .- 0.5
            demandtext = ["demand" permutedims(repeat([""],lr-1))]
            display(plot!([xpos; xpos], [zeros(lr)'; sum(demand,dims=2)'*hoursperperiod/1000], line=3, color=:black,
                            labels=demandtext))
            if lr == 21
                totelec = [sumdimdrop(annualelec[:,1:8],dims=2) sumdimdrop(annualelec[:,9:15],dims=2) sumdimdrop(annualelec[:,16:21],dims=2) annualelec[:,:TOTAL]]
                stackedbar(["EU","CAS","China","TOTAL"], collect(totelec[displayorder,:]'/1e3); labels=techlabels, left_margin=20px,
                    legend=:outertopright, size=(500,950), line=0, tickfont=14, legendfont=14, color_palette=palette, yformatter=:plain, optionlist...)
                xpos = (1:4)' .- 0.5
                totdemand = [sum(demand[1:8,:]) sum(demand[9:15,:]) sum(demand[16:21,:]) sum(demand)]
                display(plot!([xpos; xpos], [zeros(4)'; totdemand*hoursperperiod/1e3], line=3, color=:black, labels=permutedims(repeat([""],4))))
                totcost2 = [sum(Systemcost[1:8]) sum(Systemcost[9:15]) sum(Systemcost[16:21]) sum(Systemcost)]
                tothydro = [sum(existinghydro[1:8]) sum(existinghydro[9:15]) sum(existinghydro[16:21]) sum(existinghydro)]
                lcoe_tot = NamedArray(totcost2./(totdemand .- tothydro) * 1000, (["system cost (€/MWh)"], ["EU","CAS","China","TOTAL"]))
                println("\nSystem cost per MWh demand (€/MWh):")
                display(round.(lcoe_tot, digits=2))
            else
                stackedbar(["TOTAL"], collect(annualelec[displayorder,:TOTAL]')/1e3; labels=techlabels, left_margin=30px,
                    size=(350,600), line=0, tickfont=14, legendfont=14, color_palette=palette, yformatter=:plain,
                    legend = :outertopright, optionlist...)
                xpos = (1:1)' .- 0.5
                totdemand = [sum(demand)]
                display(plot!([xpos; xpos], [zeros(1)'; totdemand*hoursperperiod/1e3], line=3, color=:black,
                    xlims=(-0.2,1.2), label="demand"))
            end
            country == :BARS && return nothing
        end

        if country == :TOTAL || country == :TOT || country == :total || country == :tot
            regs = 1:length(REGION)
        elseif country == :EU || country == :eu
            regs = 1:8
        elseif country == :ASIA || country == :asia
            regs = 9:15
        elseif country == :CHINA || country == :china || country == :CH || country == :ch
            regs = 16:21
        else
            countryindex = findfirst(REGION .== country)
            countryindex == nothing && error("Region $country not in $REGION.")
            regs = [countryindex]
        end

        regelec = sumdimdrop(elec[:,:,regs], dims=3)[:,displayorder] / hoursperperiod
        regcharge = sumdimdrop(charge[:,regs], dims=2) / hoursperperiod
        regdischarge = regelec[:, findfirst(displaytechs .== :battery)]
        regdemand = sumdimdrop(demand[regs,:], dims=1)

        composite = plot(; layout=grid(2,3,widths=[.4,.2,.4,.4,.2,.4]), size=(1850,950),
                        legend=false, tickfont=16, titlefont=20, optionlist...)
        for (i,k) in enumerate([:wind, :offwind, :hydro, :pv, :pvroof, :csp])
            colors = [palette[findfirst(displaytechs .== k)]; RGB(0.9,0.9,0.9)]
            classes = (k == :offwind || k == :pvroof) ? CLASS[k][1:length(CLASS[k])÷2] : CLASS[k]
            used = [sum(Capacity[r,k,c] for r in REGION[regs]) for c in classes]
            lims = [sum(k == :hydro ? hydrocapacity[r,c] : classlimits[r,k,c] for r in REGION[regs]) for c in classes]
            stackedbar!(String.(classes), [used lims-used]; subplot=i, title=labellookup[k], tickfont=14, legendfont=14, guidefont=14,
                            ylabel="GW", line=0, color_palette=colors, bottom_margin=30px, optionlist...)
        end
        display(composite)

        if plotstoragetech != :none
            regstorage = sumdimdrop(existingstoragelevel[:,[plotstoragetech],regs], dims=3, optionlist...)
            p = plot(regstorage; size=(900,550), tickfont=16, legendfont=16, label="storage level (TWh)")
            if plotstoragetech == :battery
                plot!(regcharge/1000, label="charge (TWh/h)")
                batteryelec = sumdimdrop([Electricity[:battery,:_][r,h] for r in REGION, h in HOUR][regs,:], dims=1) / hoursperperiod
                plot!(batteryelec/1000, label="discharge (TWh/h)")
            end
            display(p)
        end

        level = StorageLevel[:battery,:_]*1000
        reglevel = sumdimdrop(level[:,regs], dims=2)
        # display(plot(HOUR,[regcharge regelec[:,12] reglevel],size=(1850,950)))

        stackedarea(HOUR, regelec; labels=techlabels, size=(900,550), line=(0.03,1,:black), tickfont=14, legendfont=14, guidefont=14,
                                    xlabel="hour of year", ylabel="GW", yformatter=:plain, color_palette=palette, optionlist...)
        plotbatterycharge && plot!(HOUR, -regcharge, color=RGB([157,87,205]/255...))
        plotbatterylevel && plot!(HOUR, reglevel, line=(:black,:dash))
        # plot!(HOUR, regdischarge, color=:green)
        display(plot!(HOUR, regdemand, c=:black, line=2, label="demand"))
        nothing
    end

    # chart(:NOR)

    # if true   # plot hydro storage & shadow prices
    #   plot(elec[:,:hydro,:NOR]/hoursperperiod/1000, size=(1850,950), tickfont=16, legendfont=16)
    #   plot!(sum(elec[:,:wind,:],2)/hoursperperiod/1000)
    #   display(plot!(vec(mean(prices,1))))
    # end

    return annualelec, capac, tcapac, chart
end

function chart_energymix_scenarios(scenarios, resultsnames, resultsfile; size=(900,550), options...)
    numscen = length(scenarios)
    scenelec, demands, hoursperperiod, displayorder, techlabels, palette = allscenarioresults(scenarios, resultsnames, resultsfile)

    println("\nShare of demand (%):")
    demshare = [techlabels "-"; round.(scenelec[displayorder,:]'./demands .* 100, digits=1)  scenarios]
    display(reverse(permutedims(demshare), dims=1))

    stackedbar(collect(scenelec[displayorder,:]')/1e3, label=techlabels, size=size, left_margin=20px,
            xticks=(1:numscen,scenarios), line=0, tickfont=12, legendfont=12, guidefont=12,
            color_palette=palette, ylabel="TWh/year", yformatter=:plain; options...)
    xpos = (1:numscen)'
    lab = fill("",(1,numscen))
    lab[1] = "demand"
    display(plot!([xpos; xpos], [zeros(numscen)'; demands'*hoursperperiod/1e3], line=3, color=:black, label=lab))
    nothing
end

function allscenarioresults(scenarios, resultsnames, resultsfile)
    numscen = length(scenarios)
    scenelec = zeros(13,numscen)
    demands = zeros(numscen)
    hoursperperiod, displayorder, techlabels, palette = nothing, nothing, nothing, nothing

    for (i,s) in enumerate(scenarios)
        println("\nLoading results: $s...")
        totalelec, totaldemand, hoursperperiod, displayorder, techlabels, palette =
                readscenariodata(resultsnames[i], resultsfile)
        scenelec[:,i] = totalelec
        demands[i] = totaldemand
    end
    return scenelec, demands, hoursperperiod, displayorder, techlabels, palette
end

function readscenariodata(resultname, resultsfile)
    println(resultname, ": ", resultsfile)
    results = loadresults(resultname, resultsfile=resultsfile)

    annualelec, capac, tcapac, chart = analyzeresults(results);
    chart(:BARS)
    println()
    display(round.(Int, tcapac))

    @unpack TECH, REGION, CLASS, HOUR = results.sets
    hoursperperiod = results.hourinfo.hoursperperiod
    totaldemand = sum(results.params[:demand])
    totalelec = [sum(sum(results.Electricity[k,c]) for c in CLASS[k]) for k in TECH]

    @unpack palette, labellookup, displaytechs = CHARTTECHS
    techlabels = [labellookup[k] for r=1:1, k in displaytechs]
    displayorder = [i for (i,k) in enumerate(TECH), d in displaytechs if d == k]

    return totalelec, totaldemand, hoursperperiod, displayorder, techlabels, palette
end

function getlines(r, tcapac)
    reg = r.sets.REGION
    ii = findall(tcapac.>0)
    tc = [round(Int, tcapac[i]*1000) for i in ii if i[1]>i[2]]
    conn = [(reg[i[1]], reg[i[2]]) for i in ii if i[1]>i[2]]
    return sortslices([tc conn], dims=1, rev=true)
end

function exportwindresults(r, name)
    classes = r.sets.CLASS[:wind][1:5]
    electricity = Dict(c => sum(r.Electricity[:wind, c], dims=1) for c in classes)
    off_electricity = Dict(c => sum(r.Electricity[:offwind, c], dims=1) for c in classes)
    inputdata = getdatafolder(Dict(:datafolder => ""))
    open(joinpath(inputdata, "windresults_$name.csv"), "w") do io
        write(io, "region,cap1,cap2,cap3,cap4,cap5,")
        write(io, "offcap1,offcap2,offcap3,offcap4,offcap5,")
        write(io, "elec1,elec2,elec3,elec4,elec5,")
        write(io, "offelec1,offelec2,offelec3,offelec4,offelec5\n")
        for (i, reg) in enumerate(r.sets.REGION)
            cap = [round(1000 * r.Capacity[reg, :wind, c], digits=3) for c in classes]
            offcap = [round(1000 * r.Capacity[reg, :offwind, c], digits=3) for c in classes]
            elec = [round(electricity[c][i], digits=3) for c in classes]
            offelec = [round(off_electricity[c][i], digits=3) for c in classes]
            write(io, string(reg), ", ")
            write(io, join(cap, ", "), ", ")
            write(io, join(offcap, ", "), ", ")
            write(io, join(elec, ", "), ", ")
            write(io, join(offelec, ", "))
            write(io, "\n")
        end
    end
end

#=
default colors
get_color_palette(:auto, default(:bgcolor), 13)

 RGBA{Float64}(0.0,0.605603,0.97868,1.0)
 RGBA{Float64}(0.888874,0.435649,0.278123,1.0)
 RGBA{Float64}(0.242224,0.643275,0.304449,1.0)
 RGBA{Float64}(0.76444,0.444112,0.824298,1.0)
 RGBA{Float64}(0.675544,0.555662,0.0942343,1.0)
 RGBA{Float64}(4.82118e-7,0.665759,0.680997,1.0)
 RGBA{Float64}(0.930767,0.367477,0.57577,1.0)
 RGBA{Float64}(0.776982,0.509743,0.146425,1.0)
 RGBA{Float64}(3.80773e-7,0.664268,0.552951,1.0)
 RGBA{Float64}(0.558465,0.593485,0.117481,1.0)
 RGBA{Float64}(5.94762e-7,0.660879,0.798179,1.0)
 RGBA{Float64}(0.609671,0.499185,0.911781,1.0)
 RGBA{Float64}(0.380002,0.551053,0.966506,1.0)

function test()
_Systemcost = AxisArray(getvalue(m.vars.Systemcost))
_CO2emissions = AxisArray(getvalue(m.vars.CO2emissions))
_FuelUse = AxisArray(getvalue(m.vars.FuelUse))
_Electricity = Supergrid.getdict(getvalue(m.vars.Electricity))
_Charging = Supergrid.getdict(getvalue(m.vars.Charging))
_StorageLevel = Supergrid.getdict(getvalue(m.vars.StorageLevel))
_Transmission = AxisArray(getvalue(m.vars.Transmission))
_TransmissionCapacity = AxisArray(getvalue(m.vars.TransmissionCapacity))
_Capacity = Supergrid.getdict(getvalue(m.vars.Capacity))
end

@time elec = [sum(getvalue(m.vars.Electricity[r,k,c,h]) for c in CLASS[k]) for h in HOUR, k in TECH, r in REGION];

@time elec = Dict(getvalue(m.vars.Electricity[r,k,c,h]) for r in REGION, k in TECH, c in CLASS[k], h in HOUR);
@time elec = Dict((k,c) => [JuMP._getValue(m.vars.Electricity[r,k,c,h]) for r in REGION, h in HOUR] for k in TECH for c in CLASS[k]);
=#

max_stars_repo_name: niclasmattsson/Supergrid

max_issues_repo_name: niclasmattsson/Supergrid

max_forks_repo_name: niclasmattsson/Supergrid


text:
# Load docstring from markdown files:
for (name, path) in [:test_ordered => joinpath(@__DIR__, "test_ordered.md")]
    try
        include_dependency(path)
        str = read(path, String)
        @eval @doc $str $name
    catch err
        @error "Failed to import docstring for $name" exception = (err, catch_backtrace())
    end
end

function getlabel(x)
    i, example = x
    if example isa NamedTuple
        return example.label
    else
        return "$i"
    end
end

function getdata(x)
    i, example = x
    if example isa NamedTuple
        return example.data
    else
        return example
    end
end

const RECURSION_LIMIT = Ref(1000)

struct RecursionLimitError <: Exception
    f
    args::Tuple
    kwargs::NamedTuple
end

RecursionLimitError(f, args) = RecursionLimitError(f, args, NamedTuple())

function Base.showerror(io::IO, err::RecursionLimitError)
    println(io, "RecursionLimitError")
    print(io, "f = ", err.f)
    for (i, a) in enumerate(err.args)
        println(io)
        print(io, "args[$i] = ")
        show(io, "text/plain", a)
    end
    for (k, v) in pairs(err.kwargs)
        println(io)
        print(io, "kargs.$k = ")
        show(io, "text/plain", v)
    end
end

"""
    recursive_vcat(splittable_iterator, [amount]; recursion_limit)

Recursively call `halve` and `vcat`.  This should be equivalent to
`collect(splittable_iterator)` for ordered iterators.

The second argument is used for computing the length of an iterator.
By default [`amount`](@ref) is used.  If `length` is defined, passing
`length` should also work.
"""
recursive_vcat(data, _len = amount; recursion_limit = RECURSION_LIMIT[]) = recursive_vcat(
    data,
    _len,
    0,
    recursion_limit,
    RecursionLimitError(recursive_vcat, (data, _len), (recursion_limit = recursion_limit,)),
)

function recursive_vcat(data, _len, recursions, limit, err)
    recursions += 1
    recursions > limit && throw(err)
    _len(data) < 2 && return vec(collect(data))
    left, right = halve(data)
    @debug "recursive_vcat(data, $length)" _len(left) _len(right) data left right
    return vcat(
        recursive_vcat(left, _len, recursions, limit, err),
        recursive_vcat(right, _len, recursions, limit, err),
    )
end

function test_recursive_halving(x)
    @debug "Testing _recursive halving_: $(getlabel(x))"
    @testset "recursive halving" begin
        if Base.IteratorSize(getdata(x)) isa Union{Base.HasLength,Base.HasShape}
            @test isequal(recursive_vcat(getdata(x), length), vec(collect(getdata(x))))
        end
        @test isequal(recursive_vcat(getdata(x)), vec(collect(getdata(x))))
    end
end

function test_ordered(examples)
    @testset "$(getlabel(x))" for x in enumerate(examples)
        @debug "Testing `vcat`: $(getlabel(x))"
        @testset "vcat" begin
            data = getdata(x)
            left, right = halve(getdata(x))
            @test isequal(
                vcat(vec(collect(left)), vec(collect(right))),
                vec(collect(getdata(x))),
            )
        end
        test_recursive_halving(x)
    end
end

@deprecate test test_ordered false

function countmap(xs)
    counts = Dict{Any,Int}()
    for x in xs
        counts[x] = get(counts, x, 0) + 1
    end
    return counts
end

"""
    SplittablesTesting.test_unordered(examples)

See [`test_ordered`](@ref Main.SplittablesTesting.test_ordered).
"""
function test_unordered(examples)
    @testset "$(getlabel(x))" for x in enumerate(examples)
        @testset "concatenation" begin
            data = getdata(x)
            left, right = halve(getdata(x))
            @test isequal(
                countmap(collect(data)),
                merge(+, countmap(collect(left)), countmap(collect(right))),
            )
        end
        test_recursive_halving(x)
    end
end

const _PUBLIC_API = [
    # A list of public APIs to be picked by SplittablesTesting.jl
    :test_ordered,
    :test_unordered,
]


text:
# Exports
export access_geometry, visible
export spacecraft_compute_access
export compute_access
export parallel_compute_access
export create_lookup_location_opportunity

###################
# Access Geometry #
###################

"""
Compute the view geometry from an observer to a specific request.

Arguments:
- `sat_ecef:Array{<:Real, 1}` Earth-fixed satellite position
- `request::Request` Request object

Returns:
- `look_angle::Float64` Look angle from the satellite to the request center. Equivalent to off-nadiar angle. [deg]
- `range::Float64` Range from the observing satellite to the location. [m]
"""
function access_geometry(sat_ecef::Array{<:Real, 1}, request::Request)
    # Satellite state
    r_sat = sat_ecef[1:3]

    # Geodetic sub-satellte point
    sat_geod 	 = sECEFtoGEOD(r_sat)
    sub_sat_geod = [sat_geod[1], sat_geod[2], 0.0]
    sub_sat_ecef = sGEODtoECEF(sub_sat_geod)

    # Look angle
    nadir_dir  = (sub_sat_ecef - r_sat)/norm(sub_sat_ecef - r_sat)
    location_dir = (request.ecef - r_sat)/norm(request.ecef - r_sat)
    look_angle = acos(dot(nadir_dir, location_dir))*180.0/pi

    # Distance from sub-satellite point to location along direct line of sight
    range = norm(r_sat - request.ecef) # range to location

    return look_angle, range
end

"""
Compute the view geometry from an observer to a specific station.

Arguments:
- `sat_ecef:Array{<:Real, 1}` Earth-fixed satellite position
- `station::GroundStation` Station object

Returns:
- `elevation::Float64` Elevation of satellite with respect to station [deg]
"""
function access_geometry(sat_ecef::Array{<:Real, 1}, station::GroundStation)
    # Satellite state
    r_sat     = sat_ecef[1:3]
    r_station = station.ecef

    azimuth, elevation, range = sSEZtoAZEL(sECEFtoSEZ(r_station, r_sat), use_degrees=true)

    return elevation, range
end

"""
Computes whether an request is visible from a given state and return boolean true/false
result.

Arguments:
- `sat_ecef::Array{<:Real, 1}` Satellite position in Earth fixed frame
- `request::Request` Request being viewed

Returns:
- `visible::Bool` Indication of whether request is visible or not
"""
function visible(sat_ecef::Array{<:Real, 1}, request::Request)
    # Compute access geometry
    look_angle, range = access_geometry(sat_ecef, request)

    r = norm(sat_ecef[1:3])
    # theta_max = request.look_angle_max*pi/180.0
    # max_range = -r*cos(theta_max) - sqrt(r^2*cos(theta_max)^2 - (r^2-R_EARTH^2)^2)

    # To be fixed, no idea why this is complex:
    max_range = R_EARTH
    # max_range = 1.5*sqrt(r^2 - R_EARTH^2)

    # @debug "max range: $max_range"
    if request.look_angle_min <= look_angle &&
       look_angle <= request.look_angle_max &&
       range <= max_range

        return true
    else
        return false
    end
end

"""
Computes whether an station is visible from a given state and return boolean true/false
result.

Arguments:
- `sat_ecef::Array{<:Real, 1}` Satellite position in Earth fixed frame
- `station::GroundStation` Station being communicatted with

Returns:
- `visible::Bool` Indication of whether station is visible or not
"""
function visible(sat_ecef::Array{<:Real, 1}, station::GroundStation)
    elevation, range = access_geometry(sat_ecef, station)

    r = norm(sat_ecef[1:3])

    if station.elevation_min <= elevation
        # dist = norm(sat_ecef[1:3] - station.ecef[1:3])
        return true
    else
        return false
    end
end

###################
# Location Access #
###################

"""
Find the time when a location transitions between being accessible/not-accessible.
"""
function find_access_boundary(tle::TLE, location::Location, epc0::Epoch, step::Real; tol::Real=1.0e-1)

    # println("Step size: $step")

    epc = deepcopy(epc0)

    if abs(step) < tol
        # If time step is below desired tolerance do nothing
        return epc0
    else
        x_ecef   = sECItoECEF(epc, state(tle, epc))
        initial_visibility = visible(x_ecef, location)

        # While visibility is the same as it was initially step in time and
        # recompute visibility
        while visible(x_ecef, location) == initial_visibility
            epc    = epc + step
            x_ecef = sECItoECEF(epc, state(tle, epc))
        end

        # Reverse search direction and half step size
        next_step = -sign(step)*max(abs(step)/2, tol/2)

        return find_access_boundary(tle, location, epc, next_step, tol=tol)
    end
end

function spacecraft_compute_access(problem::SatPlanningProblem, spacecraft::Spacecraft,
            locations::Array{<:Location, 1};
            macro_step::Real=60.0,
            tol::Real=0.01,
            orbit_fraction::Real=0.5,
            id_offset::Integer=0)

    # Access is computed over planning horizon of problem (t_start -> t_end)

    # Compute orbital period
    T = orbit_period(sCARTtoOSC(state(spacecraft.tle, problem.t_start), use_degrees=true)[1])

    opportunities = Array{Opportunity, 1}(undef, 0)

    next_step = macro_step

    for loc in locations
        # Perform macro adaptive stepsize search
        epc = problem.t_start
        while epc < problem.t_end

            # println("Current step: $epc")

            x_ecef = sECItoECEF(epc, state(spacecraft.tle, epc))

            if visible(x_ecef, loc)
                # println("Found instant of visibility - $(visible(x_ecef, loc)) - $epc")

                # Search for AOS (before initial guess epoch)
                window_open = find_access_boundary(spacecraft.tle, loc, epc, -macro_step, tol=tol)
                # println("Found window_open boundary: $window_open")

                # Search for LOS (after initial guess epoch)
                window_close = find_access_boundary(spacecraft.tle, loc, epc, macro_step, tol=tol)
                # println("Found window_close boundary: $window_close")

                # Convert window_open/window_close to elapsed time values
                # window_open = get_rel_time(problem, window_open)
                # window_close = get_rel_time(problem, window_close)

                # If zero-doppler collection is required updated pass-times and geometry profile
                collect_duration = 0
                if typeof(loc) == Request && loc.require_zero_doppler == true
                    epc_mid = window_open + (window_close - window_open)/2.0
                    window_open  = epc_mid - loc.collect_duration/2.0
                    window_close = epc_mid + loc.collect_duration/2.0
                end

                # Create Oppourtunity
                if (window_close - window_open) > 0.0
                    if typeof(loc) == Request
                        collect = Collect(window_open,
                            window_close,
                            spacecraft=spacecraft,
                            location=loc)

                        push!(opportunities, collect)
                    elseif typeof(loc) == GroundStation
                        contact = Contact(window_open,
                            window_close,
                            spacecraft=spacecraft,
                            location=loc)

                        push!(opportunities, contact)
                    else
                        throw(ErrorException("Unknown type of location: $loc"))
                    end
                end

                # Step a most of an orbit because there (likely) won't be another
                # Acces until at least 1 orbit later
                epc += T*orbit_fraction

            else
                # Take a macro step because there's a good chance there won't be an
                # opportunity until the next one
                epc += macro_step
            end
        end
    end

    # Sort opportunities in ascending order
    sort!(opportunities, by = x -> x.t_start)

    # Apply IDs for collects in ascending time order
    for (idx, opp) in enumerate(opportunities)
        # Override computed opportunity id
        opp.id = idx + id_offset
    end

    return opportunities
end

function compute_access(problem::SatPlanningProblem;
            macro_step::Real=60.0, tol::Real=0.01,
            orbit_fraction::Real=0.5)

    # All Opportunities
    all_opportunities = Array{Opportunity, 1}(undef, 0)

    # Compute all opportunities for each spacecraft
    for spacecraft in problem.spacecraft
        @debug "Computing access for spacecraft: $(spacecraft.id)"
        opportunities = spacecraft_compute_access(problem,
                            spacecraft,
                            problem.locations,
                            macro_step=macro_step, tol=tol,
                            orbit_fraction=orbit_fraction)

        # Add opportunities to list of all opportunities
        push!(all_opportunities, opportunities...)
    end

    # Sort opportunities in ascending order
    sort!(all_opportunities, by = x -> x.t_start)

    # Apply IDs for collects in ascending time order
    id_offset = length(problem.opportunities)
    for (idx, opp) in enumerate(all_opportunities)
        # Override computed opportunity id
        opp.id = idx
    end

    # Add opportunities to problem
    problem.opportunities = all_opportunities

    # Create separate arrays of contacts and requests
    problem.contacts = filter(x -> typeof(x) == Contact, problem.opportunities)
    problem.collects = filter(x -> typeof(x) == Collect, problem.opportunities)

    # Zero Location Opportunity Counts
    for loc in problem.locations
        problem.lt_loc_opps[loc.id] = Integer[]
    end

    # Update opportunity lookup table
    for opp in problem.opportunities
        problem.lt_opportunities[opp.id] = opp

        if typeof(opp) == Contact
            problem.lt_contacts[opp.id] = opp
        elseif typeof(opp) == Collect
            problem.lt_collects[opp.id] = opp
        end

        # Location -> Opportunity Lookup
        push!(problem.lt_loc_opps[opp.location.id], opp.id)
    end

    # Update opportunities with number of collects remaining
    for (img_id, opp_ids) in problem.lt_loc_opps
        # Update opportunities with number of remaining collect opps
        for (i, opp_id) in enumerate(reverse(opp_ids))
            problem.lt_opportunities[opp_id].nr = i-1
        end
    end

    return
end

function parallel_compute_access(problem::SatPlanningProblem;
            macro_step::Real=60.0, tol::Real=0.01,
            orbit_fraction::Real=0.5)

    # Create anonymous function to map array inputs to
    # Must be declared at start of function for julia reasons...
    # x is a tuple (spacecraft, locations)
    fn(p, a, b, c) = x -> spacecraft_compute_access(p, x[1], x[2], macro_step=a, tol=b, orbit_fraction=c)

    # Create work assignments
    nw = nworkers()        # Number of workers
    lw = length(problem.locations) # Length of work

    # @debug "Have $lw items of work and $nw workers")

    wa = floor(Int, lw/nw) # Average work per worker
    wr = lw - nw*wa        # Remaining worker

    # Construct assignments (function inputs) for Workers
    assignments = Tuple{Spacecraft, Array{<:Location, 1}}[]

    for sc in problem.spacecraft
         # @debug "Worker 1 convering $(1):$(wa+wr)")
        push!(assignments, (sc, problem.locations[1:(wa+wr)]))
        for i in 2:nw
            # @debug "Worker $i convering $(1+wr+(i-1)*wa):$(i*wa+wr)")
            push!(assignments, (sc, problem.locations[(1+wr+(i-1)*wa):(i*wa+wr)]))
        end
    end

    # Execute Work in paralle
    results = pmap(fn(problem, macro_step, tol, orbit_fraction), assignments)

    # Aggregate and process results

    # Add opportunities to problem
    problem.opportunities = vcat(results...)

    # Sort opportunities in ascending order
    sort!(problem.opportunities, by = x -> x.t_start)

    # Create separate arrays of contacts and requests
    problem.contacts = filter(x -> typeof(x) == Contact, problem.opportunities)
    problem.collects = filter(x -> typeof(x) == Collect, problem.opportunities)

    # Zero Location Opportunity Counts
    for loc in problem.locations
        problem.lt_loc_opps[loc.id] = Integer[]
    end

    # Update opportunity lookup table
    for (id, opp) in enumerate(problem.opportunities)
        # Assign global opportunity ID
        opp.id = id

        # Create lookup entry
        problem.lt_opportunities[opp.id] = opp

        if typeof(opp) == Contact
            problem.lt_contacts[opp.id] = opp
        elseif typeof(opp) == Collect
            problem.lt_collects[opp.id] = opp
        end

        # Location -> Opportunity Lookup
        push!(problem.lt_loc_opps[opp.location.id], opp.id)
    end

    return
end

###############################
# Action Space Precomputation #
###############################

export precompute_action_space
function precompute_action_space(problem::SatPlanningProblem; enable_resources::Bool=false)

    # First compute feasible transition lookup
    lt_feasible_transitions = Dict{Tuple{Integer, Integer}, Bool}()

    println("Computing Feasibility Lookup Table:")
    for start_opp_idx in 1:length(problem.opportunities)
        start_opp_opp = problem.lt_opportunities[start_opp_idx]

        if start_opp_idx % 100 == 0
            println("Current Index: $start_opp_idx")
        end

        for end_opp_idx in (start_opp_idx+1):length(problem.opportunities)
            end_opp_opp = problem.lt_opportunities[end_opp_idx]

            # Check feasibility constraints
            valid = true

            for constraint in problem.constraints
                if valid == false
                    break
                end

                valid = valid && constraint(problem.lt_opportunities[start_opp_idx], problem.lt_opportunities[end_opp_idx])
            end

            # Push
            lt_feasible_transitions[(start_opp_idx, end_opp_idx)] = valid
        end
    end

    println("Current Index: $(length(problem.opportunities))")

    # Dict mapping (last collect/downlink action) -> {all possible actions}
    sunpoint_actions = Tuple{Opportunity, Integer}[]
    lt_feasible_actions = Dict{Tuple{Integer, Integer}, AbstractVector{Opportunity}}()

    # Below might work but it doesn't solve the problem of having sunpointed actions
    # be part of the
    println("Computing Action Space Lookup:")
    for ca_idx in 1:length(problem.opportunities)
        ca_opp = problem.lt_opportunities[ca_idx]

        if ca_idx % 100 == 0
            println("Current Index: $ca_idx")
        end

        # Create sunpoint action once for time
        if enable_resources == true && ca_idx < length(problem.opportunities)
            sp_action = Sunpoint(id=length(problem.opportunities)+length(sunpoint_actions)+1, t_start=problem.opportunities[ca_idx+1].t_start)
            push!(sunpoint_actions, (sp_action, ca_idx+1)) # Record equivalent index of sunpoint action
        end

        for lcdo_idx in 1:ca_idx
            lcdo_opp = problem.lt_opportunities[lcdo_idx]

            # Create Array for (current action, last action)
            lt_feasible_actions[(lcdo_opp.id, ca_opp.id)] = Opportunity[]

            for fa_idx in (ca_idx+1):length(problem.opportunities)
                if lt_feasible_transitions[(lcdo_idx, fa_idx)] == true
                    push!(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)], problem.opportunities[fa_idx])
                end

                if problem.solve_breadth > 0 && length(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)]) > problem.solve_breadth
                    break
                end
            end

            # Add sunpoint action as feasible action
            if enable_resources == true && ca_idx < length(problem.opportunities)
                push!(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)], sp_action)
            end
        end
    end

    println("Current Index: $(length(problem.opportunities))")

    if enable_resources == true
        println("Computing Sunpoint Indices:")
        for (sp_action, ca_idx) in sunpoint_actions
            ca_opp = problem.lt_opportunities[ca_idx]

            if ca_idx % 100 == 0
                println("Current Sunpoint Index: $ca_idx")
            end

            for lcdo_idx in 1:ca_idx
                lcdo_opp = problem.lt_opportunities[lcdo_idx]

                # Copy Exisiting feasible actions
                lt_feasible_actions[(lcdo_opp.id, sp_action.id)] = copy(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)])
            end
        end

        println("Current Sunpoint Index: $(length(problem.opportunities))")
    end

    # Update Problem Values
    problem.lt_feasible_actions = lt_feasible_actions
    if enable_resources == true
        problem.actions = vcat(problem.opportunities, [x[1] for x in sunpoint_actions])
    else
        problem.actions = problem.opportunities
    end
    sort!(problem.actions, by = x -> x.t_start)

    # Create Action Lookup Table
    lt_actions = Dict{Integer, Opportunity}()
    for action in problem.actions
        lt_actions[action.id] = action
    end
    problem.lt_actions = lt_actions
end

export pas_fa_work
function pas_fa_work(problem::SatPlanningProblem, start_indices::AbstractVector{Integer})
    # Array to store feasible transitoins
    feasible_transitions = Tuple{Integer, Integer, Bool}[]

    # println("Checking in! - $(length(problem.opportunities)) - Work: $start_indices")

    # For Opportunity compute all feasible transitions
    for start_opp_idx in start_indices
        start_opp_opp = problem.lt_opportunities[start_opp_idx]

        for end_opp_idx in (start_opp_idx+1):length(problem.opportunities)
            end_opp_opp = problem.lt_opportunities[end_opp_idx]

            # Check feasibility constraints
            valid = true

            for constraint in problem.constraints
                if valid == false
                    break
                end

                valid = valid && constraint(problem.lt_opportunities[start_opp_idx], problem.lt_opportunities[end_opp_idx])
            end

            # Push
            push!(feasible_transitions, (start_opp_idx, end_opp_idx, valid))
        end
    end

    return feasible_transitions
end

export pas_cas_work
function pas_cas_work(problem::SatPlanningProblem, lt_transitions::Dict{Tuple{Integer, Integer}, Bool}, current_index::AbstractVector{Integer}; enable_resources::Bool=false)
    println("Checking in!")
    flush(stdout)
    sunpoint_actions = Tuple{Opportunity, Integer}[]
    lt_feasible_actions = Dict{Tuple{Integer, Integer}, AbstractVector{Opportunity}}()

    for ca_idx in current_index
        ca_opp = problem.lt_opportunities[ca_idx]

        # Create sunpoint action once for time
        if enable_resources == true && ca_idx < length(problem.opportunities)
            sp_action = Sunpoint(id=length(problem.opportunities)+length(sunpoint_actions)+1, t_start=problem.opportunities[ca_idx+1].t_start)
            push!(sunpoint_actions, (sp_action, ca_idx+1)) # Record equivalent index of sunpoint action
        end

        for lcdo_idx in 1:ca_idx
            lcdo_opp = problem.lt_opportunities[lcdo_idx]

            # Create Array for (current action, last action)
            lt_feasible_actions[(lcdo_opp.id, ca_opp.id)] = Opportunity[]

            for fa_idx in (ca_idx+1):length(problem.opportunities)
                if lt_transitions[(lcdo_idx, fa_idx)] == true
                    push!(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)], problem.opportunities[fa_idx])
                end

                if problem.solve_breadth > 0 && length(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)]) > problem.solve_breadth
                    break
                end
            end

            # Add sunpoint action as feasible action
            if enable_resources == true && ca_idx < length(problem.opportunities)
                push!(lt_feasible_actions[(lcdo_opp.id, ca_opp.id)], sp_action)
            end
        end
    end

    return lt_feasible_actions, sunpoint_actions
end

export parallel_precompute_action_space
function parallel_precompute_action_space(problem::SatPlanningProblem; enable_resources::Bool=false)

    ## Compute Feasible Transitions
    fn(p) = w -> pas_fa_work(p, w)

    # Create work assignments to minimize message passing
    lw = length(problem.opportunities)  # Length of work
    nw = nworkers()                     # Number of workers
    wa = floor(Int, lw/nw)              # Average work per worker
    wr = lw - nw*wa                     # Remaining worker

    # Construct Work Assignment
    work = AbstractVector{Integer}[]

    push!(work, collect(1:(wa+wr)))
    for i in 2:nw
        push!(work, collect((1+wr+(i-1)*wa):(i*wa+wr)))
    end

    # Execute work in parallel, aggregate
    feasibility_results = vcat(pmap(fn(problem), work)...)

    # Aggrecate results
    lt_feasible_transitions = Dict{Tuple{Integer, Integer}, Bool}()

    for item in feasibility_results
        lt_feasible_transitions[(item[1], item[2])] = item[3]
    end

    ## Compute Action Space

    fn(p, lt_f, er) = w -> pas_cas_work(p, lt_f, w, enable_resources=er)

    # Create work assignments to minimize message passing
    lw = length(problem.opportunities)  # Length of work
    nw = nworkers()                     # Number of workers
    wa = floor(Int, lw/nw)              # Average work per worker
    wr = lw - nw*wa                     # Remaining worker

    # Construct Work Assignment
    work = AbstractVector{Integer}[]

    cids = collect(shuffle(1:length(problem.opportunities)))
    push!(work, collect(cids[1:(wa+wr)]))
    for i in 2:nw
        push!(work, collect(cids[(1+wr+(i-1)*wa):(i*wa+wr)]))
    end

    # Execute work in parallel, aggregate
    println("Starting Action Space Compute")
    action_space_results = pmap(fn(problem, lt_feasible_transitions, enable_resources), work)

    # Aggrecate results
    sunpoint_actions = Tuple{Opportunity, Integer}[]
    lt_feasible_actions = Dict{Tuple{Integer, Integer}, AbstractVector{Opportunity}}()

    for (lt_feas, sp_actions) in action_space_results

        # Copy Feasibility Lookup
        for (k, v) in lt_feas
            lt_feasible_actions[k] = v
        end

        # Copy Sunpoint Actions and update ID
        for spa in sp_actions
            # Update ID
            spa[1].id = length(problem.opportunities)+length(sunpoint_actions)+1
            push!(sunpoint_actions, spa)
        end

    end

    return lt_feasible_actions

    ## Compute Sunpoint Actions

    # Update Problem Values
    # problem.lt_feasible_actions = lt_feasible_actions
    # if enable_resources == true
    #     problem.actions = vcat(problem.opportunities, [x[1] for x in sunpoint_actions])
    # else
    #     problem.actions = problem.opportunities
    # end
    # sort!(problem.actions, by = x -> x.t_start)

    # # Create Action Lookup Table
    # lt_actions = Dict{Integer, Opportunity}()
    # for action in problem.actions
    #     lt_actions[action.id] = action
    # end
    # problem.lt_actions = lt_actions
end


text:
# This file was generated, do not modify it. # hide
@with mysim begin
    @with trace begin
        Ey = α .+ β .* x
        r = (y - Ey) / σ
        (;Ey, r)
    end
end

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
module LifeVal

#using Distributed

#@everywhere begin
    using DrWatson
    @quickactivate "LifeVal"
#end

using JuliaDB, DelimitedFiles, DataFrames, XLSX, ProgressMeter #, Traceur
using Random, Distributions #not technically required, but added here so that all modules use same project

using Pkg, Pkg.Artifacts
#using CUDA, Distributed    #wip

using LifeValStructures
using LifeValData
using LifeValFunc
using LifeValReserve
using LifeValGenerate

function real_main()
    runsettings, spcode, names, outfile = LifeValData.Initialise(filePath, fileName)
    policyData, basisData, fData = LifeValData.Read_Data(runsettings, spcode, names, filePath)

    reserveList = LifeValFunc.Table_Length(policyData, runsettings)

    o = nrow(runsettings) #ProgressMeter
    q = Progress(o, 1) #ProgressMeter

    @time for runtype in eachrow(runsettings)
        # eventually create struct called policyList
        policyList = LifeValFunc.Table_to_Policies(policyData[runtype.runno], :prod, inputOrder, PolicyType)

        #CALC UNIT PRICES, this must be set per runno and only required up to max projTerm
        #ulList = filter!(policy.spcode > 49, policyList)
        #not sure how to remove this hardcoding!
        glob_proj_max = 421 #maximum(policyList[80001].projMax) + 1 #because permiums paid upfront, need 1 more at end
        unitpriceList = Array{Float64}(undef, nrow(runsettings), length(fData), 2, glob_proj_max)

        #for fund in fData[runtype.runno,:]
        #    unitpriceList[runtype.runno, fund.datfram.ID,1,:] = LifeValFunc.unitPFunction(fund.datfram.UP, (fund.datfram.V_IR - fund.datfram.V_IC), glob_proj_max)
        #    unitpriceList[runtype.runno, fund.datfram.ID,2,:] = LifeValFunc.unitPFunction(fund.datfram.UP / (1 + fund.datfram.V_IR .- fund.datfram.V_IC), (fund.datfram.V_IR .- fund.datfram.V_IC), glob_proj_max)
        #end
        #unitpriceList[runtype.runno, x.ID, 1, :] = map(x -> LifeValFunc.unitPFunction((x[1].datfram.UP, (x[1].datfram.V_IR .- x[1].datfram.V_IC), glob_proj_max)...), fData[runtype.runno,:])
        #unitpriceList[runtype.runno, x.ID, 2, :] = map(x -> LifeValFunc.unitPFunction((x[2].datfram.UP /(1+(x[2].datfram.V_IR .- x[2].datfram.V_IC)), (x[2].datfram.UP, (x[2].datfram.V_IR .- x[2].datfram.V_IC), glob_proj_max)...), fData[runtype.runno,:])

        for fund in fData[runtype.runno,:]
            unitpriceList[runtype.runno, fund.datfram.ID, 1, 1] = fund.datfram.UP #.* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
            unitpriceList[runtype.runno, fund.datfram.ID, 2, 1] = fund.datfram.UP
            for t in 2:glob_proj_max
                unitpriceList[runtype.runno, fund.datfram.ID, 1, t] = unitpriceList[runtype.runno, fund.datfram.ID, 1, (t-1)]  .* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
                unitpriceList[runtype.runno, fund.datfram.ID, 2, t] = unitpriceList[runtype.runno, fund.datfram.ID, 2, (t-1)]  .* (1 .+ fund.datfram.P_IR .- fund.datfram.P_IC) .^ (1/12)
            end
        end

        reserveList[:, runtype.runno] = map(x -> LifeValReserve.Calc(x, (runtype, basisData, unitpriceList)...), policyList)

        #outputtab = table(pol = policyList, res = reserveList)
        #writelm("out1.csv", outputtab, ",")
        #writedlm(filePath * "basis65-" * runtype.runno * ".csv", reserveList, ',')
        next!(q) #ProgressMeter
    end

    #combined = table(policyData,reserveList)

    writedlm(filePath * outfile, reserveList, ',')
end

#THE GPU ENVIRONMENT IS NOT ACTIVATED AT THE MOMENT DUE TO CuArray ERROR (bits??)
function real_main_gpu()

    runsettings, spcode, names, outfile = LifeValData.Initialise(filePath, fileName)
    policyData, basisData, fData = LifeValData.Read_Data(runsettings, spcode, names, filePath)

    reserveList = LifeValFunc.Table_Length(policyData, runsettings)
    policyList1 = CuArray(policyData)  #this is wip

    o = nrow(runsettings) #ProgressMeter
    q = Progress(o, 1) #ProgressMeter

    @time for runtype in eachrow(runsettings)
        # eventually create struct called policyList
        policyList = LifeValFunc.Table_to_Policies(policyData[runtype.runno], :prod, inputOrder, PolicyType)
        policyList = CuArray(policyList)  #this is an ERROR???

        #CALC UNIT PRICES, this must be set per runno and only required up to max projTerm
        #ulList = filter!(policy.spcode > 49, policyList)
        #not sure how to remove this hardcoding!
        glob_proj_max = 421 #maximum(policyList[80001].projMax) + 1 #because permiums paid upfront, need 1 more at end
        unitpriceList = Array{Float64}(undef, nrow(runsettings), length(fData), 2, glob_proj_max)
        #for fund in fData[runtype.runno,:]
        #    unitpriceList[runtype.runno, fund.datfram.ID,:] = unitPFunction(fund.datfram.UP, (fund.datfram.V_IR .- fund.datfram.V_IC), glob_proj_max)
        #end

        for fund in fData[runtype.runno,:]
            unitpriceList[runtype.runno, fund.datfram.ID, 1, 1] = fund.datfram.UP #.* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
            unitpriceList[runtype.runno, fund.datfram.ID, 2, 1] = fund.datfram.UP
            for t in 2:glob_proj_max
                unitpriceList[runtype.runno, fund.datfram.ID, 1, t] = unitpriceList[runtype.runno, fund.datfram.ID, 1, (t-1)]  .* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
                unitpriceList[runtype.runno, fund.datfram.ID, 2, t] = unitpriceList[runtype.runno, fund.datfram.ID, 2, (t-1)]  .* (1 .+ fund.datfram.P_IR .- fund.datfram.P_IC) .^ (1/12)
            end
        end

        reserveList[:, runtype.runno] = map(x -> LifeValReserve.Calc(x, (runtype, basisData, unitpriceList)...), policyList)

        #outputtab = table(pol = policyList, res = reserveList)
        #writelm("out1.csv", outputtab, ",")
        #writedlm(filePath * "basis65-" * runtype.runno * ".csv", reserveList, ',')
        next!(q) #ProgressMeter
    end

    #combined = table(policyData,reserveList)

    writedlm(filePath * outfile, reserveList, ',')
end

#MAKE SURE TO USE THE PREMS OUTPUT FOR PREMIUMS!
function real_main_prem()
    runsettings, spcode, names, outfile = LifeValData.Initialise(filePath, fileName)
    policyData, basisData, fData = LifeValData.Read_Data(runsettings, spcode, names, filePath)

    premList = LifeValFunc.Table_Length(policyData, runsettings)

    o = nrow(runsettings) #ProgressMeter
    q = Progress(o, 1) #ProgressMeter

    @time for runtype in eachrow(runsettings)
        # eventually create struct called policyList
        policyList = LifeValFunc.Table_to_Policies(policyData[runtype.runno], :prod, inputOrder, PolicyType)

        #CALC UNIT PRICES, this must be set per runno and only required up to max projTerm
        #ulList = filter!(policy.spcode > 49, policyList)
        #not sure how to remove this hardcoding!
        glob_proj_max = 1000 #maximum(policyList[80001].projMax) + 1 #because permiums paid upfront, need 1 more at end
        unitpriceList = Array{Float64}(undef, nrow(runsettings), length(fData), 2, glob_proj_max)

        for fund in fData[runtype.runno,:]
            unitpriceList[runtype.runno, fund.datfram.ID, 1, 1] = fund.datfram.UP #.* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
            unitpriceList[runtype.runno, fund.datfram.ID, 2, 1] = fund.datfram.UP
            for t in 2:glob_proj_max
                unitpriceList[runtype.runno, fund.datfram.ID, 1, t] = unitpriceList[runtype.runno, fund.datfram.ID, 1, (t-1)]  .* (1 .+ fund.datfram.V_IR .- fund.datfram.V_IC) .^ (1/12)
                unitpriceList[runtype.runno, fund.datfram.ID, 2, t] = unitpriceList[runtype.runno, fund.datfram.ID, 2, (t-1)]  .* (1 .+ fund.datfram.P_IR .- fund.datfram.P_IC) .^ (1/12)
            end
        end

        premList[:, runtype.runno] = map(x -> LifeValReserve.Prem(x, (runtype, basisData, unitpriceList)...), policyList)

        #INSTEAD OF THE MAP ABOVE, THE FOLLOWING CAN ALSO BE USED, 2X SLOWER...
        #p = Progress(length(policyList), barglyphs=BarGlyphs("[=> ]"))
        #progress_map(policyList, progress=p) do x
        #    premList[:, runtype.runno] .= LifeValReserve.Prem(x, (runtype, basisData, unitpriceList)...)
        #end

        #outputtab = table(pol = policyList, res = reserveList)
        #writelm("out1.csv", outputtab, ",")
        #writedlm(filePath * "basis65-" * runtype.runno * ".csv", reserveList, ',')
        next!(q) #ProgressMeter
    end

    writedlm(filePath * outfile, premList, ',')
end

#POLICYLIST test environment
#policyTable = policyData[1]
#productIndex = :prod
#columnsIndex = inputOrder

#test environment
#policyList = Table_to_Policies(policyData[2], :prod, inputOrder)
#policy=policyList[1]
#runtype=runsettings[1,:]

function real_main_gen()
    inputs = DataFrame(XLSX.readtable(fileName, "randominput")...)
    totpolicies = []
    totprems = []
    @time for poltype in eachrow(inputs)
        poltypevec = Vector(poltype)
        policyList, premList =  LifeValGenerate.Gen(poltypevec)
        for k in 1:length(policyList[:,1])
            polvect=policyList[k,:]
            push!(totpolicies,polvect)
        end
        for k in 1:length(premList[:,1])
            premvect=premList[k,:]
            push!(totprems,premvect)
        end
    end

    writedlm(filePath * "inputs69.csv", totpolicies, ',')
    writedlm(filePath * "prems69.csv", totprems, ',')
end

# ---------------------------------------------------------------------------- #

function julia_main(gpu::Int64)
    try
        if gpu == 1 #using CUDA
            real_main_gpu()
        elseif gpu == 3  #using Distributed (NOT ENABLED YET)
            real_main_cpu()
        elseif gpu == 4  #e.g. using CUDA and Distributed (NOT ENABLED YET)
            real_main_cpu()
        elseif gpu == 10  #generating policies
            real_main_gen()
        elseif gpu == 20
            real_main_prem()
        else   #using Threads
            real_main()
        end
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

#THE FOLLOWING IS TESTED IN PROD ENVIRONMENT, NOT RUN IN TEST ENVIRONMENT
#@show ARGS
#@show Base.PROGRAM_FILE
#@show DEPOT_PATH
#@show LOAD_PATH
#@show pwd()
#@show Base.active_project()
#@show Threads.nthreads()
#@show Sys.BINDIR
#display(Base.loaded_modules)
#println()

#fooifier_path() = joinpath(artifact"fooifier", "bin", "fooifier" * (Sys.iswindows() ? ".exe" : ""))

#FIXED VALUES FOR ALL RUNS / CURRENT SETTINGS FOR ALL
PolicyType = [LifeValStructures.LifeAssurance, LifeValStructures.TermAssurance, LifeValStructures.Annuity, LifeValStructures.EndowmentAssurance, LifeValStructures.LifeAssuranceUL, LifeValStructures.EndowmentAssuranceUL]
inputOrder = (:id, :sex, :age, :Term, :Term_if, :Premium, :Sum_Assured, :spc,
    :prem_esc, :sa_esc, :esc_m, :units_1, :units_99, :alloc_1, :alloc_99)

filePath = datadir("exp_raw") * "\\"   #current fixed input/output directory

print("Do you want to repeat last run settings (y/n) ? \n\n")
inpu = readline()
if inpu == "y" || inpu == "Y" #note not working for GEN currently
    fileName = filePath * "lastrun.xlsx"
    print("settings: 0 = val (default), 1 = gpu, 2-9 = multicore/+, 10 = generate, 20=pricing \n\n")
    print("What do you want to do ? \n\n")
    inpu = readline()
    num = parse(Int64, inpu)
else
    print("settings: 0 = val (default), 1 = gpu, 2-9 = multicore/+, 10 = generate, 20=pricing \n\n")
    print("What do you want to do ? \n\n")
    inpu = readline()
    num = parse(Int64, inpu)

    # print("Please copy path of input files ? \n\n") #ENABLE FOR CUSTOM DIR
    # inpu = readline()
    # filePath = inpu * "\\"

    print("Please copy name of runsettings (5 char minimum) ? \n\n")
    inpu = readline()
    if inpu[end-4:end] != ".xlsx"   #appends xlsx
        inpu = inpu * ".xlsx"
    end
    fileName = filePath * inpu
end

julia_main(num)

end # module

#LifeVal.julia_main(num)

#temp converting data (again, just lazy to go and find it each time)
 # filePath = datadir("exp_raw") * "\\"
 # dataFile = filePath * "inputsA1.csv"
 # policyDataT = loadtable(dataFile)
 # JuliaDB.save(policyDataT, filePath * "policyDataA")


text:
using JSON, Random, ArgParse
import EvaluationOfRLAlgs

include("algconstructors.jl")
include("envconstructors.jl")

function runsweep(id, seed, envconfig, algname, save_dir, trials, paramdir=nothing)
    num_episodes = envconfig["num_episodes"]
    num_endepisodes = envconfig["num_endepisodes"]

    rng = Random.MersenneTwister(seed)

    save_name = "$(algname)_$(lpad(seed, 5, '0')).csv"

    save_path = joinpath(save_dir, save_name)

    open(save_path, "w") do f
        for trial in 1:trials
            env = create_env(envconfig, rng)
            agent, hyps, logp = sample_alg(algname, env, rng, num_episodes)
            meanret = 0.
            meaneret = 0.

            returns = EvaluationOfRLAlgs.run_agent!(env, agent, num_episodes, rng)
            meanret += mean(returns)
            ereturns = EvaluationOfRLAlgs.eval_agent(env, agent, num_endepisodes, rng)
            meaneret += mean(ereturns)

            result = join([hyps..., logp, meanret, meaneret], ',')
            write(f, "$(result)\n")
            flush(f)
            println("$trial \t $(result)")
            flush(stdout)
        end
    end

end

function main()
    s = ArgParseSettings()
    @add_arg_table s begin
        "--alg-name", "-a"
            help = "name of the algorithm to run"
            arg_type = String
            required = true
        "--env-file", "-e"
            help = "file path to JSON file containing environment to run"
            arg_type = String
            required = true
        "--log-dir", "-l"
            help = "folder directory to store results"
            arg_type = String
            required = true
        "--id"
            help = "identifier for this experiment. Used in defining a random seed"
            arg_type = Int
            required = true
        "--seed", "-s"
            help = "random seed is seed + id"
            arg_type = Int
            required = true
        "--trials", "-t"
            help = "number of random trials to run"
            arg_type = Int
            default = 1
    end

    parsed_args = parse_args(ARGS, s)
    aname = parsed_args["alg-name"]
    efile = parsed_args["env-file"]
    println("Excuting algorithm: $aname")
    println("Loading environment config: $efile")
    println("Logging Directory: $(parsed_args["log-dir"])")
    flush(stdout)
    envconfig = JSON.parsefile(efile, dicttype=Dict, inttype=Int64)["env"]
    println("Environment config:")
    println(envconfig)
    flush(stdout)
    save_dir = parsed_args["log-dir"]
    save_dir = joinpath(save_dir, envconfig["logname"])
    mkpath(save_dir)
    println("Saving results in: $save_dir")
    trials = parsed_args["trials"]
    id = parsed_args["id"]
    seed = parsed_args["seed"]

    println("id $id seed=(id+seed) $(id + seed)")
    flush(stdout)
    seed = id + seed

    runsweep(id, seed, envconfig, aname, save_dir, trials)

end

main()

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
module PlotThemes

using Plots, Colors
import Plots: _invisible, _themes
export PlotTheme, plot_theme

_255_to_1(c::Symbol, colors) = RGBA(map(x-> x/255,colors[c])...)
expand_palette(bg, palette; kwargs...) = [convert(RGBA,c) for c in  distinguishable_colors(20, vcat(bg, palette); kwargs...)][2:end]
immutable PlotTheme
    bg_primary
    bg_secondary
    lines
    text
    palette
    gradient
end

PlotTheme(bg_primary, bg_secondary, lines, text, palette) = PlotTheme(bg_primary, bg_secondary, lines, text, palette, nothing)

function add_plots_theme(s, theme)
    add_theme(s,
    bg = theme.bg_secondary,
    bginside = theme.bg_primary,
    fg       = theme.lines,
    fgtext  = theme.text,
    fgguide = theme.text,
    fglegend = theme.text,
    palette = theme.palette)
    if !(theme.gradient == nothing)
        PlotUtils.register_gradient_colors(s, theme.gradient)
    end
end

function plot_theme(s)
    PlotUtils._default_gradient[] = s
    Plots.set_theme(s)
end

include("dark.jl")
include("ggplot2.jl")
include("solarized.jl")
include("sand.jl")
end # module


text:
abstract type HyperPlaneType end

abstract type OptimalityCut <: HyperPlaneType end
abstract type FeasibilityCut <: HyperPlaneType end
abstract type LinearConstraint <: HyperPlaneType end
abstract type Unbounded <: HyperPlaneType end

struct HyperPlane{H <: HyperPlaneType, T <: Real, A <: AbstractVector}
    δQ::A
    q::T
    id::Int

    function (::Type{HyperPlane})(δQ::AbstractVector, q::Real, id::Int, ::Type{H}) where H <: HyperPlaneType
        T = promote_type(eltype(δQ),Float32)
        δQ_ = convert(AbstractVector{T},δQ)
        new{H, T, typeof(δQ_)}(δQ_,q,id)
    end
end
OptimalityCut(δQ::AbstractVector,q::Real,id::Int) = HyperPlane(δQ,q,id,OptimalityCut)
FeasibilityCut(δQ::AbstractVector,q::Real,id::Int) = HyperPlane(δQ,q,id,FeasibilityCut)
LinearConstraint(δQ::AbstractVector,q::Real,id::Int) = HyperPlane(δQ,q,id,LinearConstraint)
Unbounded(id::Int) = HyperPlane{[],Inf,id,Unbounded}

SparseHyperPlane{T <: Real} = HyperPlane{<:HyperPlaneType,T,SparseVector{T,Int64}}
SparseOptimalityCut{T <: Real} = HyperPlane{OptimalityCut,T,SparseVector{T,Int64}}
SparseFeasibilityCut{T <: Real} = HyperPlane{FeasibilityCut,T,SparseVector{T,Int64}}
SparseLinearConstraint{T <: Real} = HyperPlane{LinearConstraint,T,SparseVector{T,Int64}}

function (hyperplane::HyperPlane{H})(x::AbstractVector) where H <: HyperPlaneType
    if length(hyperplane.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match",length(hyperplane.δQ),length(x))))
    end
    return hyperplane.δQ⋅x,hyperplane.q
end

function (cut::HyperPlane{OptimalityCut})(x::AbstractVector)
    if length(cut.δQ) != length(x)
        throw(ArgumentError(@sprintf("Dimensions of the cut (%d)) and the given optimization vector (%d) does not match",length(cut.δQ),length(x))))
    end
    return cut.q-cut.δQ⋅x
end

bounded(hyperplane::HyperPlane) = true
bounded(hyperplane::HyperPlane{Unbounded}) = false
function optimal(cut::HyperPlane{OptimalityCut},x::AbstractVector,θ::Real,τ::Real)
    Q = cut(x)
    return θ > -Inf && abs(θ-Q) <= τ*(1+abs(Q))
end
function active(hyperplane::HyperPlane,x::AbstractVector,τ::Real)
    Gval,g = hyperplane(x)
    return abs(Gval-g) <= τ*(1+abs(Gval))
end
function satisfied(hyperplane::HyperPlane,x::AbstractVector,τ::Real)
    Gval,g = hyperplane(x)
    return Gval >= g - τ*(1+abs(Gval))
end
function satisfied(cut::HyperPlane{OptimalityCut},x::AbstractVector,θ::Real,τ::Real)
    Q = cut(x)
    return θ > -Inf && θ >= Q - τ*(1+abs(Q))
end
function gap(hyperplane::HyperPlane,x::AbstractVector)
    Gval,g = hyperplane(x)
    return Gval-g
end
function gap(cut::HyperPlane{OptimalityCut},x::AbstractVector,θ::Real)
    if θ > -Inf
        return θ-cut(x)
    else
        return Inf
    end
end
function lowlevel(hyperplane::HyperPlane{H,T,SparseVector{T,Int}}) where {H <: HyperPlaneType, T <: Real}
    return hyperplane.δQ.nzind,hyperplane.δQ.nzval,hyperplane.q,Inf
end
function lowlevel(cut::HyperPlane{OptimalityCut,T,SparseVector{T,Int}}) where T <: Real
    nzind = copy(cut.δQ.nzind)
    nzval = copy(cut.δQ.nzval)
    push!(nzind,length(cut.δQ)+cut.id)
    push!(nzval,1.0)
    return nzind,nzval,cut.q,Inf
end

# Constructors #
# ======================================================================== #
function OptimalityCut(subproblem::SubProblem)
    @assert status(subproblem.solver) == :Optimal
    λ = getduals(subproblem.solver)
    π = subproblem.π
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s,(i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -π*λ[i]*coeff
    end
    δQ = sparsevec(cols,vals,length(subproblem.x))
    q = π*getobjval(subproblem.solver)+δQ⋅subproblem.x

    return OptimalityCut(δQ, q, subproblem.id)
end
ArtificialCut(val::Real,dim::Int,id::Int) = OptimalityCut(sparsevec(zeros(dim)), val, id)

function FeasibilityCut(subproblem::SubProblem)
    @assert status(subproblem.solver) == :Infeasible "Trying to generate feasibility cut from non-infeasible subproblem"
    λ = getduals(subproblem.solver)
    cols = zeros(length(subproblem.masterterms))
    vals = zeros(length(subproblem.masterterms))
    for (s,(i,j,coeff)) in enumerate(subproblem.masterterms)
        cols[s] = j
        vals[s] = -λ[i]*coeff
    end
    G = sparsevec(cols,vals,length(subproblem.x))
    g = getobjval(subproblem.solver)+G⋅subproblem.x

    # Scale to avoid numerical issues
    scaling = abs(g)
    if scaling == 0
        scaling = maximum(G)
    end
    G = G/scaling

    return FeasibilityCut(G, g, subproblem.id)
end

function LinearConstraint(constraint::JuMP.LinearConstraint,i::Integer)
    sense = JuMP.sense(constraint)
    if sense == :range
        throw(ArgumentError("Cannot handle range constraints"))
    end
    cols = map(v->v.col,constraint.terms.vars)
    vals = constraint.terms.coeffs * (sense == :(>=) ? 1 : -1)
    G = sparsevec(cols,vals,constraint.terms.vars[1].m.numCols)
    g = JuMP.rhs(constraint) * (sense == :(>=) ? 1 : -1)

    return LinearConstraint(G,g,i)
end

function linearconstraints(m::JuMP.Model)
    constraints = Vector{HyperPlane{LinearConstraint}}(length(m.linconstr))
    for (i,c) in enumerate(m.linconstr)
        constraints[i] = LinearConstraint(c,i)
    end
    return constraints
end

Unbounded(subprob::SubProblem) = Unbounded(subprob.id)
# ======================================================================== #


text:
module PregnancyCoagDataAnalysisKit

    # include my codes -
    include("Include.jl")

    # export types -
    export VLResult

    # export functions -
    export load_study_data_set
    export extract_minimum_complete_data_set
    export evaluate_classifier

    # data transform -
    export z_score_transform_data_set
    export z_score_transform_array
    export z_score_transform_vector

    # mle classifier -
    export mle_fit_logistic_model_classifier
    export mle_logistic_model_classifier_cross_validation

    # ols methods -
    export ols_fit_linear_model
    export ols_fit_linear_model_cross_validation

end # module

max_stars_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_issues_repo_name: varnerlab/PregnancyCoagDataAnalysisKit

max_forks_repo_name: varnerlab/PregnancyCoagDataAnalysisKit


text:
module Testing
using SimpleSparseArrays
export isapprox_full

function isapprox_full(f, args_sparse...; kw...)
    args_full = map(args_sparse) do arg
        if arg isa SimpleSparseArray
            Array(arg)
        else
            arg
        end
    end
    ret_full = f(args_full...)
    ret_sparse = f(args_sparse...)
    isapprox(ret_full, ret_sparse; kw...)
end

end#module