text:
# inspect effect of averaging litter input compare to fluctuating litter input

using Sesam
using ModelingToolkit, DifferentialEquations
using DataFrames, Tables
using Distributions
using Chain
using MTKHelpers

# name both systems s, so that can use same keys in parameters
ss = sesam3_revMM(;name=:s)
s = se = seam3(;name=:s)
t1 = 0; t2 = 100
@named pl = plant_face_fluct(;t1,t2)

p = pC = Dict(
    s.ϵ_tvr => 0.45,   # carbon use efficiency of microbial tvr (part by predators
        #which respire and corresponding amount of N must be mineralized)
    s.κ_E => 0.8,     ##<< amount of recycling enzyme turnover by biomass (
        # added to assimilable, i.e. uptake instead of R)
    s.a_E => 0.001*365,   ##<< C biomass allocated to enzymes 1/day /microbial biomass
    s.m => 0.005*365,    ##<< maintenance respiration rate   1/day /microbial biomass,
    s.τ => 1/60*365,  ##<< biomass turnover rate (12 days)
    s.k_L => 1.0,       ##<< 1/(x years)
    #s.k_L => 5.0,       ##<< 1/(x years)   # formerly 1 year
    s.k_R => 1/(40.0),        ##<< 1/(x years) # to demonstrate changes on short time scale
    ss.k_mN_L => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
        # /yr enzyme turnover 60 times a year
    ss.k_mN_R => 0.05 * 60, # enzyme half-saturation constant, in magnitude of enzymes *
    # /yr enzyme turnover 60 times a year
    s.k_N => 60,
    s.k_m => 0.05,
    s.ϵ => 0.5,      ##<< carbon use efficiency for growth respiration
    #i_L => t -> 1 - exp(-t),  # litter input
    pl.i_L0 => 400.0,         # g/m2 input per year (half NPP)
    #pl.β_Ni0 => 25.0,
    pl.β_Ni0 => 30.0,
    #pl.i_IN0 => 0.0,   ##<< input of mineral N,
    pl.i_IN0 => 0.7,   ##<< 7kg/ha/yr
    pl.k_Lagr => 12/2, # above ground litter turnover of 2 month
    pl.k_PlantN0 => 100.0,
    #pl.k_PlantN0 => 2, # try with a lower rate - cannot resupply
    #
    #P from plant model parameters not used in CN-Sesam soil model
    pl.β_Pi0 => Inf, #25*20, ## leaf litter N:P ~20(massratio Kang10)
    pl.i_IP0 => Inf, #0.65,   ##<< input of mineral P, weathering: Table3 mixed sedimentary rocks 0.65g/m2/yr Hartmann14 10.1016/j.chemgeo.2013.10.025
    pl.s_EP0 => Inf, # 0.5, # plant 1/20 of typical total microbial enzyme synthesis flux
    pl.u_PlantPmax0 => Inf,
    pl.k_PlantP0 => Inf,
)
pN = Dict(
    s.i_BN => 0.4, ##<< potential immobilization flux rate
    s.β_NEnz => 3.1,     # Sterner02: Protein (Fig. 2.2.), high N investment (low P) need
    s.β_NB => 11.0,
    s.l_N => 0.96,       #0.00262647*365     ##<< leaching rate of mineralN lN IN
    #s.l_N => 0.0,
    s.ν_N =>  0.9,     # microbial N use efficiency accounting for apparent
    ## minceralization of N during uptake in heterogeneous soils
    # for N uptake: take the defaults which take as much N as supplied by litter
    # pl.u_PlantNmax0 => Inf32, # only constrained by k_PlantN in min function
    # pl.k_PlantN0 => 10.57, #0.0289652*365     ##<< plant uptake rate first order of IN
)
p = merge(pC, pN)

get_litter_input_fake_system = () -> begin
    # solve the pl system and setup fake pl system that returns i_L from solutions
    # pl,p need to be defined

    #pls = structural_simplify(pl)
    #@unpack Lagr, i_L0, β_Ni0, i_IN0 = pls
    #ps = Dict(i_L0 => p[pl.i_L0], β_Ni0  => p[pl.β_Ni0], i_IN0 => p[pl.i_IN0])
    @named plse = embed_system(pl)
    prob_pl = ODEProblem(plse, [pl.Lagr => p[pl.i_L0]/2 / p[pl.k_Lagr]], (-500,200), p)
    sol_pl = solve(prob_pl, Vern7(), reltol=1e-5) # litterfall with high accuracy
    plf = plant_face_fluct_fake(;name=pl.name,sys=pl, sol=sol_pl, t1, t2);
    #Plots.plot(sol_pl, tspan=(-1.5,-1), vars=[pl.i_L, pl.i_L_annual, pl.i_Lagr])
    #
    # repeat with annually averaged litter input - here without fake
    # ps_plse = ProblemParSetter(plse, (pl.share_autumn, pl.Lagr))
    # prob_pl_ann = update_statepar(ps_plse, (0.0, 0.0), prob_pl);
    # sol_pl_ann = solve(prob_pl_ann, Vern7())
    # plf_ann = plant_face_fluct_fake(;name=pl.name,sys=pl, sol=sol_pl_ann);

    pl_ann = plf_ann = plant_face(;name=:pl, t1, t2)
    plf, plf_ann, sol_pl
end
plf, plf_ann, sol_pl = get_litter_input_fake_system()

i_inspect_integrated_litter_input = () -> begin
    # compose a system that integrates pl.i_L
    function rep_plf(pl;name, simplify=true)
        @parameters t
        D = Differential(t)
        sts = @variables x(t)
        # integrate of i_L to check annual litter inputs
        tmp = compose(ODESystem([
            D(x) ~ pl.i_L
            ],t,sts,[];name), pl)
        simplify ? structural_simplify(tmp) : tmp
    end
    @unpack x = plf_rep
    # @named plf_rep = rep_plf(pl)
    # prob_tmp = ODEProblem(plf_rep, [x => 0.0, pl.Lagr => p[pl.i_L0]/2 / p[pl.k_Lagr]], (-200,0), p)
    @named plf_rep = rep_plf(plf)
    prob_tmp = ODEProblem(plf_rep, [x => 0.0], (-200,0), p)
    #sol_tmp = solve(prob_tmp, solver);
    # without constrainingn reltol, the litter inputs are all over the place
    sol_tmp = solve(prob_tmp, solver, reltol=1e-5);
    # compute longterm litter input
    tmp = (sol_tmp[x] ./ (sol_tmp.t .- sol_tmp.t[1]))
    lines(sol_tmp.t, tmp)
    lines(last(sol_tmp.t,100), last(tmp,100))
    # indeed there is slightly biased upwards 413/400
    # -> adjusted by providing biasfac in solution of plf of 400/413
    #
    ts1 = range(-40.0,-2.0,length=200)
    ts2 = ts1 .+ 1.0
    collect(ts2 .- ts1) # one year difference each
    tmp = sol_tmp(ts2, idxs=x).u .- sol_tmp(ts1, idxs=x).u
    lines(ts2, tmp)
    # the annual litter input is fluctuating quite a lot
end

sim_u0steady = () -> begin
    u0C = Dict(
        s.B => 17,
        s.L => 100,
        s.R => 1100,
        #s.cumresp => 0.0,
        s.α_R => 0.1, # enzyme synthesis into L # TODO model by optimality
        pl.Lagr => p[pl.i_L0]/2 / p[pl.k_Lagr]
    )
    u0C[s.α_L] = 1.0 - u0C[s.α_R]
    u0C2 = Dict(
        s.E_L => u0C[s.α_L] * p[s.a_E] * u0C[s.B] / (p[s.k_N]),
        s.E_R => u0C[s.α_R] * p[s.a_E] * u0C[s.B] / (p[s.k_N]),
    )
    u0N = Dict(
        s.I_N => 0.04, ##<< inorganic pool gN/m2
        s.L_N => u0C[s.L]/p[pl.β_Ni0],
        s.R_N => u0C[s.R]/calculate_β_NR_sesam3(p,s) #p[s.β_NB],
        )
    u0p = Dict(
        pl.i_L => 0.0, pl.β_Ni => 0.0, pl.u_PlantNmax => 0.0, pl.k_PlantN => 0.0, pl.i_IN => 0.0, #pl.s_one => 1.0,
    )
    u00 = merge(u0C, u0C2, u0N, u0p)
    #u00[s.R]/u00[s.R_N] # smaller p[s.β_NB]
    tspan_spinup = (-800.0,0.0)
    #@named ssp0 = plant_sesam_system(se,plf)
    #
    @named ssp0 = plant_sesam_system(se,plf_ann)  # include Enzyemes in u0, use se
    #@named ssp0 = plant_sesam_system(ss,plf_ann)
    prob0 = ODEProblem(ssp0, u00, tspan_spinup, p)
    # pss_tmp = ProblemParSetter(ssp0, (pl.share_autumn, pl.Lagr))
    # prob0s = update_statepar(pss_tmp, (0.0, 0.0), prob0)
    #prob0s = deepcopy(prob0)
    sol = sol_sesam3s0 = solve(prob0);
    #plot(sol, vars=[s.R, s.L])
    #plot(sol, tspan = (-5,0), vars=[s.R]) # only slight increase still
    k = first(states(ssp0))
    u0s = Dict(k => sol[k,end] for k in states(ssp0))
    #print(u0s) # outputs code that can b pasted for u0
end

u0 = sim_u0steady()
#u0 = merge(u0, Dict(s.I_N => u0[s.I_N]*20))

#---------- simulate sesam3 with varying litter inputs
using CairoMakie, AlgebraOfGraphics
set_aog_theme!()  # to get the consistent colors

import IterTools
using DataFrames, Chain
variants = @chain IterTools.product((:seam, :sesam),(:seasonal, :annual)) begin
    DataFrame()
    rename!([:enzyme,:litter])
    leftjoin(DataFrame(enzyme=[:seam, :sesam], linestyle=[:dash, :solid]), on=:enzyme)
    transform(:enzyme => (x -> Makie.current_default_theme().palette.color[][1:length(x)]) => :color)
    transform(AsTable(Cols(:enzyme,:litter)) => ByRow(t -> "$(t.enzyme)_$(t.litter)") => :label)
end
variants[!,:sol] = Vector{Any}(fill(missing,nrow(variants)))
#typeof(variants.sol)

# inspect what goes wrong at instability with displaying state
#tspan = (-80.0,50.0)
tspan = (-200.0,50.0)
solver = Vern7() # implicit method
#solver = Vern9() # implicit method
#solver = Rodas5()
@named sep = plant_sesam_system(se,plf)
ps_se = ProblemParSetter(sep, states(sep))
probe = ODEProblem(sep, u0, tspan, p)
@named ssp = plant_sesam_system(ss,plf)
probs = ODEProblem(ssp, u0, tspan, p)
ps_ss = ProblemParSetter(ssp, states(ssp))

ue_tmp = copy(probe.u0)
check_unstable_e = (dt,u,p,t) -> begin
    !any(isnan,u) && return false
    @show t,label_state(ps_se,u)
    ue_tmp[:] .= u
    true
end
us_tmp = copy(probs.u0)
check_unstable_s = (dt,u,p,t) -> begin
    !any(isnan,u) && return false
    @show t,label_state(ps_ss,u)
    us_tmp[:] .= u
    true
end

# first both seam3 and sesam3 with fluctuating litter input
sol = variants[findfirst(variants.label .== "seam_seasonal"),:sol] = sol_seam3f =
    solve(probe, solver; reltol=1e-5, unstable_check=check_unstable_e);
sol = variants[findfirst(variants.label .== "sesam_seasonal"),:sol] = sol_sesam3f =
    solve(probs, solver; reltol=1e-5, unstable_check=check_unstable_s);

i_inspect_instability = () -> begin
    probi = remake(probe, u0 = sol[end])
    soli = solve(probi, solver; tspan=(sol.t[end], last(tspan)), abstol=1e-8, unstable_check=check_unstable_seas)

    plot(soli, vars=[s.α_L, s.α_R])
    label_state(ps_seas, u_tmp)
    plot(sol_seam3f, vars=[s.I_N])
    plot(sol_sesam3f, vars=[s.I_N])

    tend = maximum(sol_sesam3f.t)
    ts = tend .+ (-0.5, +0.5)
    tse = tend .+ (-0.5, 0.0)
    plot(sol_sesam3f, tspan=ts) # check infinities
    plot(sol_sesam3f, tspan=ts, ylim=(0,1)) # infinities
    plot(sol_sesam3f, tspan=(-75,maximum(sol_sesam3f.t)), vars=[s.I_N], legend=:topleft)
    Plots.plot!(sol_seam3f, tspan=(-75,maximum(sol_sesam3f.t)+1), vars=[s.I_N])

    plot(sol_sesam3f, tspan=ts, vars=[s.I_N])
    Plots.plot!(sol_sesam3f, tspan=ts, vars=[s.I_N])
    Plots.plot!(sol_sesam3f, tspan=ts, vars=[s.i_L])
    Plots.plot!(sol_seam3f, tspan=ts, vars=[s.i_L])

    plot(sol_seam3f, tspan=ts, vars=[s.α_R], xlim=ts)
    plot!(sol_sesam3f, tspan=tse, vars=[s.α_R], xlim=ts)

    plot(sol_sesam3f, tspan=tse, vars=[s.E_R])
    Plots.plot!(sol_seam3f, tspan=ts, vars=[s.E_R])

    plot(sol_seam3f, tspan=ts, vars=[s.C_synBC, s.C_synBN], xlim=ts)
    plot!(sol_sesam3f, tspan=tse, vars=[s.C_synBC, s.C_synBN], xlim=ts)

    plot(sol_seam3f, tspan=ts, vars=[s.w_C, s.w_N], xlim=ts)
    Plots.plot!(sol_sesam3f, tspan=tse, vars=[s.w_C, s.w_N], xlim=ts)

    plot(sol_seam3f, vars=[s.w_C, s.w_N])
    plot(sol_sesam3f, vars=[s.w_C, s.w_N])

    plot(sol_sesam3f, tspan=(-75,maximum(sol_sesam3f.t)), vars=[s.R])
    plot(sol_sesam3f, tspan=(-75,maximum(sol_sesam3f.t)), vars=[s.u_PlantN, s.u_PlantNmax])

    plot(sol_sesam3s, tspan=(-75,maximum(sol_sesam3f.t)), vars=[s.u_PlantN, s.u_PlantNmax])
    plot(sol_seam3s, tspan=(-75,maximum(sol_sesam3f.t)+1), vars=[s.u_PlantN, s.u_PlantNmax])
    plot(sol_sesam3s, tspan=(-10,0), vars=[s.I_N])
    plot(sol_sesam3s, tspan=(-10,0), vars=[s.leach_N])
    plot(sol_sesam3s, tspan=(-10,0), vars=[s.R/s.R_N])
    plot(sol_sesam3s, vars=[s.R])
end

