text:
# testing some special cases that can get missed in more complicated tests
# these are more checks for crashes than accuracy
# tests zeroth order mass action reactions and passing just one mass action reaction
# tests FW and tuple based mix of mass action jump / constant jump

using DiffEqBase, DiffEqJump
using Test
using StableRNGs
rng = StableRNG(12345)

doprint = false
#using Plots; plotlyjs()
doplot = false

methods = (RDirect(), RSSACR(), Direct(), DirectFW(), FRM(), FRMFW(), SortingDirect(), NRM(), RSSA(), DirectCR())

# one reaction case, mass action jump, vector of data
rate = [2.0]
rs = [[0 => 1]]
ns = [[1 => 1]]
jump = MassActionJump(rate, rs, ns)
prob = DiscreteProblem([100],(0.,100.))
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using vectors of data: last val = ", sol[end, end])
end
if doplot
    plothand = plot(sol, lab="vec, 0=>1")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, data as scalars
rate = 2.0
rs = [0 => 3]                   # stoich power should be ignored
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using scalar data: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="scalar, 0=>1")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, vector of empty data
rate = [2.0]
rs = [Vector{Pair{Int,Int}}()]
ns = [[1 => 1]]
jump = MassActionJump(rate, rs, ns)
prob = DiscreteProblem([100],(0.,100.))
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using vector of Pair{Int,Int}: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="vec, ()")
end
@test sol[end,end] > 200

# one reaction case, mass action jump, data as scalars
rate = 2.0
rs = Vector{Pair{Int,Int}}()
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
jump_prob = JumpProblem(prob, Direct(), jump; rng=rng)
sol = solve(jump_prob, SSAStepper())
if doprint
    println("mass act jump using scalar Pair{Int,Int}: last val = ", sol[end, end])
end
if doplot
    plot!(plothand, sol, lab="scalar, ()")
end
@test sol[end,end] > 200

# mix two rx types
# 0 -> A and A -> 0
rate = 600.0
rs = [0 => 3]                   # stoich power should be ignored
ns = [1 => 1]
jump = MassActionJump(rate, rs, ns)
ratefun = (u,p,t) -> 2.0*u[1]
affect! = function(integrator)
    integrator.u[1] -= 1
end
jump2 = ConstantRateJump(ratefun, affect!)
if doplot
    plothand2 = plot(reuse=false)
end

dep_graph = [
    [1, 2],
    [1, 2]
]
spec_to_dep_jumps = [[2]]
jump_to_dep_specs = [[1],[1]]
namedpars = (dep_graph=dep_graph, vartojumps_map=spec_to_dep_jumps, jumptovars_map=jump_to_dep_specs)

for method in methods
    local jump_prob = JumpProblem(prob, method, jump, jump2; rng=rng, namedpars...)
    local sol = solve(jump_prob, SSAStepper())

    if doplot
        plot!(plothand2, sol, label=("A <-> 0, " * string(method)))
    end

    if doprint
        println("Mix of constant and mass action jumps, method = ", typeof(method), ", sol[end] = ", sol[end,end])
    end
    @test sol[end,end] > 200
end

if doplot
    display(plothand)
    display(plothand2)
end

# check MassActionJump with parameters merge together correctly
rs1 = [1 => 1]
rs2 = [2=>1]
ns1 = [1 => -1, 2 => 1]
ns2 = [1=>1,2=>-1]
p  = [1.0,0.0]
maj1 = MassActionJump(rs1, ns1; param_idxs=1)
maj2 = MassActionJump(rs2, ns2; param_idxs=2)
js   = JumpSet(maj1,maj2)
maj  = MassActionJump([rs1,rs2],[ns1,ns2]; param_idxs=[1,2])
@test all(getfield(maj,fn) == getfield(js.massaction_jump,fn) for fn in [:scaled_rates,:reactant_stoch,:net_stoch])
@test all(maj.param_mapper.param_idxs .== js.massaction_jump.param_mapper.param_idxs)
maj1 = MassActionJump([rs1], [ns1]; param_idxs=[1])
maj2 = MassActionJump([rs2], [ns2]; param_idxs=[2])
js   = JumpSet(maj1,maj2)
@test all(getfield(maj,fn) == getfield(js.massaction_jump,fn) for fn in [:scaled_rates,:reactant_stoch,:net_stoch])
@test all(maj.param_mapper.param_idxs .== js.massaction_jump.param_mapper.param_idxs)


text:
# reduce() versions
@inline function reduce{FSA <: Union{FixedArray, Tuple}}(f, a::FSA)
    length(a) == 1 && return a[1]
    @inbounds begin
        red = f(a[1], a[2])
        for i=3:length(a)
            red = f(red, a[i])
        end
    end
    red
end

@inline function reduce(f::Functor{2}, a::Mat)
    length(a) == 1 && return a[1,1]
    @inbounds begin
        red = reduce(f, Tuple(a)[1])
        for i=2:size(a, 2)
            red = f(red, reduce(f, Tuple(a)[i]))
        end
    end
    red
end

#------------------------------------------------------------------------------
# map() machinery

# Get an expression indexing the collection `name` of type `T` for use in map()
index_expr{T <: Number}(::Type{T},     name, inds::Int...)    = :($name)
index_expr{T <: FixedArray}(::Type{T}, name, inds::Int...)    = :($name[$(inds...)])
index_expr{T <: AbstractArray}(::Type{T}, name, inds::Int...) = :($name[$(inds...)])

# Get expression checking size of collection `name` against `SIZE`
function sizecheck_expr{T <: Number}(::Type{T}, name, SIZE)
    :nothing
end
function sizecheck_expr{FSA<:FixedArray}(::Type{FSA}, name, SIZE)
    if size(FSA) == SIZE
        :nothing
    else
        :(throw(DimensionMismatch(string($FSA)*" is wrong size")))
    end
end
function sizecheck_expr{A<:AbstractArray}(::Type{A}, name, SIZE)
    quote
        # Note - should be marked with @boundscheck in 0.5
        size($name) == $SIZE || throw(DimensionMismatch(string($A)*" is wrong size"))
    end
end

# Type wrapper to signify that a similar FSA should be constructed as the map()
# output type via construct_similar()
immutable SimilarTo{FSA}; end

# Get expression to construct FSA from a nested tuple of storage generated by map()
constructor_expr{FSA <: FixedArray}(::Type{FSA}, tuple_expr::Expr) = :($FSA($tuple_expr))
constructor_expr{FSA <: FixedVectorNoTuple}(::Type{FSA}, tuple_expr::Expr) = :($FSA($tuple_expr...))
constructor_expr{FSA <: FixedArray}(::Type{SimilarTo{FSA}}, tuple_expr::Expr) = :(construct_similar($FSA, $tuple_expr))

# Generate an unrolled nested tuple of calls mapping funcname across the input,
# constructing an `OutFSA` to store the result.
#
# julia> FixedSizeArrays.unrolled_map_expr(:f, Vec{2,Bool}, (2,), (Vec{2,Int},Int), (:A,:b))
#
# generates, after cleaning up:
#
#    (FixedSizeArrays.Vec{2,Bool})(
#        tuple(tuple(f(A[1,1],b), f(A[2,1],b)),
#              tuple(f(A[1,2],b), f(A[2,2],b)))
#    )
function unrolled_map_expr(funcname, OutFSA, SIZE, argtypes, argnames)
    sizecheck = [sizecheck_expr(T,n,SIZE) for (T,n) in zip(argtypes,argnames)]
    tuple_expr = fill_tuples_expr(SIZE) do inds...
        Expr(:call, funcname,
            [index_expr(argtypes[i], argnames[i], inds...) for i=1:length(argtypes)]...
        )
    end
    quote
        $(Expr(:meta, :inline))
        $(sizecheck...)
        @inbounds rvalue = $(constructor_expr(OutFSA, tuple_expr))
        rvalue
    end
end

# map() comes in two flavours:
#
# 1) You can specify the output type with the first argument
#   map(func, ::Type{OutFSA}, arg1, arg2, ...)

# General N-ary version with explicit output.  (Unary and binary versions are
# written out below since this generates better code in julia-0.4.)
@generated function map{OutFSA<:FixedArray}(func, ::Type{OutFSA}, args...)
    argexprs = ntuple(i->:(args[$i]), length(args))
    unrolled_map_expr(:func, OutFSA, size(OutFSA), args, argexprs)
end
# Binary
@generated function map{OutFSA<:FixedArray}(func, ::Type{OutFSA}, arg1, arg2)
    unrolled_map_expr(:func, OutFSA, size(OutFSA), (arg1,arg2), (:arg1,:arg2))
end
# Unary
@generated function map{OutFSA<:FixedArray}(func, ::Type{OutFSA}, arg1)
    unrolled_map_expr(:func, OutFSA, size(OutFSA), (arg1,), (:arg1,))
end

# 2) You can let map() infer the output FixedArray type from `func`
#   map(func, fixed_arrays...)

# Binary, inferred output
@generated function map{F1<:FixedArray, F2<:FixedArray}(func, arg1::F1, arg2::F2)
    unrolled_map_expr(:func, SimilarTo{arg1}, size(F1), (arg1,arg2), (:arg1,:arg2))
end
@generated function map{F<:FixedArray}(func, arg1::F, arg2::Union{Number,AbstractArray})
    unrolled_map_expr(:func, SimilarTo{arg1}, size(F), (arg1,arg2), (:arg1,:arg2))
end
@generated function map{F<:FixedArray}(func, arg1::Union{Number,AbstractArray}, arg2::F)
    unrolled_map_expr(:func, SimilarTo{arg2}, size(F), (arg1,arg2), (:arg1,:arg2))
end
# Unary, inferred output
@generated function map{F<:FixedArray}(func, arg1::F)
    unrolled_map_expr(:func, SimilarTo{arg1}, size(F), (arg1,), (:arg1,))
end

# Nullary special case version.
#
# TODO: This is inconsistent, since it maps *indices* through the functor
# rather than using a functor with no arguments.
@generated function map{FSA <: FixedArray}(F, ::Type{FSA})
    tuple_expr = fill_tuples_expr((inds...) -> :(F($(inds...))), size(FSA))
    constructor_expr(FSA, tuple_expr)
end

@inline map{T}(::Type{T}, v::FixedArray{T}) = v

ext: jl


text:
# """
# OLD version: uses gaussian_process.jl directly. New version uses data.jl to prepare the data.
# Interactive exploration of the hyperparameter space using Interact.jl.
# """

using JLD2, Statistics, LinearAlgebra, Plots
using Interact, Blink

include("kernels.jl")
include("gaussian_process.jl")
include("../les/get_les_data.jl")

# file to gather data from
filename = togglebuttons(OrderedDict("general_strat_16_profiles" =>"general_strat_16_profiles.jld2",
                                     "general_strat_32_profiles" =>"general_strat_32_profiles.jld2"),
                                     label="LES")

D = 16 # gridpoints

# these variables are common to the files in this example
# -- CHANGE if using different files!! --
data = get_les_data(filename[])
t = data.t
Nt = length(t)

total_set = 1:(Nt-1)
training_set = 1:4:(Nt-1)
verification_set = setdiff(total_set, training_set)

# which variable to explore
V_name = togglebuttons(Dict("T" =>"Temperature [°C]", "wT"=>"Temperature flux [°C⋅m/s]"), label="profile")
# smooth the profile
smooth = toggle("smooth profile?")
# normalize the data (pre- / postprocessing)
normalize = toggle("normalize?")

γs = -3.0:0.1:3.0
σs = 0.0:0.1:2.0

# hyperparameter knobs
γ1 = slider(γs, label="log length scale, log₁₀(γ)")
σ1 = slider(σs, label="log signal variance, log₁₀(σ²)")
time_index = slider(1:40:(Nt-1), label="time [s]")

dist_metric = tabulator(OrderedDict("l²-norm"  =>  "l²-norm:  d(x,x') = || x - x' ||",
                                    "H¹-norm"  =>  "H¹-norm:  d(x,x') = || diff(x)./diff(z) - diff(x')./diff(z) ||",
                                    "H⁻¹-norm" =>  "H⁻¹-norm: d(x,x') = || diff(x).*diff(z) - diff(x').*diff(z) ||"
                                    ))

kern = tabulator(OrderedDict("Squared exponential"          => "Squared exponential kernel:           k(x,x') = σ * exp( - ||x-x'||² / 2γ² )",
                             "Matern 1/2"                    => "Matérn with ʋ=1/2:                    k(x,x') = σ * exp( - ||x-x'|| / γ )",
                             "Matern 3/2"                    => "Matérn with ʋ=3/2:                    k(x,x') = σ * (1+c) * exp(-√(3)*||x-x'||)/γ)",
                             "Matern 5/2"                    => "Matérn with ʋ=5/2:                    k(x,x') = σ * ( 1 + √(5)*||x-x'||)/γ + 5*||x-x'||²/(3*γ^2) ) * exp(-√(5)*||x-x'||)/γ)",
                             "Rational quadratic w/ α=1"     => "Rational quadratic kernel:            k(x,x') = σ * (1+(x-x')'*(x-x')/(2*α*(γ²))^(-α)",
                            ))

get_data(filename::String) = get_les_data(filename)

function get_V(V_name, data)
    if V_name=="Temperature [°C]"
        return data.T
    end
    if V_name=="Temperature flux [°C⋅m/s]"
        return data.wT
    else
        throw(error())
    end
end

function get_d(dist_metric)
    if dist_metric==1; return l2_norm end
    if dist_metric==2; return h1_norm end
    if dist_metric==3; return hm1_norm
    else
        throw(error())
    end
end

get_vavg(V)       = [custom_avg(V[:,j], D) for j in 1:Nt] # compress variable array to D values per time
get_x(vavg)       = vavg[1:(Nt-1)] # (v₀, v₁, ... ,v_(Nt-1)) (Nt-1)-length array of D-length inputs
get_y(vavg)       = vavg[2:(Nt)]   # (v₁, v₂, ... ,v_Nt    ) (Nt-1)-length array of D-length targets
get_x_train(x)    = x[training_set]
get_y_train(y)    = y[training_set]
get_z(data)       = custom_avg(data.z, D)

function get_kernel(k::Int64, γ, σ)
    # convert from log10 scale
    γ = 10^γ
    σ = 10^σ
  if k==1; return SquaredExponentialKernelI(γ, σ) end
  if k==2; return Matern12I(γ, σ) end
  if k==3; return Matern32I(γ, σ) end
  if k==4; return Matern52I(γ, σ) end
  if k==5; return RationalQuadraticI(γ, σ, 1.0)
  else; throw(error()) end
end

function plot_kernel(kernel::Kernel, d, z)
    kmat = [kernel_function(kernel; d=d, z=z)(i,j) for i in 1:10:Nt, j in 1:10:Nt]# fill kernel mx with values
    return heatmap(kmat, title = "Covariance Matrix", xaxis=(:false), yaxis=(:flip, :false), clims=(0.0,100), legend=true)
end

function get_gp(x_train, y_train, kernel::Kernel, normalize, d, z)
    𝒢 = construct_gpr(x_train, y_train, kernel; distance_fn=d, z=z, normalize=normalize);
    return 𝒢
end

function get_gpr_pred(gp::GP, x, vavg)
    # predict temperature profile from start to finish without the training data
    gpr_prediction = similar(vavg[1:Nt-1])
    starting = x[1]
    gpr_prediction[1] = starting
    for i in 1:(Nt-2)
        gpr_prediction[i+1] = prediction([gpr_prediction[i]], gp)
    end
    return gpr_prediction
end

function plot_profile(gp::GP,V,V_name,time_index,gpr_prediction,z)
    exact = V[:,time_index+1]
    day_string = string(floor(Int, t[time_index]/86400))
    p = scatter(gpr_prediction[time_index+1], z, label = "GP")
    plot!(exact, z, legend = :topleft, label = "LES", xlabel = "$(V_name)", ylabel = "depth", title = "day " * day_string)
    return p
end

function plot_log_error(gp::GP, time_index, x, y, vavg)
    # gpr_error = collect(verification_set)*1.0
    # # greedy check
    # for j in eachindex(verification_set)
    #     test_index = verification_set[j]
    #     y_prediction = prediction([x[test_index]], gp)
    #     error = norm(y_prediction - y[test_index])
    #     gpr_error[j] = error
    # end
    # mean_error = sum(gpr_error)/length(gpr_error)
    # error_plot_log = histogram(log.(gpr_error), title = "log(Error)", xlabel="log(Error)", ylabel="Frequency",xlims=(-20,0),ylims=(0,250), label="frequency")
    # vline!([log(mean_error)], line = (4, :dash, 0.8), label="mean error")
    # vline!([log(gpr_error[time_index])], line = (1, :solid, 0.6), label="error at t=$(time_index)")

    #compute error for true check
    gpr_error = zeros(Nt-2)
    gpr_prediction = get_gpr_pred(gp, x, vavg)
    for i in 1:Nt-2
        exact    = y[i+1]
        predi    = gpr_prediction[i+1]
        gpr_error[i] = l2_norm(exact, predi) # euclidean distance
    end
    mean_error = sum(gpr_error)/(Nt-2)

    error_plot_log = histogram(log.(gpr_error), title = "log(error) at each timestep of the full evolution", xlabel="log(Error)", ylabel="Frequency",xlims=(-20,0),ylims=(0,250), label="frequency")
    vline!([log(mean_error)], line = (4, :dash, 0.8), label="mean error")
    vline!([log(gpr_error[time_index])], line = (1, :solid, 0.6), label="error at t=$(time_index)")

end

function plot_hyp_landscp(k::Int64, x_train, y_train, x, y, vavg, normalize, d, z)

    mlls = zeros(length(γs)) # mean log marginal likelihood
    mes  = zeros(length(γs)) # mean error (greedy check)
    mets  = zeros(length(γs)) # mean error (true check)

    for i in 1:length(γs)
        kernel = get_kernel(k, γs[i], 0.0)
        𝒢 = construct_gpr(x_train, y_train, kernel; distance_fn=d, z=z, normalize=normalize);
        mlls[i] = -1*mean_log_marginal_loss(y_train, 𝒢::GP, add_constant=false)

        #compute mean error for greedy check (same as in plot log error)
        total_error = 0.0
        # greedy check
        for j in eachindex(verification_set)
            test_index = verification_set[j]
            y_prediction = prediction([x[test_index]], 𝒢)
            error = l2_norm(y_prediction, y[test_index])
            total_error += error
        end
        mes[i] = total_error/length(verification_set)

        total_error = 0.0
        #compute mean error for true check
        gpr_prediction = get_gpr_pred(𝒢, x, vavg)
        for i in 1:Nt-2
            exact    = y[i+1]
            predi    = gpr_prediction[i+1]
            total_error += l2_norm(exact, predi) # euclidean distance
        end
        mets[i] = total_error/(Nt-2)

    end

    mll_plot = plot(γs, mlls, xlabel="log(γ)", title="negative mean log marginal likelihood, P(y|X)", legend=false, yscale=:log10) # 1D plot: mean log marginal loss vs. γ
    # me_plot  = plot(γs, mets,  xlabel="log(γ)", title="mean error on full evolution ('true check'), min = $(minimum(mets))", legend=false, yscale=:log10)  # 1D plot: mean error vs. γ
    me_plot  = plot(γs, mets,  xlabel="log(γ)", title="min = $(minimum(mets))", legend=false, yscale=:log10)  # 1D plot: mean error vs. γ

    return plot(mll_plot, me_plot, layout = @layout [a ; b])
end

#updating variables
#output               function          args
data            = map(get_data,         filename)
d               = map(get_d,            dist_metric)
V               = map(get_V,            V_name, data)
vavg            = map(get_vavg,         V)
x               = map(get_x,            vavg)
y               = map(get_y,            vavg)
z               = map(get_z,            data)
x_train         = map(get_x_train,      x)
y_train         = map(get_y_train,      y)
k               = map(get_kernel,       kern, γ1, σ1)
k_plot          = map(plot_kernel,      k, d, z)
gp              = map(get_gp,           x_train, y_train, k, normalize, d, z)
gpr_prediction  = map(get_gpr_pred,     gp, x_train, vavg)
profile_plot    = map(plot_profile,     gp, V, V_name, time_index, gpr_prediction, z)
log_error_plot  = map(plot_log_error,   gp, time_index, x, y, vavg)
hyp_landscape   = map(plot_hyp_landscp, kern, x_train, y_train, x, y, vavg, normalize, d, z)

println(normalize[])

# layout
top    = vbox(hbox(filename, V_name, vbox(smooth, normalize)), hbox(kern, dist_metric), hbox(k_plot, hyp_landscape))
middle = vbox(γ1, σ1, time_index)
bottom = hbox(profile_plot, log_error_plot)
ui = vbox(top, middle, bottom) # aligns vertically

# Blink GUI
window = Window()
body!(window, ui)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using JuMP
using Base.Test
using ReverseDiffSparse2

function compare_jump(m::Model)

    nl = JuMPNLPEvaluator(m)
    nvar = MathProgBase.numvar(m)
    x = rand(nvar)
    println("JuMP init time:")
    @time MathProgBase.initialize(nl, [:Grad,:ExprGraph])

    objexpr = MathProgBase.obj_expr(nl)

    println("RDS2 init time:")

    @time nd,const_values = expr_to_nodedata(objexpr)
    @time adj = adjmat(nd)

    storage = zeros(length(nd))
println("RDS2 Feval time")
    forward_eval(storage,nd,adj,const_values,x)
    @time fval = forward_eval(storage,nd,adj,const_values,x)

    @show fval
    MathProgBase.eval_f(nl,x)
println("JuMP Feval time")
    @time fval = MathProgBase.eval_f(nl,x)
    @show fval

    reverse_storage = zeros(length(nd))
    grad = zeros(nvar)
    ∇f = zeros(nvar)
    MathProgBase.eval_grad_f(nl, ∇f, x)
    reverse_eval(grad,reverse_storage,storage,nd,adj,const_values,x)
    @test isapprox(grad,∇f)
    println("RDS2 grad time")
    @time reverse_eval(grad,reverse_storage,storage,nd,adj,const_values,x)
    @time reverse_eval(grad,reverse_storage,storage,nd,adj,const_values,x)
    println("JuMP grad time")
    @time MathProgBase.eval_grad_f(nl, ∇f, x)
    @time MathProgBase.eval_grad_f(nl, ∇f, x)
    #for i in 1:10
    #    @profile reverse_eval(grad,reverse_storage,storage,nd,adj,const_values,x)
    #end
end

function mod1(N)

    m = Model()
    @defVar(m, x[1:N])
    c = 10

    @setNLObjective(m, Min, c*sum{x[1]*sum{sin(x[i]);iseven(i)},i=1:N})

    return m
end

#m = mod1(100000)
#compare_jump(m)


text:
using Pkg; Pkg.activate("GraphFusedElasticNet.jl")
using Distributed

# only set to true if running interactive
# not if running from batch file
debug_mode = false  # Base.isinteractive()
if debug_mode
    walltime = 120.0
    ngens = 4
    gensize = 16
    cvsplits = 5
    data_targets = readdir("productivity_splits")
    num_procs = 6
else
    walltime = parse(Float64, ARGS[1])
    ngens = parse(Int, ARGS[2])
    gensize = parse(Int, ARGS[3])
    cvsplits = parse(Int, ARGS[4])
    data_targets = ARGS[5:end]
    num_procs = 1
end
if num_procs > 1
    addprocs(num_procs - 1, exeflags="--project=GraphFusedElasticNet.jl")
end

# relevant path (ADAPT TO CLUSTER IF NECESSARY)
@everywhere begin
    # instantiate and load libraries in every process, define globals
    using JSON
    using DataFrames
    using Dates
    using Distributions
    using Random
    using CSV
    using Base.Threads
    using DelimitedFiles
    using StatsFuns
    using GraphFusedElasticNet

    workdir() = pwd()
    traildata_relpath() = "processed_data/trails.json"
    betas_relpath() = "best_betas/"
    fitmetrics_relpath() = "modelfit_metrics/"
    data_relpath() = "productivity_splits/"
end

# data reader
@everywhere function loadtrails()
    fname = joinpath(workdir(), traildata_relpath())
    trails = JSON.parsefile(fname)
    ptr = convert(Vector{Int}, trails["ptr"])
    brks = convert(Vector{Int}, trails["brks"])
    wts = convert(Vector{Float64}, trails["wts"])
    istemp = convert(Vector{Bool}, trails["istemp"])
    num_nodes = Int(trails["num_nodes"])
    return ptr, brks, wts, istemp, num_nodes
end

# generates sets of nodes to take out during cv
@everywhere function generate_cvsets(y, cvsplits)
    # make the cv splits
    # the only trick is not to count the integers that
    # have missing data
    attempts = y[:, 1]
    N = length(attempts)
    cvsets = [Set{Int}() for i in 1:cvsplits]
    iobs = [i for (i, a) in enumerate(attempts) if a >= 1.0]
    shuffle!(iobs)
    Nobs = length(iobs)
    splitsize = Nobs ÷ cvsplits
    for k in 1:cvsplits
        for i in ((k - 1) * splitsize + 1):(k * splitsize)
            push!(cvsets[k], iobs[i])
        end
    end
    return cvsets
end

# prepare training and test data fro a cvset
@everywhere function get_train_data(y, test_idx)
    train = copy(y)
    for i in 1:size(y, 1)
        if i ∈ test_idx
            train[i, :] = [0.0, 0.0]
        end
    end
    train
end

# fits a binomial GFEN model for given params and data
@everywhere function fit_model(
    ytrain, ptr, brks, wts, istemp, λ1, λ2, η1, η2, modelopts, fitopts; ϵ=1e-5
)
    N = size(ytrain, 1)
    lambdasl1 = [w * (t ? η1 : λ1) for (t, w) in zip(istemp, wts)]
    lambdasl2 = [w * (t ? η2 : λ2) for (t, w) in zip(istemp, wts)]
    # define and train model
    model = BinomialGFEN(ptr, brks, lambdasl1, lambdasl2; modelopts...)
    succ = ytrain[:, 1] .+ ϵ
    att = ytrain[:, 2] .+ 2ϵ
    fit!(model, succ, att; fitopts...)
    return model
end

# finds best hyperparams using cross validation
@everywhere function cv_eval(
    y, ptr, brks, wts, istemp, pars, cvsets, modelopts, fitopts, usethreads
)
    # for each cv split get the mse error
    N = size(y, 1)
    cvsplits = length(cvsets)
    npars = length(pars)
    cv_logll = zeros(npars)
    thrdid = zeros(Int, npars)
    avtime = zeros(npars)
    avalpha = zeros(npars)
    # prepare the tree structure and the bint
    Nobs = sum(y[:, 2])
    Threads.@threads for k = 1:npars
        λ1, λ2, η1, η2 = pars[k]
        for i = 1:cvsplits
            thrdid[k] = threadid()
            # get the cv vector with missing data
            ytrain = get_train_data(y, cvsets[i])
            runningtime = @elapsed begin
                model = fit_model(
                    ytrain, ptr, brks, wts, istemp,
                    λ1, λ2, η1, η2, modelopts, fitopts)
            end
            beta = model.beta
            avtime[k] += runningtime
            avalpha[k] += model.admm_penalty
            # compute the out-of-sample likelihood
            ll = 0.0
            for j in cvsets[i]
                s, N = y[j, :]
                β = beta[j]
                if beta[j] >= 0
                    ll += - N * log(1.0 + exp(-β)) - (N - s) * β
                else
                    ll += s * β - N * log(1.0 + exp(β))
                end
            end
            cv_logll[k] = ll
        end
    end
    cv_logll /= Nobs
    avtime /= cvsplits
    avalpha /= cvsplits
    cv_logll, thrdid, avtime, avalpha
end