# next parameterize plant model for non-fluctuating litter input
# this is achieved by setting the share of autumn input to zero
# pse_s = ProblemParSetter(sep, (pl.share_autumn, pl.Lagr))
# probe_ann = update_statepar(pse_s, (0.0, 0.0), probe);
# pss_s = ProblemParSetter(ssp, (pl.share_autumn, pl.Lagr))
# probs_ann = update_statepar(pss_s, (0.0, 0.0), probs);

@named sep_ann = plant_sesam_system(se,plf_ann)
probe_ann = ODEProblem(sep_ann, u0, tspan, p)
@named ssp_ann = plant_sesam_system(ss,plf_ann)
probs_ann = ODEProblem(ssp_ann, u0, tspan, p)

sol = variants[findfirst(variants.label .== "seam_annual"),:sol] = sol_seam3s =
    solve(probe_ann, solver; abstol=1e-8);
sol = variants[findfirst(variants.label .== "sesam_annual"),:sol] = sol_sesam3s =
    solve(probs_ann, solver; abstol=1e-8);

function plot_vars(vars, tspan=tspan; kwargs...)
    pl = plot(sol_seam3f; tspan=ts, vars=vars, label="seam", linestyle=:dash, kwargs...)
    plot!(pl, sol_seam3s; tspan=ts, vars=vars, label="seam annual", linestyle=:dash, kwargs...)
    plot!(pl, sol_sesam3s; tspan=ts, vars=vars, label="sesam annual", kwargs...)
    plot!(pl, sol_sesam3f; tspan=ts, vars=vars, label="sesam", xlab="time (yr)", kwargs...)
end

cm2inch(x) = x/2.54
cm2inch.((8.3,7))
figpath = "tmp"
#figpath = "inst/22paper_upscaling" # interactivly set to override paper figures

i_plot_Plots = () -> begin
    #depr: now using Makie instead of Plots
    #using Plots
    pla = (
        palette = :Dark2_6,
        dpi = dpi_pub,
        size = inch2px.(cm2inch.((8.3,7)), dpi_pub),
    )
    ts = (-5.0,min(5.0,maximum(sol_sesam3f.t)))
    plot(sol_seam3f; tspan=ts, vars=[pl.i_L, pl.i_L_annual], xlab="time (yr)",
        ylab="litter input (g/m2/yr)", pla...)
    savefig(joinpath(figpath,"fluct_litterinput.pdf"))

    plot_vars([s.E_R]; ylab="enzyme pool E_R (g/m2)", tspan=(-2,0), legend=:topleft, pla...)
    savefig(joinpath(figpath,"fluct_E_R.pdf"))

    plot_vars([s.leach_N]; ylab="N leaching (g/m2/yr)", tspan=ts, legend=:topleft, pla...)
    savefig(joinpath(figpath,"fluct_Nleach.pdf"))
end

sol = sol_sesam3f
ts = (-5.0,min(5.0,maximum(sol.t)))
tsi = (1:length(sol.t))[first(ts) .<= sol.t .<= last(ts)]

# vars=[pl.i_L, pl.i_L_annual]
using CairoMakie, AlgebraOfGraphics
set_aog_theme!()
legend_theme = Theme(
    Legend=(
        bgcolor=(:white,0.8),
        framevisible=true,
        rowgap=0.2,
        framewidth=0.2,
        framecolor=(:lightgrey,0.8),
        padding=(5.0,5.0,2.0,2.0),
        patchsize=(20.0,10.0),
        margin=(2,2,2,2),
))
update_theme!(legend_theme)
update_theme!(
    figure_padding = 1,
    Lines = (linewidth = 0.8,),
    Series = (linewidth=0.8,),
    )
axis_theme = Theme(
    Axis=(
        #xlabelsize=20,
        xgridvisible=true,
        ygridvisible=true,
        #xgridstyle=:dash, ygridstyle=:dash,
        xtickalign=1.0, ytickalign=1.0, yticksize=3, xticksize=3,
        xticklabelpad=0,
        xlabelpadding=0,
        #ylabelpadding=0,
        #leftspinevisible = false,
        rightspinevisible = true,
        # bottomspinevisible = false,
        topspinevisible = true,
        topspinecolor = :darkgray,
        rightspinecolor = :darkgray,
    )
)
update_theme!(axis_theme)

#Makie.current_default_theme().palette.color[]
#Makie.current_default_theme().Lines
#Makie.current_default_theme().Legend
#Makie.current_default_theme().Axis

function plotm_vars!(ax, vars, tspan=tspan; variants = variants, legend_position=:rr, kwargs...)
    #row = first(Tables.namedtupleiterator(select(variants, :sol, :label, :linestyle, :color)));
    for row in Tables.namedtupleiterator(select(variants, :sol, :label, :linestyle, :color))
        #@show typeof(row)
        series_sol!(ax, row.sol, vars; tspan, labels=[row.label],
        linestyle=row.linestyle,
        solid_color=row.color,
        )
    end
    axislegend(ax, unique=true, position = legend_position)
    display(fig)
    fig
end
#save("tmp.pdf", fig, pt_per_unit = 1)
float_to_month = x -> ["Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"][Int32(round(mod(x,1)*12)+1)]

# include("cairo_makie_util.jl") # moved to MTKHelpers

size7 = cm2inch.((8.3,7.0))

fig, ax = pdf_figure_axis(size7, xlabel = "Time", ylabel="enzyme pool E_R (g/m2)");
plotm_vars!(ax, [s.E_R], (-2,0); variants = variants[[1,2],:], legend_position=:lt)
ax.xticks = [-2.0, -1.5, -1, -0.5, 0.0]
ax.xtickformat = xs -> [float_to_month(x) for x in xs]
display(fig)
save(joinpath(figpath,"fluct_E_R.pdf"), fig, pt_per_unit = 1)

fig, ax = pdf_figure_axis(size7, xlabel = "Time (yr)", ylabel="litter input (g/m2/yr)");
ts = (-5.0,min(5.0,maximum(sol_seam3f.t)))
series_sol!(ax, sol_seam3f, [pl.i_L, pl.i_L_annual], tspan=ts, labels=["seasonal","annual"], linewidth=0.8)
axislegend(ax, unique=true, valign = :top, halign=:left, margin=(2,2,2,2))
display(fig)
save(joinpath(figpath,"fluct_litterinput.pdf"), fig, pt_per_unit = 1)

fig, ax = pdf_figure_axis(size7, xlabel = "Time (yr)", ylabel="N leaching (g/m2/yr)");
plotm_vars!(ax, [s.leach_N], (-5,5); variants = variants[[1,3,4],:], legend_position=:lt)
#plotm_vars!(ax, [s.leach_N], (-5,15); variants = variants[[1,3,4],:], legend_position=:lt)
save(joinpath(figpath,"fluct_Nleach.pdf"), fig, pt_per_unit = 1)

# fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="N limitation (g/g)");
# plotm_vars!(ax, [s.lim_N], (-5,15); variants = variants[[1,3,4],:], legend_position=:lt)
# #display(fig)
# save(joinpath(figpath,"fluct_Nlim.pdf"), fig, pt_per_unit = 1)

fig, ax = pdf_figure_axis(size7, xlabel = "Time", ylabel="normalized value");
ts = (-2.3, -1.3)
cols = Makie.current_default_theme().palette.color[];
max_i_Lagr = maximum(sol_pl[pl.i_Lagr][sol_pl.t .< 0])
series_sol!(ax, sol_pl, [pl.i_Lagr/max_i_Lagr], tspan=ts, labels=["i_Lagr"], linewidth=0.8, solid_color=cols[4])
max_i_L = maximum(sol_sesam3f[pl.i_L][sol_sesam3f.t .< 0])
# series_sol!(ax, sol_sesam3f, [pl.i_L/max_i_L], tspan=ts, labels=["normalized soil litter input"], linewidth=0.8, solid_color=cols[2])
max_L = maximum(sol_sesam3f[s.L][sol_sesam3f.t .< 0])
max_E_L = maximum(sol_sesam3f[s.E_L][sol_sesam3f.t .< 0])
min_alpha, max_alpha = extrema(sol_sesam3f[s.α_L][ts[1] .<= sol_sesam3f.t .<= ts[2]])
r_alpha = max_alpha - min_alpha
#
series_sol!(ax, sol_sesam3f, [pl.i_L/max_i_L, s.L/max_L, (s.α_L - min_alpha+ r_alpha/2)/r_alpha/2, s.E_L/max_E_L], tspan=ts, labels=["i_L","L","α_L","E_L"], linewidth=0.8)
axislegend(ax, unique=true, valign = :bottom, halign=:right, margin=(2,2,2,2))
#CairoMakie.ylims!(ax, 0, 1.7) # separate
ax.xticks = collect(-26:2:-14) ./ 12
ax.xtickformat = xs -> [float_to_month(x) for x in xs]
display(fig)
save(joinpath(figpath,"fluct_litter_delay_composition.pdf"), fig, pt_per_unit = 1)

i_plot = () -> begin
    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="R (g/m2)");
    plotm_vars!(ax, [s.R], (first(tspan),0); variants = variants[[1,2,3,4],:], legend_position=:lt)
    # solved: after 200 years still increasing?
    #    slightly higher annual integrated litter input -> avoid steep slopes + reltol in solver

    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="R (g/m2)");
    plotm_vars!(ax, [s.R], (-5,0); variants = variants[[1,2,3,4],:], legend_position=:lt)

    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="I_N (g/m2)");
    plotm_vars!(ax, [s.I_N], (-5,0); variants = variants[[1,2,3,4],:], legend_position=:lt)
    # consistently higher than with constant litter input

    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="α_L (1/1)");
    plotm_vars!(ax, [s.α_L], (-1,0); variants = variants[[1,2,3],:], legend_position=:lb)
    # annual cycle: after winter shiftring towards R and in autum rapidly shifting towards L

    # plant uptake matches maximum
    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="plant uptake (g/m2/yr)");
    series_sol!(ax, sol_sesam3f, [s.u_PlantNmax, s.u_PlantN], tspan=(-2,2), linewidth=0.8)
    axislegend(ax, unique=true, valign = :top, halign=:left, margin=(2,2,2,2))
    display(fig)

    # also the same N pool
    fig, ax = pdf_figure_axis(xlabel = "Time (yr)", ylabel="CN of R pool (g/g)");
    plotm_vars!(ax, [s.R/s.R_N], (-1,0); variants = variants[[2,4],:], legend_position=:lb)

    # using Plots
    plot(sol, vars=[s.R])
    ts = tspan
    ts = (-5.0,min(5.0,maximum(sol_sesam3f.t)))
    plot(sol_seam3f; tspan=ts, vars=[pl.i_L, pl.i_L_annual], xlab="time (yr)",
        ylab="litter input (g/m2/yr)", pla...)

    #plot(sol, tspan=ts, vars=[pl.Lagr])

    plot_vars([s.R]; ylab="residue pool R (g/m2)", tspan=ts, pla...)
    #plot_vars([s.dR]; ylab="change residue pool dR (g/m2)", tspan=ts)
    plot_vars([s.L]; ylab="litter pool L (g/m2)", tspan=ts, pla...)
    plot_vars([s.L+s.R]; ylab="SOM stocks (L+R) (g/m2)", tspan=ts, pla...)
    plot_vars([s.r_tot]; ylab="respiration (g/m2/yr)", tspan=ts, pla...)
    plot_vars([s.leach_N]; ylab="N leaching (g/m2/yr)", tspan=ts, legend=:bottomright, pla...)
    plot_vars([s.α_R]; ylab="enzyme allocation to E_R (g/g)", tspan=ts, pla...)
    plot_vars([s.E_R]; ylab="residue enzyme pool E_R (g/m2)", tspan=ts, pla...)

    plot(sol_sesam3f, vars = [pl.i_L, pl.i_L_annual]; ylab="litter input (g/m2/yr)",
        tspan=ts, xlab="time (yr)")

    ts = (-2,0)
    # almost no difference between steady enzyme and explicit enzyme
    plot_vars([s.E_R]; ylab="residue enzyme pool E_R (g/m2)", tspan = ts, ylim=(0.065,0.075))

    ts = (-2,5)
    ts = (-1,1)
    ts = (0,10) .+ tspan[1]
    plot(sol, tspan=ts, vars=[pl.i_L, pl.i_L_annual])
    plot(sol, tspan=ts, vars=[pl.i_Lagr, pl.dec_Lagr, pl.d_Lagr])

    x = (9.5/12:1/(12*10):11.5/12) .+ 1
    scatter(x, first.(sol.(x)))
    plot(sol, tspan=extrema(x), vars=[pl.i_Lagr, pl.dec_Lagr])

    plot(sol, vars=[pl.i_L / pl.β_Ni])
end


text:
### Create a list of flights based on the Flight id, this output should include the passenger Id,
# relevant IATA/FAA codes, the departure time, the arrival time
# (times to be converted toHH:MM:SS format), and the flight times.
reload("nmr")

nmr.runjob(nmr.NMR(2, ["AComp_Passenger_data.csv"], nmr.mapper_parserecordacomp, "acompdata.json"))
c = nmr.runcombiner("acompdata.json", nmr.combiner_parsejson)

@test typeof(c) == Vector{String}
@test length(c) == 264
@test JSON.parse(c[1])["totalflighttime"] == 1049

ext: jl

max_stars_repo_name: deyandyankov/nmr

max_issues_repo_name: deyandyankov/nmr

max_forks_repo_name: deyandyankov/nmr


text:
vertices(array :: AbstractArray{Bool}) = sum(array)

function ortvec(N, i)
    vec = zeros(Int, N)
    vec[i] = 1

    return CartesianIndex(vec...)
end

function segments(array :: AbstractArray{Bool, N}) where N
    δ = [ortvec(N, i) for i in 1:N]

    indices    = CartesianIndices(array)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(fidx)

    mapreduce(+, indices) do idx
        slice = array[idx:min(idx + uidx, lidx)]

        mapreduce(+, 1:N) do i
            idx2 = uidx + δ[i]
            (checkbounds(Bool, slice, idx2) && slice[uidx] && slice[idx2]) ? 1 : 0
        end
    end
end

function checkface(array :: AbstractArray{Bool, N},
                   δ1    :: CartesianIndex{N},
                   δ2    :: CartesianIndex{N}) where N
    uidx = array |> CartesianIndices |> first |> oneunit

    idx1 = uidx
    idx2 = uidx + δ1
    idx3 = uidx + δ2
    idx4 = uidx + δ1 + δ2

    return checkbounds(Bool, array, idx2) &&
        checkbounds(Bool, array, idx3) &&
        checkbounds(Bool, array, idx4) &&
        array[idx1] && array[idx2] && array[idx3] && array[idx4]
end

function faces(array :: AbstractArray{Bool, 3})
    indices    = CartesianIndices(array)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(fidx)

    mapreduce(+, indices) do idx
        slice = array[idx:min(idx + uidx, lidx)]

        checkface(slice, CartesianIndex(1, 0, 0), CartesianIndex(0, 1, 0)) +
            checkface(slice, CartesianIndex(1, 0, 0), CartesianIndex(0, 0, 1)) +
            checkface(slice, CartesianIndex(0, 1, 0), CartesianIndex(0, 0, 1))
    end
end