# main program
@everywhere function fit_split(filename, walltime, cvsplits, ngens, gensize, eval_corners=true)
    # read trail data (it is then copied to each worker by pmap)
    ptr, brks, wts, istemp, num_nodes = loadtrails()

    # set up gaussian process with smoothing parameters
    # (7 ^ 4) = 2401 parameter space size
    # Kernel matrix is 2401 x 2401 (196,882 entries)
    # for much larger spaces consider more efficient
    # implementations of gaussian processes
    sl1 = [1e-3,  0.25, 0.5, 0.75, 1.25, 2.5, 5.0]
    sl2 = [1e-3,  0.5, 1.0, 2.5, 5.0, 10.0, 20.0]
    tl1 = [1e-3,  0.25, 0.5, 0.75, 1.25, 2.5, 5.0]
    tl2 = [1e-3,  0.5, 1.0, 2.5, 5.0, 10.0, 20.0]

    hparams = hcat([[a, b, c, d]
                    for a in sl1
                    for b in sl2
                    for c in tl1
                    for d in tl2]...)

    nl = size(hparams, 2)
    a, σ, b = 0.1, 0.00001, 0.0001^2
    X = log.(hparams)
    gpsampler = GaussianProcessSampler(X, a=a, σ=σ, b=b)
    gp_offset = 0.0  # empirically assigned to running obs mean

    # model parameters
    modelopts = Dict{Symbol, Any}(
        :admm_balance_every => 10,
        :admm_init_penalty => 0.1,
        :admm_residual_balancing => true,
        :admm_adaptive_inflation => true,
        :reltol => 1e-3,
        :admm_min_penalty => 0.1,
        :admm_max_penalty => 5.0,
        :abstol => 0.,
        :save_norms => true,
        :save_loss => true)
    fitopts = Dict{Symbol, Any}(
        :walltime => walltime,
        :parallel => false)

    # read data
    fname = joinpath(workdir(), data_relpath(), filename)
    y = readdlm(fname, ',')
    cvsets = generate_cvsets(y, cvsplits)

    # df = Vector{DataFrame}(undef, ngens)
    results = []
    corners = Set{Int}()
    for gen in 1:ngens
        # sample lambdas from generation
        if gen == 1
            idx = sample(1:nl, gensize, replace=false)
            prev_pred = fill(-Inf, gensize)
            prev_sd = fill(0.0, gensize)
        else
            idx, _, _, prev_pred, prev_sd  = gpsample(gpsampler, gensize)
        end

        # override to eval first corners
        if eval_corners && length(corners) < 16
            # override to evaluate first in all the corners!
            overriden = 0
            for sl1_ in [minimum(sl1), maximum(sl1)]
                for sl2_ in [minimum(sl2), maximum(sl2)]
                    for tl1_ in [minimum(tl1), maximum(tl1)]
                        for tl2_ in [minimum(tl2), maximum(tl2)]
                            for j in 1:nl
                                if (
                                    overriden < gensize
                                    && !(j in corners)
                                    && hparams[1, j] == sl1_
                                    && hparams[2, j] == sl2_
                                    && hparams[3, j] == tl1_
                                    && hparams[4, j] == tl2_
                                )
                                    overriden += 1
                                    idx[overriden] = j
                                    push!(corners, j)
        end end end end end end end

        # run for every lambda in thread
        fns = fill(filename, gensize)
        pids = fill(getpid(), gensize)
        hosts = fill(gethostname()[1:8], gensize)
        gens = fill(gen, gensize)

        # obtain cv losses and update gaussian process
        usethreads = true
        pars = [hparams[:, i] for i in idx]
        cv_logll, thrds, avtime, avalpha = cv_eval(
            y, ptr, brks, wts, istemp,
            pars, cvsets,
            modelopts, fitopts, usethreads)
        addobs!(gpsampler, idx, cv_logll)
        # update the offset helps with better estimates
        # at far regions from observed data
        # it's like using an empirical prior, it doesn't
        # affect much regions with densely observed data
        N0, N1 = (gen - 1) * gensize, gen * gensize
        gp_offset = (gp_offset *  N0 +  sum(cv_logll)) / N1
        gpsampler.offset = gp_offset

        slambdasl1 = hparams[1, idx]
        slambdasl2 = hparams[2, idx]
        tlambdasl1 = hparams[3, idx]
        tlambdasl2 = hparams[4, idx]
        df = DataFrame(
            fn=fns, pid=pids, host=hosts, thread=thrds,
            λsl1=slambdasl1, λsl2=slambdasl2, λtl1=tlambdasl1, λtl2=tlambdasl2,
            cv_logll=cv_logll, time=avtime, alpha=avalpha, gen=gens,
            paramid=idx, prev_pred=prev_pred, prev_sd=prev_sd)
        println(df)
        push!(results, df)
    end

    # now obtain the best hyperperameters and train the best model
    results = vcat(results...)
    pred, band = gpeval(gpsampler)
    results.final_pred = pred[results.paramid]
    results.final_sd = band[results.paramid]
    #
    bestidx = argmax(results.final_pred)
    best_predicted = results.final_pred[bestidx]
    λ1, λ2, η1, η2 = hparams[:, results.paramid[bestidx]]

    gpdf = DataFrame(
        slambdasl1=hparams[1, :],
        slambdasl2=hparams[2, :],
        tlambdasl1=hparams[3, :],
        tlambdasl2=hparams[4, :],
        pred=pred,
        band=band,
        tested=gpsampler.tested)

    modelopts[:reltol] /= 10.0
    fitopts = Dict{Symbol, Any}(
        :parallel => true,
        :walltime => 16.0 * walltime)
    runtime_final = @elapsed begin
        model = fit_model(
            y, ptr, brks, wts, istemp,
            λ1, λ2, η1, η2, modelopts, fitopts)
    end
    betasfile = joinpath(
        workdir(),
        betas_relpath(),
        "betas_" * filename
    )
    resultsfile = joinpath(
        workdir(),
        fitmetrics_relpath(),
        "cvloss_" * filename
    )
    gpfile = joinpath(
        workdir(),
        fitmetrics_relpath(),
        "gp_" * filename
    )

    # write results
    CSV.write(resultsfile, results)
    CSV.write(gpfile, gpdf)
    open(betasfile, "w") do f
        writedlm(f, model.beta, ',')
    end

    println("Best model:")
    println("   split:     ", filename)
    println("   params:    ", (λ1, λ2, η1, η2))
    println("   predicted: ", best_predicted)
    println("   runtime:   ", runtime_final)
    println("   generation:   ", results.gen[bestidx])
    steps = model.steps
    conv = model.converged
    t = "$(Dates.Time(Dates.now()))"[1:5]
    pnorm = model.prim_norms[end]
    dnorm = model.dual_norms[end]
    loss = model.loss[end]
    admm_penalty = model.admm_penalty
    println("steps=$steps, conv=$conv, t=$t, loss=$loss, pnorm=$pnorm, dnorm=$dnorm, alpha=$admm_penalty")
end

# define and call main routine using distributed computing
function main(data_targets, walltime, ngens, gensize, cvsplits)
    # if in slurm cluster environment (TACC)
    local_rank = haskey(ENV, "SLURM_PROCID") ? parse(Int, ENV["SLURM_PROCID"]) : -1

    if local_rank >= 0
        # running in MPI style of parallel program, each program has
        #   a rank and selects one target
        data_targets = [data_targets[local_rank + 1]]
    end

    if num_procs > 1
        @sync @distributed for filename in data_targets
            # running in in-program distributed mode, the program will
            #   schedule a file for each subprocess
            m = "Processing $(filename) in process $(getpid()) with rank $(local_rank)..."
            println(m)
            fit_split(filename, walltime, cvsplits, ngens, gensize)
        end
    else
        for filename in data_targets
            m = "Processing $(filename) in process $(getpid()) with rank $(local_rank)..."
            println(m)
            fit_split(filename, walltime, cvsplits, ngens, gensize)
        end
    end
end

main(data_targets, walltime, ngens, gensize, cvsplits)

max_stars_repo_name: mauriciogtec/gfen-reproduce

max_issues_repo_name: mauriciogtec/gfen-reproduce

max_forks_repo_name: mauriciogtec/gfen-reproduce


text:
"""rainfall noises"""
function rainfall_noises(sp, data::Dict, params::Dict, t::Int)
        SDDP.parameterize(sp, collect(1:size(data["hydro"]["scenario_probabilities"],2)), data["hydro"]["scenario_probabilities"][cidx(t,data["hydro"]["size_inflow"][1]),:]) do ω
            # if JuMP.MathOptInterface.get(sp,JuMP.MathOptInterface.VariablePrimalStart(), JuMP.all_variables(sp)[end]) == nothing
            #     JuMP.MathOptInterface.set.(sp,JuMP.MathOptInterface.VariablePrimalStart(), JuMP.all_variables(sp), NaN)
            # end
            nostart = findall(x-> x == nothing, JuMP.MathOptInterface.get.(sp,JuMP.MathOptInterface.VariablePrimalStart(), JuMP.all_variables(sp)))
            for theta in nostart
                JuMP.MathOptInterface.set(sp,JuMP.MathOptInterface.VariablePrimalStart(), JuMP.all_variables(sp)[theta], sp.ext[:lower_bound])
            end
            for i in 1:data["hydro"]["nHyd"]
                JuMP.fix(sp[:inflow][i], data["hydro"]["Hydrogenerators"][i]["inflow"][t,ω]; force=true)
            end
        end
    return nothing
end

"""creates hydro balance constraint"""
function constraint_hydro_balance(sp, data::Dict, params::Dict)
    k = 0.0036 # conversion factor from m3/s to hm3
    @constraint(sp, hydro_balance[i=1:data["hydro"]["nHyd"]],   sp[:reservoir][i].out == sp[:reservoir][i].in + (sp[:inflow][i] - sp[:outflow][i])*(k*params["stage_hours"]) - sp[:spill][i] +
                                                    sum(sp[:spill][j] for j in data["hydro"]["Hydrogenerators"][i]["upstrem_hydros_spill"]) +
                                                    sum(sp[:outflow][j] for j in data["hydro"]["Hydrogenerators"][i]["upstrem_hydros_turn"])*(k*params["stage_hours"])
    )
    return nothing
end

"""creates energy constraints which bind the discharge with the active energy injected to the grid"""
function constraint_hydro_generation(sp, data::Dict, pm::AbstractPowerModel)
    @constraint(sp, turbine_energy[i=1:data["hydro"]["nHyd"]; data["hydro"]["Hydrogenerators"][i]["index_grid"] != nothing], PowerModels.var(pm, 0, 1, :pg)[data["hydro"]["Hydrogenerators"][i]["i_grid"]]*data["powersystem"]["baseMVA"] == sp[:outflow][i]*data["hydro"]["Hydrogenerators"][i]["production_factor"]
    )
    return nothing
end

"""add deficit variables"""
function constraint_mod_deficit(sp, data::Dict, pm::AbstractPowerModel)
    for i=1:length(PowerModels.con(pm, 0, 1, :kcl_p))
        set_coefficient(PowerModels.con(pm, 0, 1, :kcl_p)[i], sp[:deficit][i], -1)
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dccad7c8-62fb-11eb-226f-393c73301bcb
begin
	using Dates
	using Pkg
end

# ╔═╡ 5938a9a6-6099-11eb-1ac9-1fb12d9c9237
using ReinforcementLearning

# ╔═╡ ad1302aa-609a-11eb-0a22-892b5c80750d
using Random

# ╔═╡ 06013422-609b-11eb-24e0-790eb0272183
using Plots

# ╔═╡ ac7e107c-609b-11eb-2d50-4b50517e1840
using Flux:InvDecay

# ╔═╡ bed48ef6-62fb-11eb-327a-ed29787dda3c
md"""
# How to Write a Customized Environment in ReinforcementLearning.jl?
"""

# ╔═╡ e262295a-62fb-11eb-32c0-051c8995d14a
md"""
- Last Update: $(now())
- Julia Version: $VERSION
- ReinforcementLearning.jl Version: $([v.version for (k,v) in Pkg.dependencies() if v.name=="ReinforcementLearning"][1])
"""

# ╔═╡ 03852aa2-6099-11eb-39db-ffba5eeade98
md"""
The first step to apply algorithms in ReinforcementLearning.jl is to define the problem you want to solve in a recognizable way. Here we'll demonstrate how to write many different kinds of environments based on interfaces defined in `ReinforcementLearningBase.jl`

The most commonly used interfaces to describe reinforcement learning tasks is [OpenAI/Gym](https://gym.openai.com/). Inspired by it, we expand those interfaces a little to utilize the multiple-dispatch in Julia and to cover multi-agent environments.

## The Minimal Interfaces to Implement

Many interfaces in [ReinforcementLearningBase.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearningBase.jl) have a default implementation. So in most cases, you only need to implement the following functions to define a customized environment:

```julia
action_space(env::YourEnv)
state(env::YourEnv)
state_space(env::YourEnv)
reward(env::YourEnv)
is_terminated(env::YourEnv)
reset!(env::YourEnv)
(env::YourEnv)(action)
```

## An Example: The LotteryEnv

Here we use an example introduced in [Monte Carlo Tree Search: A Tutorial](https://www.informs-sim.org/wsc18papers/includes/files/021.pdf) to demonstrate how to write a simple environment.

The game is defined like this: assume you have \$10 in your pocket, and you are faced with the following three choices:

1. Buy a PowerRich lottery ticket (win \$100M w.p. 0.01; nothing otherwise);
1. Buy a MegaHaul lottery ticket (win \$1M w.p. 0.05; nothing otherwise);
1. Do not buy a lottery ticket.

This game is a one-shot game. It terminates immediately after taking an action
and a reward is received. First we define a concrete subtype of `AbstractEnv`
named `LotteryEnv`:
"""

# ╔═╡ d4d6a2e0-6099-11eb-10ad-dd67c607ea0a
Base.@kwdef mutable struct LotteryEnv <: AbstractEnv
    reward::Union{Nothing, Int} = nothing
end

# ╔═╡ dee89ee4-6099-11eb-2c5a-9d05e6c4eb86
md"""
`LotteryEnv` has only one field named `reward`, by default it is initialized with `nothing`. Now let's implement the necessary interfaces:
"""

# ╔═╡ ef142626-6099-11eb-1d25-e1278a8c393e
RLBase.action_space(env::LotteryEnv) = (:PowerRich, :MegaHaul, nothing)

# ╔═╡ f02be148-6099-11eb-0b1b-95bd6aab04cb
md"""
Here `RLBase` is just an alias for `ReinforcementLearningBase`.
"""

# ╔═╡ 0f0ec27e-609a-11eb-3557-031e6004e78a
begin
	RLBase.reward(env::LotteryEnv) = env.reward
	RLBase.state(env::LotteryEnv) = !isnothing(env.reward)
	RLBase.state_space(env::LotteryEnv) = [false, true]
	RLBase.is_terminated(env::LotteryEnv) = !isnothing(env.reward)
	RLBase.reset!(env::LotteryEnv) = env.reward = nothing
end

# ╔═╡ 1ba95468-609a-11eb-1678-8bd0f3906606
md"""
Because the lottery game is just a simple one-shot game. If the `reward` is
`nothing` then the game is not started yet and we say the game is in state
`false`, otherwise the game is terminated and the state is `true`. So the result
of `state_space(env)` describes the possible states of this environment. By `reset!`
the game, we simply assign the reward with `nothing`, meaning that it's in the
initial state again.

The only left one is to implement the game logic:
"""

# ╔═╡ 52fd1580-609a-11eb-3de2-ed1e5ad95005
function (x::LotteryEnv)(action)
    if action == :PowerRich
        x.reward = rand() < 0.01 ? 100_000_000 : -10
    elseif action == :MegaHaul
        x.reward = rand() < 0.05 ? 1_000_000 : -10
    elseif isnothing(action)
        x.reward = 0
    else
        @error "unknown action of $action"
    end
end

# ╔═╡ 57f86f42-609a-11eb-19af-2d1051ad81dd
md"""

## Test Your Environment

A method named `RLBase.test_runnable!` is provided to rollout several
simulations and see whether the environment we defined is functional.
"""

# ╔═╡ 5d8d1c48-609a-11eb-1fb4-a979b8aca3dc
env = LotteryEnv()

# ╔═╡ 6307f9e0-609a-11eb-1f92-c594d4aabeae
RLBase.test_runnable!(env)

# ╔═╡ 894620fa-609a-11eb-2665-93f5393ba35e
md"""
It is a simple smell test which works like this:

```
for _ in 1:n_episode
    reset!(env)
    while !is_terminated(env)
        env |> action_space |> rand |> env
    end
end
```
"""

# ╔═╡ 9a6f9636-609a-11eb-2357-9d10f0a52d7a
md"""
One step further is to test that other components in ReinforcementLearning.jl
also work. Similar to the test above, let's try the `RandomPolicy` first:
"""

# ╔═╡ b67979f0-609a-11eb-1efd-bb1407e6affc
run(RandomPolicy(action_space(env)), env, StopAfterEpisode(1_000))

# ╔═╡ fc508c3e-609a-11eb-2b19-df3956c2fc7d
md"""
If no error shows up, then it means our environment at least works with the
`RandomPolicy` 🎉🎉🎉. Next, we can add a hook to collect the reward in each
episode to see the performance of the `RandomPolicy`.
"""

# ╔═╡ 04299018-609b-11eb-185e-83b941aa7462
begin
	hook = TotalRewardPerEpisode()
	run(RandomPolicy(action_space(env)), env, StopAfterEpisode(1_000), hook)
	plot(hook.rewards)
end

# ╔═╡ a07eed5a-609b-11eb-2bbb-f318c9a796b5
md"""
A random policy is usually not very meaningful. Here we'll use a tabular based
monte carlo method to estimate the state-action value. (You may choose
appropriate algorithms based on the problem you're dealing with.)
"""

# ╔═╡ b443632a-609b-11eb-3976-5d199171c779
p = QBasedPolicy(
    learner = MonteCarloLearner(;
            approximator=TabularQApproximator(
                ;n_state = length(state_space(env)),
                n_action = length(action_space(env)),
                opt = InvDecay(1.0)
            )
        ),
    explorer = EpsilonGreedyExplorer(0.1)
)

# ╔═╡ b9ebd5aa-609b-11eb-11cc-739730e6208e
p(env)

# ╔═╡ be7c6ac6-609b-11eb-2b32-8dd208db1c57
md"""
Oops, we get an error here. So what does it mean?

Before answering this
question, let's spend some time on understanding the policy we defined above.
A
[`QBasedPolicy`](https://juliareinforcementlearning.org/ReinforcementLearning.jl/latest/rl_core/#ReinforcementLearningCore.QBasedPolicy)
contains two parts: a `learner` and an `explorer`. The `learner` *learn* the
state-action value function (aka *Q* function) duiring interactions with the
`env`. The `explorer` is used to select an action based on the Q value returned
by the `learner`. Here the [`EpsilonGreedyExplorer(0.1)`](https://juliareinforcementlearning.org/ReinforcementLearning.jl/latest/rl_core/#ReinforcementLearningCore.EpsilonGreedyExplorer) will select the action of the
largest value with probability `0.9` and select a random one with probability `0.1`. Inside of the `MonteCarloLearner`, a `TabularQApproximator` is
used to estimate the Q value.

That's the problem! A `TabularQApproximator` only
accepts states of type `Int`.
"""

# ╔═╡ e408057c-609b-11eb-27d5-ff2c082e98c6
p.learner.approximator(1, 1)  # Q(s, a)

# ╔═╡ ea06f580-609b-11eb-1828-ddc0e5744cc5
p.learner.approximator(1)     # [Q(s, a) for a in action_space(env)]

# ╔═╡ ed393fee-609b-11eb-1706-830287c2fbc7
p.learner.approximator(false)

# ╔═╡ f0eabad0-609b-11eb-0422-25937841f8e6
md"""
OK, now we know where the problem is. But how to fix it?

A initial idea is to
rewrite the `RLBase.state(env::LotteryEnv)` function to force it return an
`Int`. That's workable. But in some cases, we may be using environments written
by others and it's not very easy to modify the code directly. Fortunatelly, some
[built-in wrappers](https://github.com/JuliaReinforcementLearning/ReinforcementLearningEnvironments.jl/tree/master/src/environments/wrappers) are provided to help us transform the environment.
"""

# ╔═╡ 03eebcc6-609c-11eb-352c-398708e60ec3
wrapped_env = ActionTransformedEnv(
    StateOverriddenEnv(
        env,
        s -> s ? 1 : 2
    ),
    action_space_mapping = _ -> Base.OneTo(3),
    action_mapping = i -> action_space(env)[i]
)

# ╔═╡ 13535078-609c-11eb-0a3a-873d4d82a3b4
p(wrapped_env)

# ╔═╡ 3ad3704c-609c-11eb-15e5-0709c3246396
md"""
Nice job! Now we are ready to run the experiment:
"""

# ╔═╡ 4108c73c-609c-11eb-19ff-af4e1641eb8a
begin
	h = TotalRewardPerEpisode()
	run(p, wrapped_env, StopAfterEpisode(1_000), h)
	plot(h.rewards)
end

# ╔═╡ 59f80b90-609c-11eb-30ef-776f1154d897
md"""
If you are observant enough, you'll find that our policy is not updating at all!!!
"""

# ╔═╡ ab96c892-609c-11eb-08b4-e3f1d1d70049
p.learner.approximator.table

# ╔═╡ b63d5192-609c-11eb-1bdc-21dd1f56cf22
md"""
Well, actually the policy is running in the **evaluation** mode here. We'll explain it in another blog. For now, you only need to know that we can wrap the policy in an `Agent` to train the policy.
"""

# ╔═╡ c4e97848-609d-11eb-0e9d-dda4843007bd
agent = Agent(;
	policy=p,
	trajectory=VectorSARTTrajectory()
)

# ╔═╡ d6575e4c-609d-11eb-01be-63f485339089
new_hook = TotalRewardPerEpisode()

# ╔═╡ a4a1e976-609d-11eb-01aa-b32b0af8e34d
run(agent, wrapped_env, StopAfterStep(100_000), new_hook)

# ╔═╡ 6c707936-609e-11eb-02dd-f14440c00c2f
p.learner.approximator.table

# ╔═╡ 642b61fc-609c-11eb-23a4-25002f4a0bd3
md"""
!!! note
	Always remember that each algorithm usually only works in some specific environments, just like the `QBasedPolicy` above. So choose the right tool wisely 😉.
"""

# ╔═╡ c3b1dde0-60a6-11eb-214d-b7c966400eae
md"""
## More Complicated Environments

The above `LotteryEnv` is quite simple. Many environments we are interested in fall in the same category. Beyond that, there're still many other kinds of environments. You may take a glimpse at the [table](https://github.com/JuliaReinforcementLearning/ReinforcementLearningEnvironments.jl#built-in-environments) to see how many different types of environments are supported in [ReinforcementLearningZoo.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearningEnvironments.jl).

To distinguish different kinds of environments, some common traits are defined in [ReinforcementLearningBase.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearningBase.jl). Now we'll explain them one-by-one.

### `StateStyle`

In the above `LotteryEnv`, `state(env::LotteryEnv)` simply returns a `true` or `false`. But in some other environments, the function name `state` may be kind of vague. People with different background often talk about the same thing with different names. You may be interested in this discussion: [What is the difference between an observation and a state in reinforcement learning?](https://ai.stackexchange.com/questions/5970/what-is-the-difference-between-an-observation-and-a-state-in-reinforcement-learn) To avoid confusion when executing `state(env)`, the environment designer can explicitly define `state(::AbstractStateStyle, env::YourEnv)`. So that users can fetch necessary information on demand. Following are some built-in state styles:
"""

# ╔═╡ bd9cf900-60a9-11eb-1b11-8bc74ae5da60
subtypes(RLBase.AbstractStateStyle)

# ╔═╡ acf33c12-60af-11eb-00a8-b3e6a35dd256
md"""
Note that every state style may have different representations, `String`, `Array`, `Graph` and so on. All the above state styles can accept a data type as parameter. For example:
"""

# ╔═╡ 12fa94ec-60b0-11eb-0bd7-a5da4368db56
RLBase.state(::Observation{String}, env::LotteryEnv) = is_terminated(env) ? "Game Over" : "Game Start"

# ╔═╡ 43d13376-60b0-11eb-3f05-8f50e494764b
md"""
For environments which support many different kinds of states, developers should specify all the supported state styles. For example:
"""

# ╔═╡ e06eccf4-60b5-11eb-1511-31e0795dc87b
tp = TigerProblemEnv();

# ╔═╡ ccffa9cc-60b5-11eb-140c-676065b1c34f
StateStyle(tp)

# ╔═╡ e7564254-60b5-11eb-24e4-af25d66f7eba
state(tp, Observation{Int64}())

# ╔═╡ 4ef7543e-60b6-11eb-2883-735a555b0244
state(tp, InternalState{Int64}())

# ╔═╡ 56a92700-60b6-11eb-1af1-6500aa59fee2
state(tp)

# ╔═╡ 5cfd5da8-60b6-11eb-014d-e14fec02fd5d
DefaultStateStyle(tp)

# ╔═╡ 78fddc1c-60b6-11eb-1266-75360b57e742
md"""
### `DefaultStateStyle`

The `DefaultStateStyle` trait returns the first element in the result of `StateStyle` by default.

For algorithm developers, they usually don't care about the state style. They can assume that the default state style is always well defined and simply call `state(env)` to get the right representation. So for environments of many different representations, `state(env)` will be dispatched to `state(DefaultStateStyle(env), env)`. And we can use the [`DefaultStateStyleEnv`](https://juliareinforcementlearning.org/ReinforcementLearning.jl/latest/rl_envs/#ReinforcementLearningEnvironments.DefaultStateStyleEnv-Union{Tuple{E},%20Tuple{S}}%20where%20E%20where%20S) wrapper to override the pre-defined `DefaultStateStyle(::YourEnv)`.
"""

# ╔═╡ 8ddb2a40-60b6-11eb-1396-d9880e304217
md"""
### `RewardStyle`

For games like Chess, Go or many card game, we only get the reward at the end of an game. We say this kind of games is of `TerminalReward`, otherwise we define it as `StepReward`. Actually the `TerminalReward` is a special case of `StepReward` (for non-terminal steps, the reward is `0`). The reason we still want to distinguish these two cases is that, for some algorithms there may be a more efficient implementation for `TerminalReward` style games.

"""

# ╔═╡ 97cba458-60b6-11eb-2b9e-b5ca1f9a5c73
RewardStyle(tp)

# ╔═╡ b2dab6da-60b6-11eb-275b-57d0a3e65da8
RewardStyle(MontyHallEnv())

# ╔═╡ b8e097f8-60b5-11eb-27ab-7f73e2d7c502
md"""
### `ActionStyle`

For some environments, the valid actions in each step may be different. We call this kind of environments are of `FullActionSet`. Otherwise, we say the environment is of `MinimalActionSet`. A typical built-in environment with `FullActionSet` is the `TicTacToeEnv`. Two extra methods must be implemented:
"""

# ╔═╡ 81614ebc-60b5-11eb-15b9-897cdcbde04f
ttt = TicTacToeEnv();

# ╔═╡ a1c73c42-60b3-11eb-3db7-d71785bb4370
ActionStyle(ttt)

# ╔═╡ 6127b88e-60b5-11eb-2bbc-4d078ce6b44e
legal_action_space(ttt)

# ╔═╡ 8a0e6db0-60b5-11eb-3263-295fc6a951a2
legal_action_space_mask(ttt)

# ╔═╡ 9558df16-60b5-11eb-3fbd-77c38bbe6ca4
md"""
### `NumAgentStyle`

In the above `LotteryEnv`, only one player is involved in the environment. In many board games, usually multiple players are engaged.
"""

# ╔═╡ 4060756c-60b7-11eb-3e2f-a59544eae5c8
NumAgentStyle(env)

# ╔═╡ 4d36f23e-60b7-11eb-39be-494f4cb0fa24
NumAgentStyle(ttt)

# ╔═╡ 566c4a16-60b7-11eb-2ced-0b04dbc2953d
md"""
For multi-agent environments, some new APIs are introduced. The meaning of some APIs we've seen are also extended.

First, multi-agent environment developers must implement `players` to distinguish different players.
"""

# ╔═╡ 17cbc880-60b8-11eb-0b27-d3cbdadb4093
players(ttt)

# ╔═╡ 4ade9f26-60b9-11eb-3a64-497b0ab2c25f
current_player(ttt)