function volumes(array :: AbstractArray{Bool})
    indices    = CartesianIndices(array)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(fidx)

    mapreduce(+, fidx:(lidx - uidx)) do idx
        slice = array[idx:(idx + uidx)]
        all(isone, slice)
    end
end

"""
    euler_characteristic(array :: AbstractArray{Bool})

Calculate Euler characteristic for a set of points `S` on a regular
cubic or square grid. The points are defined by a binary array `array`
in such a manner that if `array[idx] == true` then `idx ∈ S`.
"""
function euler_characteristic end

euler_characteristic(array :: AbstractArray{Bool, 2}) =
    vertices(array) - segments(array) + volumes(array)

euler_characteristic(array :: AbstractArray{Bool, 3}) =
    vertices(array) - segments(array) + faces(array) - volumes(array)


text:
using Test
using GraphPlot # For plotting tests
using Neo4j
using DistributedFactorGraphs
using Pkg

## To run the IIF tests, you need a local Neo4j with user/pass neo4j:test
# To run a Docker image
# Install: docker pull neo4j
# Run: docker run -d --publish=7474:7474 --publish=7687:7687 --env NEO4J_AUTH=neo4j/test neo4j
##

# If you want to enable debugging logging (very verbose!)
# logger = SimpleLogger(stdout, Logging.Debug)
# global_logger(logger)

# Test each interface
# Still test LightDFG and MetaGraphsDFG for the moment until we remove in 0.4.2
apis = [
    GraphsDFG,
    DistributedFactorGraphs.MetaGraphsDFG,
    DistributedFactorGraphs.SymbolDFG,
    LightDFG]
for api in apis
    @testset "Testing Driver: $(api)" begin
        @info "Testing Driver: $(api)"
        global testDFGAPI = api
        include("interfaceTests.jl")
    end
end

# Test that we don't export LightDFG and MetaGraphsDFG
@testset "Deprecated Drivers Test" begin
    @test_throws UndefVarError SymbolDFG{NoSolverParams}()
    @test_throws UndefVarError MetaGraphsDFG{NoSolverParams}()
end

# Test special cases

@testset "Plotting Tests" begin
    include("plottingTest.jl")
end

@testset "LightDFG subtype tests" begin
    for type in [(var=DFGVariableSummary, fac=DFGFactorSummary), (var=SkeletonDFGVariable,fac=SkeletonDFGFactor)]
        @testset "$(type.var) and $(type.fac) tests" begin
            @info "Testing $(type.var) and $(type.fac)"
            global VARTYPE = type.var
            global FACTYPE = type.fac
            include("LightDFGSummaryTypes.jl")
        end
    end
end

if get(ENV, "IIF_TEST", "") == "true"

    Pkg.add("IncrementalInference")
    # TODO: Remove this once we move to v0.5.0
    Pkg.add(PackageSpec(name="IncrementalInference", rev="enhancement/compare_move_dfg"))
    @info "------------------------------------------------------------------------"
    @info "These tests are using IncrementalInference to do additional driver tests"
    @info "------------------------------------------------------------------------"

    using IncrementalInference

    apis = [
        GraphsDFG{NoSolverParams}(),
        LightDFG{NoSolverParams}(),
        # DistributedFactorGraphs.MetaGraphsDFG{NoSolverParams}(),
        # DistributedFactorGraphs.SymbolDFG{NoSolverParams}(),
        CloudGraphsDFG{SolverParams}("localhost", 7474, "neo4j", "test",
                                    "testUser", "testRobot", "testSession",
                                    nothing,
                                    nothing,
                                    IncrementalInference.decodePackedType,
                                    IncrementalInference.rebuildFactorMetadata!,
                                    solverParams=SolverParams())
            ]
    for api in apis
        @testset "Testing Driver: $(typeof(api))" begin
            @info "Testing Driver: $(api)"
            global dfg = deepcopy(api)
            include("iifInterfaceTests.jl")
        end

        @testset "FileDFG Testing Driver: $(typeof(api))" begin
            @info "FileDFG Testing Driver: $(typeof(api))"
            global dfg = deepcopy(api)
            include("fileDFGTests.jl")
        end
    end

    @testset "CGStructure Tests for CGDFG" begin
        # Run the CGStructure tests
        include("CGStructureTests.jl")
    end
else
    @warn "Skipping IncrementalInference driver tests"
end


text:
using ADCME
using PyPlot
using PyCall
mpl = pyimport("tikzplotlib")
db = Database("noise.db")
res = execute(db, """
select noise, kappa from tau
""")|>collect

noise = [x[1] for x in res]
noise = round.(noise, sigdigits = 2)
kappa = [map(z->parse(Float64, z), split(x[2],',')) for x in res]

close("all")
for k = 1:length(noise)
    plot(1:100:40000, kappa[k][1:100:40000], label = "\$\\sigma_n = $(noise[k])\$")
end
hlines(0.06, 0, 40000, linestyle = "--", color = "k")
legend()
xlabel("Iterations")
ylabel("\$\\tau\$")
grid("on")
minorticks_on()
mpl.save("with_noise.tex")
savefig("with_noise.png")


text:
struct AggressivenessBeliefMDP <: MDP{AggressivenessBelief, MLAction}
    up::AggressivenessUpdater
end

function generate_sr(p::AggressivenessBeliefMDP, b_old::AggressivenessBelief, a::MLAction, rng::AbstractRNG)
    up = p.up
    pomdp = get(up.problem)
    s = rand(rng, b_old)
    sp = generate_s(pomdp, s, a, rng)
    o = generate_o(pomdp, s, a, sp)

    b_new = AggressivenessBelief(CorrelatedIDMMOBIL(
                                 get(up.problem).dmodel.behaviors), o,
                                 Vector{Vector{Float64}}(length(o.cars)),
                                 Vector{Vector{Float64}}(length(o.cars)))

    rsum = 0.0
    particles = Vector{MLState}(up.nb_sims)
    stds = max.(agg_stds(b_old), 0.01)
    for i in 1:up.nb_sims
        if rand(up.rng) < up.p_resample_noise
            s = rand(up.rng, b_old, up.resample_noise_factor*stds)
        else
            s = rand(up.rng, b_old)
        end
        particles[i], r = generate_sr(get(up.problem), s, a, up.rng)
        rsum += r
    end

    cweights_from_particles!(b_new, get(up.problem), o, particles, up.params)

    for i in 1:length(o.cars)
        if isempty(b_new.cweights[i])
            # println("car $i has empty weights")
            b_new.particles[i] = rand(up.rng, up.nb_sims)
            b_new.cweights[i] = 1.0:1.0:up.nb_sims
        end
    end

    return b_new, rsum/up.nb_sims
end

# actions(p::AggressivenessBeliefMDP) = actions(get(p.up.problem))
actions(p::AggressivenessBeliefMDP, b::AggressivenessBelief) = actions(get(p.up.problem), b.physical)
discount(p::AggressivenessBeliefMDP) = discount(get(p.up.problem))

struct ABMDPSolver <: Solver
    solver
    updater
end

function solve(sol::ABMDPSolver, pomdp)
    up = deepcopy(sol.updater)
    set_problem!(up, pomdp)
    return solve(sol.solver, AggressivenessBeliefMDP(up))
end

struct BehaviorBeliefMDP{G} <: MDP{BehaviorParticleBelief{G}, MLAction}
    up::BehaviorParticleUpdater
end

BehaviorBeliefMDP(up) = BehaviorBeliefMDP{typeof(get(up.problem).dmodel.behaviors)}(up)

function generate_sr(p::BehaviorBeliefMDP, b_old::BehaviorParticleBelief, a::MLAction, rng::AbstractRNG)
    up = p.up
    pomdp = get(up.problem)
    s = rand(rng, b_old)
    sp = generate_s(pomdp, s, a, rng)
    o = generate_o(pomdp, s, a, sp)

    b_new::BehaviorParticleBelief=BehaviorParticleBelief(get(up.problem).dmodel.behaviors, o,
                                        Vector{Vector{IDMMOBILBehavior}}(length(o.cars)),
                                        Vector{Vector{Float64}}(length(o.cars)))

    rsum = 0.0
    particles = Vector{MLState}(up.nb_sims)
    samples = lv_resample(b_old, up)
    for i in 1:up.nb_sims
        particles[i], r = generate_sr(get(up.problem), samples[i], a, up.rng)
        rsum += r
    end

    cweights_from_particles!(b_new, get(up.problem), o, particles, up.params)

    for i in 1:length(o.cars)
        if isempty(b_new.cweights[i])
            b_new.particles[i] = [rand(up.rng, b_new.gen) for i in 1:up.nb_sims]
            b_new.cweights[i] = collect(1.0:convert(Float64, up.nb_sims))./up.nb_sims
        end
    end

    return b_new, rsum/up.nb_sims
end

actions(p::BehaviorBeliefMDP, b::BehaviorParticleBelief) = actions(get(p.up.problem), b.physical)
discount(p::BehaviorBeliefMDP) = discount(get(p.up.problem))

struct BBMDPSolver <: Solver
    solver
    updater
end

function solve(sol::BBMDPSolver, pomdp)
    up = deepcopy(sol.updater)
    set_problem!(up, pomdp)
    return solve(sol.solver, BehaviorBeliefMDP(up))
end


text:
mutable struct Permute_Tuner
    inst::Permute_Instruct
    velocities::Matrix{Float64}
    successes::BitMatrix #(memoryxfunc)
    tabular_display::DataFrame
    time_history::Vector{Float64}
    override::Bool
end

"""
    Permute_Tuner(instruction)

Generate a Permute_Tuner for the given instruction.
"""
function Permute_Tuner(instruction::Permute_Instruct)
    nfuncs=length(instruction.funcs)
    funcnames=Vector{String}()
    vels=ones(TUNING_MEMORY*instruction.func_limit,nfuncs)
    succs=trues(TUNING_MEMORY*instruction.func_limit,nfuncs)
    for (idx,func) in enumerate(instruction.funcs)
        length(instruction.args[idx])>0 ? (kwstr="(+$(length(instruction.args[idx]))kwa)") : (kwstr="")
        push!(funcnames, string(nameof(func),kwstr))
    end

    tabular_display=DataFrame("Function"=>funcnames, "Succeed"=>zeros(Int64,nfuncs), "Fail"=>zeros(Int64, nfuncs),"Velocity"=>ones(nfuncs), "Weights"=>instruction.weights)

    instruction.override_time>0. ? (override=true) : (override=false)
    return Permute_Tuner(instruction,vels,succs,tabular_display,zeros(CONVERGENCE_MEMORY),override)
end

"""
    tune_weights!(tuner, call_report)

Given a call_report from permute_IPM(), adjust tuner's Permute_Instruct weights for function success rate and likelihood surface velocity.
"""
function tune_weights!(tuner::Permute_Tuner, call_report::Vector{Tuple{Int64,Float64,Float64}})
    for call in call_report
        funcidx,time,distance=call
        distance!==-Inf && (tuner.velocities[:,funcidx]=update_velocity!(tuner.velocities[:,funcidx],time,distance)) #do not push velocity to array if it has -Inf probability (usu no new model found)
        if call===call_report[end]
            tuner.successes[:,funcidx]=update_sucvec!(tuner.successes[:,funcidx],true)
            tuner.tabular_display[funcidx,"Succeed"]+=1
        else
            tuner.successes[:,funcidx]=update_sucvec!(tuner.successes[:,funcidx],false)
            tuner.tabular_display[funcidx,"Fail"]+=1
        end
    end
    tuner.override && mean(tuner.time_history)>tuner.inst.override_time ? (tuner.inst.weights=tuner.inst.override_weights; tuner.tabular_display[!,"Weights"]=tuner.inst.override_weights) :
        update_weights!(tuner)
end

function update_velocity!(velvec,time,distance)
    popfirst!(velvec) #remove first value
    vel = distance - log(time)
    push!(velvec,distance-log(time))
end

function update_sucvec!(sucvec, bool)
    popfirst!(sucvec)
    push!(sucvec, bool)
end

function update_weights!(t::Permute_Tuner)
    mvels=[mean(t.velocities[:,n]) for n in 1:length(t.inst.funcs)]
    t.tabular_display[!,"Velocity"]=copy(mvels)

    any(i->i<0,mvels) && (mvels.+=-minimum(mvels)+1.) #to calculate weights, scale negative values into >1.
    pvec=[mvels[n]*(sum(t.successes[:,n])/length(t.successes[:,n])) for n in 1:length(t.inst.funcs)]
    pvec./=sum(pvec)
    (any(pvec.<t.inst.min_clmps) || any(pvec.>t.inst.max_clmps)) && clamp_pvec!(pvec,t.inst.min_clmps,t.inst.max_clmps)

    @assert isprobvec(pvec)
    t.inst.weights=pvec; t.tabular_display[!,"Weights"]=pvec
end

"""
    clamp_pvec(pvec, tuner)

Clamp the values of a probability vector between the minimums and maximums provided by a Permute_Tuner.
"""
            function clamp_pvec!(pvec, min_clmps, max_clmps)
                #logic- first accumulate on low values, then distribute excess from high values
                any(pvec.<min_clmps) ? (low_clamped=false) : (low_clamped=true)
                while !low_clamped
                    vals_to_accumulate=pvec.<min_clmps
                    vals_to_deplete=pvec.>min_clmps

                    depletion=sum(min_clmps[vals_to_accumulate].-pvec[vals_to_accumulate])
                    pvec[vals_to_accumulate].=min_clmps[vals_to_accumulate]
                    pvec[vals_to_deplete].-=depletion/sum(vals_to_deplete)

                    !any(pvec.<min_clmps)&&(low_clamped=true)
                end

                any(pvec.>max_clmps) ? (high_clamped=false) : (high_clamped=true)
                while !high_clamped
                    vals_to_deplete=pvec.>max_clmps
                    vals_to_accumulate=pvec.<max_clmps

                    depletion=sum(pvec[vals_to_deplete].-max_clmps[vals_to_deplete])
                    pvec[vals_to_deplete].=max_clmps[vals_to_deplete]
                    pvec[vals_to_accumulate].+=depletion/sum(vals_to_accumulate)

                    !any(pvec.>max_clmps)&&(high_clamped=true)
                end
            end

function Base.show(io::IO, tuner::Permute_Tuner; progress=false)
    show(io, tuner.tabular_display, rowlabel=:I, summary=false)
    progress && return(size(tuner.tabular_display,1)+4)
end


text:
using SACBP
using LinearAlgebra

function gradientControlUpdate(model::UKFPosRange,
                               s::AugState{PhysPos{T},VecBelMvNormal{T}},
                               UArray::Vector{<:PosControl},
                               u_param_min::Vector{<:Real},
                               u_param_max::Vector{<:Real},
                               Q::Matrix{T},
                               dto::Float64) where T <: Real
    gradientPolicy = GradientMultiTargetLocalizationPolicy();
    u_val = vec(control_nominal(gradientPolicy,model,s,u_param_min,u_param_max,dto,Q));
    UArray_new = copy(UArray)
    for ii = 1:length(UArray_new)
        if UArray_new[ii].t < round(s.t + dto,digits=5)
            UArray_new[ii] = PosControl(UArray_new[ii].t,u_val);
        end
    end
    return UArray_new
end;
#=
function gradientControlVal(model::UKFPosRange,
                            s::AugState{PhysPos{T},VecBelMvNormal{T}},
                            u_param_min::Vector{<:Real},
                            u_param_max::Vector{<:Real},
                            Q::Matrix{T},
                            dto::Float64) where T <: Real
    gradientPolicy = GradientMultiTargetLocalizationPolicy();
    u_val = vec(control_nominal(gradientPolicy,model,s,u_param_min,u_param_max,dto,Q));
    return u_val
end;
=#


text:
using Images, FileIO

function main(h::Integer, w::Integer, side::Bool=false)
    W0 = w >> 1
    H0 = h >> 1
    @inline function motecolor(x::Integer, y::Integer)
        h = clamp(180 * (atan2(y - H0, x - W0) / π + 1.0), 0.0, 360.0)
        return HSV(h, 0.5, 0.5)
    end

    img = zeros(RGB{N0f8}, h, w)
    img[H0, W0] = RGB(1, 1, 1)
    free = trues(h, w)
    free[H0, W0] = false
    for i in eachindex(img)
        x = rand(1:h)
        y = rand(1:w)
        free[x, y] || continue
        mc = motecolor(x, y)
        for j in 1:1000
            xp = x + rand(-1:1)
            yp = y + rand(-1:1)
            contained = checkbounds(Bool, img, xp, yp)
            if contained && free[xp, yp]
                x, y = xp, yp
                continue
            else
                if side || (contained && !free[xp, yp])
                    free[x, y] = false
                    img[x, y] = mc
                end
                break
            end
        end
    end
    return img
end

imgnoside = main(256, 256)
imgwtside = main(256, 256, true)
save("data/browniantree_noside.jpg", imgnoside)
save("data/browniantree_wtside.jpg", imgwtside)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 0f2b2e38-f2bb-11eb-2f1a-eb1c29869308
begin
	using Pkg
	Pkg.add(url="https://github.com/samurai-kant/Fairness.jl.git")
	Pkg.add("Flux")
	Pkg.add("Downloads")
	Pkg.add("CSV")
	Pkg.add("DataFrames")
	Pkg.add("Statistics")
	Pkg.add("LinearAlgebra")
end

# ╔═╡ 8ca89ec2-f2bb-11eb-3ffe-0130194a8f43
begin
	using Downloads
	using Fairness
	using Flux
	using CSV
	using DataFrames
	using Statistics
	using LinearAlgebra
end

# ╔═╡ 4d38534e-f2bc-11eb-2b66-2f694e3e2410
begin
	const MODULE_DIR = "/Users/agango/Desktop/Datasets"
	const DATA_DIR = joinpath(MODULE_DIR)
	function ensure_download(url::String, file::String)
	    cd(DATA_DIR) # This is to ensue that the dataset is not downloaded to /tmp instead of ./data
	    fpath = joinpath(DATA_DIR, file)
	    if !isfile(fpath)
	        Downloads.download(url, file)
	    end
	end
end

# ╔═╡ be81eefa-f2bb-11eb-368d-775fef1a0a89
begin
	url = "http://archive.ics.uci.edu/ml/machine-learning-databases/communities/communities.data"
	fname = "communities_crime.data"
	cols = ["state", "county", "community", "communityname", "fold", "population", "householdsize", "racepctblack", "racePctWhite", "racePctAsian", "racePctHisp", "agePct12t21", "agePct12t29", "agePct16t24", "agePct65up", "numbUrban", "pctUrban", "medIncome", "pctWWage", "pctWFarmSelf", "pctWInvInc", "pctWSocSec", "pctWPubAsst", "pctWRetire", "medFamInc", "perCapInc", "whitePerCap", "blackPerCap", "indianPerCap", "AsianPerCap", "OtherPerCap", "HispPerCap", "NumUnderPov", "PctPopUnderPov", "PctLess9thGrade", "PctNotHSGrad", "PctBSorMore", "PctUnemployed", "PctEmploy", "PctEmplManu", "PctEmplProfServ", "PctOccupManu", "PctOccupMgmtProf", "MalePctDivorce", "MalePctNevMarr", "FemalePctDiv", "TotalPctDiv", "PersPerFam", "PctFam2Par", "PctKids2Par", "PctYoungKids2Par", "PctTeen2Par", "PctWorkMomYoungKids", "PctWorkMom", "NumIlleg", "PctIlleg", "NumImmig", "PctImmigRecent", "PctImmigRec5", "PctImmigRec8", "PctImmigRec10", "PctRecentImmig", "PctRecImmig5", "PctRecImmig8", "PctRecImmig10", "PctSpeakEnglOnly", "PctNotSpeakEnglWell", "PctLargHouseFam", "PctLargHouseOccup", "PersPerOccupHous", "PersPerOwnOccHous", "PersPerRentOccHous", "PctPersOwnOccup", "PctPersDenseHous", "PctHousLess3BR", "MedNumBR", "HousVacant", "PctHousOccup", "PctHousOwnOcc", "PctVacantBoarded", "PctVacMore6Mos", "MedYrHousBuilt", "PctHousNoPhone", "PctWOFullPlumb", "OwnOccLowQuart", "OwnOccMedVal", "OwnOccHiQuart", "RentLowQ", "RentMedian", "RentHighQ", "MedRent", "MedRentPctHousInc", "MedOwnCostPctInc", "MedOwnCostPctIncNoMtg", "NumInShelters", "NumStreet", "PctForeignBorn", "PctBornSameState", "PctSameHouse85", "PctSameCity85", "PctSameState85", "LemasSwornFT", "LemasSwFTPerPop", "LemasSwFTFieldOps", "LemasSwFTFieldPerPop", "LemasTotalReq", "LemasTotReqPerPop", "PolicReqPerOffic", "PolicPerPop", "RacialMatchCommPol", "PctPolicWhite", "PctPolicBlack", "PctPolicHisp", "PctPolicAsian", "PctPolicMinor", "OfficAssgnDrugUnits", "NumKindsDrugsSeiz", "PolicAveOTWorked", "LandArea", "PopDens", "PctUsePubTrans", "PolicCars", "PolicOperBudg", "LemasPctPolicOnPatr", "LemasGangUnitDeploy", "LemasPctOfficDrugUn", "PolicBudgPerPop", "ViolentCrimesPerPop"]
	ensure_download(url, fname)
	fpath = joinpath(DATA_DIR, fname)
	df = DataFrame(CSV.File(fname, header=cols, silencewarnings=true, delim=","); copycols = false)
	df = dropmissing(df, names(df))
	X1 = df[!, names(df)[1:127]]
	y1 = df[!, names(df)[128]] .> 0.7
	y1= categorical(y1)
	(X1, y1)
end

# ╔═╡ ebaea712-f2bc-11eb-01a9-f1f933afac79
begin
	df_cleaned = df[:,eltype.(eachcol(df)).!=String]
	df_cleaned = select!(df_cleaned, Not([:"state",:"fold"]))
end

# ╔═╡ c7151616-f2bc-11eb-148d-2578d6582076
begin
	crime_rate_70_percentile = quantile!(df[:,128], 0.7)
	race_black_70_percentile = quantile!(df[:,"racepctblack"],0.7)

	crimes_df = df[:,"ViolentCrimesPerPop"] .> crime_rate_70_percentile
	groups_df = df[:,"racepctblack"] .> race_black_70_percentile
end

# ╔═╡ 820c7cf6-f2be-11eb-338e-c9d4758793f4
begin
	df_cleaned.mergekey = zeros(1994)
	df_pos = df_cleaned[crimes_df,:]
	df_neg = df_cleaned[crimes_df.==false,:]
end

# ╔═╡ e6c0f142-f2cb-11eb-2bfb-7943a7231cd7
queries = rand(1:10,
822613)

# ╔═╡ 9a530fc2-f2bf-11eb-0d6f-4d4cf253221e
begin
	pair_features = outerjoin(df_pos, df_neg, on=:"mergekey",makeunique=true)
	pair_features = select!(pair_features, Not([:"mergekey"]))
	pos_groups = pair_features[:,"racepctblack"] .> race_black_70_percentile
	neg_groups = pair_features[:,"racepctblack_1"] .> race_black_70_percentile
	pair_features_pos = Matrix(pair_features[:,1:convert(Int64,size(pair_features)[2]/2)])
	pair_features_neg = Matrix(pair_features[:,convert(Int64,(size(pair_features)[2]/2+1)):convert(Int64,size(pair_features)[2])])
	pair_features_pos = convert(Array{Float64,2}, pair_features_pos)
	pair_features_neg = convert(Array{Float64,2}, pair_features_neg)
end

# ╔═╡ e18a1a3c-f2c1-11eb-3b91-f104bd2247e8
begin
	W = zeros(1,100)
	predict(x) = W*x'
	function loss(x1, x2)
		sum(max.(0,(predict(x1) .- predict(x2))))
	end
	function lossbin(x1, x2)
		sum(predict(x1) .> predict(x2))
	end
end

# ╔═╡ bf04688e-f2cc-11eb-0978-3d287ac97dfa
select_all = ones(
822613)

# ╔═╡ e621bd38-f2d1-11eb-07a5-87ed8cdc286d
.~neg_groups

# ╔═╡ 15e9fac6-f2cd-11eb-3cd5-1d2b9c9349ff
begin
	mask = queries.==5
	mask_pos = mask .* (pos_groups .* .~neg_groups)
	mask_neg = mask .* (.~pos_groups .* neg_groups)
end

# ╔═╡ c4f94664-f301-11eb-2fb8-d9b9b7838e93
ϵ = 0.01

# ╔═╡ b2768480-f34b-11eb-1f9b-fd342bb6269b
M = ones(4,4)*0.25

# ╔═╡ db96504c-f356-11eb-2698-3bc73a8e0e57
M

# ╔═╡ 109ef703-768e-46c7-be37-18f9c7b68601
λ2 = eigvecs(M)[findmax(abs.(eigvals(M)))[2],:]

# ╔═╡ 33f47df7-988a-4e7d-90f8-05a8aff5019f
findmax(λ2)[1]

# ╔═╡ 9cc956ee-9527-48e2-b47c-11f32dcc075f
opt = Descent(0.1)

# ╔═╡ f2944635-f285-4b87-b5ee-26be6d5d035e
λ=randn(4)

# ╔═╡ b717b880-f2c2-11eb-35ab-1b1df15ac2bb
grads = gradient(()->λ[1]*loss(pair_features_pos[mask,:], pair_features_neg[mask,:])+λ[2]*loss(pair_features_pos[mask_pos,:], pair_features_neg[mask_pos,:])+λ[3]*loss(pair_features_pos[mask_neg,:], pair_features_neg[mask_neg,:])+λ[4], Flux.params(W))

# ╔═╡ a88a06c2-21a4-4090-a4c8-e18589376e02
grads[W]

# ╔═╡ aaa1de28-f2fe-11eb-17c4-bfffadf1e873
grads_λ = gradient(()->λ[2]*(lossbin(pair_features_pos[mask_pos,:], pair_features_neg[mask_pos,:])-ϵ)+λ[3]*(lossbin(pair_features_pos[mask_neg,:], pair_features_neg[mask_neg,:])-ϵ)+λ[4], Flux.params(λ))

# ╔═╡ 00583ace-f2c3-11eb-0c44-0915bf9d683b
M_update = grads_λ[λ]*λ'

# ╔═╡ 155049ba-f351-11eb-28b2-47f133b8fae5
M_new = M.^(M_update*0.001)

# ╔═╡ ea3e4730-f354-11eb-3399-dff8b0275932
begin
	M[:,1] = M_new[:,1]/sum(abs.(M_new[:,1]))
	M[:,2] = M_new[:,2]/sum(abs.(M_new[:,2]))
	M[:,3] = M_new[:,3]/sum(abs.(M_new[:,3]))
	M[:,4] = M_new[:,4]/sum(abs.(M_new[:,4]))
end

# ╔═╡ 61697053-7bf9-4a47-8f29-b9367488b123
grads_λ[λ]

# ╔═╡ ecce427a-a17e-4278-80dd-5ffb9f3700e9
eigvecs(M)[:,4]

# ╔═╡ b4cb96ea-ab14-4c07-8c0f-f26fc4fb988b
M*eigvecs(M)[:,4]

# ╔═╡ 84a755fb-671f-41ca-baa7-b9e05acb3104
λ

# ╔═╡ 58b24dd7-6198-4895-8ead-3e888b5d656a
begin
	for j in 1:10
		λ = eigvecs(M)[:,findmax(abs.(eigvals(M)))[2]]
		grads = gradient(()->λ[1]*loss(pair_features_pos[mask,:], pair_features_neg[mask,:])+λ[2]*loss(pair_features_pos[mask_pos,:], pair_features_neg[mask_pos,:])+λ[3]*loss(pair_features_pos[mask_neg,:], pair_features_neg[mask_neg,:])+λ[4], Flux.params(W))
		grads_λ = gradient(()->λ[2]*(lossbin(pair_features_pos[mask_pos,:], pair_features_neg[mask_pos,:])-ϵ)+λ[3]*(lossbin(pair_features_pos[mask_neg,:], pair_features_neg[mask_neg,:])-ϵ)+λ[4], Flux.params(λ))
		Flux.update!(opt, W, grads[W])
		M_update = grads_λ[λ]*λ'
		M_new = M.*(ℯ.^(M_update*0.001))
		M[:,1] = M_new[:,1]/sum(abs.(M_new[:,1]))
		M[:,2] = M_new[:,2]/sum(abs.(M_new[:,2]))
		M[:,3] = M_new[:,3]/sum(abs.(M_new[:,3]))
		M[:,4] = M_new[:,4]/sum(abs.(M_new[:,4]))
	end
end

# ╔═╡ 301fae52-f5c2-4fa2-a719-b300a47b897d
M

# ╔═╡ 2de75817-f6bc-48aa-9ca1-5fa5f25b38d2
W

# ╔═╡ b33a8bc8-a5a0-4f74-9beb-aab808be324c
M

# ╔═╡ e163f362-c0a4-4174-81f7-c4efc701b558
M*eigvecs(M)[:,findmax(abs.(eigvals(M)))[2]]