# ╔═╡ 34e6633a-60b8-11eb-375f-cd0093097caa
md"""

| Single Agent | Multi-Agent |
| ------------:| -----------:|
| `state(env)` | `state(env, player)`|
| `reward(env)`| `reward(env, player)`|
| `env(action)`| `env(action, player)`|
| `action_space(env)`| `action_space(env, player)`|
| `state_space(env)`| `state_space(env, player)`|
| `is_terminated(env)` | `is_terminated(env, player)`|

Note that the APIs in single agent is still valid, only that they all fall back to the perspective from the `current_player(env)`.
"""

# ╔═╡ 61f6fd5c-60b9-11eb-1062-f1aff5da84af
md"""
#### `UtilityStyle`

In multi-agent environments, sometimes the sum of rewards from all players are always `0`. We call the `UtilityStyle` of these environments `ZeroSum`. `ZeroSum` is a special case of `ConstantSum`. In cooperational games, the reward of each player are the same. In this case, they are called `IdenticalUtility`. Other cases fall back to `GeneralSum`.
"""

# ╔═╡ 44cad7c0-60ba-11eb-2570-b1da3d16c973
md"""
#### `InformationStyle`

If all players can see the same state, then we say the `InformationStyle` of these environments are of `PerfectInformation`. They are a special case of `ImperfectInformation` environments.
"""

# ╔═╡ 80f3f7a4-60ba-11eb-301f-71b73ece720f
md"""
#### `DynamicStyle`

All the environments we've seen so far were of `Sequential` style, meaning that at each step, only **ONE** player was allowed to take an action. Alternatively there are `Simultaneous` environments, where all the players take actions simultaneously without seeing each other's action in advance. Simultaneous environments must take a collection of actions from different players as input.
"""

# ╔═╡ 6a58f464-60bb-11eb-12ab-592202a61d9f
rps = RockPaperScissorsEnv();

# ╔═╡ 5cf93fca-60bb-11eb-1376-37b9322258d4
action_space(rps)

# ╔═╡ 668165fe-60bb-11eb-0b0b-bf0fc2f2e36a
rps(rand(action_space(rps)))

# ╔═╡ 7e056cde-60bb-11eb-17e6-3bcefeeae0b0
md"""
#### `ChanceStyle`

If there's no `rng` in the environment, everything is deterministic afer taking each action, then we call the `ChanceStyle` of these environments are of `Deterministic`. Otherwise, we call them `Stochastic`. One special case is that, in [Extensive Form Games](https://en.wikipedia.org/wiki/Extensive-form_game), a chance node is envolved. And the action probability of this special player is known. For these environments, we need to have the following methods defined:
"""

# ╔═╡ 2b203eb2-60bc-11eb-30e6-9beaf61fe7a9
kp = KuhnPokerEnv();

# ╔═╡ 48707c20-60bc-11eb-2ac4-0da804dd9381
chance_player(kp)

# ╔═╡ 426a9bb2-60bc-11eb-1efd-f385e77e9c22
prob(kp, chance_player(kp))

# ╔═╡ 51cdadce-60bc-11eb-01f7-e1e9b6362ad3
chance_player(kp) in players(kp)

# ╔═╡ 6b9e0564-60bc-11eb-196e-09b9b8f557ba
md"""
## Examples

Finally we've gone through all the details you need to know for how to write a customized environment. You're encouraged to take a look at the examples provided in [ReinforcementLearningEnvironments.jl](https://github.com/JuliaReinforcementLearning/ReinforcementLearningEnvironments.jl). Feel free to create an issue there if you're still not sure how to describe your problem with the interfaces defined in this package.
"""

# ╔═╡ Cell order:
# ╟─bed48ef6-62fb-11eb-327a-ed29787dda3c
# ╟─dccad7c8-62fb-11eb-226f-393c73301bcb
# ╟─e262295a-62fb-11eb-32c0-051c8995d14a
# ╟─03852aa2-6099-11eb-39db-ffba5eeade98
# ╠═5938a9a6-6099-11eb-1ac9-1fb12d9c9237
# ╠═d4d6a2e0-6099-11eb-10ad-dd67c607ea0a
# ╟─dee89ee4-6099-11eb-2c5a-9d05e6c4eb86
# ╠═ef142626-6099-11eb-1d25-e1278a8c393e
# ╟─f02be148-6099-11eb-0b1b-95bd6aab04cb
# ╠═0f0ec27e-609a-11eb-3557-031e6004e78a
# ╟─1ba95468-609a-11eb-1678-8bd0f3906606
# ╠═52fd1580-609a-11eb-3de2-ed1e5ad95005
# ╟─57f86f42-609a-11eb-19af-2d1051ad81dd
# ╠═5d8d1c48-609a-11eb-1fb4-a979b8aca3dc
# ╠═6307f9e0-609a-11eb-1f92-c594d4aabeae
# ╟─894620fa-609a-11eb-2665-93f5393ba35e
# ╟─9a6f9636-609a-11eb-2357-9d10f0a52d7a
# ╠═ad1302aa-609a-11eb-0a22-892b5c80750d
# ╠═b67979f0-609a-11eb-1efd-bb1407e6affc
# ╟─fc508c3e-609a-11eb-2b19-df3956c2fc7d
# ╠═06013422-609b-11eb-24e0-790eb0272183
# ╠═04299018-609b-11eb-185e-83b941aa7462
# ╟─a07eed5a-609b-11eb-2bbb-f318c9a796b5
# ╠═ac7e107c-609b-11eb-2d50-4b50517e1840
# ╠═b443632a-609b-11eb-3976-5d199171c779
# ╠═b9ebd5aa-609b-11eb-11cc-739730e6208e
# ╟─be7c6ac6-609b-11eb-2b32-8dd208db1c57
# ╠═e408057c-609b-11eb-27d5-ff2c082e98c6
# ╠═ea06f580-609b-11eb-1828-ddc0e5744cc5
# ╠═ed393fee-609b-11eb-1706-830287c2fbc7
# ╟─f0eabad0-609b-11eb-0422-25937841f8e6
# ╠═03eebcc6-609c-11eb-352c-398708e60ec3
# ╠═13535078-609c-11eb-0a3a-873d4d82a3b4
# ╟─3ad3704c-609c-11eb-15e5-0709c3246396
# ╠═4108c73c-609c-11eb-19ff-af4e1641eb8a
# ╟─59f80b90-609c-11eb-30ef-776f1154d897
# ╠═ab96c892-609c-11eb-08b4-e3f1d1d70049
# ╟─b63d5192-609c-11eb-1bdc-21dd1f56cf22
# ╠═c4e97848-609d-11eb-0e9d-dda4843007bd
# ╠═d6575e4c-609d-11eb-01be-63f485339089
# ╠═a4a1e976-609d-11eb-01aa-b32b0af8e34d
# ╠═6c707936-609e-11eb-02dd-f14440c00c2f
# ╟─642b61fc-609c-11eb-23a4-25002f4a0bd3
# ╟─c3b1dde0-60a6-11eb-214d-b7c966400eae
# ╠═bd9cf900-60a9-11eb-1b11-8bc74ae5da60
# ╟─acf33c12-60af-11eb-00a8-b3e6a35dd256
# ╠═12fa94ec-60b0-11eb-0bd7-a5da4368db56
# ╟─43d13376-60b0-11eb-3f05-8f50e494764b
# ╠═e06eccf4-60b5-11eb-1511-31e0795dc87b
# ╠═ccffa9cc-60b5-11eb-140c-676065b1c34f
# ╠═e7564254-60b5-11eb-24e4-af25d66f7eba
# ╠═4ef7543e-60b6-11eb-2883-735a555b0244
# ╠═56a92700-60b6-11eb-1af1-6500aa59fee2
# ╠═5cfd5da8-60b6-11eb-014d-e14fec02fd5d
# ╟─78fddc1c-60b6-11eb-1266-75360b57e742
# ╟─8ddb2a40-60b6-11eb-1396-d9880e304217
# ╠═97cba458-60b6-11eb-2b9e-b5ca1f9a5c73
# ╠═b2dab6da-60b6-11eb-275b-57d0a3e65da8
# ╟─b8e097f8-60b5-11eb-27ab-7f73e2d7c502
# ╠═81614ebc-60b5-11eb-15b9-897cdcbde04f
# ╠═a1c73c42-60b3-11eb-3db7-d71785bb4370
# ╠═6127b88e-60b5-11eb-2bbc-4d078ce6b44e
# ╠═8a0e6db0-60b5-11eb-3263-295fc6a951a2
# ╟─9558df16-60b5-11eb-3fbd-77c38bbe6ca4
# ╠═4060756c-60b7-11eb-3e2f-a59544eae5c8
# ╠═4d36f23e-60b7-11eb-39be-494f4cb0fa24
# ╟─566c4a16-60b7-11eb-2ced-0b04dbc2953d
# ╠═17cbc880-60b8-11eb-0b27-d3cbdadb4093
# ╠═4ade9f26-60b9-11eb-3a64-497b0ab2c25f
# ╟─34e6633a-60b8-11eb-375f-cd0093097caa
# ╟─61f6fd5c-60b9-11eb-1062-f1aff5da84af
# ╟─44cad7c0-60ba-11eb-2570-b1da3d16c973
# ╟─80f3f7a4-60ba-11eb-301f-71b73ece720f
# ╠═6a58f464-60bb-11eb-12ab-592202a61d9f
# ╠═5cf93fca-60bb-11eb-1376-37b9322258d4
# ╠═668165fe-60bb-11eb-0b0b-bf0fc2f2e36a
# ╟─7e056cde-60bb-11eb-17e6-3bcefeeae0b0
# ╠═2b203eb2-60bc-11eb-30e6-9beaf61fe7a9
# ╠═48707c20-60bc-11eb-2ac4-0da804dd9381
# ╠═426a9bb2-60bc-11eb-1efd-f385e77e9c22
# ╠═51cdadce-60bc-11eb-01f7-e1e9b6362ad3
# ╟─6b9e0564-60bc-11eb-196e-09b9b8f557ba


text:
"""
Original publication:
A. Ghasemi, M. T. Manzuri, H. R. Rabiee, M. H.
Rohban, and S. Haghiri. Active one-class learning by
kernel density estimation. In 2011 IEEE International
Workshop on Machine Learning for Signal Processing,
pages 1–6, Sept 2011.
"""
struct MinimumLossQs <: DataBasedQs
    p_x::Array{Float64, 1}
    bw_method::Union{String, U} where U <: Real
    p_inlier::Float64
    function MinimumLossQs(x::Array{T, 2}; bw_method="scott", p_inlier=nothing) where T <: Real
        ((p_inlier === nothing) || !(0 <= p_inlier <= 1)) && throw(ArgumentError("Invalid inlier probability $(p_inlier)."))
        return new(multi_kde(x, bw_method)(x), bw_method, p_inlier)
    end
end

function qs_score(qs::MinimumLossQs, x::Array{T, 2}, labels::Dict{Symbol, Array{Int, 1}})::Array{Float64, 1} where T <: Real
    haskey(labels, :Lin) || throw(MissingLabelTypeException(:Lin))
    s_t_a(u) = leave_out_one_cv_kde(hcat(x[:, labels[:Lin]], u), qs.bw_method)
    s_t_b(u) = haskey(labels, :Lout) ? mean(multi_kde(hcat(x[:, labels[:Lin]], u), qs.bw_method)(x[:, labels[:Lout]])) : 0
    s_o_a = leave_out_one_cv_kde(x[:, labels[:Lin]], qs.bw_method)
    s_o_b(u) = haskey(labels, :Lout) ? mean(multi_kde(x[:, labels[:Lin]], qs.bw_method)(hcat(x[:, labels[:Lout]], u))) : 0
    s_t(u) = s_t_a(u) - s_t_b(u)
    s_o(u) = s_o_a - s_o_b(u)
    p_inlier = length(labels[:Lin]) / size(x, 2)
    s(u) = qs.p_inlier * s_t(u) + (1 - qs.p_inlier) * s_o(u)
    return [s(x[:, i]) for i in 1:size(x, 2)]
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule Minuit2_Julia_Wrapper_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("Minuit2_Julia_Wrapper")
JLLWrappers.@generate_main_file("Minuit2_Julia_Wrapper", UUID("65cff058-67fb-5034-9600-a8f6a0feb90b"))
end  # module Minuit2_Julia_Wrapper_jll


text:
module PolyhedralRelaxations

using DataStructures
using JuMP

import ForwardDiff

import Logging
import LoggingExtras

# Setup Logging
include("logging.jl")
function __init__()
    global _DEFAULT_LOGGER = Logging.current_logger()
    global _LOGGER = Logging.ConsoleLogger(;
        meta_formatter = PolyhedralRelaxations._pr_metafmt,
    )

    return Logging.global_logger(_LOGGER)
end

const EPS = 1e-6
const INF = 1e12

include("common.jl")
include("univariate_lp_relaxation.jl")
include("univariate_milp_relaxation.jl")
include("bilinear_relaxation.jl")
include("api.jl")

end


text:
module Hello

using ArgParse
using PyCall
using PyPlot

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--opt1"
            help = "an option with an argument"
        "--opt2", "-o"
            help = "another option with an argument"
            arg_type = Int
            default = 0
        "--flag1"
            help = "an option without argument, i.e. a flag"
            action = :store_true
        "arg1"
            help = "a positional argument"
            required = true
    end

    return parse_args(s)
end

function main()
    println(pyimport("sys")[:executable])
    plot([1,2,3])
    savefig("saved.png")
    parsed_args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in parsed_args
        println("  $arg  =>  $val")
    end
end

Base.@ccallable function julia_main(ARGS::Vector{String})::Cint
    main()
    return 0
end

end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
using DataStructures
using Base.Order
using ResumableFunctions
import Base.length

mutable struct HierarchicalSoftmaxNode
	parent::Int32
	branch::Bool
end

struct HierarchicalOutput
	code::Array{Int8}
	path::Array{Int}
end

length(out::HierarchicalOutput) = length(out.path)

function HierarchicalSoftmaxNode()
	return HierarchicalSoftmaxNode(Int32(0), false)
end