# ╔═╡ Cell order:
# ╠═0f2b2e38-f2bb-11eb-2f1a-eb1c29869308
# ╠═8ca89ec2-f2bb-11eb-3ffe-0130194a8f43
# ╠═4d38534e-f2bc-11eb-2b66-2f694e3e2410
# ╠═be81eefa-f2bb-11eb-368d-775fef1a0a89
# ╠═ebaea712-f2bc-11eb-01a9-f1f933afac79
# ╠═c7151616-f2bc-11eb-148d-2578d6582076
# ╠═820c7cf6-f2be-11eb-338e-c9d4758793f4
# ╠═e6c0f142-f2cb-11eb-2bfb-7943a7231cd7
# ╠═9a530fc2-f2bf-11eb-0d6f-4d4cf253221e
# ╠═e18a1a3c-f2c1-11eb-3b91-f104bd2247e8
# ╠═bf04688e-f2cc-11eb-0978-3d287ac97dfa
# ╠═e621bd38-f2d1-11eb-07a5-87ed8cdc286d
# ╠═15e9fac6-f2cd-11eb-3cd5-1d2b9c9349ff
# ╠═c4f94664-f301-11eb-2fb8-d9b9b7838e93
# ╠═b2768480-f34b-11eb-1f9b-fd342bb6269b
# ╠═b717b880-f2c2-11eb-35ab-1b1df15ac2bb
# ╠═aaa1de28-f2fe-11eb-17c4-bfffadf1e873
# ╠═00583ace-f2c3-11eb-0c44-0915bf9d683b
# ╠═155049ba-f351-11eb-28b2-47f133b8fae5
# ╠═ea3e4730-f354-11eb-3399-dff8b0275932
# ╠═db96504c-f356-11eb-2698-3bc73a8e0e57
# ╠═109ef703-768e-46c7-be37-18f9c7b68601
# ╠═33f47df7-988a-4e7d-90f8-05a8aff5019f
# ╠═a88a06c2-21a4-4090-a4c8-e18589376e02
# ╠═61697053-7bf9-4a47-8f29-b9367488b123
# ╠═9cc956ee-9527-48e2-b47c-11f32dcc075f
# ╠═f2944635-f285-4b87-b5ee-26be6d5d035e
# ╠═ecce427a-a17e-4278-80dd-5ffb9f3700e9
# ╠═b4cb96ea-ab14-4c07-8c0f-f26fc4fb988b
# ╠═84a755fb-671f-41ca-baa7-b9e05acb3104
# ╠═58b24dd7-6198-4895-8ead-3e888b5d656a
# ╠═301fae52-f5c2-4fa2-a719-b300a47b897d
# ╠═2de75817-f6bc-48aa-9ca1-5fa5f25b38d2
# ╠═b33a8bc8-a5a0-4f74-9beb-aab808be324c
# ╠═e163f362-c0a4-4174-81f7-c4efc701b558


text:
using CairoMakie
using Images, FileIO
using ImageTransformations
using WindMixing
using JLD2
using WindMixing: animate_profiles_fluxes_final

# convective_adjustment = parse(Bool, ARGS[1])
convective_adjustment = false
num = parse(Int, ARGS[1])

if convective_adjustment
    CA_str = "CA"
else
    CA_str = "noCA"
end

u_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "u.png"))))
v_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "v.png"))))
T_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "T.png"))))
uw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "uw.png"))))
vw_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "vw.png"))))
wT_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "wT.png"))))
Ri_img = rotr90(load(assetpath(joinpath(pwd(), "images_plot", "Ri.png"))))
z_img = load(assetpath(joinpath(pwd(), "images_plot", "z.png")))

axis_images = (
    u = u_img,
    v = v_img,
    T = T_img,
    z = z_img,
    uw = uw_img,
    vw = vw_img,
    wT = wT_img,
    Ri = Ri_img,
)

EXTRACTED_FILE_DIR = joinpath(pwd(), "extracted_training_output")

EXTRACTED_FILE_PATH = joinpath(EXTRACTED_FILE_DIR, "NDE_3sim_diurnal_18simBFGST0.8nograd_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8_extracted.jld2")
FILE_DIR = joinpath(pwd(), "final_data_2", "NDE_3sim_diurnal_18simBFGST0.8nograd_divide1f5_gradient_smallNN_leakyrelu_rate_2e-4_T0.8_$(CA_str)")

# if !ispath(FILE_DIR)
#     mkdir(FILE_DIR)
# end

file = jldopen(EXTRACTED_FILE_PATH)
train_files = file["training_info/train_files"]
loss_scalings = file["training_info/loss_scalings"]
train_parameters = file["training_info/parameters"]
close(file)

test_files_all = [
    "wind_-5e-4_diurnal_5e-8",
    "wind_-3.5e-4_diurnal_5e-8",
    "wind_-2e-4_diurnal_5e-8",

    "wind_-4.5e-4_diurnal_4e-8",
    "wind_-4.5e-4_diurnal_2e-8",
    "wind_-3e-4_diurnal_4e-8",
    "wind_-3e-4_diurnal_2e-8",

    "wind_-5.5e-4_diurnal_5.5e-8",
    "wind_-1.5e-4_diurnal_5.5e-8",

    "wind_-5.5e-4_new",

    "wind_-5.5e-4_heating_-3.5e-8",
    "wind_-1.5e-4_heating_-3.5e-8",
    "wind_-5.5e-4_cooling_3.5e-8",
    "wind_-1.5e-4_cooling_3.5e-8",
]

test_files = [test_files_all[num]]

ν₀ = train_parameters["ν₀"]
ν₋ = train_parameters["ν₋"]
ΔRi = train_parameters["ΔRi"]
Riᶜ = train_parameters["Riᶜ"]
Pr = train_parameters["Pr"]

# solve_oceananigans_modified_pacanowski_philander_nn(test_files, EXTRACTED_FILE_PATH, FILE_DIR,
#                                                         timestep=1, convective_adjustment=convective_adjustment)

diurnal = occursin("diurnal", test_files[1])

if num <= 3
    DIR_NAME = "train_$(test_files[1])"
    animation_type = "Training"
elseif num >= 4 && num <= 7
    DIR_NAME = "test_$(test_files[1])"
    animation_type = "Interpolating"
else
    DIR_NAME = "test_$(test_files[1])"
    animation_type = "Extrapolating"
end

n_trainings = length(train_files)

plot_data = NDE_profile_oceananigans(joinpath(FILE_DIR, DIR_NAME), train_files, test_files,
                                  ν₀=ν₀, ν₋=ν₋, ΔRi=ΔRi, Riᶜ=Riᶜ, Pr=Pr,
                                  loss_scalings=loss_scalings,
                                  OUTPUT_PATH=joinpath(FILE_DIR, DIR_NAME, "profiles_fluxes_oceananigans.jld2"))

animate_profiles_fluxes_final(plot_data, axis_images,
                    joinpath(FILE_DIR, DIR_NAME, "$(test_files[1])_60fps"),
                    animation_type=animation_type, n_trainings=n_trainings, training_types="Wind + Diurnal Heating & Cooling", fps=60)


text:
mutable struct DelayMNRMJumpAggregation{T,S,F1,F2,RNG,DG,PQ} <: AbstractDSSAJumpAggregator
    next_jump::Int
    prev_jump::Int
    next_jump_time::T
    end_time::T
    cur_rates::Vector{T}
    sum_rate::T
    ma_jumps::S
    rates::F1
    affects!::F2
    save_positions::Tuple{Bool,Bool}
    rng::RNG
    dep_gr::DG
    pq::PQ
    next_delay::Union{Nothing,Vector{Int}}
    num_next_delay::Union{Nothing,Vector{Int}}
    time_to_next_jump::T
    dt_delay::T
    dep_gr_delay::Union{Nothing,Dict{Int,Vector{Int}}}
end

function DelayMNRMJumpAggregation(nj::Int, njt::T, et::T, crs::Vector{T}, sr::T, maj::S, rs::F1, affs!::F2, sps::Tuple{Bool,Bool}, rng::RNG; num_specs, dep_graph=nothing, kwargs...) where {T,S,F1,F2,RNG}

    # a dependency graph is needed and must be provided if there are constant rate jumps
    if dep_graph === nothing
        if (get_num_majumps(maj) == 0) || !isempty(rs)
            error("To use ConstantRateJumps with the Delay Modified Next Reaction Method (Delay MNRM) algorithm a dependency graph must be supplied.")
        else
            dg = make_dependency_graph(num_specs, maj)
        end
    else
        dg = dep_graph
        # make sure each jump depends on itself
        add_self_dependencies!(dg)
    end

    pq = MutableBinaryMinHeap{T}()

    nd = nothing
    nnd = nothing
    ttnj = zero(et)
    dt_delay = zero(et)
    dg_delay = nothing
    DelayMNRMJumpAggregation{T,S,F1,F2,RNG,typeof(dg),typeof(pq)}(nj, nj, njt, et, crs, sr, maj, rs, affs!, sps, rng, dg, pq, nd, nnd, ttnj, dt_delay, dg_delay)
end

############################# Required Functions ##############################
# creating the JumpAggregation structure (function wrapper-based constant jumps)
function aggregate(aggregator::DelayMNRM, u, p, t, end_time, constant_jumps,
    ma_jumps, save_positions, rng; kwargs...)

    # handle constant jumps using function wrappers
    rates, affects! = get_jump_info_fwrappers(u, p, t, constant_jumps)
    build_jump_aggregation(DelayMNRMJumpAggregation, u, p, t, end_time, ma_jumps,
        rates, affects!, save_positions, rng; num_specs=length(u), kwargs...)
end

# set up a new simulation and calculate the first jump / jump time
function initialize!(p::DelayMNRMJumpAggregation, integrator, u, params, t)
    fill_rates_and_get_times!(p, u, params, t)
    p.dep_gr_delay = dep_gr_delay(p, integrator)
    find_next_delay_dt!(p, integrator)
    generate_jumps!(p, integrator, u, params, t)
    nothing
end

# execute one jump, changing the system state
function execute_jumps!(p::DelayMNRMJumpAggregation, integrator, u, params, t)
    # execute jump
    update_state_delay!(p, integrator, u, t)
    # update current jump rates and times
    update_dependent_rates_delay!(p, integrator, u, params, t)
    nothing
end

# calculate the next jump / jump time
# just the top of the priority queue
function generate_jumps!(p::DelayMNRMJumpAggregation, integrator, u, params, t)
    next_jump_time_reaction, next_jump = top_with_handle(p.pq)
    dt_reaction = next_jump_time_reaction - t
    dt_delay_generation!(p, integrator)
    compare_delay!(p, integrator.de_chan, p.dt_delay, dt_reaction, t)
    if p.next_delay != nothing
        p.next_jump = 0
    else
        p.next_jump = next_jump
    end
    nothing
end

# recalculate jump rates for jumps that depend on the just executed jump (p.next_jump)
function update_dependent_rates_delay!(p::DelayMNRMJumpAggregation, integrator, u, params, t)
    if p.next_delay == nothing  # if next reaction is not delay reaction
        @inbounds dep_rxs = p.dep_gr[p.next_jump]
    else
        # find the dep_rxs w.r.t next_delay vectors
        dep_rxs_ = [p.dep_gr_delay[p.next_delay[i]] for i in eachindex(p.next_delay)]
        dep_rxs = reduce(vcat, dep_rxs_)
    end
    @unpack cur_rates, rates, ma_jumps = p
    num_majumps = get_num_majumps(ma_jumps)
    @inbounds for rx in dep_rxs
        oldrate = cur_rates[rx]

        # update the jump rate
        @inbounds cur_rates[rx] = calculate_jump_rate(ma_jumps, num_majumps, rates, u, params, t, rx)

        # calculate new jump times for dependent jumps
        if rx != p.next_jump && oldrate > zero(oldrate)
            if cur_rates[rx] > zero(eltype(cur_rates))
                update!(p.pq, rx, t + oldrate / cur_rates[rx] * (p.pq[rx] - t))
            else
                update!(p.pq, rx, typemax(t))
            end
        else
            if cur_rates[rx] > zero(eltype(cur_rates))
                update!(p.pq, rx, t + randexp(p.rng) / cur_rates[rx])
            else
                update!(p.pq, rx, typemax(t))
            end
        end
    end
    nothing
end

# reevaulate all rates, recalculate all jump times, and reinit the priority queue
function fill_rates_and_get_times!(p::DelayMNRMJumpAggregation, u, params, t)

    # mass action jumps
    majumps = p.ma_jumps
    cur_rates = p.cur_rates
    # num_de_chan = length(integrator.de_chan) # num_de_chan
    pqdata = Vector{typeof(t)}(undef, length(cur_rates))
    @inbounds for i in 1:get_num_majumps(majumps)
        cur_rates[i] = evalrxrate(u, i, majumps)
        pqdata[i] = t + randexp(p.rng) / cur_rates[i]
    end

    # constant rates
    rates = p.rates
    idx = get_num_majumps(majumps) + 1
    @inbounds for rate in rates
        cur_rates[idx] = rate(u, params, t)
        pqdata[idx] = t + randexp(p.rng) / cur_rates[idx]
        idx += 1
    end

    # setup a new indexed priority queue to storing rx times
    p.pq = MutableBinaryMinHeap(pqdata)
    nothing
end


text:
struct StrictlyGreaterToStrictlyLessBridge{
    T,
    F<:MOI.AbstractScalarFunction,
    G<:MOI.AbstractScalarFunction,
} <: MOIBC.FlipSignBridge{T,Strictly{T,MOI.GreaterThan{T}},Strictly{T,MOI.LessThan{T}},F,G}
    constraint::CI{F,Strictly{T,MOI.LessThan{T}}}
end
function MOIBC.map_set(::Type{<:StrictlyGreaterToStrictlyLessBridge}, set::Strictly{T,MOI.GreaterThan{T}}) where T
    return Strictly(MOI.LessThan(-set.set.lower))
end
function MOIBC.inverse_map_set(::Type{<:StrictlyGreaterToStrictlyLessBridge}, set::Strictly{T,MOI.LessThan{T}}) where T
    return Strictly(MOI.GreaterThan(-set.set.upper))
end
function MOIBC.concrete_bridge_type(
    ::Type{<:StrictlyGreaterToStrictlyLessBridge{T}},
    G::Type{<:MOI.AbstractScalarFunction},
    ::Type{Strictly{T,MOI.GreaterThan{T}}},
) where {T}
    F = MOIU.promote_operation(-, T, G)
    return StrictlyGreaterToStrictlyLessBridge{T,F,G}
end


text:
# https://github.com/Inazuma110/HCBE.jl/issues/7
# Deprecated

"""
  UnionFind{T <: Integer}

  **Deprecated!!**
  In the future, using [DataStructures.jl](https://github.com/JuliaCollections/DataStructures.jl)'s DisjointSet.
"""
mutable struct UnionFind{T <: Integer}
  parent:: Vector{T}  # parent[root] is the negative of the size

  function UnionFind{T}(nodes::T) where T<:Integer
    if nodes <= 0
      throw(ArgumentError("invalid argument for nodes: $nodes"))
    end

    parent = -ones(T, nodes)
    new{T}(parent)
  end
end

UnionFind(nodes::Integer) = UnionFind{typeof(nodes)}(nodes)

function root(uf::UnionFind{T}, x::T)::T where T<:Integer
  if uf.parent[x] < 0
    return x
  else
    # uf.parent[x] = root{T}(uf, uf.parent[x])
    # return uf.parent[x]
    return uf.parent[x] = root(uf, uf.parent[x])
  end
end

function issame(uf::UnionFind{T}, x::T, y::T)::Bool where T<:Integer
  return root(uf, x) == root(uf, y)
end

function size(uf::UnionFind{T}, x::T)::T where T<:Integer
  return -uf.parent[root(uf, x)]
end