@resumable function softmax_path(nodes::Array{HierarchicalSoftmaxNode},
		V::Integer, id::Integer) :: Tuple{Int32, Int8}
	while true
		node = nodes[id]
		if node.parent == 0 break; end
		@assert node.parent > V
		@yield (Int32(node.parent - V), Int8(node.branch))
		id = node.parent
	end
end

function build_huffman_tree(freqs::Array{Tf}) where {Tf <: Number}
	V = length(freqs)
	nodes = Array{HierarchicalSoftmaxNode}(undef, V)
	for v in 1:V
		nodes[v] = HierarchicalSoftmaxNode()
	end

	freq_ord = By(wf -> wf[2])
	heap = heapify!([(nodes[v], freqs[v]) for v in 1:V], freq_ord)

	function pop_initialize!(parent::Int, branch::Bool)
		node = heappop!(heap, freq_ord)
		node[1].parent = Int32(parent)
		node[1].branch = branch
		return node[2]
	end

	L = V
	while length(heap) > 1
		L += 1
		node = HierarchicalSoftmaxNode()
		push!(nodes, node)

		freq = 1
		freq = pop_initialize!(L, true) + pop_initialize!(L, false)
		heappush!(heap, (node, freq), freq_ord)
	end

	@assert(length(heap) == 1, string(heap))

	return nodes
end

function convert_huffman_tree(nodes::Array{HierarchicalSoftmaxNode}, V::Integer)
	outputs = Array{HierarchicalOutput, 1}(undef, V)
	for v in 1:V
		code = Array{Int8, 1}()
		path = Array{Int, 1}()

		for (n, branch) in softmax_path(nodes, V, v)
			push!(code, branch)
			push!(path, n)
		end

		outputs[v] = HierarchicalOutput(code, path)
	end

	return outputs
end

export HierarchicalSoftmaxNode, softmax_path


text:
include("util.jl")

recipe1 = LinkedListNode(3, nothing, nothing)
recipe2 = LinkedListNode(7, recipe1, nothing)
recipe1.next = recipe2
recipes = LinkedList(recipe1, recipe2, 2)

function get_digits(x::Int64)::Array{Int64}
    digits = string(x)
    return [parse(Int64, c) for c in digits]
end

elf1 = recipe1
elf2 = recipe2
function dostep(elf1, elf2)
    combined_score = elf1.data+elf2.data
    new_recipes = get_digits(combined_score)
    for digit in new_recipes
        push!(recipes, digit)
    end
    elf1_movement = elf1.data+1
    elf2_movement = elf2.data+1
    for i in 1:(elf1_movement)
        if elf1.next != nothing
            elf1 = elf1.next
        else
            elf1 = recipes.head
        end
    end
    for i in 1:(elf2_movement)
        if elf2.next != nothing
            elf2 = elf2.next
        else
            elf2 = recipes.head
        end
    end
    return elf1,elf2
end

part_1 = 765071
while recipes.length < (part_1+10)
    global elf1, elf2
    elf1,elf2 = dostep(elf1, elf2)
end

a = recipes[part_1+1:part_1+10]

function stringatize(l)
    return join([string(c) for c in l])
end

function check_for(recipes, part_2)
    l = length(part_2)
    if recipes.length <= l
        return false
    end
    last_l = stringatize(recipes[-l:-1])
    next_last_l = stringatize(recipes[-l-1:-2])
    if last_l == part_2
        println(recipes.length-l)
        return true
    elseif next_last_l == part_2
        println(recipes.length-l-1)
        return true
    else
        return false
    end
end

part_2 = "765071"
c = 0
while !(check_for(recipes, part_2))
    global elf1, elf2, c
    c += 1
    if c % 100000 == 0
        println(c)
    end
    elf1,elf2 = dostep(elf1, elf2)
end


text:
if !("." in LOAD_PATH)
    push!(LOAD_PATH, ".")
end
include("common.jl")
import BenchmarkTools

sim = init_sim("Bomex")
tendencies = copy(sim.state.prog)
update_n(sim, tendencies, 1) # compile first
trial = BenchmarkTools.@benchmark update_n($sim, $tendencies, 1)
show(stdout, MIME("text/plain"), trial)
println()

sim = init_sim("TRMM_LBA")
tendencies = copy(sim.state.prog)
update_n(sim, tendencies, 1) # compile first
trial = BenchmarkTools.@benchmark update_n($sim, $tendencies, 1)
show(stdout, MIME("text/plain"), trial)
println()

sim = init_sim("Rico")
tendencies = copy(sim.state.prog)
update_n(sim, tendencies, 1) # compile first
trial = BenchmarkTools.@benchmark update_n($sim, $tendencies, 1)
show(stdout, MIME("text/plain"), trial)
println()


text:
module InterpretMe

include("pdp.jl")
export partial_dependence
end # module


text:
struct Env
    paths::Dict
end

mutable struct Project
    name::String
    dat::Dict{String,Any}
    tileRegister::DataFrame
    Project() = new()
end

mutable struct OsmFeatures
    OsmId::Int64
    OSG_UPRNs::Array{Int64}
    BlockId_HA::String
    featureType::String
    nodes::Array{Int64}
    LatLon::Array{Float64}
    centre::Array{Float64}
    Orientation::Float64
    w::Float16
    d::Float16
    tp::Float16
    tol::Float16
    tor::Float16
    area::Float16
    flag::String
    arch::String
    shape::String
end

mutable struct Osm
    region::String
    ways
    nodeList::Array{Int64}
    nodes::Array{XMLElement,1}
    ind::Array{Int64}
    frame::Array{Float64}
    Osm() = new()
end

mutable struct Toggle
    bldgs::Widget{:toggle,Bool}
    roads::Widget{:toggle,Bool}
    coast::Widget{:toggle,Bool}
    postcodes::Widget{:toggle,Bool}
    Toggle() = new()
end
mutable struct Dropdown
    markerfield::Widget{:dropdown,Any}
    zone::Widget{:dropdown,Any}
    orient::Widget{:dropdown,Any}
    Dropdown() = new()
end
mutable struct Radio
    spare1::Widget{:radiobuttons,Any}
    Radio() = new()
end
mutable struct Button
    clr::Widget{:button,Int64}
    generateArch::Widget{:button,Int64}
    #selectArch::Widget{:button,Int64}
    generateIdf::Widget{:button,Int64}
    Button() = new()
end
mutable struct Spinbox
    dimensionalRounding::Widget{:spinbox,Union{Nothing, Int64}}
    maxPlanDepth::Widget{:spinbox,Union{Nothing, Int64}}
    Spinbox() = new()
end
mutable struct Check
    census::Widget{:checkbox,Bool}
    epc::Widget{:checkbox,Bool}
    homeAnalytics::Widget{:checkbox,Bool}
    acorn::Widget{:checkbox,Bool}
    Check() = new()
end

mutable struct Wgt
    toggle::Toggle
    dropdown::Dropdown
    radio::Radio
    button::Button
    check::Check
    spinbox::Spinbox
end

mutable struct UI
    title::HTML{String}
    titlePane::Node{WebIO.DOM}
    mainWindow::Window
    subPanelLeft1::Node{WebIO.DOM}
    panelLeft::Node{WebIO.DOM}
    panelRight::Node{WebIO.DOM}
    p0::PlotlyJS.SyncPlot
    plt_macro::PlotlyJS.SyncPlot
    macro_point::Tuple{Float64,Float64}
    wgt::Wgt
    UI() = new()
end

struct UI_launch
    title::HTML{String}
    window::Window
    wgt_dropdown::Widget{:dropdown,Any}
    wgt_button_go::Widget{:button,Int64}
    wgt_button_newProject::Widget{:button,Int64}
    wgt_textbox_newProject::Widget{:textbox,String}
end
mutable struct GIS
    datum
    origin
    transf_LLAfromENU
    transf_ENUfromLLA
    GIS() = new()
end

mutable struct Data
    Census::DataFrame
    EST::DataFrame
    ByBlock::DataFrame
    OSM::DataFrame
    Acorn::DataFrame
    Data() = new()
end

mutable struct Recentre2
    fix
    tmp_PostcodeGroups
    yscale
    xscale
    iOsmLLA
    BlockId
    focalPt
    refocusId
    frame
    angleSetEST
    angleSetOSM
    BlockSubsetEST
    BlockSubsetOSM
    iNearestOSM
    refocusNudge
    frameExpansion
    refocusLLA
    iOsmFocus
    Recentre2() = new()
end


text:
riskregionnames = ["Machakos/Muranga" "Mandera" "Baringo/Nakuru" "Nairobi" "Uasin Gishu" "Marsabit" "Garissa" "Nakuru/Narok" "Turkana" "Taita Taveta" "Kitui/Meru" "Kilifi/Mombasa" "Kericho/Kisumu" "Kilifi/Lamu" "Kakamega/Kisumu" "Wajir" "Kajiado/Kisumu" "Homa bay/Migori" "Samburu/Laikipia" "Kilifi/Kwale" "Total"]
age_cats = ["0-4","5-9","10-14","15-19","20-24","25-29","30-34","35-39","40-44","45-49","50-54","55-59","60-64","65-69","70-74","75-79","80+"]
rr2county_mat = matread("data/conversion_matrix.mat")["conversion_matrix"]
T = readtable("data/2019-population_census-report-per-county.csv")
county_names = convert.(String,T.County)

using LinearAlgebra
function plot_total_incidence(results_group,treatments::Tuple{Float64,Real},i)
    τ,ϵ_D = treatments
    inc_D1 = results_group[1][i][1]
    inc_D2 = results_group[2][i][1]
    inc_A2 = results_group[2][i][2]
    plt = plot(1:365,inc_D1[21,:,1].+1,
                fontfamily="Helvetica",
                lw = 3,
                legend = :topright,
                ribbon =(inc_D1[21,:,2],inc_D1[21,:,3]),
                # fillalpha = 0.15,
                # xlims = (0.,100),
                yscale = :log10,
                xlabel = "Days",
                yticks = ([1,10,100,1000,10000],["0" "10" "100" "1000" "10000"]),
                ylabel = "Daily incidence")
        # plot!(plt,1:365,inc_D2[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario",
        #             ribbon =(inc_D2[21,:,2],inc_D2[21,:,3]),
        #             fillalpha = 0.15)
        #             # yscale = :log10,
        #             # xlabel = "Days",
        #             # ylabel = "Daily incidence")
        # plot!(plt,1:365,inc_A2[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario: undetected",
        #             ls = :dot,
        #             # ribbon =(inc_D2[21,:,2],inc_D2[21,:,3]),
        #             fillalpha = 0.15)
        #             # yscale = :log10,
        #             # xlabel = "Days",
        #             # ylabel = "Daily incidence")

        # plot!(plt,1:365,inc_D3[21,:,1].+inc_A3[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario, delta = 0.8",
        #             # ribbon =(inc_D3[21,:,2],inc_D3[21,:,3]),
        #             fillalpha = 0.15,
        #             yscale = :log10,
        #             xlabel = "Days",
        #             ylabel = "Daily incidence")
    return plt
end

function plot_total_incidence_group(scenario_group,treatment_group,treatment_num,rel_transmission_perc)
        plt = plot()
        tracing_rate,e_D = treatment_group[treatment_num]
        for (i,scenarioresults) in enumerate(scenario_group)
                inc_D = scenarioresults[treatment_num][1]
                plot!(plt,inc_D[21,:,1].+1,
                fontfamily="Helvetica",
                lw = 3,
                lab= "rel. infect. of subclinical: $(rel_transmission_perc[i])%",
                legend = :topright,
                ribbon =(inc_D[21,:,2],inc_D[21,:,3]),
                fillalpha = 0.15,
                yticks = ([1,10,100,1000,10000],["0" "10" "100" "1000" "10000"]),
                yscale = :log10,
                xlabel = "Days after detecting established CoV transmission",
                ylabel = "Daily incidence")
        end
        return plt
end

function plot_total_incidence_group(plt,scenario_group,treatment_group,treatment_num,rel_transmission_perc)
        tracing_rate,e_D = treatment_group[treatment_num]
        for (i,scenarioresults) in enumerate(scenario_group)
                inc_D = scenarioresults[treatment_num][1]
                plot!(plt,inc_D[21,:,1].+1,
                fontfamily="Helvetica",
                lw = 3,
                lab= "rel. infect. of subclinical: $(rel_transmission_perc[i])%",
                legend = :topright,
                ribbon =(inc_D[21,:,2],inc_D[21,:,3]),
                fillalpha = 0.15,
                yticks = ([1,10,100,1000,10000],["0" "10" "100" "1000" "10000"]),
                yscale = :log10,
                xlabel = "Days after detecting established CoV transmission",
                ylabel = "Daily incidence")
        end
        return plt
end

function plot_total_incidence(results)
    inc_D1 = results[1][1]

    plt = plot(1:365,inc_D1[21,:,1].+1,
                fontfamily="Helvetica",
                lw = 3,
                legend = :topright,
                ribbon =(inc_D1[21,:,2],inc_D1[21,:,3]),
                # fillalpha = 0.15,
                # xlims = (0.,100),
                yscale = :log10,
                xlabel = "Days",
                ylabel = "Daily incidence + 1")
        # plot!(plt,1:365,inc_D2[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario",
        #             ribbon =(inc_D2[21,:,2],inc_D2[21,:,3]),
        #             fillalpha = 0.15)
        #             # yscale = :log10,
        #             # xlabel = "Days",
        #             # ylabel = "Daily incidence")
        # plot!(plt,1:365,inc_A2[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario: undetected",
        #             ls = :dot,
        #             # ribbon =(inc_D2[21,:,2],inc_D2[21,:,3]),
        #             fillalpha = 0.15)
        #             # yscale = :log10,
        #             # xlabel = "Days",
        #             # ylabel = "Daily incidence")

        # plot!(plt,1:365,inc_D3[21,:,1].+inc_A3[21,:,1].+1,
        #             lw = 3,
        #             lab="MERS-like scenario, delta = 0.8",
        #             # ribbon =(inc_D3[21,:,2],inc_D3[21,:,3]),
        #             fillalpha = 0.15,
        #             yscale = :log10,
        #             xlabel = "Days",
        #             ylabel = "Daily incidence")
    return plt
end

function plot_total_incidence(results,treatment_rates)
    inc_D = results[1][1]
    plt = plot(1:365,inc_D[21,:,1].+1,
                fontfamily="Helvetica",
                lw = 3,
                lab="No isolation",
                legend = :topleft,
                ribbon =(inc_D[21,:,2],inc_D[21,:,3]),
                fillalpha = 0.15,
                xlims = (30.,100),
                yscale = :log10,
                xlabel = "Days",
                ylabel = "Daily incidence")
    # for i = 2:3
    #     τ = treatment_rates[i]
    #     inc_D = results[i][1]
    #     plot!(plt,1:365,inc_D[21,:,1].+1,
    #                 lw = 3,
    #                 lab="Isolation in mean $(round(1/τ,digits = 1)) days",
    #                 ribbon =(inc_D[21,:,2],inc_D[21,:,3]),
    #                 fillalpha = 0.15,
    #                 xlims = (30.,100),
    #                 # yscale = :log10,
    #                 xlabel = "Days",
    #                 ylabel = "Daily incidence")
    # end
    return plt
end

function plot_incidence_spatial(results,treatment_rates,i)
    ordering = sortperm(results[i][1][1:20,40,1],rev = true)
    τ = treatment_rates[i]
    inc_D = results[i][1]
    plt = plot()
    for i in ordering
        plot!(plt,inc_D[i,:,1].+1,
                    yticks = ([1,10,100,1000],["0" "10" "100" "1000"]),
                    fillalpha = 0.1,
                    ribbon =(inc_D[i,:,2],inc_D[i,:,3]),
                    lab = riskregionnames[i],
                    lw=2,
                    xlims = (0.,100.),
                    yscale = :log10,
                    legend = :outertopright,
                    legendfontsize = 8.9);
    end
    xlabel!(plt,"Days after detecting established CoV transmission")
    ylabel!(plt,"Daily incidence")

    return plt
end

function plot_incidence_spatial(results,treatment_rates,i,rr2county_mat)
    inc_D_county = rr2county_mat*results[i][1][1:20,:,1]
    lb = rr2county_mat*results[i][1][1:20,:,2]
    ub = rr2county_mat*results[i][1][1:20,:,3]

    ordering = sortperm(inc_D_county[:,40],rev = true)
    τ = treatment_rates[i]
    plt = plot()
    for k in ordering
        plot!(plt,inc_D_county[k,:].+1,
                    yticks = ([1,10,100,1000],["0" "10" "100" "1000"]),
                    fillalpha = 0.05,
                    ribbon =(lb[k,:],ub[k,:]),
                    lab = county_names[k],
                    lw=2,
                    xlims = (0.,100.),
                    yscale = :log10,
                    legend = :outertopright,
                    legendfontsize = 8.9);
    end
    xlabel!(plt,"Days after detecting established CoV transmission")
    ylabel!(plt,"Daily incidence")

    return plt
end

function plot_total_incidence_by_age(scenario_group,treatment_rates,rel_transmission_perc,i)
    d1,d2,d3 = size(scenario_group[1][i][4][:,:,:])
    n = length(scenario_group)
    cases_age = zeros(d2,n,3)
    label_scen = "$(rel_transmission_perc[1])% rel. infectiousness"
    for l in 2:length(rel_transmission_perc)
        label_scen = hcat(label_scen,"$(rel_transmission_perc[l])% rel. infectiousness")
    end
    for (n,results) in enumerate(scenario_group)
        for a = 1:d2
            cases_age[a,n,1] = median([sum(results[i][4][:,a,k]) for k =1:1000])
            cases_age[a,n,2] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.025)
            cases_age[a,n,3] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.975)
        end
    end
    plt = groupedbar(cases_age[:,:,1]./1e3,
                    fontfamily="Helvetica",
                    lab =label_scen,
                    xticks = (1:17,age_cats))

    return plt

end
function plot_total_incidence_by_age(scenario_group,i,labels)
    d1,d2,d3 = size(scenario_group[1][i][4][:,:,:])
    n = length(scenario_group)
    cases_age = zeros(d2,n,3)
    for (n,results) in enumerate(scenario_group)
        for a = 1:d2
            cases_age[a,n,1] = median([sum(results[i][4][:,a,k]) for k =1:1000])
            cases_age[a,n,2] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.025)
            cases_age[a,n,3] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.975)
        end
    end
    plt = groupedbar(cases_age[:,:,1]./1e3,
                    fontfamily="Helvetica",
                    lab =labels,
                    xticks = (1:17,age_cats))

    return plt

end

function plot_total_incidence_by_age_normed(scenario_group,treatment_rates,rel_transmission_perc,i)
    d1,d2,d3 = size(scenario_group[1][i][4][:,:,:])
    n = length(scenario_group)
    cases_age = zeros(d2,n,3)
    label_scen = "$(rel_transmission_perc[1])% rel. infectiousness"
    for l in 2:length(rel_transmission_perc)
        label_scen = hcat(label_scen,"$(rel_transmission_perc[l])% rel. infectiousness")
    end
    for (n,results) in enumerate(scenario_group)
        for a = 1:d2
            cases_age[a,n,1] = median([sum(results[i][4][:,a,k]) for k =1:1000])
            cases_age[a,n,2] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.025)
            cases_age[a,n,3] = quantile([sum(results[i][4][:,a,k]) for k =1:1000],0.975)
        end
    end
    for l = 1:n
            cases_age[:,l,1] .= LinearAlgebra.normalize(cases_age[:,l,1],1)
    end
    plt = groupedbar(cases_age[:,:,1],
                    fontfamily="Helvetica",
                    lab =label_scen,
                    xticks = (1:17,age_cats))

    return plt

end

function plot_total_incidence_by_treatment(results,treatment_rates)
    cases= [sum(results[i][4][21,:,1]) for i = 1:4]
    case_L = [sqrt(sum(results[i][4][21,:,2].^2)) for i = 1:4]
    case_U = [sqrt(sum(results[i][4][21,:,3].^2)) for i = 1:4]

    plt = scatter(collect(1:4)*2,cases,
            # yerror = (case_L,case_U),
            lab = "",
            xticks = (collect(1:4)*2,["n/t","3 wks","2 wks","1 wk"]),
            ms = 10)
    title!(plt,"Scenario A: Summed median total cases by treatment")
    ylabel!(plt,"Total cases")
    xlabel!(plt,"Av. time to treatment")
    ylims!(plt,(0.,2.5e7))
    return plt
end


text:
module Charts

using Genie, Stipple

import DataFrames
import Genie.Renderer.Html: HTMLString, normal_element, register_normal_element
import Stipple: Undefined, UNDEFINED

export PlotOptions, PlotData, PlotSeries, plot

const DEFAULT_WRAPPER = Genie.Renderer.Html.template

register_normal_element("apexchart", context = @__MODULE__)

"""
  PlotOptions
A data structure representing the plot options.
"""
Base.@kwdef mutable struct PlotOptions
  chart_animations_enabled::Bool = false
  chart_animations_easing::Union{String,Symbol} = :easeinout
  chart_animations_speed::Int = 800
  chart_background::String = "#fff"
  chart_font_family::String = "Helvetica, Arial, sans-serif"
  chart_fore_color::String = "#373d3f"
  chart_height::Union{Int,String} = "auto"
  chart_offset_x::Int = 0
  chart_offset_y::Int = 0
  chart_selection_enabled::Bool = true
  chart_selection_type::Union{String,Symbol} = :xy
  chart_sparkline_enabled::Bool = false
  chart_stacked::Bool = false
  chart_stack_type::String = "normal"
  chart_toolbar_show::Bool = false
  chart_type::Union{String,Symbol} = :line
  chart_width::Union{Int,String} = "100%"
  chart_zoom_enabled::Bool = true
  chart_zoom_type::Union{String,Symbol} = :xy

  colors::Union{Vector{String},Undefined} = ["#2E93fA", "#66DA26", "#0A557A", "#E91E63", "#FF9800"]

  data_labels_enabled::Bool = false

  fill_opacity::Union{Int,Float64} = 1.0

  grid_row_colors::Union{Vector{String},Undefined} = ["#EEEEEE", "transparent"]
  grid_row_opacity::Union{Int,Float64} = 1.0
  grid_show::Bool = true
  grid_xaxis_lines_show::Bool = false
  grid_yaxis_lines_show::Bool = false

  labels::Vector{String} = String[]
  legend_show::Bool = true
  legend_position::Union{String,Symbol} = :bottom
  legend_font_size::String = "14px"
  legend_font_family::String = "Helvetica, Arial"

  no_data_text::String = ""

  plot_options_area_fill_to::Union{String,Symbol} = :origin

  plot_options_bar_border_radius::Int = 0
  plot_options_bar_column_width::String = "100%"
  plot_options_bar_data_labels_position::Union{String,Symbol} = :center
  plot_options_bar_ending_shape::Union{String,Symbol} = :flat
  plot_options_bar_horizontal::Bool = false

  plot_options_bubble_min_bubble_radius::Union{Int,Undefined} = UNDEFINED
  plot_options_bubble_max_bubble_radius::Union{Int,Undefined} = UNDEFINED

  plot_options_pie_size::Union{Int,Undefined} = UNDEFINED
  plot_options_pie_start_angle::Int = 0
  plot_options_pie_end_angle::Int = 360
  plot_options_pie_expand_on_click::Bool = true
  plot_options_pie_offset_x::Int = 0
  plot_options_pie_offset_y::Int = 0
  plot_options_pie_custom_scale::Int = 1
  plot_options_pie_data_labels_offset::Int = 0
  plot_options_pie_data_labels_min_angle_to_show_label::Int = 10
  plot_options_pie_donut_size::Union{Int,String} = 65
  plot_options_pie_donut_background::String = "transparent"
  plot_options_pie_donut_labels_show::Bool = false
  plot_options_pie_donut_labels_name::Dict = Dict()
  plot_options_pie_donut_labels_value::Dict = Dict()
  plot_options_pie_donut_labels_total::Dict = Dict()

  stroke_curve::Union{String,Symbol,Vector{String},Vector{Symbol}} = :smooth
  stroke_show::Bool = true
  stroke_width::Union{Int,Vector{Int}} = 2
  stroke_colors::Vector{String} = colors

  subtitle_text::String = ""
  subtitle_align::Union{String,Symbol} = :left
  subtitle_style_font_size::String = "12px"

  theme_mode::Union{String,Symbol} = :light
  theme_palette::Union{String,Symbol} = :palette1

  title_text::String = ""
  title_align::Union{String,Symbol} = :left
  title_margin::Int = 10
  title_style_font_family::Union{String,Undefined} = UNDEFINED
  title_style_font_size::String = "14px"
  title_style_font_weight::Union{Int,Symbol,String} = :bold
  title_style_color::String = "#263238"

  tooltip_enable::Bool = true

  xaxis_categories::Union{Vector{String},Vector{Float64}} = String[]
  xaxis_decimals_in_float::Union{Int,Undefined} = UNDEFINED
  xaxis_labels_show::Bool = true
  xaxis_max::Union{Int,Float64,String,Undefined} = UNDEFINED
  xaxis_min::Union{Int,Float64,String,Undefined} = UNDEFINED
  xaxis_tick_amount::Union{Int,Float64,String,Undefined} = UNDEFINED
  xaxis_tick_placement::Union{String,Symbol} = :between # :on
  xaxis_type::Union{String,Symbol} = :category

  yaxis_decimals_in_float::Union{Int,Undefined} = UNDEFINED
  yaxis_labels_show::Bool = true
  yaxis_max::Union{Int,Float64,String,Undefined} = UNDEFINED
  yaxis_min::Union{Int,Float64,String,Undefined} = UNDEFINED
  yaxis_show::Bool = true
  yaxis_tick_amount::Union{Int,Float64,String,Undefined} = UNDEFINED

  extra_properties::Dict = Dict()
  extra_options::Dict = Dict()
end

"""
  PlotData
A data structure storing data for plotting.
"""
Base.@kwdef mutable struct PlotData{T<:Vector}
  data::T = T[]
  key::Union{Symbol,Nothing} = :data
end

function PlotData(data::DataFrames.DataFrame; key=:data)
  PlotData(key = key, data = [Array(r) for r in DataFrames.eachrow(data)])
end

"""
  PlotData(data)
Transform initial data in PlotData struct.

 # Examples
 ```jldoctest
 julia> vector = [1, 2, 3, 4];

 julia> PlotData(vector)
PlotData{Vector{Int64}}([1, 2, 3, 4], :data)
```
"""
function PlotData(data::Vector{T}) where {T}
  PlotData(data = data)
end

function PlotData(data...) where {T}
  PlotData([data...])
end

"""
  PlotSeries
A data structure to store PlotData struct with name.
"""
Base.@kwdef mutable struct PlotSeries
  name::String = ""
  plotdata::PlotData = PlotData(Any[])
end

function plot(fieldname::Union{Symbol,String};
              options::Union{Symbol,Nothing} = nothing,
              wrap::Function = DEFAULT_WRAPPER,
              args...) :: String

  k = (Symbol(":series"),)
  v = Any["$fieldname"]

  if options !== nothing
    k = (k..., Symbol("options!"))
    push!(v, options)
  end

  wrap() do
    apexchart(; args..., NamedTuple{k}(v)...)
  end
end

#===#

function Base.parse(::Type{Vector{PlotSeries}}, d::Vector{Any})
  [PlotSeries(name = x["name"], data = x["data"]) for x in d]
end

function Base.parse(::Type{PlotSeries}, x::Dict{String,Any})
  PlotSeries(name = x["name"], data = x["data"])
end

function Base.convert(::Type{Vector{PlotSeries}}, a::Vector{Dict{String,Any}})
  [parse(PlotSeries, x) for x in a]
end

function Base.parse(::Type{Dict{String,Any}}, ::PlotSeries)
  # TODO: not implemented
end

function Base.parse(::Type{PlotData}, x::Vector{Any})
  PlotData(x)
end

function Base.parse(::Type{PlotOptions}, ::Dict{String,Any})
  # error("Not implemented") # todo: add parser
end

#===#