function unite!(uf::UnionFind{T}, x::T, y::T)::Bool where T<:Integer
  x = root(uf, x)
  y = root(uf, y)
  if x == y
    return false
  end
  if uf.parent[x] > uf.parent[y]
    x, y = y, x
  end
  # unite smaller tree(y) to bigger one(x)
  uf.parent[x] += uf.parent[y]
  uf.parent[y] = x
  return true
end


text:
function make_time_series_cache(
    ::Type{T},
    component,
    name,
    initial_time,
    len::Int;
    ignore_scaling_factors=true,
) where {T <: PSY.StaticTimeSeries}
    return IS.StaticTimeSeriesCache(
        T,
        component,
        name,
        start_time=initial_time,
        ignore_scaling_factors=ignore_scaling_factors,
    )
end

function make_time_series_cache(
    ::Type{T},
    component,
    name,
    initial_time,
    horizon::Int;
    ignore_scaling_factors=true,
) where {T <: PSY.AbstractDeterministic}
    return IS.ForecastCache(
        T,
        component,
        name,
        start_time=initial_time,
        horizon=horizon,
        ignore_scaling_factors=ignore_scaling_factors,
    )
end

function make_time_series_cache(
    ::Type{PSY.Probabilistic},
    component,
    name,
    initial_time,
    horizon::Int;
    ignore_scaling_factors=true,
)
    return IS.ForecastCache(
        PSY.Probabilistic,
        component,
        name,
        start_time=initial_time,
        horizon=horizon,
        ignore_scaling_factors=ignore_scaling_factors,
    )
end

function get_time_series_values!(
    time_series_type::Type{T},
    model::DecisionModel,
    component,
    name,
    multiplier_id::Int,
    initial_time,
    horizon::Int;
    ignore_scaling_factors=true,
) where {T <: PSY.Forecast}
    if !use_time_series_cache(get_settings(model))
        return IS.get_time_series_values(
            T,
            component,
            name,
            start_time=initial_time,
            len=horizon,
            ignore_scaling_factors=ignore_scaling_factors,
        )
    end

    cache = get_time_series_cache(model)
    key = TimeSeriesCacheKey(IS.get_uuid(component), T, name, multiplier_id)
    if haskey(cache, key)
        ts_cache = cache[key]
    else
        ts_cache = make_time_series_cache(
            time_series_type,
            component,
            name,
            initial_time,
            horizon,
            ignore_scaling_factors=ignore_scaling_factors,
        )
        cache[key] = ts_cache
    end

    ts = IS.get_time_series_array!(ts_cache, initial_time)
    return TimeSeries.values(ts)
end

function get_time_series_values!(
    ::Type{T},
    model::EmulationModel,
    component::U,
    name,
    multiplier_id::Int,
    initial_time,
    len::Int=1;
    ignore_scaling_factors=true,
) where {T <: PSY.StaticTimeSeries, U <: PSY.Component}
    if !use_time_series_cache(get_settings(model))
        return IS.get_time_series_values(
            T,
            component,
            name,
            start_time=initial_time,
            len=len,
            ignore_scaling_factors=ignore_scaling_factors,
        )
    end

    cache = get_time_series_cache(model)
    key = TimeSeriesCacheKey(IS.get_uuid(component), T, name, multiplier_id)
    if haskey(cache, key)
        ts_cache = cache[key]
    else
        ts_cache = make_time_series_cache(
            T,
            component,
            name,
            initial_time,
            len,
            ignore_scaling_factors=ignore_scaling_factors,
        )
        cache[key] = ts_cache
    end

    ts = IS.get_time_series_array!(ts_cache, initial_time)
    return TimeSeries.values(ts)
end


text:
# Improves performance on my machine.
BLAS.set_num_threads(1)

using Ensemble
using CARMA
using HDF5

usage = "julia run_carma.jl DATAFILE P Q [NLIVE]"

if size(ARGS, 1) < 3
    println(usage)
    exit(1)
end

datafile = ARGS[1]
p = parse(Int, ARGS[2])
q = parse(Int, ARGS[3])

nlive = 1024
if size(ARGS, 1) == 4
    nlive = int(ARGS[4])
end

nmcmc = 128

outfile = "state-$(p)-$(q).hdf5"
ckpt_file = "state-$(p)-$(q).ckpt"

data = readdlm(datafile)
data = data[sortperm(data[:,1]),:]

ts = Float64[data[1,1]]
ys = Float64[data[1,2]]
dys = Float64[data[1,3]]

for i in 2:size(data,1)
    t = data[i,1]
    y = data[i,2]
    dy = data[i,3]

    if t == ts[end]
        dy2 = dy*dy
        dys2 = dys[end]*dys[end]

        yy = (y*dys2 + ys[end]*dy2)/(dy2 + dys2)
        dyy = dy*dys[end]/sqrt(dys2 + dy2)

        ys[end] = yy
        dys[end] = dyy
    else
        push!(ts, t)
        push!(ys, y)
        push!(dys, dy)
    end
end

# Fix any zeros in dy---set to minimum positive dy
dys[dys.<=0] = minimum(dys[dys.>0])

post = Kalman.CARMAKalmanPosterior(ts, ys, dys, p, q)

function logl(x)
    Kalman.log_likelihood(post, x)
end
function logp(x)
    Kalman.log_prior(post, x)
end

if ispath(ckpt_file)
    nest_state = h5open(ckpt_file, "r") do f
        EnsembleNest.NestState(f, logl=logl, logp=logp)
    end
else
    nest_state = EnsembleNest.NestState(logl, logp, Kalman.init(post, nlive), nmcmc)
end

EnsembleNest.run!(nest_state, 0.01, ckpt_file=ckpt_file)

h5open(outfile, "w") do f
    ng = g_create(f, "nest_state")
    write(ng, nest_state)

    pg = g_create(f, "carma_data")
    pg["ts", "compress", 3, "shuffle", ()] = ts
    pg["ys", "compress", 3, "shuffle", ()] = ys
    pg["dys", "compress", 3, "shuffle", ()] = dys
    pg["p"] = p
    pg["q"] = q
end

if ispath(ckpt_file)
    rm(ckpt_file)
end


text:
## -------  Ben-Porath

"""
    $(SIGNATURES)

For a single college (not a ModelObject)
"""
mutable struct HcProdFct <: AbstractHcProdFct
    tfp  ::  Double
    timeExp  ::  Double
    hExp  ::  Double
    # Depreciation rate
    deltaH :: Double
    # Ability scale
    aScale  ::  Double
    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double
end

# -----  All colleges

Base.@kwdef mutable struct HcProdSwitches <: AbstractHcProdSwitches
    # Same exponents on time and h. If yes, ignore `hExp`.
    sameExponents :: Bool = true
    timeExp :: Double = 0.6
    calTimeExp :: Bool = true
    hExp :: Double = 0.6
    hExpLb :: Double = -1.0
    calHExp :: Bool = false
    deltaH :: Double = 0.0
    calDeltaH :: Bool = false
    aScale :: Double = 0.2
    calAScale :: Bool = true
end

"""
    $(SIGNATURES)

ModelObject that specifies human capital production functions for all colleges. Ben-Porath functional form.

Since all colleges share some parameters, we need a model object that keeps
track of parameters that are common or differ by college.
"""
mutable struct HcProdFctSet <: AbstractHcProdSet
    objId :: ObjectId
    switches :: HcProdSwitches
    nc :: CollInt

    # Calibrated parameters
    tfpV  ::  IncreasingVector
    timeExp  ::  Double
    hExp  ::  Double
    deltaH :: Double
    # Ability scale
    aScale  ::  Double

    # Fixed time cost per course
    timePerCourse :: Double
    # Study time per course minimum (this is assigned when study time very low)
    minTimePerCourse :: Double

    pvec :: ParamVector
end

## -------------- Ben-Porath: For one college

function make_test_hprod()
    deltaH = 0.05;
    return HcProdFct(1.2, 0.3, 0.45, deltaH, 1.2, 0.02, 0.005)
end

function validate_hprod(hS :: HcProdFct)
    isValid = (tfp(hS) > 0.01)  &&  (0.0 < time_exp(hS) ≤ 1.0);
    return isValid
end

"""
    $(SIGNATURES)

H produced (before shock is realized). Nonnegative.

# Arguments
- nTriedV
    number of courses attempted this period.
- h0V
    h endowments, so that `hV - h0V` is learning.
"""
function dh(hS :: HcProdFct, abilV,  hV, h0V, timeV, nTriedV)
    sTimeV = study_time_per_course(hS, timeV, nTriedV);
    return nTriedV .* (hV .^ h_exp(hS)) .* (sTimeV .^ hS.timeExp) .*
        tfp(hS) .* exp.(hS.aScale .* abilV);
end

function show_string(hS :: HcProdFct)
    tfpStr = round(tfp(hS); digits = 2);
    hExpStr = round(h_exp(hS); digits = 2);
    tExpStr = round(time_exp(hS); digits = 2);
    aScaleStr = round(hS.aScale; digits = 2);
    return "dh = $tfpStr  h ^ $hExpStr  t ^ $tExpStr  exp($aScaleStr a)"
end

function Base.show(io :: IO, hS :: HcProdFct)
    print(io,  "H prod fct:  $(show_string(hS))")
end

## ---------------------  For all colleges

function settings_table(h :: HcProdSwitches)
    same_exponents(h)  ?  expStr = "the same"  :  expStr = "different";
    ddh = delta_h(h);
    cal_delta_h(h)  ?  deprecStr = "fixed at $ddh"  :  deprecStr = "calibrated";
    return [
        "H production function" "Cobb Douglas";
        "Exponents on time and h"  expStr;
        "Depreciation"  deprecStr
    ]
end

function settings_list(h :: HcProdSwitches, st)
    eqnHChange = ["H production function", "eqnHChange",  eqn_hchange(h)];
    return [eqnHChange]
end

function validate_hprod(s :: HcProdSwitches)
    isValid = true;
    if s.sameExponents && s.calHExp
        @warn "Cannot calibrate hExp if it has the same value as timeExp"
        isValid = false;
    end
    return isValid
end

function validate_hprod_set(h :: HcProdFctSet)
    isValid = (h.nc > 1)  &&  (h.timeExp > 0.0)  &&
        (h.aScale > 0.0)  &&  (h.timePerCourse > 0.0);
    isValid = isValid  &&  (1.0 > delta_h(h) >= 0.0);
    return isValid
end

function eqn_hchange(h :: HcProdSwitches)
    "\\hTfp \\sTimePerCourse^{\\hTimeExp}  \\hc^{\\hHExp} e^{\\hAScale \\abil}"
end

# Initialize with defaults
function make_hc_prod_set(objId :: ObjectId, nc :: Integer,
    switches :: HcProdSwitches)

    @assert validate_hprod(switches);

    pTimePerCourse = init_time_per_course();
    pTfp = init_tfp(objId, nc);

    timeExp = switches.timeExp;
    pTimeExp = Param(:timeExp, ldescription(:hTimeExp), lsymbol(:hTimeExp),
        timeExp, timeExp, 0.2, 0.9, switches.calTimeExp);

    pDeltaH = init_h_depreciation(switches);
    if switches.sameExponents
        hExp = switches.timeExp;
        calHExp = false;
    else
        hExp = switches.hExp;
        calHExp = switches.calHExp;
    end
    # Allowing negative hExp for concavity.
    pHExp = Param(:hExp, ldescription(:hHExp), lsymbol(:hHExp),
        hExp, hExp, switches.hExpLb, 0.9, calHExp);

    aScale = switches.aScale;
    pAScale = Param(:aScale, ldescription(:hAScale), lsymbol(:hAScale),
        aScale, aScale, 0.02, 2.0, switches.calAScale);

    pvec = ParamVector(objId,  [pTimeExp, pHExp, pDeltaH, pAScale, pTimePerCourse]);
    # Min study time required per course. Should never bind.
    minTimePerCourse =
        hours_per_week_to_mtu(0.1 / data_to_model_courses(1));

    return HcProdFctSet(objId, switches, nc, pTfp, timeExp, hExp,
        ModelParams.value(pDeltaH), aScale,
        pTimePerCourse.value, minTimePerCourse, pvec)
end

# TFP by college.
function init_tfp(objId :: ObjectId, nc :: Integer;
    tfp0 = 0.2, tfp0Lb = 0.02, tfp0Ub = 4.0)
    # First college
    # tfp0 = 0.2;
    pTfp0 = Param(:x0, ldescription(:hTfpOne), lsymbol(:hTfpOne),
        tfp0, tfp0, tfp0Lb, tfp0Ub, true);

    # Increments for other colleges.
    dTfpV = fill(0.5, nc - 1);
    pTfpGrad = Param(:dxV, ldescription(:hTfpGrad), lsymbol(:hTfpGrad),
        dTfpV, dTfpV, dTfpV .* 0.1, dTfpV .* 10, true);

    ownId = make_child_id(objId, :tfpV);
    pvec = ParamVector(ownId,  [pTfp0, pTfpGrad]);
    return IncreasingVector(ownId, pvec,  tfp0, dTfpV)
end

make_test_hc_prod_set() =
    make_hc_prod_set(ObjectId(:HProd), 4, HcProdSwitches(deltaH = 0.05));

# Make h production function for one college
function make_h_prod(hs :: HcProdFctSet, iCollege :: Integer)
    return HcProdFct(tfp(hs, iCollege), time_exp(hs), h_exp(hs),
        delta_h(hs), hs.aScale,
        hs.timePerCourse,  hs.minTimePerCourse);
end

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


text:
"""
    BadStatesProperty{N<:Real} <: Property

Type that represents a safety property characterized by a set of bad states.
The property is satisfied by a given set of states if the intersection with the
set of bad states is empty.

### Fields

- `bad`     -- convex set representing the bad states
- `witness` -- witness point (empty vector if not set)

### Notes

The following formula characterizes whether a set ``X`` satisfies a safety
property characterized by a set of bad states 𝑃:

```math
    X \\models 𝑃 \\iff X ∩ 𝑃.\\texttt{bad} = ∅
```
"""
mutable struct BadStatesProperty{N<:Real} <: Property
    bad::LazySet
    witness::Vector{N}

    BadStatesProperty{N}(bad::LazySet) where {N<:Real} = new(bad, N[])
end

# type-less convenience constructor
BadStatesProperty(bad::LazySet{N}) where {N<:Real} =
    BadStatesProperty{N}(bad)

"""
    dim(𝑃::BadStatesProperty)::Int

Return the dimension of a property with bad states.

### Input

- `𝑃` -- safety property with bad states

### Output

The dimension of the bad states.
"""
function dim(𝑃::BadStatesProperty)::Int
    return dim(𝑃.bad)
end

"""
    check(𝑃::BadStatesProperty, X::LazySet; witness::Bool=false)

Checks whether a convex set is disjoint from the set of bad states.

### Input

- `𝑃`       -- safety property with bad states
- `X`       -- convex set
- `witness` -- (optional, default: `false`) flag for returning a counterexample
               if the property is violated

### Output

Let ``Y`` be the bad states represented by 𝑃.
* If `witness` option is deactivated: `true` iff ``X ∩ Y = ∅``
* If `witness` option is activated:
  * `(true, [])` iff ``X ∩ Y = ∅``
  * `(false, v)` iff ``X ∩ Y ≠ ∅`` and ``v ∈ X ∩ Y``
"""
@inline function check(𝑃::BadStatesProperty, X::LazySet; witness::Bool=false)
    return isdisjoint(X, 𝑃.bad, witness)
end


text:
"""
    Origin(indices...)
    Origin(origin::Tuple)
    Origin(origin::CartesianIndex)

A helper type to construct OffsetArray with given origin.

The `origin` of an array is defined as the index of its first element, i.e., `first.(axes(A))`.

# Example

```jldoctest; setup=:(using OffsetArrays)
julia> a = [1 2; 3 4];

julia> OffsetArray(a, OffsetArrays.Origin(0, 1))
2×2 OffsetArray(::$(Array{Int,2}), 0:1, 1:2) with eltype $Int with indices 0:1×1:2:
 1  2
 3  4

julia> OffsetArray(a, OffsetArrays.Origin(0)) # short notation for `Origin(0, 0)`
2×2 OffsetArray(::$(Array{Int, 2}), 0:1, 0:1) with eltype $Int with indices 0:1×0:1:
 1  2
 3  4
```
"""
struct Origin{T<:Union{Tuple{Vararg{Int}}, Int}}
    index::T
end
Origin(I::Tuple{Vararg{Int}}) = Origin{typeof(I)}(I)
Origin(I::Tuple{Vararg{Number}}) = Origin(map(Int, I))
Origin(I::CartesianIndex) = Origin(Tuple(I))
Origin(I::Number...) = Origin(I)
# Origin(0) != Origin((0, )) but they work the same with broadcasting
Origin(n::Number) = Origin{Int}(Int(n))

(o::Origin)(A::AbstractArray) = o.index .- first.(axes(A))

Base.Broadcast.broadcastable(o::Origin) = Ref(o)


text:
# copied and tweaked from https://github.com/staticfloat/LLVMBuilder
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
using BinaryBuilder

# Collection of sources required to build LLVM
llvm_ver = "6.0.1"
sources = [
    "http://releases.llvm.org/$(llvm_ver)/llvm-$(llvm_ver).src.tar.xz" =>
    "b6d6c324f9c71494c0ccaf3dac1f16236d970002b42bb24a6c9e1634f7d0f4e2",
    "http://releases.llvm.org/$(llvm_ver)/cfe-$(llvm_ver).src.tar.xz" =>
    "7c243f1485bddfdfedada3cd402ff4792ea82362ff91fbdac2dae67c6026b667",
    "http://releases.llvm.org/$(llvm_ver)/compiler-rt-$(llvm_ver).src.tar.xz" =>
    "f4cd1e15e7d5cb708f9931d4844524e4904867240c306b06a4287b22ac1c99b9",
    "http://releases.llvm.org/$(llvm_ver)/libcxx-$(llvm_ver).src.tar.xz" =>
    "7654fbc810a03860e6f01a54c2297a0b9efb04c0b9aa0409251d9bdb3726fc67",
    "http://releases.llvm.org/$(llvm_ver)/libcxxabi-$(llvm_ver).src.tar.xz" =>
    "209f2ec244a8945c891f722e9eda7c54a5a7048401abd62c62199f3064db385f",
    "http://releases.llvm.org/$(llvm_ver)/polly-$(llvm_ver).src.tar.xz" =>
    "e7765fdf6c8c102b9996dbb46e8b3abc41396032ae2315550610cf5a1ecf4ecc",
    "http://releases.llvm.org/$(llvm_ver)/libunwind-$(llvm_ver).src.tar.xz" =>
    "a8186c76a16298a0b7b051004d0162032b9b111b857fbd939d71b0930fd91b96",
    "http://releases.llvm.org/$(llvm_ver)/lld-$(llvm_ver).src.tar.xz" =>
    "e706745806921cea5c45700e13ebe16d834b5e3c0b7ad83bf6da1f28b0634e11",
    "patches",
    # Julia source
    "JuliaSource",
    # Add pre-built Julia tarballs
    "JuliaBinary",
    # Add pre-built LLVM tarballs
    "LLVM",
    # Add libcxxffi source
    "libcxxffi"
]

# Since we kind of do this LLVM setup twice, this is the shared setup start:
script_setup = raw"""
# We want to exit the program if errors occur.
set -o errexit

cd $WORKSPACE/srcdir/

# First, move our other projects into llvm/projects
for f in *.src; do
    # Don't symlink llvm itself into llvm/projects...
    if [[ ${f} == llvm-*.src ]]; then
        continue
    fi

    # clang lives in tools/clang and not projects/cfe
    if [[ ${f} == cfe-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/clang
    elif [[ ${f} == polly-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/polly
    elif [[ ${f} == lld-*.src ]]; then
        mv $(pwd)/${f} $(echo llvm-*.src)/tools/lld
    else
        mv $(pwd)/${f} $(echo llvm-*.src)/projects/${f%-*}
    fi
done

# Next, boogie on down to llvm town
cd llvm-*.src

# Update config.guess/config.sub stuff
# update_configure_scripts

# Apply all our patches
for f in $WORKSPACE/srcdir/llvm_patches/*.patch; do
    echo "Applying patch ${f}"
    patch -p1 < ${f}
done
"""

# Next, we will Bash recipe for building across all platforms
script = script_setup * raw"""
# build libcxxffi
cd $WORKSPACE/srcdir/

mkdir LLVMBinary
tar -C LLVMBinary --strip-components=1 -xf *LLVM*.tar.gz

cd $WORKSPACE/srcdir/
mkdir build && cd build
CMAKE_FLAGS="-DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain"
CMAKE_FLAGS="${CMAKE_FLAGS} -DJULIA_SOURCE_PREFIX=$WORKSPACE/srcdir/julia"
CMAKE_FLAGS="${CMAKE_FLAGS} -DJULIA_BINARY_PREFIX=$WORKSPACE/srcdir/juliabin"
CMAKE_FLAGS="${CMAKE_FLAGS} -DLLVMBUILDER_PREFIX=$WORKSPACE/srcdir/LLVMBinary"
cmake .. ${CMAKE_FLAGS}
make -j${nproc} VERBOSE=1
make install VERBOSE=1

mkdir -p ${prefix}/src/clang-6.0.1/include
mkdir -p ${prefix}/src/llvm-6.0.1/include
mkdir -p ${prefix}/build/clang-6.0.1/lib/clang/6.0.1/include
mkdir -p ${prefix}/build/clang-6.0.1/include
mkdir -p ${prefix}/build/llvm-6.0.1/include
cp -r $WORKSPACE/srcdir/llvm-6.0.1.src/tools/clang/include/* ${prefix}/src/clang-6.0.1/include/
cp -r $WORKSPACE/srcdir/llvm-6.0.1.src/include/* ${prefix}/src/llvm-6.0.1/include/
cp -r $WORKSPACE/srcdir/LLVMBinary/include/* ${prefix}/build/llvm-6.0.1/include/
cp -r $WORKSPACE/srcdir/LLVMBinary/lib/* ${prefix}/lib/
cp -r $WORKSPACE/srcdir/LLVMBinary/lib/clang/6.0.1/include/* ${prefix}/build/clang-6.0.1/lib/clang/6.0.1/include/
cp -r ${prefix}/build/llvm-6.0.1/include/* ${prefix}/build/clang-6.0.1/include/

cd $WORKSPACE/srcdir
make -f GenerateConstants.Makefile BASE_LLVM_BIN=$WORKSPACE/srcdir/LLVMBinary BASE_JULIA_BIN=$WORKSPACE/srcdir/juliabin BASE_JULIA_SRC=$WORKSPACE/srcdir/julia LLVM_VERSION=6.0.1
cp $WORKSPACE/srcdir/clang_constants.jl ${prefix}/build/

if [[ ${target} == *mingw32* ]] && [[ ${nbits} == 64 ]]; then
    cp -r $WORKSPACE/srcdir/LLVMBinary/bin/* ${prefix}/bin
    mkdir -p ${prefix}/mingw/include
    mkdir -p ${prefix}/mingw/sys-root/include
    cp -r /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/include/* ${prefix}/mingw/include/
    cp -r /opt/x86_64-w64-mingw32/x86_64-w64-mingw32/sys-root/include/* ${prefix}/mingw/sys-root/include/
fi

if [[ ${target} == *mingw32* ]] && [[ ${nbits} == 32 ]]; then
    cp -r $WORKSPACE/srcdir/LLVMBinary/bin/* ${prefix}/bin
    mkdir -p ${prefix}/mingw/include
    mkdir -p ${prefix}/mingw/sys-root/include
    cp -r /opt/i686-w64-mingw32/i686-w64-mingw32/include/* ${prefix}/mingw/include/
    cp -r /opt/i686-w64-mingw32/i686-w64-mingw32/sys-root/include/* ${prefix}/mingw/sys-root/include/
fi

"""

platforms = [
        # BinaryProvider.Linux(:i686; libc=:glibc, compiler_abi=CompilerABI(:gcc7)),
         BinaryProvider.Linux(:x86_64; libc=:glibc, compiler_abi=CompilerABI(:gcc7)),
        # BinaryProvider.Linux(:aarch64; libc=:glibc, compiler_abi=CompilerABI(:gcc7)),
        # BinaryProvider.Linux(:armv7l; libc=:glibc, compiler_abi=CompilerABI(:gcc7)),
        # BinaryProvider.MacOS(:x86_64; compiler_abi=CompilerABI(:gcc7)),
        # BinaryProvider.Windows(:i686; compiler_abi=CompilerABI(:gcc7)),
        # BinaryProvider.Windows(:x86_64; compiler_abi=CompilerABI(:gcc7))
    ]

# The products that we will ensure are always built
products(prefix) = [
    # libraries
    LibraryProduct(prefix, "libcxxffi",  :libcxxffi)
]

build_tarballs(ARGS, "libcxxffi", v"0.0.0", sources, script, platforms, products, [])

max_stars_repo_name: MaximeBouton/CxxBuilder

max_issues_repo_name: MaximeBouton/CxxBuilder

max_forks_repo_name: MaximeBouton/CxxBuilder


text:
function JuMP.function_string(::Type{MIME("text/plain")}, p::ParameterRef)
    par_name = name(p)
    if !isempty(par_name)
        return par_name
    else
        return "noname_param"
    end
end

function JuMP.function_string(::Type{MIME("text/latex")}, p::ParameterRef)
    par_name = name(p)
    if !isempty(par_name)
        # TODO: This is wrong if parameter name constains extra "]"
        return replace(replace(par_name, "[" => "_{", count = 1), "]" => "}")
    else
        return "nonameParam"
    end
end

function JuMP.function_string(mode, a::ParameterJuMP.PAE, show_constant=true)
    ret = ""
    str_v = JuMP.function_string(mode, a.v, false)
    if str_v != "0"
        ret = str_v
    end
    str_p = JuMP.function_string(mode, a.p, false)
    if str_p != "0"
        if ret == ""
            ret = str_p
        else
            if str_p[1] == '-'
                ret = ret * " - " * str_p[2:end]
            else
                ret = ret * " + " * str_p
            end
        end
    end
    if !JuMP._is_zero_for_printing(a.v.constant) && show_constant
        ret = string(ret, JuMP._sign_string(a.v.constant),
                     JuMP._string_round(abs(a.v.constant)))
    end
    return ret
end

max_stars_repo_name: joaquimg/ParameterJuMP

max_issues_repo_name: joaquimg/ParameterJuMP

max_forks_repo_name: joaquimg/ParameterJuMP


text:
function Restarts.restart(exception::Exception, restarts::Pair...)
    for (cvar, cont) in restarts
        f = cvar[]
        f === nothing && continue
        y = f(exception)::Union{Some,Nothing}
        y === nothing && continue
        return cont(something(y))
    end
    throw(exception)
end

macro define_restart(name::Symbol)
    m = @__MODULE__
    quote
        $m.@contextvar $name::$Any = nothing
    end |> esc
end

# TODO: type check?
Restarts.with(f, restarts::Pair...) = with_context(f, restarts...)


text:
using Supervisors
using Test

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


text:
# This file was generated, do not modify it.

using RDatasets
using DataFrames

boston = dataset("MASS", "Boston");

typeof(boston)

names(boston)

first(boston, 4)

boston.Crim[1:5]

boston[3, 5]

boston[1:5, [:Crim, :Zn]]

boston[1:5, 1:2]

b1 = select(boston, [:Crim, :Zn, :Indus])
first(b1, 2)

b2 = select(boston, Not(:NOx))
first(b2, 2)

select!(b1, Not(:Crim))
first(b1, 2)

using StatsBase
describe(boston, :min, :max, :mean, :median, :std)

foo(x) = sum(abs.(x)) / length(x)
d = describe(boston, :mean, :median, :foo => foo)
first(d, 3)

using Statistics

mat = convert(Matrix, boston)
mat[1:3, 1:3]

boston.Crim_x_Zn = boston.Crim .* boston.Zn;

mao = dataset("gap", "mao")
describe(mao, :nmissing)

std(mao.Age)

std(skipmissing(mao.Age))

iris = dataset("datasets", "iris")
first(iris, 3)

unique(iris.Species)

gdf = groupby(iris, :Species);

subdf_setosa = gdf[1]
describe(subdf_setosa, :min, :mean, :max)

df = DataFrame(a=1:3, b=4:6)
combine(df, :a => sum, nrow)

foo(v) = v[1:2]
combine(df, :a => maximum, :b => foo)

bar(v) = v[end-1:end]
combine(df, :a => foo, :b => bar)

combine(groupby(iris, :Species), :PetalLength => mean)

gdf = groupby(iris, :Species)
combine(gdf, :PetalLength => mean => :MPL, :PetalLength => std => :SPL)

combine(gdf, names(iris, Not(:Species)) .=> std)

names(iris, Not(:Species))


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for graph generators

############################################################################################################################

@testset "Graph Generation" begin
   ls = map(string, 1:10)

   g1 = Graph(10)
   @test nv(g1) == 10
   @test ne(g1) == 0
   @test isempty(edges(g1))

   g2 = Graph(10, 20)
   @test nv(g2) == 10
   @test ne(g2) > 0

   g3 = Graph(10, ls)
   @test nv(g3) == 10
   @test ne(g3) == 0
   @test encode(g3) == ls

   g4 = Graph(10, ls, 20)
   @test nv(g4) == 10
   @test ne(g4) > 0
   @test encode(g4) == ls

   @test emptygraph(10) == g1

   g = randgraph(10, 20)
   @test nv(g) == 10
   @test ne(g) > 0

   g = randgraph(10)
   @test nv(g) == 10
   @test ne(g) > 0

   g = completegraph(10)
   @test nv(g) == 10
   @test ne(g) == 90
end