function Stipple.watch(vue_app_name::String, fieldtype::R{PlotSeries}, fieldname::Symbol, channel::String, model::M)::String where {M<:ReactiveModel}
  string(vue_app_name, raw".\$watch('", fieldname, "', function(newVal, oldVal){

  });\n\n")
end

#===#

function Stipple.render(pd::PlotData{T}, fieldname::Union{Symbol,Nothing} = nothing) where {T<:Vector}
  pd.key === nothing ? pd.data : [Dict(pd.key => pd.data)]
end

function Stipple.render(pdv::Vector{PlotData{T}}, fieldname::Union{Symbol,Nothing} = nothing) where {T<:Vector}
  [Dict(pd.key => pd.data) for pd in pdv]
end

function Stipple.render(ps::PlotSeries, fieldname::Union{Symbol,Nothing} = nothing)
  Dict(:name => ps.name, ps.plotdata.key => ps.plotdata.data)
end

function Stipple.render(psv::Vector{PlotSeries}, fieldname::Union{Symbol,Nothing} = nothing)
  [Dict(:name => ps.name, ps.plotdata.key => ps.plotdata.data) for ps in psv]
end

function Stipple.render(po::PlotOptions, fieldname::Union{Symbol,Nothing} = nothing)
  val = Dict(
    # chart
    :chart => Dict(
      :animations => Dict(
        :enabled => po.chart_animations_enabled,
        :easing => po.chart_animations_easing,
        :speed => po.chart_animations_speed
      ),
      :background => po.chart_background,
      :fontFamily => po.chart_font_family,
      :foreColor => po.chart_fore_color,
      :height => po.chart_height,
      :offsetX => po.chart_offset_x,
      :offsetY => po.chart_offset_y,
      :selection => Dict(
        :enabled => po.chart_selection_enabled,
        :type => po.chart_selection_type
      ),
      :sparkline => Dict(
        :enabled => po.chart_sparkline_enabled
      ),
      :stacked => po.chart_stacked,
      :stackType => po.chart_stack_type,
      :toolbar => Dict(
        :show => po.chart_toolbar_show
      ),
      :type => po.chart_type,
      :width => po.chart_width,
      :zoom => Dict(
        :enabled => po.chart_zoom_enabled,
        :type => po.chart_zoom_type
      )
    ),

    :colors => po.colors,

    :dataLabels => Dict(
      :enabled => po.data_labels_enabled
    ),

    :fill => Dict(
      :opacity => po.fill_opacity
    ),

    # grid
    :grid => Dict(
      :show => po.grid_show,
      :row => Dict(
        :colors => po.grid_row_colors,
        :opacity => po.grid_row_opacity
      ),
      :xaxis => Dict(
        :lines => Dict(
          :show => po.grid_xaxis_lines_show
        )
      ),
      :yaxis => Dict(
        :lines => Dict(
          :show => po.grid_yaxis_lines_show
        )
      )
    ),

    :labels => po.labels,

    :legend => Dict(
      :position => po.legend_position,
      :fontFamily => po.legend_font_family,
      :fontSize => po.legend_font_size,
      :show => po.legend_show
    ),

    :noData => Dict(
      :text => po.no_data_text
    ),

    :stroke => Dict(
      :curve => po.stroke_curve,
      :show => po.stroke_show,
      :width => po.stroke_width,
      :colors => po.stroke_colors
    ),

    :subtitle => Dict(
      :align => po.subtitle_align,
      :text => po.subtitle_text,
      :style => Dict(
        :fontSize => po.subtitle_style_font_size
      )
    ),

    :theme => Dict(
      :mode => po.theme_mode,
      :palette => po.theme_palette
    ),

    :title => Dict(
      :text => po.title_text,
      :align => po.title_align,
      :margin => po.title_margin,
      :style => Dict(
        :color => po.title_style_color,
        :fontFamily => po.title_style_font_family,
        :fontSize => po.title_style_font_size,
        :fontWeight => po.title_style_font_weight
      )
    ),

    :tooltip => Dict(
      :enable => po.tooltip_enable
    ),

    :xaxis => Dict(
      :categories => po.xaxis_categories,
      :decimalsInFloat => po.xaxis_decimals_in_float,
      :tickAmount => po.xaxis_tick_amount,
      :tickPlacement => po.xaxis_tick_placement,
      :type => po.xaxis_type,
      :max => po.xaxis_max,
      :min => po.xaxis_min,
      :labels => Dict(
        :show => po.xaxis_labels_show
      )
    ),

    :yaxis => Dict(
      :decimalsInFloat => po.yaxis_decimals_in_float,
      :labels => Dict(
        :show => po.yaxis_labels_show
      ),
      :max => po.yaxis_max,
      :min => po.yaxis_min,
      :show => po.yaxis_show,
      :tickAmount => po.yaxis_tick_amount,
    )
  )

  isempty(po.extra_properties) || (val = recursive_merge(val, po.extra_properties))

  plot_options = if po.chart_type == :area
    Dict(
      :area => Dict(
        :fillTo => po.plot_options_area_fill_to
      )
    )
  elseif po.chart_type == :bar
    Dict(
      :bar => Dict(
        :borderRadius => po.plot_options_bar_border_radius,
        :columnWidth => po.plot_options_bar_column_width,
        :dataLabels => Dict(
          :position => po.plot_options_bar_data_labels_position
        ),
        :endingShape => po.plot_options_bar_ending_shape,
        :horizontal => po.plot_options_bar_horizontal
      )
    )
  elseif po.chart_type == :bubble
    Dict(
      :bubble => Dict(
        :minBubbleRadius => po.plot_options_bubble_min_bubble_radius,
        :maxBubbleRadius => po.plot_options_bubble_max_bubble_radius
      )
    )
  elseif po.chart_type == :pie
    Dict(
      :pie => Dict(
        :size  => po.plot_options_pie_size,
        :startAngle => po.plot_options_pie_start_angle,
        :endAngle   => po.plot_options_pie_end_angle,
        :expandOnClick => po.plot_options_pie_expand_on_click,
        :offsetX    => po.plot_options_pie_offset_x,
        :offsetY    => po.plot_options_pie_offset_y,
        :customScale  => po.plot_options_pie_custom_scale,
        :dataLabels => Dict(
          :offset   => po.plot_options_pie_data_labels_offset,
          :minAngleToShowLabel => po.plot_options_pie_data_labels_min_angle_to_show_label
        ),
        :donut => Dict(
          :size => endswith(string(po.plot_options_pie_donut_size) |> strip, '%') ?
                    string(po.plot_options_pie_donut_size) |> strip :
                    (string(po.plot_options_pie_donut_size) |> strip) * "%",
          :background => po.plot_options_pie_donut_background,
          :labels => Dict(
            :show => po.plot_options_pie_donut_labels_show,
            :name => po.plot_options_pie_donut_labels_name,
            :value => po.plot_options_pie_donut_labels_value,
            :total => po.plot_options_pie_donut_labels_total
          )
        ),
      )
    )
  else
    Dict()
  end

  plot_options = recursive_merge(plot_options, get(po.extra_properties, :plotOptions, Dict()))
  isempty(po.extra_options) || (plot_options = recursive_merge(plot_options, po.extra_options))
  isempty(plot_options) || (val[:plotOptions] = plot_options)

  val
end

# #===#

recursive_merge(x::AbstractDict...) = merge(recursive_merge, x...)
recursive_merge(x...) = x[end]

end


text:
function samplefile(infname::AbstractString, outfname::AbstractString, pct::Float64=0.0005, sorted::Bool=true)
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    exists(dfs, infname) || (return 0)
    fptr = open(dfs, infname, "r")
    wptr = open(dfs, outfname, "w")
    nsamp = 0
    try
        nsamp = samplefile(fptr, wptr, pct, sorted)
    finally
        close(fptr)
        close(wptr)
    end
    nsamp
end

function samplefile(fptr::IO, wptr::IO, pct::Float64, sorted::Bool=true)
    len = filesize(fptr)
    data = Array(UInt8, len)
    read!(fptr, data)
    datastr = bytestring(data)
    dataarr = split(datastr)

    nsamp = round(Int, length(dataarr) * pct + 1)
    dataarr = sample(dataarr, nsamp)
    sorted && sort!(dataarr)

    for s in dataarr
        write(wptr, s, '\n')
    end
    nsamp
end

function samplerandfile(fname::AbstractString)
    println("samplerandfile starting at $(myid())")
    locfname = "$(fname)-$(myid())"
    d,f = splitdir(locfname)
    sampfname = joinpath(TeraSort.SAMPDIR, f)
    nsamp = samplefile(locfname, sampfname)
    println("samplerandfile sampled $nsamp records at $(myid())")
    nsamp
end

function samplefilenames(fname::AbstractString, wrkrs)
    d,f = splitdir(fname)
    sampfname = joinpath(TeraSort.SAMPDIR, f)
    ["$sampfname-$id" for id in wrkrs]
end

function getpartitionpoints(infname::AbstractString, outfname::AbstractString, nparts::Int)
    println("determining partition points...")
    dfs = HDFSClient(TeraSort.DFSHOST, TeraSort.DFSPORT, UserGroupInformation(TeraSort.DFSUSER))
    fptr = open(dfs, infname, "r")
    wptr = open(dfs, outfname, "w")
    splpts = splitpoints(Int(filesize(fptr)), nparts)
    for splpt in splpts
        seek(fptr, max(0, splpt-TeraSort.RECSZ))
        s = readline(fptr)
        (length(s) == TeraSort.RECSZ) || (s = readline(fptr))
        write(wptr, s)
    end
    close(fptr)
    close(wptr)
    println("wrote $(length(splpts)) partition points for $nparts partitions")
    nothing
end

function distsamplefile()
    wrkrs = workers()
    nwrkrs = length(wrkrs)
    refs = Array(RemoteRef, nwrkrs)
    idx = 1
    fname = joinpath(INDIR, "randstr")
    for nid in wrkrs
        refs[idx] = @spawnat nid TeraSort.samplerandfile(fname)
        idx += 1
    end
    nrecs = 0
    for idx in 1:nwrkrs
        nrecs += fetch(refs[idx])
    end
    sampfnames = samplefilenames(fname, wrkrs)
    mergesorted("$fname-sample", sampfnames...)
    getpartitionpoints("$fname-sample", "$fname-partitions", 5)
    nrecs
end

function terasample()
    nrecs = t = b = g = 0
    nrecs, t, b, g = @timed distsamplefile()
    println("sampled $nrecs records")
    println("time:$t, memused:$b, gctime:$g")
end


text:
#Copyright Nathan Parrish 2018
#University of Colorado, Boulder
#Colorado Center for Astrodynamics Research
#Celestial and Spaceflight Mechanics Laboratory
#
#
#This module, "LowThrustOpt", is designed for low-thrust spacecraft trajectory
#optimization, especially in the Earth-Moon circular restricted 3-body problem
#dynamics. The algorithms implemented here are also valid in any other
#dynamical environment.
#
#Direct multiple shooting
#Indirect multiple shooting

#To add the LowThrustOpt directory to the path for a single session:
# push!(LOAD_PATH, "C:/.../path/to/containing/folder")

module LowThrustOpt

#Earth-Moon system:
const global MU = 0.012150585609624037 # MU = mu_moon/(mu_moon + mu_planet)
const global DU = 384747.96285603708 #Earth-Moon distance (km) (km per distance unit)
const global TU = 375699.81732246041 #time units (seconds per time unit)
const global r_moon = 1737. #km
const global r_earth = 6378. #km
const global day = 86400. #seconds
##Calculate mu_moon based on MU for consistency with other code:
const global mu_planet = 398600.4415
const global mu_moon = (MU * mu_planet) / (1 - MU)

#External packages required:
using DifferentialEquations
using Ipopt
using JuMP
using Interpolations
using ForwardDiff
using GeneralCode
using Plots

export multiShoot_CRTBP_direct
export plotTrajPlotly_direct
export meshRefine_direct
export multiShoot_CRTBP_indirect
export plotTrajPlotly_indirect
export controlLaw_cart
export CRTBP_stateCostate_deriv!
export CRTBP_prop_EP_NNControl_deriv!
export jacobiConstant
export interpInitialStates
export find_τ
export densify
export reduceFuel_indirect
export addTimeFinal

include("multiShoot_CRTBP_direct.jl")
include("multiShoot_CRTBP_indirect.jl")
include("CRTBP_stateCostate_deriv.jl")
include("CRTBP_prop_EP_deriv.jl")
include("HelperFunctions.jl")

end

max_stars_repo_name: travelingspaceman/LowThrustOpt

max_issues_repo_name: travelingspaceman/LowThrustOpt

max_forks_repo_name: travelingspaceman/LowThrustOpt


text:
### Utilities

# take a column

refcolumn(x::Matrix, j::Integer) = (n = size(x,1); pointer_to_array(pointer(x, (j-1) * n + 1), n))

# Functions to construct q-matrix

function qmatrix!(q::Array{Float64}, labels::Vector{Int})
    @check_argdims size(q, 1) == length(labels)
    fill!(q, 0.0)
    for i = 1 : size(q, 1)
        q[i, labels[i]] = 1.0
    end
    q
end

function qmatrix!(q::Array{Float64}; p=0.0) # n x k matrix, each row sum to 1
    if !(0. <= p <= 1.)
        throw(ArgumentError("The value of p must be in [0, 1]"))
    end

    n = size(q, 1)
    K = size(q, 2)

    if p == 0.0
        rand!(q)
        bmultiply!(q, 1.0 ./ sum(q, 2), 1)
    elseif p == 1.0
        fill!(q, 0.0)
        rg = 1:K
        for i = 1:n
            k = rand(rg)
            q[i, k] = 1.0
        end
    else
        rand!(q)
        bmultiply!(q, (1 - p) ./ sum(q, 2), 1)
        rg = 1:K
        for i = 1:n
            k = rand(rg)
            q[i,k] += p
        end
    end
    q
end

qmatrix(labels::Vector{Int}, K::Int) = qmatrix!(Array(Float64, length(labels), K), labels)

qmatrix(n::Int, K::Int; p=0.0) = qmatrix!(Array(Float64, n, K); p=p)


text:
using DashiBoard
using Test

include("preprocessing.jl")