text:
"""
```
impulse_responses(β, Σ, n_obs_shock, horizon, shock_size = 1;
    method = :cholesky, flip_shocks = false, use_intercept = true,
    frequency_band = (2π/32, 2π/6)) where {S<:Real}
```
computes the impulse responses of a VAR system represented in the form

```
yₜ = Xₜβ + ϵₜ,
```
where `Xₜ` stacks the lags of yₜ (with dimensions n_observables x n_regressors), and

```
ϵₜ ∼ 𝒩 (0, Σ).
```

### Inputs
* `β::AbstractMatrix{S}`: coefficient matrix
* `Σ::AbstractMatrix{S}`: innovations variance-covariance matrix
* `n_obs_shock::Int`: index of the observable corresponding to the orthogonalized shock
    causing the impulse response.
* `shock_size::S`: number of standard deviations of the shock

### Keywords
* `method::Symbol`: type of impulse response to compute. The available options are
    `:cholesky` (default), `:maximum_business_cycle_variance` or `:maxBC`, and
    `:cholesky_long_run` or `:choleskyLR`. See `?cholesky_shock`, `?maxBC_shock`,
    and `?cholesky_long_run_shock`.
* `flip_shocks::Bool`: by default, we compute the impulse responses to a negative shock.
    Set `flip_shocks = true` to obtain a positive shock.
* `use_intercept::Bool`: `impulse_responses` assumes `β` has constant term(s). If there
    are no such terms, then `use_intercept` must be set to `false`.
* `frequency_band::Tuple{S,S}`: See `?maxBC_shock`.

### Outputs
* `Y::AbstractMatrix`: Impulse response matrix with dimensions horizons x n_observables
"""
function impulse_responses(β::AbstractMatrix{S}, Σ::AbstractMatrix{S}, n_obs_shock::Int,
                           horizon::Int, shock_size::S = one(S);
                           method::Symbol = :cholesky,
                           flip_shocks::Bool = false,
                           use_intercept::Bool = true,
                           frequency_band::Tuple{S,S} =
                           (2*π/32, 2*π/6)) where {S<:Real}

    # Compute dimensions
    n = size(β, 2)
    lags = convert(Int, use_intercept ? (size(β, 1) - 1) / n : size(β, 1) / n)

    # Compute impact based on IRF type
    Y = zeros(S, n, lags + horizon)
    Y[:, lags + 1] = if method == :cholesky
        cholesky_shock(Σ, n, n_obs_shock, shock_size;
                       flip_shocks = flip_shocks)
    elseif method == :maximum_business_cycle_variance || method == :maxBC
        maxBC_shock(β, Σ, n, n_obs_shock, shock_size, lags, frequency_band;
                    flip_shocks = flip_shocks)
    elseif method == :choleskyLR || method == :cholesky_long_run
        cholesky_long_run_shock(β, Σ, n_obs_shock, n, lags, shock_size;
                                flip_shocks = flip_shocks)
    else
        error("IRF method $(string(method)) has not been implemented.")
    end

    # For efficiency
    if use_intercept
        β = @views β[2:end, :]
    end

    # Compute impulse response
    for t = 2:horizon
        xT = vec(Y[:, lags + t - 1:-1:lags + t - lags])
        Y[:, lags + t] = vec(xT' * β)
    end

    return Y[:, lags + 1:end]
end

"""
```
cholesky_shock(Σ, n, n_obs_shock, shock_size, flip_shocks = false) where {S<:Real}
```
computes a Cholesky-identified shock to the specified observable.

Consider a VAR
```
yₜ = Xₜ β + uₜ,
```
where `Xₜ` stacks the `p` lags of `yₜ` and `uₜ ∼ 𝒩 (0, Ω)`.
We identify orthogonal shocks `ϵₜ` with identity covariance by computing
```
uₜ = cholesky(Ω).L * ϵₜ.
```

### Inputs
* `Σ::AbstractMatrix{S}`: innovations variance-covariance matrix
* `n::Int`: number of observables
* `n_obs_shock::Int`: index of the observable corresponding to the orthogonalized shock
    causing the impulse response.
* `shock_size::S`: number of standard deviations of the shock

### Keywords
* `flip_shocks::Bool`: by default, we compute the impulse responses to a negative shock.
    Set `flip_shocks = true` to obtain a positive shock.
"""
function cholesky_shock(Σ::Matrix{S}, n::Int, n_obs_shock::Int,
                        shock_size::S; flip_shocks::Bool = false) where {S<:Real}
    cholmat = cholesky((Σ + Σ') ./ 2).L
    vec_shock = zeros(n)
    vec_shock[n_obs_shock] = flip_shocks ? shock_size : -shock_size # negative by DSGE convention
    return (cholmat * vec_shock)'
end

"""
```
maxBC_shock(β, Σ, n, n_obs_shock, shock_size, lags, frequency_band,
    flip_shocks = false) where {S<:Real}
```
computes the shock which maximizes the cyclical variance
explained for a chosen observable.
The observable's index is specified by `n_obs_shock`, and the
frequencies of the cycle are specified by `frequency_band`.

A typical business cycle's frequencies are given by (2π / 32, 2π / 6).

### Inputs
* `β::AbstractMatrix{S}`: coefficient matrix
* `Σ::AbstractMatrix{S}`: innovations variance-covariance matrix
* `n::Int`: number of observables
* `n_obs_shock::Int`: index of the observable for which we compute
    the shock that maximizes the observable's explained variance
* `shock_size::S`: number of standard deviations of the shock
* `lags::Int`: number of lags in VAR system
* `frequency_band::Tuple{S,S}`: the frequencies between which the variance of
    the observable specified by `n_obs_shock` will be maximized.

### Keywords
* `flip_shocks::Bool`: by default, we compute the impulse responses to a negative shock.
    Set `flip_shocks = true` to obtain a positive shock.
"""
function maxBC_shock(β::Matrix{S}, Σ::Matrix{S}, n::Int, n_obs_shock::Int, shock_size::S,
                     lags::Int, frequency_band::Tuple{S,S};
                     flip_shocks::Bool = false) where {S<:Real}
    if lags * n < size(β,1)
        β = @views β[2:end, :]
    end

    cholmat = cholesky((Σ + Σ') ./ 2).L
    increment = abs(frequency_band[1] - frequency_band[2]) / 200.
    V = zeros(S, n, n) # variance
    eminusif = zeros(Complex{S}, 1, 1, lags)
    for f = frequency_band[1]:increment:round(frequency_band[2], digits=10) # not rounding sometimes leads to one fewer loop than desired
        eminusif[1, 1, :] = exp.(-im .* f .* collect(1:lags))
        sumB = dropdims(sum(reshape(β', n, n, lags) .*
                           repeat(eminusif, n, n, 1); dims = 3), dims = 3)
        invA = (Matrix{Complex{S}}(I, n, n) - sumB) \ cholmat
        V += reshape(real.(kron(conj(invA[n_obs_shock, :]), invA[n_obs_shock, :])), n, n) .*
            increment ./ abs(frequency_band[1] - frequency_band[2])
    end
    eigout = eigen(V)
    q = eigout.vectors[:, argmax(eigout.values)]
    q .*= sign(q[n_obs_shock])
    q .*= flip_shocks ? shock_size : -shock_size # negative by DSGE convention

    return (cholmat * q)'
end

"""
```
cholesky_long_run_shock(β, Σ, n, n_obs_shock, shock_size, lags, frequency_band,
    flip_shocks = false) where {S<:Real}
```
computes the long-run Cholesky-identified shock to the observable
specified by `n_obs_shock`.

Given a VAR system
```
yₜ = B₁yₜ₋₁ + ... + B₁yₜ₋ₚ + Γϵₜ,      ϵₜ ∼ 𝒩 (0, Σ),
```

the long-run covariance matrix is
```
S̃ = B̃⁻¹ Σ (B̃⁻¹)'
```

and the Cholesky identification is given by
```
ΓΓ' = Σ ⇒ Γ = B̃ * cholesky(S̃).
```

### Inputs
* `β::AbstractMatrix{S}`: coefficient matrix
* `Σ::AbstractMatrix{S}`: innovations variance-covariance matrix
* `n::Int`: number of observables
* `n_obs_shock::Int`: index of the observable corresponding to the orthogonalized shock
    causing the impulse response.
* `shock_size::S`: number of standard deviations of the shock
* `lags::Int`: number of lags in VAR system
* `frequency_band::Tuple{S,S}`: the frequencies between which the variance of
    the observable specified by `n_obs_shock` will be maximized.

### Keywords
* `flip_shocks::Bool`: by default, we compute the impulse responses to a negative shock.
    Set `flip_shocks = true` to obtain a positive shock.
"""
function cholesky_long_run_shock(β::Matrix{S}, Σ::Matrix{S}, n_obs_shock::Int, n::Int,
                                 lags::Int, shock_size::S;
                                 flip_shocks::Bool = false) where {S<:Real}
    if n * lags < size(β, 1)
        β = β[2:end,:] # don't need the constant
    end

    # Compute decomposition
    B̃ = Matrix{S}(I, n, n) - dropdims(sum(reshape(β', n, n, lags), dims = 3), dims = 3)
    S̃ = B̃ \ (Σ * inv(B̃)')             # LR covariance = S̃ = B̃⁻¹ * Σ * B̃⁻¹' =>
    Γ = B̃ * cholesky((S̃ + S̃') ./ 2).L # S = B̃ \ (Σ * B̃⁻¹')

    # Compute shock
    vec_shock = zeros(n)
    vec_shock[n_obs_shock] = flip_shocks ? shock_size : -shock_size # negative by DSGE convention
    return (Γ * vec_shock)'
end

"""
```
impulse_responses(β, Σ, coint_mat, n_obs_shock, horizon, shock_size = 1;
    method = :cholesky, flip_shocks = false, use_intercept = true,
    frequency_band = (2π/32, 2π/6)) where {S<:Real}
```
computes the impulse responses of a VECM system represented in the form

```
Δyₜ = eₜ₋₁ βₑ + Xₜ βᵥ + ϵₜ,
```
where `βₑ` are the coefficients for the error correction terms;
`eₜ` are the error correction terms specifying the cointegrating relationships;
`βᵥ` are the coefficients for the VAR terms (including the intercept);
`Xₜ` are the lags of observables in period `t`, i.e. `yₜ₋₁, yₜ₋2, ..., yₜ₋ₚ`;
and `ϵₜ ∼ 𝒩 (0, Σ)`. We assume that `β = [βₑ; βᵥ]`.

### Inputs
* `β::AbstractMatrix{S}`: coefficient matrix
* `Σ::AbstractMatrix{S}`: innovations variance-covariance matrix
* `coint_mat::AbstractMatrix{S}`: matrix specifying the cointegrating relationships.
    Multiplying `coint_mat * data`, where `data` is an `n_observables × T` matrix, should yield
    an `n_coint × T` matrix, where `n_coint` are the number of cointegrating
    relationships and `T` are the number of periods of data.
* `n_obs_shock::Int`: index of the observable corresponding to the orthogonalized shock
    causing the impulse response.
* `shock_size::S`: number of standard deviations of the shock

### Keywords
* `method::Symbol`: type of impulse response to compute. The available option is
    `:cholesky` (default).
* `flip_shocks::Bool`: by default, we compute the impulse responses to a negative shock.
    Set `flip_shocks = true` to obtain a positive shock.
* `use_intercept::Bool`: `impulse_responses` assumes `β` has constant term(s). If there
    are no such terms, then `use_intercept` must be set to `false`.
* `frequency_band::Tuple{S,S}`: See `?maxBC_shock`.

### Outputs
* `Y::AbstractMatrix`: Impulse response matrix with dimensions horizons x n_observables
"""
function impulse_responses(β::AbstractMatrix{S}, Σ::AbstractMatrix{S},
                           coint_mat::AbstractMatrix{S}, n_obs_shock::Int,
                           horizon::Int, shock_size::S = one(S);
                           method::Symbol = :cholesky,
                           flip_shocks::Bool = false,
                           use_intercept::Bool = true,
                           frequency_band::Tuple{S,S} =
                           (2*π/32, 2*π/6)) where {S<:Real}

    # Compute dimensions
    n = size(Σ, 1)
    n_coint = size(coint_mat, 1)
    lags = convert(Int, use_intercept ? (size(β, 1) - 1 - n_coint) / n : (size(β, 1) - n_coint) / n)

    # Compute impact based on IRF type
    Y = zeros(S, n, lags + horizon)
    Y[:, lags + 1] = if method == :cholesky
        cholesky_shock(Σ, n, n_obs_shock, shock_size;
                       flip_shocks = flip_shocks)
    # TODO
    # elseif method == :maximum_business_cycle_variance || method == :maxBC
    #     maxBC_shock(β, Σ, n, n_obs_shock, shock_size, lags, frequency_band;
    #                 flip_shocks = flip_shocks)
    # elseif method == :choleskyLR || method == :cholesky_long_run
    #     cholesky_long_run_shock(β, Σ, n_obs_shock, n, lags, shock_size;
    #                             flip_shocks = flip_shocks)
    else
        error("IRF method $(string(method)) has not been implemented.")
    end

    # For efficiency
    if use_intercept
        β = β[vcat(1:n_coint, n_coint + 2:n_coint + 1 + n * lags), :]
    end

    # Compute impulse response
    addcoint = zeros(S, n_coint) # assume deviations in cointegrating relationships also start at zero
    for t = 2:horizon
        addcoint += coint_mat * Y[:, lags + t - 1]
        xT = vcat(addcoint, vec(Y[:, lags + t - 1:-1:t])) # stacks addcoint w/ lag yₜ₋₁, ..., yₜ₋ₚ
        Y[:, lags + t] = vec(xT' * β)
    end

    return Y[:, lags + 1:end]
end


text:
import SMC
import Distributions
import JSON
using DataFrames

@everywhere begin
    using SMC
    using Distributions
    include("dpmm.jl")
    include("smc_samplers.jl")
    include("../aide.jl")
end

include("estimation.jl")

data_dir = "data"

# load subsampled dataset
# (a random subset of 100 elements from the original dataset)
data = readcsv("$data_dir/galaxy_velocities_subset.txt")[:]

# zero mean and unit variance
data = data - mean(data)
data = data / std(data)
num_data = length(data)

# set prior
m_prior = Normal(0, 1) # mean cluster center
r_prior = Gamma(1, 1) # r * data-precision is the precision of cluster centers
nu_prior = Gamma(1, 1)# expected data-precision of cluster is nu/s
s_prior = Gamma(1, 1)
nign_params_prior = NIGNParamsPrior(m_prior, r_prior, nu_prior, s_prior)
alpha_prior = Gamma(1, 1)
hypers = Hyperparameters(nign_params_prior, alpha_prior)

num_sweeps_list = [0, 1, 2, 4]
num_particles_list= [1, 2, 4, 10, 20, 40, 100, 200, 400]

# the number of runs of inference to use for diagnostic is twice the number of
# AIDE replicates, because each AIDE replicate involves running the target
# inference sampler twice
num_aide_replicates = 62
num_diagnostic_replicates = num_aide_replicates * 2
println("generating sampler data...")
gold_standard_sampler = make_optimal_proposal_smc_sampler(data, hypers, 4, 400)
data = generate_sampler_data(data, hypers, num_aide_replicates, num_diagnostic_replicates,
                             gold_standard_sampler)
writetable("$data_dir/galaxy_sampler_data.csv", data)
println("done!")