text:
"""
```
function scalar_reduce(args...)
```
Each individual iteration returns n scalars. The output is reduced to n vectors,
where the i-th vector contains all of the i-th scalars from each iteration.

The return type of reduce functions must be the same type as the tuple of
arguments passed in. If args is a tuple of Vector{Float64}, then the return
argument will be a Vector{Float64}.

e.g.
a, b = @parallel (scalar_reduce) for i in 1:10000
           [[1], [2]]
       end
a = [1, 1, 1, ...]
b = [2, 2, 2, ...]

Input/Output type: Vector{Vector{Float64}}
"""
function scalar_reduce(args...)
    return_arg = args[1]
    for (i, arg) in enumerate(args[2:end])
        for (j, el) in enumerate(arg)
            append!(return_arg[j], el)
        end
    end
    return return_arg
end

"""
```
function vector_reduce(args...)
```
Each individual iteration returns n Vector types; we vector-reduce to n matrices, where
the i-th column of that matrix corresponds to the i-th vector from an individual iteration.

Input/Output type: Vector{Matrix{Float64}}
"""
function vector_reduce(args...)
    nargs1 = length(args)    # The number of times the loop is run
    nargs2 = length(args[1]) # The number of variables output by a single run

    return_arg = args[1]
    for i in 1:nargs2
        for j in 2:nargs1
            return_arg[i] = hcat(return_arg[i], args[j][i])
        end
    end
    return return_arg
end

"""
```
function scalar_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function scalar_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Vector{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = arg
    end
    return return_arg
end

"""
```
function vector_reshape(args...)
```
Function ensures type conformity of the return arguments.
"""
function vector_reshape(args...)
    n_args = length(args)
    return_arg = Vector{Matrix{Float64}}(undef, n_args)
    for i in 1:n_args
        arg = typeof(args[i]) <: Vector ? args[i] : [args[i]]
        return_arg[i] = reshape(arg, length(arg), 1)
    end
    return return_arg
end

"""
```
function generate_free_blocks(n_free_para, n_blocks)
```

Return a Vector{Vector{Int64}} where each internal Vector{Int64} contains a subset of the range 1:n_free_para of randomly permuted indices. This is used to index out random blocks of free parameters from the covariance matrix for the mutation step.
"""
function generate_free_blocks(n_free_para::Int64, n_blocks::Int64)
    rand_inds = shuffle(1:n_free_para)

    subset_length     = cld(n_free_para, n_blocks) # ceiling division
    last_block_length = n_free_para - subset_length*(n_blocks - 1)

    blocks_free = Vector{Vector{Int64}}(undef, n_blocks)
    for i in 1:n_blocks
        if i < n_blocks
            blocks_free[i] = rand_inds[((i-1)*subset_length + 1):(i*subset_length)]
        else
            # To account for the fact that the last block may be smaller than the others
            blocks_free[i] = rand_inds[end-last_block_length+1:end]
        end
    end
    return blocks_free
end

"""
```
function generate_all_blocks(blocks_free, free_para_inds)
```

Return a Vector{Vector{Int64}} where each internal Vector{Int64} contains indices corresponding to those in `blocks_free` but mapping to `1:n_para` (as opposed to `1:n_free_para`). These blocks are used to reconstruct the particle vector by inserting the mutated free parameters into the size `n_para,` particle vector, which also contains fixed parameters.
"""
function generate_all_blocks(blocks_free::Vector{Vector{Int64}}, free_para_inds::Vector{Int64})
    n_free_para = length(free_para_inds)
    ind_mappings = Dict{Int64, Int64}()

    for (k, v) in zip(1:n_free_para, free_para_inds)
        ind_mappings[k] = v
    end

    blocks_all = similar(blocks_free)
    for (i, block) in enumerate(blocks_free)
        blocks_all[i] = similar(block)
        for (j, b) in enumerate(block)
            blocks_all[i][j] = ind_mappings[b]
        end
    end
    return blocks_all
end

function init_stage_print(cloud::ParticleCloud;
                          verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    if use_fixed_schedule
        println("--------------------------")
            println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
    else
        println("--------------------------")
            println("Iteration = $(cloud.stage_index)")
    end
	println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
	println("--------------------------")
        println("c = $(cloud.c)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
	println("--------------------------")
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(cloud.particles[1])
            println("$(cloud.particles[1].keys[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
	    end
    end
end

function init_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
                          verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    if use_fixed_schedule
        println("--------------------------")
            println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
    else
        println("--------------------------")
            println("Iteration = $(cloud.stage_index)")
    end
	println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
	println("--------------------------")
        println("c = $(cloud.c)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
	println("--------------------------")
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(para_symbols)
            println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
	    end
    end
end

function end_stage_print(cloud::ParticleCloud;
                         verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    total_sampling_time_minutes = cloud.total_sampling_time/60
    if use_fixed_schedule
        expected_time_remaining_sec = (cloud.total_sampling_time/cloud.stage_index)*(cloud.n_Φ - cloud.stage_index)
        expected_time_remaining_minutes = expected_time_remaining_sec/60
    end

    println("--------------------------")
    if use_fixed_schedule
        println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
        println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
        println("estimated time remaining: $(round(expected_time_remaining_minutes, digits = 4)) minutes")
    else
        println("Iteration = $(cloud.stage_index)")
        println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
    end
    println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
    println("--------------------------")
        println("c = $(cloud.c)")
        println("accept = $(cloud.accept)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
    println("--------------------------")
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(cloud.particles[1])
            println("$(cloud.particles[1].keys[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
        end
    end
end
function end_stage_print(cloud::Cloud, para_symbols::Vector{Symbol};
                         verbose::Symbol=:low, use_fixed_schedule::Bool = true)
    total_sampling_time_minutes = cloud.total_sampling_time/60
    if use_fixed_schedule
        expected_time_remaining_sec = (cloud.total_sampling_time/cloud.stage_index)*(cloud.n_Φ - cloud.stage_index)
        expected_time_remaining_minutes = expected_time_remaining_sec/60
    end

    println("--------------------------")
    if use_fixed_schedule
        println("Iteration = $(cloud.stage_index) / $(cloud.n_Φ)")
        println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
        println("estimated time remaining: $(round(expected_time_remaining_minutes, digits = 4)) minutes")
    else
        println("Iteration = $(cloud.stage_index)")
        println("time elapsed: $(round(total_sampling_time_minutes, digits = 4)) minutes")
    end
    println("--------------------------")
        println("phi = $(cloud.tempering_schedule[cloud.stage_index])")
    println("--------------------------")
        println("c = $(cloud.c)")
        println("accept = $(cloud.accept)")
        println("ESS = $(cloud.ESS[cloud.stage_index])   ($(cloud.resamples) total resamples.)")
    println("--------------------------")
    if VERBOSITY[verbose] >= VERBOSITY[:high]
        μ = weighted_mean(cloud)
        σ = weighted_std(cloud)
        for n=1:length(para_symbols)
            println("$(para_symbols[n]) = $(round(μ[n], digits = 5)), $(round(σ[n], digits = 5))")
        end
    end
end


text:
"""
    output_gt(target, src...)
---
This is a temporary function to output the 012 genotypes and ID
names in two separate files for other programs.

The result files are `target.id` and `target.gt`.
"""
function output_gt(target, src...)
    title("Output G with $src")
    item("Check if files available and qualified")
    nlc = Int[]
    for i in src
        if !isfile("$i.bim")
            warning("$i not available")
            return
        end
        push!(nlc, countlines("$i.bim"))
    end
    for n in nlc[2:end]
        if n ≠ nlc[1]
            warning("Files are different")
            return
        end
    end
    pathto, out = splitdir(target)
    if length(pathto) == 0
        warning("Writing G matrix in the current directory")
    else
        isdir(pathto) || mkdir(pathto)
    end
    done()

    item("Merging $src")
    tmp = mktempdir(".")
    write("$tmp/merge.lst", join(src, '\n'), '\n')
    run(`$plink --sheep --merge-list $tmp/merge.lst --out $tmp/plink`)
    open("$target.id", "w") do io
        for line in eachline("$tmp/plink.fam")
            println(io, split(line)[2])
        end
    end
    run(`$plink --sheep --bfile $tmp/plink --recode A --out $tmp/plink`)
    open("$tmp/plink.raw", "r") do io
        _  = readline(io)       # skip the line of SNP names
        buffer = IOBuffer()
        nid = 0
        for line in eachline(io)
            nid += 1
            t = 1
            for _ in 1:6
                t = findnext(' ', line, t) + 1
            end
            write(buffer, replace(line[t:end], ' '=>""))
        end
        open("$target.gt", "w") do gt
            GT = reshape(take!(buffer), :, nid)
            for i in 1:size(GT)[2]
                println(gt, String(GT[:, i]))
            end
        end
    end
    rm(tmp, recursive=true, force=true)
    done()
end

ext: jl


text:
using Documenter, BART

makedocs(;
    modules=[BART],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/brandondbutcher/BART.jl/blob/{commit}{path}#L{line}",
    sitename="BART.jl",
    authors="Brandon Butcher, University of Iowa",
    assets=String[],
)

deploydocs(;
    repo="github.com/brandondbutcher/BART.jl",
)


text:
using Flux
using Flux: logitcrossentropy, onehotbatch, onecold, @epochs
using Flux.Data: DataLoader
using MLDatasets: MNIST
using Random
using Zygote

function get_dataloaders(batch_size::Int, shuffle::Bool)
    train_x, train_y = MNIST.traindata(Float32)
    test_x, test_y = MNIST.testdata(Float32)

    train_y, test_y = onehotbatch(train_y, 0:9), onehotbatch(test_y, 0:9)

    train_loader = DataLoader(train_x, train_y, batchsize=batch_size, shuffle=shuffle)
    test_loader = DataLoader(test_x, test_y, batchsize=batch_size, shuffle=shuffle)

    return train_loader, test_loader
end

function accuracy(data_loader, model)
    acc_correct = 0
    for (x_batch, y_batch) in data_loader
        batch_size = size(x_batch)[end]
        acc_correct += sum(onecold(model(x_batch)) .== onecold(y_batch)) / batch_size
    end
    return acc_correct / length(data_loader)
end

function create_model(input_dim, dropout_ratio, hidden_dim, classes)
    return Chain(
        Flux.flatten,
        Dense(input_dim, hidden_dim, relu),
        Dropout(dropout_ratio),
        Dense(hidden_dim, classes)
    )
end

function main(num_epochs, batch_size, shuffle, η)
    train_loader, test_loader = get_dataloaders(batch_size, shuffle)

    model = create_model(28*28, 0.2, 128, 10)
    trainable_params = Flux.params(model)

    optimiser = ADAM(η)
    loss(x,y) = logitcrossentropy(model(x), y)

    @epochs num_epochs Flux.train!(loss, trainable_params, train_loader, optimiser)

    testmode!(model)
    @show accuracy(train_loader, model)
    @show accuracy(test_loader, model)
    println("Complete!")
end

if abspath(PROGRAM_FILE) == @__FILE__
    batch_size = 64
    shuffle_data = true
    η = 0.0001
    num_epochs = 1
    main(num_epochs, batch_size, shuffle_data, η)
end

max_stars_repo_name: alecokas/ml-framework-comparison

max_issues_repo_name: alecokas/ml-framework-comparison

max_forks_repo_name: alecokas/ml-framework-comparison


text:
function P_F_equation(cur_reactor::AbstractReactor, cur_value::Number)
  cur_G_T = cur_value

  cur_G_T /= K_F(cur_reactor)

  cur_G_T /= K_n(cur_reactor) ^ 2

  cur_G_T /= cur_reactor.sigma_v

  cur_equation = Equation(-1, 0, 2, cur_G_T)

  cur_equation
end


text:
"""
is_crash(scene::Scene)
return true if the ego car is in collision in the given scene, do not check for collisions between
other participants
"""
function is_crash(scene::Scene)
ego = scene[findfirst(1, scene)]
@assert ego.id == 1
if ego.state.v ≈ 0
    return false
end
ego_safety_distance_side = Vehicle(ego.state, VehicleDef(2, ego.def.length, ego.def.width+1.5), ego.id)
for veh in scene
    if veh.id != 1
        if AutomotivePOMDPs.is_colliding(ego_safety_distance_side, veh)
            #println(veh)
            println("-----------------> Collision <----------------------")
            return true
        end
    end
end
return false
end

function get_observations_state_space(model::PedestrianAvoidancePOMDPFrenet, ego::Vehicle, sensor_observations::Vector{Vehicle} )

    o = Dict{Int64, SingleOCFObs}()
    for object in sensor_observations

        object_posF = Frenet(proj(object.state.posG, get_lane(model.env.roadway, ego.state),
                                    model.env.roadway, move_along_curves=false), model.env.roadway)

        delta_s = object_posF.s - ego.state.posF.s - ego.def.length/2 - object.def.width/2
        delta_t = object_posF.t - ego.state.posF.t
        delta_theta = object_posF.ϕ - ego.state.posF.ϕ

         if ( AutomotivePOMDPs.is_observable_fixed(ego.state, VehicleState(VecSE2(ego.state.posG.x+delta_s, ego.state.posG.y+delta_t, 0.), 0.0), model.pomdp.env) == false )
            occluded = true
        else
            occluded = false
        end

        if ( occluded == false && delta_s < model.pomdp.S_MAX && delta_s > model.pomdp.S_MIN && delta_t < model.pomdp.T_MAX  &&  delta_t > model.pomdp.T_MIN )
            o[object.id] = SingleOCFState(ego.state.posF.t, ego.state.v, delta_s, delta_t, delta_theta, object.state.v)
            #println(o[object.id])
        end

       # println("PED: t: ", object_posF.t, " / ego: t: ", ego.state.posF.t, " s: ", ego.state.posF.s)
       # println("delta_s: ", delta_s, " delta_t: ", delta_t)
    end
    o[PEDESTRIAN_OFF_KEY] = get_state_absent(model.pomdp, ego.state.posF.t, ego.state.v)
    return o
end

function get_state_absent(pomdp::SingleOCFPOMDP, ego_y::Float64, ego_v::Float64)

    (ego_y_state_space,ego_v_state_space) = getEgoDataInStateSpace(pomdp, ego_y, ego_v)
    return SingleOCFState(ego_y_state_space, ego_v_state_space, -10.0, -10.0, pomdp.PED_THETA_RANGE[1] ,0.0)

end

function get_state_absent(pomdp::SingleOCFPOMDP, s::SingleOCFState)

    (ego_y_state_space,ego_v_state_space) = getEgoDataInStateSpace(pomdp, s.ego_y, s.ego_v)
    return SingleOCFState(ego_y_state_space, ego_v_state_space, -10.0, -10.0, pomdp.PED_THETA_RANGE[1] ,0.0)

end

function is_state_absent(pomdp::SingleOCFPOMDP, s::SingleOCFState)

    if (s.ped_s == -10.0 && s.ped_T == -10.0 )
        return true;
    else
        return false;
    end

end

function is_observation_absent(pomdp::SingleOCFPOMDP, o::SingleOCFObs)

    if (o.ped_s > pomdp.S_MAX || o.ped_s < pomdp.S_MIN || o.ped_T > pomdp.T_MAX ||  o.ped_T < pomdp.T_MIN )
        return true;
    else
        return false;
    end

end

function getEgoDataInStateSpace(pomdp::SingleOCFPOMDP, y::Float64, v::Float64)
    y_grid = RectangleGrid(pomdp.EGO_Y_RANGE)
    (id, weight) = interpolants(y_grid, [y])

    ego_y = ind2x(y_grid, findmax(id)[1])[1]

    v_grid = RectangleGrid(pomdp.EGO_V_RANGE)
    (id, weight) = interpolants(v_grid, [v])
    ego_v = ind2x(v_grid, findmax(id)[1])[1]
    return ego_y, ego_v
end

function getObstructionCorner(pomdp::SingleOCFPOMDP, obstacle::ConvexPolygon)

    x = Vector{Float64}(undef, obstacle.npts)
    y = Vector{Float64}(undef, obstacle.npts)
    for i = 1:obstacle.npts
        x[i] = obstacle.pts[i].x
        y[i] = obstacle.pts[i].y
    end

    ego_pos = VecE2(pomdp.ego_vehicle.state.posG.x, pomdp.ego_vehicle.state.posG.y)
    delta_s = maximum(x) - ego_pos.x - pomdp.ego_vehicle.def.length/2
    if delta_s > 0
        right_side = true
        if ( ego_pos.y > mean(y) )
            delta_t = -(ego_pos.y -  maximum(y))
            right_side = true
        else
            delta_t = minimum(y) - ego_pos.y
            right_side = false
        end
        return (delta_s, delta_t, right_side)
    else
        return (1000., 1000., false)
    end
end

function calulateHiddenPositionsRightSide(pomdp::SingleOCFPOMDP, obst_s::Float64, obst_T::Float64)

    sT_pos = []
    if ( obst_s < pomdp.S_MAX)
        idx = findfirst(x -> x >= obst_s, pomdp.S_RANGE)
        s_grid = pomdp.S_RANGE[idx:end]

        idx = findlast(x -> x < obst_T, pomdp.T_RANGE)
        T_grid = pomdp.T_RANGE[2:idx]

        thetha = atan(obst_T, obst_s)
        for s in s_grid
            dT = tan(thetha)*(s-obst_s)
            for T in T_grid
                if T < obst_T+dT
                    push!(sT_pos, [s, T])
                end
            end
        end
    end

    return sT_pos
end

function calulateHiddenPositionsLeftSide(pomdp::SingleOCFPOMDP, obst_s::Float64, obst_T::Float64)

    sT_pos = []
    if ( obst_s < pomdp.S_MAX)
        idx = findfirst(x -> x >= obst_s, pomdp.S_RANGE)
        s_grid = pomdp.S_RANGE[idx:end]

        idx = findfirst(x -> x >= obst_T, pomdp.T_RANGE)
        T_grid = pomdp.T_RANGE[idx:end-1]

        thetha = atan(obst_T, obst_s)
        for s in s_grid
            dT = tan(thetha)*(s-obst_s)
            for T in T_grid
                if T > obst_T+dT
                    push!(sT_pos, [s, T])
                end
            end
        end
    end

    return sT_pos
end

function objects_time_delay(model::PedestrianAvoidancePOMDPFrenet)

    objects_tracked = Vehicle[]

    for object in model.sensor_observations

	    if haskey(model.objects_tracked, object.id) == false   # new object
            model.objects_tracked[object.id] = model.t_current + 0.2
            #println(model.objects_tracked)
        else
	        if ( model.t_current >= model.objects_tracked[object.id] )   # how long is the object already visible
                model.objects_tracked[object.id] = model.t_current
                push!(objects_tracked, object)
            end
        end
    end

    # delete objects if they are older than 0.2s, because if object occurs again, the tracker needs again some time
    for oid in keys(model.objects_tracked)
        if ( model.objects_tracked[oid] < model.t_current - 0.2 )
            delete!(model.objects_tracked, oid)
        end
    end

    model.sensor_observations = deepcopy(objects_tracked)
    return model.sensor_observations
end


text:
function kernel(x::REALMAT, x₀::REALMAT, kernel_struct::Kernel)
    @assert size(x, 2) == size(x₀, 2) "The number of columns of x and x₀ do not match."
    # Dispatch to kernel
    K = kernel_calculation(x, x₀, kernel_struct::Kernel)
    return K
end
kernel(x::REALMAT, kernel_struct::Kernel) = kernel(x, x, kernel_struct)

# Linear
function kernel_calculation(x, x₀, k::Linear)
    K = [dot(x[i1, :], x[i2, :]) for i1 in 1:size(x, 1), i2 in 1:size(x₀, 1)]
    return K * k.σ
end

# Radial Basis Function (Squared Exponential)
function kernel_calculation(x, x₀, k::RBF)
    K = zeros(size(x, 1), size(x₀, 1))
    for i1 in 1:size(x, 1)
        for i2 in 1:size(x₀, 1)
            K[i1, i2] = k.σ * exp(-sum((x[i1, :] .- x₀[i2, :]).^2)/(2*k.l^2))
        end
    end
    return K
end

# Rational Quadratic
function kernel_calculation(x, x₀, k::RationalQuadratic)
    K = zeros(size(x, 1), size(x₀, 1))
    for i1 in 1:size(x, 1)
        for i2 in 1:size(x₀, 1)
            K[i1, i2] = k.σ * (1.0 + sum((x[i1, :] .- x₀[i2, :])^2)/(2 * k.a * k.l^2))^-k.a
        end
    end
    return K
end

# Periodic
function kernel_calculation(x, x₀, k::Periodic)
    K = zeros(size(x, 1), size(x₀, 1))
    for i1 in 1:size(x, 1)
        for i2 in 1:size(x₀, 1)
            K[i1, i2] = k.σ * exp(-2 * sin(sum(abs.(x[i1, :] .- x₀[i2, :])))^2/k.l^2)
        end
    end
    return K
end

# Locally Periodic
function kernel_calculation(x, x₀, k::LocallyPeriodic)
    K = zeros(size(x, 1), size(x₀, 1))
    for i1 in 1:size(x, 1)
        for i2 in 1:size(x₀, 1)
            K[i1, i2] = k.σ * exp(-2 * sin(sum(abs.(x[i1, :] .- x₀[i2, :])))^2/k.l^2) * exp(-sum((x[i1, :] .- x[i2, :].^2)/(2*k.l^2)))
        end
    end
    return K
end


text:
module MTD

## distributed packages
using Combinatorics
using Dates
using Distributions
using StatsBase

## personal packages
using SparseProbVec

include("mtd_inference.jl")
include("mtdg_inference.jl")
include("mmtd_inference.jl")
include("gof.jl")
include("general.jl")
include("reduce.jl")

end # module

max_stars_repo_name: mheiner/MTD

max_issues_repo_name: mheiner/MTD

max_forks_repo_name: mheiner/MTD


text:
module TabularTMazeEvalSet

using Curiosity
using MinimalRLCore
const TTMU = Curiosity.TabularTMazeUtils
include("../experiments/tabular_tmaze.jl")
using Random
using StatsBase
using GVFHordes
using Statistics
using JLD2

StatsBase.sample(p::GVFHordes.GVFParamFuncs.FunctionalPolicy, s, actions) =
    sample(Weights([p.func(s, a) for a in actions]))

function gen_dataset()
    parsed = TabularTMazeExperiment.default_args()
    parsed["cumulant_schedule"] = "Constant"
    parsed["cumulant"] = 1.0

    horde = TabularTMazeExperiment.get_horde(parsed, 21, 4, nothing)

    cumulant_schedule = TTMU.get_cumulant_schedule(parsed)

    exploring_starts = parsed["exploring_starts"]
    env = TabularTMaze(exploring_starts, cumulant_schedule)
    start_states = Curiosity.valid_state_mask()

    function del!(arr, ind_to_del)
        arr = arr[1:end .!= ind_to_del]
        return arr
    end
    #NOTE: This lets the agent start in terminal states. This is not normally possible.
    # Work backwards with index so you dont have to reindex due to shifting
    # ind_to_delete = [21, 17, 5, 1]
    # for i in ind_to_delete
    #     start_states = del!(start_states, i)
    # end
    observations = []
    actions = []
    total_start_states = []
    for s in start_states
        MinimalRLCore.reset!(env,s)
        for a in get_actions(env)
            push!(observations, MinimalRLCore.get_state(env))
            push!(actions, a)
            push!(total_start_states, s)
        end
    end
    num_returns = 2
    γ_thresh=1e-6
    horde_rets = zeros(length(horde.gvfs), length(observations))
    for (gvf_i,gvf) in enumerate(horde.gvfs)
        rets = monte_carlo_returns(env, gvf, total_start_states, actions, num_returns, γ_thresh)
        rets_avg = mean(hcat(rets...),dims=1) # stack horizontally and then average over runs
        horde_rets[gvf_i,:] = rets_avg
    end

    # For tabular tmaze, the goal indices are implemented in such a way
    #the MC rollout doesnot trigger end episode until 1 step later.
    ind_to_zero = []
    for ind in 1:length(observations)
        if sum(observations[ind][2:end]) != 0.0
            push!(ind_to_zero, ind)
        end
    end
    for ind in ind_to_zero
        horde_rets[:,ind] .= 0.0
    end
    return horde_rets, observations, actions
end

function save_data(rets,obs,actions)
    TTMazeUniformEvalSet = Dict()
    TTMazeUniformEvalSet["ests"] = rets
    TTMazeUniformEvalSet["actions"] = actions
    TTMazeUniformEvalSet["states"] = obs
    @save "./src/data/TTMazeUniformEvalSet.jld2" TTMazeUniformEvalSet
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
using NCDatasets
using Distributed
using SharedArrays
using Formatting

const jobs    = Channel{Tuple}(32)
const results = Channel{Tuple}(32)

var = ARGS[1]
ens = parse(Int, ARGS[2])

varname = var
s_map_file = format("clim_b.e11.B1850C5CN.f09_g16.{:03d}.pop.h.{}.050001-059912.nc", ens, varname)
wgt_file = "wgt_gx1v6_to_gx3v7.nc"
d_file = format("clim_LENS_B1850C5CN_{:03d}_gx3v7_{}.nc", ens, varname)

ds_w = Dataset(wgt_file, "r")
ds_s = Dataset(s_map_file, "r")
ds_d = Dataset(d_file, "c")

missing_value = 1e20

NN_idx = convert(Array{Integer}, nomissing(ds_w["NN_idx"][:], 0))
NNN = size(NN_idx)[1]

d_Nt = ds_s.dim["time"]
d_Nx = ds_w.dim["d_Nx"]
d_Ny = ds_w.dim["d_Ny"]
d_Nz = ds_s.dim["z_t"]

d_N = d_Nx * d_Ny

defDim(ds_d, "Nx", d_Nx)
defDim(ds_d, "Ny", d_Ny)
defDim(ds_d, "Nz", d_Nz)
defDim(ds_d, "time", Inf)

# Write horizontal grid information
for (varname, vardata, dims) in (
    ("lat", ds_w["d_lat"][:], ("Nx", "Ny",)),
    ("lon", ds_w["d_lon"][:], ("Nx", "Ny",)),
    ("z_w_top", ds_s["z_w_top"][:], ("Nz",)),
    ("z_w_bot", ds_s["z_w_bot"][:], ("Nz",)),
)
    println("varname: ", varname)
    v = defVar(ds_d, varname, Float64, dims)
    v.attrib["_FillValue"] = missing_value
    v[:] = vardata
end

println("NNN: ", NNN)

# ===== [BEGIN] parallel computation code =====

jobs_total = d_Nt * d_Nz
jobs_left  = jobs_total

function make_jobs()

    for k = 1:d_Nz, t = 1:d_Nt
        put!(jobs, (k, t))
    end

end

function do_work()

    #println("This is my pid: ", myid())
    for (k, t) in jobs
        d_data = zeros(Float64, d_Nx, d_Ny) # do it layer by layer
        s_data = reshape(nomissing(ds_s[varname][:, :, k, t], NaN), :)

        for i = 1 : length(d_data)
            if NN_idx[1, i] == 0
                d_data[i] = missing_value
                continue
            end

            members = 0.0
            for j = 1:NNN
                if isfinite(NN_idx[j, i])
                    d_data[i] += s_data[NN_idx[j, i]]
                    members += 1.0
                end
            end

            if members == 0.0
                println("Warning: no valid points to interpolate this ocean grid. Index: (", mod(i, d_Nx), ", ", floor(i/d_Nx), ")")
            end

            d_data[i] /= members

        end
        put!(results, (k, t, d_data))
    end
end

# ===== [END] parallel computation code =====

println("Workers: ", workers())

@async make_jobs()

for p in workers()
    @async do_work()
end

d_var = defVar(ds_d, varname, Float64, ("Nx", "Ny", "Nz", "time"))
d_var.attrib["_FillValue"] = missing_value

@elapsed while jobs_left > 0

    completed_jobs = jobs_total - jobs_left

    print(format("\rProgress: {:.2f}% ({:d} / {:d})",  completed_jobs / jobs_total * 100.0, completed_jobs, d_Nt))

    k, t, d_data = take!(results)

    d_var[:, :, k, t] = d_data

    global jobs_left -= 1

    if jobs_left == 0
        println()
    end
end

final_data = nomissing(d_var[:], NaN)

println("Total data count: ", length(final_data))
println("Total valid data count: ", sum(isfinite.(final_data)))
println("Total data holes count: ", sum(isnan.(final_data)))

close(ds_w)
close(ds_s)
close(ds_d)

println("Program ends.")

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
# Piracy to avoid code caching and reduce risks of silly
# errors and hard-to-understand error messages.
function serve(path="page")
    isdir(path) || error("Couldn't find path '$path'.")
    bk = pwd()
    cd(path)
    try
        F.serve(clear=true)
    catch e
        @show e
    finally
        cd(bk)
    end
    return nothing
end

# Should only be called in the deploy (on github-action)
function optimize(; input="page", output="")
    isdir(input) || error("Couldn't find the folder '$input'.")
    occursin("/", output) && error("No depth allowed in output.")
    bk = pwd()
    cd(input)
    F.optimize(minify=false)
    # Purge CSS to decrease bootstrap size massively
    io = IOBuffer()
    run(pipeline(`$(NodeJS.npm_cmd()) root`, stdout=io))
    nodepath = String(take!(io))
    run(`$(strip(nodepath))/purgecss/bin/purgecss --css __site/css/bootstrap.min.css --content __site/index.html --output __site/css/bootstrap.min.css`)

    isempty(output) && (cd(bk); return nothing)

    # copy the content of `__site` to `__site/$output`, it's
    # the user's responsibility to check this is valid (i.e. that
    # there isn't already a folder $output...)
    outp = mkpath(joinpath("__site", output))
    for obj in readdir("__site")
        obj == output && continue
        src = joinpath("__site", obj)
        dst = joinpath(outp, obj)
        mv(src, dst; force=true)
    end
    cd(bk)
    return nothing
end

"""
    newpage()
"""
function newpage(; path="page", overwrite=false)
    if isdir(path)
        if !overwrite
            error("Path '$path' already exists, use `overwrite=true` " *
                  "if you wish to remove that folder and start again.")
        else
            rm(path, recursive=true)
        end
    end
    mkdir(path)
    store = joinpath(dirname(pathof(PkgPage)), "web")
    for obj in readdir(store)
        obj == "DeployPage.yml" && continue
        src = joinpath(store, obj)
        dst = joinpath(path, obj)
        cp(src, dst)
    end
    # Pkg makes files read-only but we don't want that
    for (root, _, files) ∈ walkdir(path)
        for file in files
            chmod(joinpath(root, file), 0o644)
        end
    end
    # Try placing the `DeployPage.yml`
    name = "DeployPage.yml"
    gapath = mkpath(joinpath(".github", "workflows"))
    src = joinpath(store, name)
    dst = joinpath(gapath, name)
    cp(src, dst, force=false)
    chmod(dst, 0o644)
    return nothing
end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Maur좩cio Carvalho Maciel, Paula B. Staudt, Rafael P. Soares
#* $Id$
#*--------------------------------------------------------------------
type splitter_n
	splitter_n()=begin
		new(
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Outlet Streams",
				:Lower=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet stream",
				:PosX=>0,
				:PosY=>0.5001,
				:Symbol=>"_{in}"
			)),
			fill(stream(Dict{Symbol,Any}(
				:Brief=>"Outlet streams",
				:PosX=>1,
				:PosY=>0.5,
				:Symbol=>"_{out}"
			)),(NOutlet)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Distribution of Outlets",
				:Default=>0.5,
				:Symbol=>"\\phi"
			)),(NOutlet)),
			[
				:(sum(frac) = 1),
				:(Outlet([1:NOutlet]).F = Inlet.F*frac([1:NOutlet])),
				:(Outlet([1:NOutlet]).z = Inlet.z),
				:(Outlet([1:NOutlet]).P = Inlet.P),
				:(Outlet([1:NOutlet]).h = Inlet.h),
				:(Outlet([1:NOutlet]).T = Inlet.T),
				:(Outlet([1:NOutlet]).v = Inlet.v),
			],
			[
				"","Flow","Composition","Pressure","Enthalpy","Temperature","Vapourisation Fraction",
			],
			[:NOutlet,],
			[:Inlet,:Outlet,:frac,]
		)
	end
	NOutlet::DanaInteger
	Inlet::stream
	Outlet::Array{stream}
	frac::Array{fraction}
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export splitter_n
function setEquationFlow(in::splitter_n)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end
function atributes(in::splitter_n,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=false
	fields[:Icon]="icon/splitter_n"
	fields[:Brief]="Model of a splitter (NOT Handled by the GUI)"
	fields[:Info]="== Assumptions ==
* thermodynamics equilibrium
* adiabatic

== Specify ==
* the inlet stream
* (Noutlet - 1) fraction of split of the outlet streams:

	frac(i) = (Mole Flow of the outlet stream i /
				Mole Flow of the inlet stream)
				where i = 1, 2,...,Noutlet
"
	drive!(fields,_)
	return fields
end
splitter_n(_::Dict{Symbol,Any})=begin
	newModel=splitter_n()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# All pole digital filter for speech synthesis from LPC

poledf_delay(order::Int) = zeros(order)

type AllPoleDF <: LinearPredictionVariantsSynthesisFilter
    delay::Vector{Float64}

    function AllPoleDF(order::Int)
        new(poledf_delay(order))
    end
end

function filt!(f::AllPoleDF, x::Float64, a::Vector{Float64})
    d = delay(f)
    order = length(a) - 1

    y = x
    @inbounds for m=order:-1:2
        y -= a[m+1] * d[m]
        d[m] = d[m-1]
    end

    y -= a[2] * d[1]
    d[1] = y

    return y
end

function to_filtcoef(f::AllPoleDF, state::SpectralParamState{LinearPredictionCoef})
    has_loggain(state) && return copy(rawdata(state))

    b = copy(rawdata(state))
    for i=1:size(b, 2)
        b[1,i] = log(b[1,i])
    end
    b
end


text:
################################
# Coded by Julien Pascal
# Last modification : 30/04/2016
# Replication of the the following paper:

####################################################################################
# "On the Dynamics of Unemployment and Wage Distributions"
# Paper from Jean Marc Robin (2011)
# Available here: https://sites.google.com/site/jmarcrobin/research/publications
####################################################################################

##########################################################################################
# This version of the objective function returns the matrix variance-covariance of the
# moments of the simulated data. To be used to build an efficent MSM estimators.
# The length of each sub-sample is set to be equal to the length of the empirical data.
# In our case, this represents 255 quarters
#
# As an input, it has to the parameters estimated in a first minimization step
###########################################################################################

function objective_function_3(parameters::Vector)

	length_simulated_period = 255 #equal to the length of the empirical data
	number_repetitions = 200 #number of simulated economies
	# discard the first 10th moments to not take into account initial conditions:
	number_moments_to_keep = round(Int,0.9*number_repetitions)
	simulated_moments = zeros(number_moments_to_keep,7); #each column is a different moment, each row is an observation: here 50 observations and 7 moments

    ##################
    # Function that returns the inverse of lognormal distribution
    function logninv(p,mu,sigma)
        #output: scalar x
        logx0 = -sqrt(2)*erfcinv(2*p);
        x = exp(sigma*logx0 + mu);
        return x
    end
    #################

    ##########################
    #define match productivity:
    function yi_m(y_index,x_index, ygrid, xgrid)
        match = ygrid[y_index,1]*xgrid[x_index,1]
        return match
    end
    ##########################

    ##########################
    # define opportunity cost of employment:
    function zi_m(y_index,x_index, ygrid, xgrid, z0, alpha)
        op_cost = z0 + alpha*(yi_m(y_index,x_index, ygrid, xgrid) - z0) #call the function defined above
        return op_cost
    end

    ##########################
    # Exit from unemployment rate: ft
    function ft_f(St_m,ut_m,ut,lm,lambda0,M)
    #St_m a vector of dimension M times 1

        ft = 0;
        for i=1:M
            #calculate the sum
            if St_m[1,i]> 0
                ft = ft + ut_m[1,i]*lm[i,1];
            end
        end
        #multiply by lamba0 and divide by the overall level of unemployment ut:
        ft = lambda0*ft/ut;
        return ft
    end

    #function that calculates the job destruction rate
    function st_f(St_m,ut_m,ut,lm, delta, M)

        st = 0;
        for i=1:M
        #calculate the integral
            if St_m[1,i] <= 0
                st = st + (1 - ut_m[1,i])*lm[i,1];
            end
        end

        #finish the calculation:
        st = delta + (1 - delta)*st/(1 - ut);
        return st
    end

    #Calculates the next period's unemployment rate
    function ut_plus1_f(St_m, ut_m, delta, lambda0, M)

        #initialization:
        ut_plus1_m = zeros(1,M);

        for i=1:M
            #if the surplus is negative, unemployment is 1:
            if St_m[1,i] <= 0
                ut_plus1_m[1,i]= 1;
            else
                ut_plus1_m[1,i]= ut_m[1,i] + delta*(1 - ut_m[1,i]) - lambda0*ut_m[1,i];
            end
        end
        return ut_plus1_m
    end

    ############
    # Parameters
    ############

    #parameters to estimate:
    """
    z0 = parameters[1];
    sigma = parameters[2];
    pho = parameters[3];
    lambda0 = parameters[4];
    eta =  parameters[5];
    mu =  parameters[6];

    x_lower_bound = parameters[7]
    delta = parameters[8]

    k = 0.12;
    lambda1  = k*lambda0;
    #x_lower_bound = 0.73
    alpha = 0.64
    tau = 0.5
    #delta = 0.046  #exogenous layoff rate
    r = 0.05/4 # interest rate
    discount = (1 - delta)/(1 + r)
    epsilon = 0.002
    """

    z0 = parameters[1];
    sigma = parameters[2];
    lambda0 = parameters[3];
    eta =  parameters[4];
    mu =  parameters[5];
    x_lower_bound = parameters[6] #lower bound of the ability level

    pho = 0.94
    k = 0.12;
    lambda1  = k*lambda0;
    #x_lower_bound = 0.73
    alpha = 0.64
    tau = 0.5
    delta = 0.042 # exogenous layoff rate

    r = 0.05/4 # interest rate
    discount = (1 - delta)/(1 + r)
    epsilon = 0.002

    # grid parameters
    N = 100 #number of states
    M = 500 #number of ability levels

    # Define the grid for ability x
    xgrid_space = linspace(x_lower_bound,x_lower_bound+1,M); # column vector with M rows

    xgrid = Array{Float64}(M,1);
    for i=1:M
    xgrid[i,1] = xgrid_space[i];
    end

    # Define the grid for match ability yi:
        # 1. a:
        agrid = linspace((0+epsilon),(1-epsilon),N);

        # 2. yi:
        #intialization:
        ygrid = Array{Float64}(N,1);

        for i=1:N
            ygrid[i] = logninv(agrid[i],0,sigma); #calculate the inverse of the lognormal
        end

    #Markov transition matrix:
    Markov = zeros(N,N);

    # number of dimension
    # AR1 parameter:

    cop = Copulas.NormalCopula(2,pho)

    #Gaussian copula pdf:
    for j=1:N #move along the column:
        for i=1:N #move along the lines:
            Markov[i,j] = Copulas.dnormCopula([agrid[i] agrid[j]], cop)[];
        end
    end

    #normalize so that each row sum to 1:
    for i = 1:N
        Markov[i,:] = Markov[i,:]./sum(Markov[i,:]);
    end

    # distribution of workers on the grid x
    lm = Array{Float64}(M,1);

    # load the beta pdf
    b = Beta(eta, mu)

    for i=1:M
        lm[i,1] = pdf(b, xgrid[i]-x_lower_bound); #evaluation of the beta pdf at a given point
    end

    #normalization so that mass of workers sum up to 1:
    lm[:,1] = lm[:,1]./sum(lm[:,1]);

    # Pre calculate the opportunity cost along the grid of ability:
    # N times M matrix
    # rows = aggregate staes
    # columns = ability level
    zi_m_grid = zeros(N,M);

    for i=1:N
        for m=1:M
            zi_m_grid[i,m] = zi_m(i,m, ygrid, xgrid, z0, alpha);
        end
    end

    ##########################################################
    # Calculate the match surplus by value function iteration:
    # N times M matrices
    # N = aggregate level
    # M = ability level

    G = zeros(N,M); #intialization
    p = zeros(N,M);

    # Pre calculate values:
    for i = 1:N
        for m = 1:M
            G[i,m] = yi_m(i,m, ygrid, xgrid) - zi_m_grid[i,m];
            p[i,m] = yi_m(i,m, ygrid, xgrid);
        end
    end

    # N times M matrix
    # column index indicates the workers' type
    # row index indicates the shock index

        Si_m = zeros(N,M);

        #tol = 0.0001;
        tol = 0.01;
        maxits = 300;
        dif = tol+tol;
        its = 1;

        #initialization:
        up = ones(N,M);
        up_plus1 = zeros(N,M);
        compare = zeros(N,M);

        while dif>tol

             up_plus1 =  G + discount*Markov*max(up,compare)

              dif = norm(up_plus1 - up);

            up = up_plus1;

            its = its + 1;

            if its > maxits
                break
            end

        end

        Si_m = up;

###########################
# Simulation of the economy
###########################
number_periods = length_simulated_period *number_repetitions; # one period = one quarter

discard = round(Int,floor(number_periods/10)); #get rid of the first 10th observations

y_index_r = Array{Int16}(number_periods,1); #store the indexes of the shock
y_r = Array{Float64}(number_periods,1) #store value of the shock

#initial shock:
y_index_r[1,1] = 50; #initial shock index
y_r[1,1] = ygrid[y_index_r[1,1], 1]; #initial shock value

#Initialization:
ut_m_r = ones(number_periods,M); #row for the period, column for ability
St_m_r = Array{Float64}(number_periods,M); #row for the period, column for ability
ut_r = Array{Float64}(number_periods,1); #column vector of different observations
ft_r = Array{Float64}(number_periods,1); #column vector of different observations
qt_r = Array{Float64}(number_periods,1); #column vector of different observations
st_r = Array{Float64}(number_periods,1); #column vector of different observations

########################
# Loop over the economy:
for t=1:(number_periods-1)

    # Calculate the aggregate unemployment
    ut_r[t,1] = dot(ut_m_r[t,:],lm[:,1]); #dot product.

    # Calculate the surplus given the actual value of the shock
    St_m_r[t,:] = Si_m[y_index_r[t,1],:];

    #Remark: everytime I select a row, Julia treats it as a column vector
    #That is why I have to transpose my inputs in the next function
    # Exit rate from unemployment:
    ft_r[t,1] = ft_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, lambda0, M);

    #Job destruction rate:
    st_r[t,1] = st_f(St_m_r[t,:]', ut_m_r[t,:]', ut_r[t,1], lm, delta, M);

    # Law of motion of unemployment:
    ut_m_r[t+1,:] = ut_plus1_f(St_m_r[t,:]', ut_m_r[t,:]', delta, lambda0, M);

    # New shock from the markov transition matrix:
    r = rand(); #draw a number from the uniform distribution on 0, 1

    # I use the Markov transition matrix previously defined:
    prob = Markov[y_index_r[t,1],:];

    #I stock the index in the Markov matrix, as well as the value of the
    #shock y
    y_index_r[t+1,1] = sum(r .>= cumsum([0; prob]));

    y_r[t+1,1] = ygrid[y_index_r[t+1,1], 1];
end

#Check for odd values. These entities cannot be negative or zero:
for i = discard:(number_periods-1)
    if (isless(ut_r[i,1],0))
        error("negative unemployment rate")
    elseif (isless(ft_r[i,1],0))
        error("negative job finding rate")
    elseif (isless(st_r[i,1],0))
        error("negative job separation rate")
    end
end

moment_index = 1 #to index the observation of simulated moments
for t=1:number_periods
    #Stack the moments into the matrix "simulated moments"
    if t>discard #not take into account the fist 10% of the economy
    	if (t % length_simulated_period) == 0 #if the period is a multiple of length_simulated_period

    		mean_prod = mean(y_r[(t-254):t,1]) # Production yt
			std_prod = std(log(y_r[(t-254):t,1]))

			mean_unemployment = mean(ut_r[(t-254):t,1]) #Unemployment ut
			std_unemployment = std(log(ut_r[(t-254):t,1]))

			mean_exit_rate = mean(ft_r[(t-254):t,1]) #exit rate from unemployment ft
			std_exit_rate = std(log(ft_r[discard:(number_periods-1),1]))

            mean_job_destruction_rate = mean(st_r[(t-254):t,1]) # job separation rate st
            std_job_destruction_rate = std(log(st_r[discard:(number_periods-1),1]))

            #kurtosis_unemployment = kurtosis(log(ut_r[(t-254):t,1])) + 3 #returns excess kurtosis, so I add 3

            #b1 = [mean_prod; std_prod; mean_unemployment; std_unemployment; mean_exit_rate; std_exit_rate, mean_job_destruction_rate; std_job_destruction_rate] #simulated moments
            #b0 = [1.0; 0.023;  0.059;  0.210; 0.761; 0.085; 0.046; 0.156]

			simulated_moments[moment_index,1] = mean_prod #productivity yt
			simulated_moments[moment_index,2] = std_prod

			simulated_moments[moment_index,3] = mean_unemployment #unemployment rate ut
			simulated_moments[moment_index,4] = std_unemployment

            simulated_moments[moment_index,5] = mean_exit_rate #job finding rate ft
			simulated_moments[moment_index,6] = std_exit_rate

            simulated_moments[moment_index,7] = mean_job_destruction_rate #job separation rate st
            #simulated_moments[moment_index,8] = std_job_destruction_rate

            #simulated_moments[moment_index,5] = kurtosis_unemployment
            #simulated_moments[moment_index,5] = kurtosis_unemployment

			moment_index = moment_index + 1
    	end
    end
end

#Calculate the variance-covariance matrix of the simulated moments:
St = cov(simulated_moments) #by default, considers that column = variable and row = one observation

#Calculate the inverse of this matrix:
inverse_St = inv(St)

return Dict("St"=> St, "inverse_St" => inverse_St, "simulated_moments"=> simulated_moments)
end


text:
const PAR_MAGIC = "PAR1"
const SZ_PAR_MAGIC = length(PAR_MAGIC)
const SZ_FOOTER = 4
const SZ_VALID_PAR = 2*SZ_PAR_MAGIC + SZ_FOOTER

# page is the unit of compression
mutable struct Page
    colchunk::ColumnChunk
    hdr::PageHeader
    pos::Int
    uncompressed_data::Vector{UInt8}
    nextpos::Int64
end

"""
Keeps a cache of pages read from a file.
Pages are kept as weak refs, so that they can be collected when there's memory pressure.
"""
struct PageLRU
    refs::Dict{Tuple{ColumnChunk,Int64},WeakRef}
    lck::ReentrantLock
    function PageLRU()
        new(Dict{Tuple{ColumnChunk,Int64},WeakRef}(), ReentrantLock())
    end
end

function cacheget(fetcher, lru::PageLRU, chunk::ColumnChunk, startpos::Int64)
    key = (chunk,startpos)
    lock(lru.lck) do
        page = haskey(lru.refs, key) ? lru.refs[key].value : nothing
        if page === nothing
            page = fetcher()::Page
            lru.refs[key] = WeakRef(page)
        end
        return page
    end
end

"""
    Parquet.File(path; map_logical_types) => Parquet.File

Represents a Parquet file at `path` open for reading. Options to map logical types can be provided via `map_logical_types`.

`map_logical_types` can be one of:

- `false`: no mapping is done (default)
- `true`: default mappings are attempted on all columns (bytearray => String, int96 => DateTime)
- A user supplied dict mapping column names to a tuple of type and a converter function

Returns a `Parquet.File` type that keeps a handle to the open file and the file metadata and also holds a LRU cache of raw bytes of the pages read.
"""
mutable struct File
    path::String
    handle::IOStream
    meta::FileMetaData
    schema::Schema
    page_cache::PageLRU
end

function File(path::AbstractString; map_logical_types::Dict=TLogicalTypeMap())
    f = open(path)
    try
        return File(path, f; map_logical_types=map_logical_types)
    catch ex
        close(f)
        rethrow(ex)
    end
end

function File(path::AbstractString, handle::IOStream; map_logical_types::Dict=TLogicalTypeMap())
    is_par_file(handle) || error("Not a parquet format file: $path")
    meta_len = metadata_length(handle)
    meta = metadata(handle, path, meta_len)
    typemap = merge!(TLogicalTypeMap(), map_logical_types)
    File(String(path), handle, meta, Schema(meta.schema, typemap), PageLRU())
end

function close(par::Parquet.File)
    empty!(par.page_cache.refs)
    close(par.handle)
end

schema(par::Parquet.File) = par.schema

colname(par::Parquet.File, col::ColumnChunk) = colname(metadata(par,col))
colname(col::ColumnMetaData) = col.path_in_schema
function colnames(par::Parquet.File)
    names = Vector{Vector{String}}()
    cs = Int[]
    ns = String[]
    for x in par.schema.schema[2:end]
        if Parquet.num_children(x) > 0
            push!(cs, x.num_children)
            push!(ns, x.name)
        else
            if !isempty(cs)
                push!(names, [ns; x.name])
                cs[end] -= 1
                if cs[end] == 0
                    pop!(cs)
                    pop!(ns)
                end
            else
                push!(names, [x.name])
            end
        end
    end
    names
end

ncols(par::Parquet.File) = length(colnames(par))
nrows(par::Parquet.File) = par.meta.num_rows

coltype(par::Parquet.File, col::ColumnChunk) = coltype(metadata(par,col))
coltype(col::ColumnMetaData) = col._type

# return all rowgroups in the par file
rowgroups(par::Parquet.File) = par.meta.row_groups

function rowgroup_row_positions(par::Parquet.File)
    rgs = rowgroups(par)
    positions = Array{Int64}(undef, length(rgs)+1)
    idx = 1
    positions[idx] = Int64(1)
    for rg in rgs
        positions[idx+=1] = rg.num_rows
    end
    cumsum!(positions, positions)
end

columns(par::Parquet.File, rowgroupidx) = columns(par, rowgroups(par)[rowgroupidx])
columns(par::Parquet.File, rowgroup::RowGroup) = rowgroup.columns
columns(par::Parquet.File, rowgroup::RowGroup, colname::Vector{String}) = columns(par, rowgroup, [colname])
function columns(par::Parquet.File, rowgroup::RowGroup, cnames::Vector{Vector{String}})
    R = ColumnChunk[]
    for col in columns(par, rowgroup)
        (colname(par,col) in cnames) && push!(R, col)
    end
    R
end

##
# Iterator for pages in a column chunk
mutable struct ColumnChunkPages
    par::Parquet.File
    col::ColumnChunk
    startpos::Int64
    endpos::Int64

    function ColumnChunkPages(par::Parquet.File, col::ColumnChunk)
        startpos = page_offset(par, col)
        endpos = end_offset(par, col)
        new(par, col, startpos, endpos)
    end
end
eltype(::Type{ColumnChunkPages}) = Page
Base.iterate(ccp::ColumnChunkPages) = iterate(ccp, ccp.startpos)
function Base.iterate(ccp::ColumnChunkPages, startpos::Int64)
    if startpos >= ccp.endpos
        return nothing
    end

    page = cacheget(ccp.par.page_cache, ccp.col, startpos) do
        par = ccp.par
        io = par.handle
        seek(io, startpos)
        pagehdr = read_thrift(io, PageHeader)

        page_data_pos = position(io)
        pagesz = page_size(pagehdr)
        data = _use_mmap ? Mmap.mmap(io, Vector{UInt8}, (pagesz,), page_data_pos; grow=false, shared=false) : read!(io, Array{UInt8}(undef, pagesz))
        codec = metadata(ccp.par, ccp.col).codec

        if (codec != CompressionCodec.UNCOMPRESSED)
            uncompressed_sz = pagehdr.uncompressed_page_size
            #uncompressed_data = _use_mmap ? Mmap.mmap(Mmap.Anonymous(), Vector{UInt8}, (uncompressed_sz,), 0) : Array{UInt8}(undef, uncompressed_sz)
            uncompressed_data = Array{UInt8}(undef, uncompressed_sz)
            if codec == CompressionCodec.SNAPPY
                Snappy.snappy_uncompress(data, uncompressed_data)
            elseif codec == CompressionCodec.GZIP
                readbytes!(GzipDecompressorStream(IOBuffer(data)), uncompressed_data)
            elseif codec == CompressionCodec.ZSTD
                readbytes!(ZstdDecompressorStream(IOBuffer(data)), uncompressed_data)
            else
                error("Unknown compression codec for column chunk: $codec")
            end
            (length(uncompressed_data) == uncompressed_sz) || error("failed to uncompress page. expected $(uncompressed_sz), got $(length(uncompressed_data)) bytes")
        else
            uncompressed_data = data
        end

        nextpos = page_data_pos + pagesz
        page = Page(ccp.col, pagehdr, page_data_pos, uncompressed_data, nextpos)
    end

    page, page.nextpos
end

##
# Iterator for page values in a column chunk
mutable struct ColumnChunkPageValues{T}
    ccp::ColumnChunkPages
    max_repn::Int64
    max_defn::Int64
    has_repn_levels::Bool
    has_defn_levels::Bool
    repn_out::OutputState{Int32}
    defn_out::OutputState{Int32}
    valdict_out::OutputState{T}
    vals_out::OutputState{T}
    converter_fn::Function
end

function ColumnChunkPageValues(par::Parquet.File, col::ColumnChunk, ::Type{T}, converter_fn::Function=identity) where {T}
    cname = colname(par, col)

    max_repn = max_repetition_level(par.schema, cname)
    max_defn = max_definition_level(par.schema, cname)
    has_repn_levels = ((length(cname) > 1) && (max_repn > 0))
    has_defn_levels = !isrequired(par.schema, cname)

    repn_out = OutputState(Int32, 0)
    defn_out = OutputState(Int32, 0)
    valdict_out = OutputState(T, 0)
    vals_out = OutputState(T, 0)

    ccp = ColumnChunkPages(par, col)
    ColumnChunkPageValues{T}(ccp, max_repn, max_defn, has_repn_levels, has_defn_levels, repn_out, defn_out, valdict_out, vals_out, converter_fn)
end

function eltype(::Type{ColumnChunkPageValues{T}}) where {T}
   NamedTuple{(:value,:repn_level,:defn_level),Tuple{OutputState{T},OutputState{Int32},OutputState{Int32}}}
end

function Base.iterate(ccpv::ColumnChunkPageValues{T}) where {T}
    iterate(ccpv, ccpv.ccp.startpos)
end

function map_dict_vals(valdict::OutputState{T1}, vals::OutputState{T1}, map_vals::Vector{T2}) where {T1, T2}
    if !isempty(valdict.data) && (valdict.offset > 0)
        num_values = length(map_vals)
        @inbounds for idx in 1:num_values
            vals.data[idx] = valdict.data[map_vals[idx]+1]
        end
        vals.offset += num_values
    end
end

function Base.iterate(ccpv::ColumnChunkPageValues{T}, startpos::Int64) where {T}
    if startpos >= ccpv.ccp.endpos
        return nothing
    end

    read_data_page = false
    nextpos = startpos

    while !read_data_page
        page, nextpos = iterate(ccpv.ccp, nextpos)

        pagetype = page.hdr._type
        num_values = page_num_values(page)
        inp = InputState(page.uncompressed_data, 0)

        if (pagetype === PageType.DATA_PAGE) || (pagetype === PageType.DATA_PAGE_V2)
            read_data_page = true
            ccpv.has_repn_levels && reset_to_size(ccpv.repn_out, num_values)
            ccpv.has_defn_levels && reset_to_size(ccpv.defn_out, num_values)

            @debug("reading a data page for columnchunk")
            enc, defn_enc, repn_enc = page_encodings(page)
            nmissing = read_levels_and_nmissing(inp, ccpv.defn_out, ccpv.repn_out, defn_enc, repn_enc, Int(ccpv.max_defn), Int(ccpv.max_repn), num_values)
            nnonmissing = num_values - nmissing
            reset_to_size(ccpv.vals_out, nnonmissing)

            if enc === Encoding.PLAIN_DICTIONARY || enc === Encoding.RLE_DICTIONARY
                map_vals = read_data_dict(inp, nnonmissing)
                map_dict_vals(ccpv.valdict_out, ccpv.vals_out, map_vals)
            else
                if ccpv.converter_fn === identity
                    read_plain_values(inp, ccpv.vals_out, nnonmissing)
                else
                    read_plain_values(inp, ccpv.vals_out, nnonmissing, ccpv.converter_fn, ccpv.ccp.col.meta_data._type)
                end
            end
        elseif pagetype === PageType.DICTIONARY_PAGE
            ensure_additional_size(ccpv.valdict_out, num_values)
            if ccpv.converter_fn === identity
                read_plain_values(inp, ccpv.valdict_out, num_values)
            else
                read_plain_values(inp, ccpv.valdict_out, num_values, ccpv.converter_fn, ccpv.ccp.col.meta_data._type)
            end
        else
            error("unsupported page type $typ")
        end
    end
    iterator_result = NamedTuple{(:value,:repn_level,:defn_level),Tuple{OutputState{T},OutputState{Int32},OutputState{Int32}}}((ccpv.vals_out, ccpv.repn_out, ccpv.defn_out))
    iterator_result, nextpos
end

##
# layer 2 access
# can access decoded values from pages
function read_levels(inp::InputState, out::OutputState{Int32}, max_val::Int, enc::Int32, num_values::Int32) # levels are always 32 bits
    bit_width = UInt8(@bitwidth(max_val))
    @assert(bit_width !== 0)
    #@debug("reading levels. enc:$enc ($(Thrift.enumstr(Encoding,enc))), max_val:$max_val, num_values:$num_values")

    if enc === Encoding.RLE
        byte_width = @bit2bytewidth(bit_width)
        read_hybrid(inp, out, num_values, bit_width, byte_width)
    elseif enc === Encoding.BIT_PACKED
        read_bitpacked_run_old(inp, out, num_values, bit_width)
    else
        error("unsupported encoding $enc ($(Thrift.enumstr(Encoding,enc))) for levels")
    end
end

function read_levels_and_nmissing(inp::InputState, defn_out::OutputState{Int32}, repn_out::OutputState{Int32}, defn_enc::Int32, repn_enc::Int32, max_defn::Int, max_repn::Int, num_values::Int32)
    # read repetition levels. skipped if all columns are at 1st level
    if !isempty(repn_out.data)
        read_levels(inp, repn_out, max_repn, repn_enc, num_values)
    end

    # read definition levels. skipped if column is required
    nmissing = Int32(0)
    if !isempty(defn_out.data)
        defn_levels = defn_out.data
        defn_offset = defn_out.offset
        read_levels(inp, defn_out, max_defn, defn_enc, num_values)
        @inbounds for idx in 1:num_values
            (defn_levels[idx+defn_offset] === Int32(0)) && (nmissing += Int32(1))
        end
    end

    nmissing
end

# column and page metadata
open(par::Parquet.File, col::ColumnChunk) = open(par.handle, par.path, col)
close(par::Parquet.File, col::ColumnChunk, io) = (par.handle == io) || close(io)
function open(io, path::AbstractString, col::ColumnChunk)
    if hasproperty(col, :file_path)
        @debug("opening file to read column metadata", file=col.file_path, offset=col.file_offset)
        open(col.file_path)
    else
        if io === nothing
            @debug("opening file to read column metadata", file=path, offset=col.file_offset)
            open(path)
        else
            @debug("reading column metadata", offset=col.file_offset)
            io
        end
    end
end

function metadata(io, path::AbstractString, col::ColumnChunk)
    fio = open(io, path, col)
    seek(fio, col.file_offset)
    meta = read_thrift(fio, ColumnMetaData)
    (fio !== io) && close(fio)
    meta
end

function page_offset(par::Parquet.File, col::ColumnChunk)
    colmeta = metadata(par, col)
    offset = colmeta.data_page_offset
    hasproperty(colmeta, :index_page_offset) && (offset = min(offset, colmeta.index_page_offset))
    hasproperty(colmeta, :dictionary_page_offset) && (offset = min(offset, colmeta.dictionary_page_offset))
    offset
end
end_offset(par::Parquet.File, col::ColumnChunk) = page_offset(par, col) + metadata(par,col).total_compressed_size

page_size(page::PageHeader) = hasproperty(page, :compressed_page_size) ? page.compressed_page_size : page.uncompressed_page_size

const INVALID_ENC = Int32(-1)
page_encodings(page::Page) = page_encodings(page.hdr)
function page_encodings(page::PageHeader)
    hasproperty(page, :data_page_header) ? page_encodings(page.data_page_header) :
    hasproperty(page, :data_page_header_v2) ? page_encodings(page.data_page_header_v2) :
    hasproperty(page, :dictionary_page_header) ? page_encodings(page.dictionary_page_header) :
    (INVALID_ENC,INVALID_ENC,INVALID_ENC)
end
page_encodings(page::DictionaryPageHeader) = (page.encoding,INVALID_ENC,INVALID_ENC)
page_encodings(page::DataPageHeader) = (page.encoding, page.definition_level_encoding, page.repetition_level_encoding)
page_encodings(page::DataPageHeaderV2) = (page.encoding, Encoding.RLE, Encoding.RLE)

page_num_values(page::Page) = page_num_values(page.hdr)
function page_num_values(page::PageHeader)
    hasproperty(page, :data_page_header) ? page_num_values(page.data_page_header) :
    hasproperty(page, :data_page_header_v2) ? page_num_values(page.data_page_header_v2) :
    hasproperty(page, :dictionary_page_header) ? page_num_values(page.dictionary_page_header) : Int32(0)
end
page_num_values(page::Union{DataPageHeader,DataPageHeaderV2,DictionaryPageHeader}) = page.num_values

# file metadata
read_thrift(buff::Array{UInt8}, ::Type{T}) where {T} = read(TCompactProtocol(TMemoryTransport(buff)), T)
read_thrift(io::IO, ::Type{T}) where {T} = read(TCompactProtocol(TFileTransport(io)), T)
read_thrift(t::TR, ::Type{T}) where {TR<:TTransport,T} = read(TCompactProtocol(t), T)

function metadata_length(io)
    sz = filesize(io)
    seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER)

    # read footer size as little endian signed Int32
    read_fixed(io, Int32)
end

function metadata(io, path::AbstractString, len::Integer=metadata_length(io))
    @debug("reading file metadata", len)
    sz = filesize(io)
    seek(io, sz - SZ_PAR_MAGIC - SZ_FOOTER - len)
    meta = read_thrift(io, FileMetaData)
    meta
end

metadata(par::Parquet.File) = par.meta

#=
function fill_column_metadata(par::Parquet.File)
    meta = par.meta
    # go through all column chunks and read metadata from file offsets if required
    for grp in meta.row_groups
        for col in grp.columns
            metadata(par, col)
        end
    end
end
=#

function metadata(par::Parquet.File, col::ColumnChunk)
    if !hasproperty(col, :meta_data)
        col.meta_data = metadata(par.handle, par.path, col)
    end
    col.meta_data
end

# file format verification
function is_par_file(fname::AbstractString)
    open(fname) do io
        return is_par_file(io)
    end
end

function is_par_file(io)
    sz = filesize(io)
    (sz > SZ_VALID_PAR) || return false

    seekstart(io)
    magic = Array{UInt8}(undef, 4)
    read!(io, magic)
    (String(magic) == PAR_MAGIC) || return false

    seek(io, sz - SZ_PAR_MAGIC)
    magic = Array{UInt8}(undef, 4)
    read!(io, magic)
    (String(magic) == PAR_MAGIC) || return false

    true
end


text:
mutable struct Cell
    objs::Vector{DefectObj}
end

mutable struct CellContainer
    cells::Array{Cell}
    cellCutoff::Float64
    dCellIndexs::Vector{Vector{Int64}}
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mcke

max_forks_repo_name: permissionx/mcke


text:
#!/usr/bin/julia

# Copyright (C) 2022 Matheus Fernandes Bigolin <mfrdrbigolin@disroot.org>
# SPDX-License-Identifier: MIT

"Day Twenty-One, Dirac Dice."

include("Utils.jl")

using .Utils

"All the possible rolls of the deterministic die."
const DET_DIE = @. mod((10:-1:1) - 4, 10)

"""Using the deterministic dice and provided the position, `(P₁, P₂)`, return
the product of the losing player’s score with the number of dice rolls."""
function solve1((P₁, P₂))
    S₁ = S₂ = 0
    rolled = 0

    while S₁ < 1000 && S₂ < 1000
        for i in eachindex(DET_DIE)
            if S₁ >= 1000 || S₂ >= 1000
                break
            end

            rolled += 1

            if !iseven(i)
                P₁ += DET_DIE[i]
                S₁ += mod1(P₁, 10)
            else
                P₂ += DET_DIE[i]
                S₂ += mod1(P₂, 10)
            end
        end
    end

    return min(S₁, S₂)*3rolled
end

"The distribution of the rolls of a dirac die in a single round."
const DIRAC_DIE = [3 => 1, 4 => 3, 5 => 6, 6 => 7, 7 => 6, 8 => 3, 9 => 1]

"Storage of game states."
GMS = Dict{Tuple{Tuple{Int, Int}, Tuple{Int, Int}}, Tuple{Int, Int}}()

"""Using the dirac dice and given the initial position, `(P₁, P₂)`, return the
number of victories of each player."""
function solve2((P₁, P₂), (S₁, S₂) = (0, 0))
    key = ((P₁, P₂), (S₁, S₂))

    if S₁ >= 21
        return (1, 0)
    elseif S₂ >= 21
        return (0, 1)
    elseif haskey(GMS, key)
        return GMS[key]
    end

    vic = (0, 0)

    for (roll, n) in DIRAC_DIE
        P_new = mod1(P₁ + roll, 10)
        S_new = S₁ + P_new

        v₂, v₁ = solve2((P₂, P_new), (S₂, S_new)) .* n
        vic = vic .+ (v₁, v₂)
    end

    GMS[key] = vic
    return vic
end

const PLAYER_REX = r"Player \d starting position: (\d+)"

function main()
    usage_and_exit(length(ARGS) != 1)

    data = readlines(first(ARGS))
    players = [parse(Int, pl.captures[1]) for pl in match.(PLAYER_REX, data)]

    println(solve1(players))
    println(first(solve2(players)))
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
# function k2(u,v)
#
# 	û = [1-u[1], u[2]]
# 	v̂ = [1-v[1], v[2]]
#
# 	n1 = neighborhood(Testchart, û)
# 	n2 = neighborhood(Sourcechart, v̂)
# 	x = cartesian(n1)
# 	y = cartesian(n2)
# 	return(Kernel(x,y)*jacobian(n1)*jacobian(n2))
#
# end
#
#
#
#
# function k2p(u,v)
#
# 	û = [1-u[1], u[2]]
# 	v̂ = [1-v[1], v[2]]
#
# 	return(Kernel(û,v̂))
#
# end


text:
using Oceananigans
using CUDA
using Glob
using Test

include("utils_for_runtests.jl")

archs = test_architectures()

#####
##### Checkpointer tests
#####

function test_model_equality(test_model, true_model)
    CUDA.@allowscalar begin
        @test all(test_model.velocities.u.data     .≈ true_model.velocities.u.data)
        @test all(test_model.velocities.v.data     .≈ true_model.velocities.v.data)
        @test all(test_model.velocities.w.data     .≈ true_model.velocities.w.data)
        @test all(test_model.tracers.T.data        .≈ true_model.tracers.T.data)
        @test all(test_model.tracers.S.data        .≈ true_model.tracers.S.data)
        @test all(test_model.timestepper.Gⁿ.u.data .≈ true_model.timestepper.Gⁿ.u.data)
        @test all(test_model.timestepper.Gⁿ.v.data .≈ true_model.timestepper.Gⁿ.v.data)
        @test all(test_model.timestepper.Gⁿ.w.data .≈ true_model.timestepper.Gⁿ.w.data)
        @test all(test_model.timestepper.Gⁿ.T.data .≈ true_model.timestepper.Gⁿ.T.data)
        @test all(test_model.timestepper.Gⁿ.S.data .≈ true_model.timestepper.Gⁿ.S.data)
        @test all(test_model.timestepper.G⁻.u.data .≈ true_model.timestepper.G⁻.u.data)
        @test all(test_model.timestepper.G⁻.v.data .≈ true_model.timestepper.G⁻.v.data)
        @test all(test_model.timestepper.G⁻.w.data .≈ true_model.timestepper.G⁻.w.data)
        @test all(test_model.timestepper.G⁻.T.data .≈ true_model.timestepper.G⁻.T.data)
        @test all(test_model.timestepper.G⁻.S.data .≈ true_model.timestepper.G⁻.S.data)
    end
    return nothing
end

"""
Run two coarse rising thermal bubble simulations and make sure

1. When restarting from a checkpoint, the restarted model matches the non-restarted
   model to machine precision.

2. When using set!(test_model) to a checkpoint, the new model matches the non-restarted
   simulation to machine precision.

3. run!(test_model, pickup) works as expected
"""
function test_thermal_bubble_checkpointer_output(arch)

    #####
    ##### Create and run "true model"
    #####

    Nx, Ny, Nz = 16, 16, 16
    Lx, Ly, Lz = 100, 100, 100
    Δt = 6

    grid = RectilinearGrid(size=(Nx, Ny, Nz), extent=(Lx, Ly, Lz))
    closure = IsotropicDiffusivity(ν=4e-2, κ=4e-2)
    true_model = NonhydrostaticModel(architecture=arch, grid=grid, closure=closure,
                                     buoyancy=SeawaterBuoyancy(), tracers=(:T, :S))
    test_model = deepcopy(true_model)

    # Add a cube-shaped warm temperature anomaly that takes up the middle 50%
    # of the domain volume.
    i1, i2 = round(Int, Nx/4), round(Int, 3Nx/4)
    j1, j2 = round(Int, Ny/4), round(Int, 3Ny/4)
    k1, k2 = round(Int, Nz/4), round(Int, 3Nz/4)
    CUDA.@allowscalar true_model.tracers.T.data[i1:i2, j1:j2, k1:k2] .+= 0.01

    true_simulation = Simulation(true_model, Δt=Δt, stop_iteration=5)

    checkpointer = Checkpointer(true_model, schedule=IterationInterval(5), force=true)
    push!(true_simulation.output_writers, checkpointer)

    run!(true_simulation) # for 5 iterations

    checkpointed_model = deepcopy(true_simulation.model)

    true_simulation.stop_iteration = 9
    run!(true_simulation) # for 4 more iterations

    #####
    ##### Test `set!(model, checkpoint_file)`
    #####

    set!(test_model, "checkpoint_iteration5.jld2")

    @test test_model.clock.iteration == checkpointed_model.clock.iteration
    @test test_model.clock.time == checkpointed_model.clock.time
    test_model_equality(test_model, checkpointed_model)

    # This only applies to QuasiAdamsBashforthTimeStepper:
    @test test_model.timestepper.previous_Δt == checkpointed_model.timestepper.previous_Δt

    #####
    ##### Test pickup from explicit checkpoint path
    #####

    test_simulation = Simulation(test_model, Δt=Δt, stop_iteration=9)

    # Pickup from explicit checkpoint path
    run!(test_simulation, pickup="checkpoint_iteration0.jld2")

    @info "Testing model equality when running with pickup=checkpoint_iteration0.jld2."
    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
    @test test_simulation.model.clock.time == true_simulation.model.clock.time
    test_model_equality(test_model, true_model)

    run!(test_simulation, pickup="checkpoint_iteration5.jld2")
    @info "Testing model equality when running with pickup=checkpoint_iteration5.jld2."

    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
    @test test_simulation.model.clock.time == true_simulation.model.clock.time
    test_model_equality(test_model, true_model)

    #####
    ##### Test `run!(sim, pickup=true)
    #####

    # Pickup using existing checkpointer
    test_simulation.output_writers[:checkpointer] =
        Checkpointer(test_model, schedule=IterationInterval(5), force=true)

    run!(test_simulation, pickup=true)
    @info "    Testing model equality when running with pickup=true."

    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
    @test test_simulation.model.clock.time == true_simulation.model.clock.time
    test_model_equality(test_model, true_model)

    run!(test_simulation, pickup=0)
    @info "    Testing model equality when running with pickup=0."

    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
    @test test_simulation.model.clock.time == true_simulation.model.clock.time
    test_model_equality(test_model, true_model)

    run!(test_simulation, pickup=5)
    @info "    Testing model equality when running with pickup=5."

    @test test_simulation.model.clock.iteration == true_simulation.model.clock.iteration
    @test test_simulation.model.clock.time == true_simulation.model.clock.time
    test_model_equality(test_model, true_model)

    rm("checkpoint_iteration0.jld2", force=true)
    rm("checkpoint_iteration5.jld2", force=true)

    return nothing
end

function run_checkpointer_cleanup_tests(arch)
    grid = RectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))
    model = NonhydrostaticModel(architecture=arch, grid=grid,
                                buoyancy=SeawaterBuoyancy(), tracers=(:T, :S)
                                )

    simulation = Simulation(model, Δt=0.2, stop_iteration=10)

    simulation.output_writers[:checkpointer] = Checkpointer(model, schedule=IterationInterval(3), cleanup=true)
    run!(simulation)

    [@test !isfile("checkpoint_iteration$i.jld2") for i in 1:10 if i != 9]
    @test isfile("checkpoint_iteration9.jld2")

    rm("checkpoint_iteration9.jld2", force=true)

    return nothing
end

for arch in archs
    @testset "Checkpointer [$(typeof(arch))]" begin
        @info "  Testing Checkpointer [$(typeof(arch))]..."
        test_thermal_bubble_checkpointer_output(arch)
        run_checkpointer_cleanup_tests(arch)
    end
end


text:
# module Bukdu

module Tag

export form_for, label, text_input, select, textarea, submit

import ....Bukdu
import ....Bukdu.Octo: Changeset, change

function tag_id(f, field)
    string(lowercase(string(typeof(f))), '_', field)
end

function tag_name(f, field)
    string(lowercase(string(typeof(f))), '[', field, ']')
end

function field_id(name, field)
end

function field_name(name, field)
end

function field_value(form, field; selected=false)
end

function build(tag::String, changeset::Union{Void,Changeset}, field::Symbol, opts; body=nothing, LF=false)
    result = "<$tag"
    for (k,v) in opts
        result *= string(' ', k, "=\"", isa(v, Function) ? v(changeset.model, field) : v, '"')
    end
    linefeed = LF ? "\n" : ""
    string(result, isa(body, Void) ? " />" : ">$linefeed$body</$tag>")
end

function form_for(block::Function, changeset::Changeset, args...; kw...)
    list = Vector{Tuple{Symbol,Union{Function,String}}}(vcat(args..., kw...))
    opts = Dict(list)
    if haskey(opts, :method)
        verb = opts[:method]
    else
        verb = get
        push!(list, (:method, get))
    end
    sym_accept_charset = Symbol("accept-charset")
    if !haskey(opts, sym_accept_charset)
        push!(list, (sym_accept_charset, "utf-8"))
    end

    function form_for_args(action, list::Vector)
        vec = Vector()
        for (k,v) in list
            if :action == k
                push!(vec, k=>action)
            elseif :method == k
                method = string(Base.function_name(v))
                push!(vec, k=>method)
            else
                push!(vec, k=>v)
            end
        end
        vec
    end

    build_form = (action) -> build("form", changeset, :form, form_for_args(action, list); body=block(changeset), LF=true)

    action = opts[:action]
    if isa(action, Function)
        for route in Bukdu.RouterRoute.routes
            if route.action == action && route.verb == verb
                return build_form(route.path)
            end
        end
        action_name = Base.function_name(action)
        method = string(Base.function_name(verb))
        throw(Bukdu.NoRouteError("not found a route for $action_name $method"))
    else
        return build_form(action)
    end
end

function value_from_changeset(changeset::Changeset, field::Symbol, value=nothing)
    if isa(value, Void)
        if haskey(changeset.changes, field)
            return changeset.changes[field]
        elseif field in fieldnames(typeof(changeset.model))
            return getfield(changeset.model, field)
        end
    end
    return nothing
end

function label(changeset::Changeset, field::Symbol, body=nothing; kw...)
    build("label", changeset, field, (:for=>tag_id, kw...); body=body)
end

function text_input(changeset::Changeset, field::Symbol, value=nothing; kw...)
    value = value_from_changeset(changeset, field, value)
    build("input", changeset, field, (:id=>tag_id, :name=>tag_name, :type=>"text", :value=>value, kw...))
end

function select_option(changeset::Changeset, field::Symbol, options, value=nothing)
    # broadcast #
    # string(join(string.("    <option value=\"", options, "\">", options, "</option>"), '\n'), '\n')
    selected_value = value_from_changeset(changeset, field, value)
    string(join(map(x->
        string("    <option value=\"", x, '"', (selected_value==x ? " selected" : ""), ">", x, "</option>"), options), '\n'), '\n')
end

function select(changeset::Changeset, field::Symbol, options, value=nothing; kw...)
    build("select", changeset, field, (:id=>tag_id, :name=>tag_name, kw...); body=select_option(changeset, field, options, value), LF=true)
end

function textarea(changeset::Changeset, field::Symbol, value=""; kw...)
    build("textarea", changeset, field, (:id=>tag_id, :name=>tag_name, kw...); body=value, LF=true)
end

function submit(value; kw...)
    build("input", nothing, :submit, (:type=>"submit", :value=>value, kw...))
end

end # module Bukdu.Tag


text:
using KUnet
KUnet.atype(Array)
KUnet.ftype(Float64)
nc = 3
nd = 5
x = nothing
z = nothing
net = [Mmul(zeros(nc,nd)),PercLoss()]
for i=1:2
    x = floor(10*rand(nd,1))
    z = zeros(nc,1); z[rand(1:nc)]=1;
    @show x'
    @show z'
    @show forw(net, x)'
    @show back(net, z)
    update!(net)
    @show net[1].w.data
    @show net[1].w.diff
end


text:
module FittingWrapper

using JavaCall

const AbstractRealDistribution = @jimport org.hipparchus.distribution.continuous.AbstractRealDistribution
const Class = @jimport java.lang.Class
const EmpiricalDistribution = @jimport org.hipparchus.stat.fitting.EmpiricalDistribution
const File = @jimport java.io.File
const JString = @jimport java.lang.String
const List = @jimport java.util.List
const MixtureMultivariateNormalDistribution = @jimport org.hipparchus.distribution.multivariate.MixtureMultivariateNormalDistribution
const MultivariateNormalMixtureExpectationMaximization = @jimport org.hipparchus.stat.fitting.MultivariateNormalMixtureExpectationMaximization
const Object = @jimport java.lang.Object
const RandomGenerator = @jimport org.hipparchus.random.RandomGenerator
const StatisticalSummary = @jimport org.hipparchus.stat.descriptive.StatisticalSummary
const URL = @jimport java.net.URL

export EmpiricalDistribution
export MultivariateNormalMixtureExpectationMaximization
export cumulative_probability
export density
export equals
export estimate
export fit
export get_bin_count
export get_bin_stats
export get_class
export get_fitted_model
export get_generator_upper_bounds
export get_log_likelihood
export get_next_value
export get_numerical_mean
export get_numerical_variance
export get_sample_stats
export get_support_lower_bound
export get_support_upper_bound
export get_upper_bounds
export hash_code
export inverse_cumulative_probability
export is_loaded
export is_support_connected
export load
export log_density
export notify
export notify_all
export probability
export re_seed
export reseed_random_generator
export to_string
export wait

include("empirical_distribution.jl")
include("multivariate_normal_mixture_expectation_maximization.jl")

end


text:
f(x, y...) = y
@assert f(1, 2, 3) == (2, 3)

g(x; y...) = y
@assert g(1, a=2, b=5) == [(:a,2),(:b,5)]


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

using Markdown
using InteractiveUtils

# ╔═╡ bb28a1be-25cb-11ec-3a42-a3577c700396
begin
	import Pkg;Pkg.add( url="https://github.com/grahamstark/ScottishTaxBenefitModel.jl")
	using ScottishTaxBenefitModel

	using Plots
end

# ╔═╡ b9c7ccae-afbd-4e15-b508-44bc0803bdab
begin
	using Test
    using .ModelHousehold
	using .STBParameters
	using .STBIncomes
	using .Definitions
	using .GeneralTaxComponents
	using .SingleHouseholdCalculations
	using .RunSettings
	using .Utils
end

# ╔═╡ 9b44570c-bd0f-427a-bdbb-497737ede8cc
begin
	# sys21_22 = load_file( "../params/sys_2021_22.jl" )
	using DataFrames,CSV, Dates
	#
# full dataset is available ..
#
const IS_LOCAL = isdir("/mnt/data/frs/")

# pids for example people
# see make_pid
const RUK_PERSON = 320190010101

const SCOT_HEAD = 320190010201
const SCOT_SPOUSE = 320190010202

const DEFAULT_NUM_TYPE = Float64

const TEST_BASE_DATE = Date( 2021, 06, 14 ) # All the tests ran on this date; `TEST_BASE_DATE` is used
   # to adjust some ages e.g. social_security_ages_tests.jl so tests also pass next year.

function get_default_it_system(
   ;
  year     :: Integer=2019,
  scotland :: Bool = true,
  weekly   :: Bool = true )::Union{Nothing,IncomeTaxSys}
  it = nothing
  if year == 2019
     it = IncomeTaxSys{DEFAULT_NUM_TYPE}()
     if ! scotland
        it.non_savings_rates = [20.0,40.0,45.0]
        it.non_savings_thresholds = [37_500, 150_000.0]
        it.non_savings_basic_rate = 1
     end
     if weekly
        weeklyise!( it )
     end
  end
  it
end

function get_default_uc( ;
   year :: Integer = 2019,
   weekly :: Bool = true )
   uc = UniversalCreditSys{DEFAULT_NUM_TYPE}()
   if weekly
      weeklyise!( uc )
   else
      uc.taper /= 100
   end
   return uc
end

function to_nearest_p( x, y :: Real, ps :: Real = 1 ) :: Bool
   diff = abs(x-y)
   return diff <= 0.01*ps
   # round(x, digits=2) == round(y, digits=2)
end

function init_data(; reset :: Bool = false )
   nhh = get_num_households()
   num_people = -1
   if( nhh == 0 ) || reset
      @time nhh, num_people,nhh2 = initialise( Settings() )
   end
   (nhh,num_people)
end

function get_system(; scotland::Bool ) :: TaxBenefitSystem
   sys = TaxBenefitSystem{DEFAULT_NUM_TYPE}()
   weeklyise!(sys)
   # overwrite IT to get RuK system as needed
   # println( itn )
   sys.it = get_default_it_system( year=2019, scotland=scotland, weekly=true )
   #
   # So, assuming that the default has these things set ...
   #
   @assert sys.scottish_child_payment.amount > 0
   if ! scotland
      sys.scottish_child_payment.amount = 0.0
      sys.scottish_child_payment.maximum_age = 0
      sys.scottish_child_payment.qualifying_benefits = []
      sys.nmt_bens.carers.scottish_supplement = 0.0
   end
   return sys
end

@enum SS_Examples cpl_w_2_children_hh single_parent_hh single_hh childless_couple_hh mbu

function get_ss_examples()::Dict{SS_Examples, ModelHousehold.Household}
    d = Dict{SS_Examples, Household}()
    @time names = ExampleHouseholdGetter.initialise(settings=settings)
    d[cpl_w_2_children_hh] = ExampleHouseholdGetter.get_household( "example_hh1" )
    d[single_parent_hh] = ExampleHouseholdGetter.get_household( "single_parent_1" )
    d[single_hh] = ExampleHouseholdGetter.get_household( "example_hh2" )
    d[childless_couple_hh] = ExampleHouseholdGetter.get_household("mel_c2_scot")
    d[mbu] =  ExampleHouseholdGetter.get_household("mbu_example")
    return d
end

const EXAMPLES = get_ss_examples()
const SPARE_CHILD = EXAMPLES[cpl_w_2_children_hh].people[320190000104]
const SPARE_ADULT = get_head( EXAMPLES[single_hh])

function unemploy!( pers::Person )
   pers.usual_hours_worked = 0
   pers.actual_hours_worked = 0
   pers.employment_status = Unemployed
   delete!( pers.income, wages )
   delete!( pers.income, self_employment_income )

end

function employ!( pers::Person, wage=600.00 )
   pers.usual_hours_worked = 40
   pers.actual_hours_worked = 40
   pers.employment_status = Full_time_Employee
   pers.income[wages] = wage
end

function disable_slightly!( pers::Person )
   pers.employment_status = Permanently_sick_or_disabled
   pers.health_status = Bad
   pers.has_long_standing_illness = true
   pers.adls_are_reduced = reduced_a_little
   pers.how_long_adls_reduced = v_12_months_or_more
   pers.disabilities[mobility] = true
   pers.disabilities[stamina] = true
end

function disable_seriously!( pers::Person )
   pers.employment_status = Permanently_sick_or_disabled
   pers.health_status = Very_Bad
   pers.has_long_standing_illness = true
   pers.adls_are_reduced = reduced_a_lot
   pers.how_long_adls_reduced = v_12_months_or_more
   pers.disabilities[mobility] = true
   pers.disabilities[stamina] = true
end

function enable!( pers::Person )
   pers.dla_mobility_type = missing_lmh
   pers.dla_self_care_type = missing_lmh
   pers.pip_mobility_type = no_pip
   pers.pip_daily_living_type = no_pip
   pers.health_status = Good
   pers.has_long_standing_illness = false
   pers.adls_are_reduced = not_reduced
   pers.how_long_adls_reduced = Missing_Illness_Length
   pers.disabilities = Disability_Dict{Bool}()
end

function blind!( pers :: Person )
   pers.disabilities[vision ] = true
   pers.registered_blind = true
end

function unblind!( pers :: Person )
   delete!(pers.disabilities, vision )
   pers.registered_blind = false
end

function deafen!( pers :: Person )
   pers.disabilities[ hearing ] = true
   pers.registered_deaf = true
end

function undeafen!( pers :: Person )
   delete!(pers.disabilities, hearing )
   pers.registered_deaf = false
end

function carer!( pers :: Person )
   pers.income[carers_allowance] = 100.0
   pers.is_informal_carer = true
   pers.hours_of_care_given = 10
   pers.employment_status = Looking_after_family_or_home
end

function uncarer!( pers :: Person )
   delete!(pers.income,carers_allowance)
   pers.is_informal_carer = false
   pers.hours_of_care_given = 0
end

function retire!( pers :: Person )
   pers.usual_hours_worked = 0
   pers.employment_status = Retired
end

# FIXME relationships fixup
"""
Add a child to the 1st benefit unit
"""
function add_child!( hh :: ModelHousehold.Household, age :: Integer, sex :: Sex )::BigInt
   head = get_head(hh)
   np = deepcopy( SPARE_CHILD )
   empty!( np.income )
   np.relationships[head.pid] = Son_or_daughter_incl_adopted
   # TODO fill in other relationships
   np.default_benefit_unit = head.default_benefit_unit
   np.pid = maximum( keys( hh.people ))+1
   np.age = age
   np.sex = sex
   hh.people[ np.pid ] = np
   head.relationships[np.pid] = Parent
   spouse = get_spouse( hh )
   if spouse !== nothing
      spouse.relationships[np.pid] = Parent
   end
   # FIXME other adults in othe BUS
   make_eq_scales!( hh )
   return np.pid
end

function add_non_dependent!(
   hh  :: Household,
   age :: Integer,
   sex :: Sex ) :: BigInt

   head = get_head(hh)
   np = deepcopy( SPARE_ADULT )
   bus = get_benefit_units( hh )
   nbus = size(bus)[1]
   np.pid = maximum( keys( hh.people ))+1
   np.relationships[head.pid] = Son_or_daughter_incl_adopted
   # TODO fill in other relationships
   np.age = age
   np.sex = sex
   np.default_benefit_unit = nbus + 1
   hh.people[ np.pid ] = np
   bus = get_benefit_units( hh )
   nnbus = size(bus)[1]
   @assert nnbus == nbus + 1
   @assert get_head( bus[nnbus] ) == np
   make_eq_scales!( hh )
   return np.pid
end

function delete_person!( hh :: Household, pid :: BigInt )
   delete!( hh.people, pid )
   make_eq_scales!( hh )
end

function delete_child!( hh :: Household )
   chpids = child_pids( hh )
   if size(chpids)[1] > 0
      delete_person!( hh, chpids[1])
   end
   make_eq_scales!( hh )
end

"""
 if the test was written at TEST_BASE_DATE, what age would we have to make somebody
 to be sure the test will still work in some later year?
"""
function age_now( age :: Int ) :: Int
   yd = (Date(now()) - TEST_BASE_DATE).value ÷ 365 # leap years; no function for this
   return age + Int(yd)
end

function set_childrens_ages!( hh :: Household, ages ... )
   nc = length(ages)[1]
   nset = 0
   for (pid, pers) in hh.people
      if is_child(pers)
         nset = min( nc, nset+1)
         pers.age = ages[nset]
      end
   end
end

function init_benefit_unit_result( bu :: BenefitUnit ) :: BenefitUnitResult
   return init_benefit_unit_result( DEFAULT_NUM_TYPE, bu )
end

function make_hh(
   ;
   adults   :: Int = 1,
   children :: Int = 0,
   earnings :: Real = -1,
   rent     :: Real = -1,
   rooms    :: Int  = 4,
   age      :: Int = -1,
   spouse_age :: Int = -1,
   tenure   :: Tenure_Type = Private_Rented_Furnished ) :: Household
   hh = nothing
   if adults == 2
      if children > 0
         hh = deepcopy( EXAMPLES[cpl_w_2_children_hh])
      else
         hh = deepcopy( EXAMPLES[childless_couple_hh])
      end
   elseif adults == 1
      if children > 0
         hh = deepcopy( EXAMPLES[single_parent_hh])
      else
         hh = deepcopy( EXAMPLES[single_hh])
      end
   else
      error("can't do $adults adults yet")
   end
   hh.tenure = tenure
   hh.ct_band = Band_B
   num_kids = num_children( hh )
   if num_kids < children
      for i in (num_kids+1):children
         sex = (i % 2) == 0 ? Male : Female
         add_child!( hh, i, sex ) # use the counter as age
      end
   elseif num_kids > children
      delete_child!( hh )
   end
   nc = num_children( hh )
   na = num_adults( hh )
   @assert nc == children "num_childen=$nc but requested=$children"
   @assert na == adults "num_adults=$na but requested=$adults"
   head = get_head( hh )
   if age != -1
      head.age = age
   end
   if spouse_age != -1
      @assert na == 2 "need 2 adults for spouse age to be meaningful"
      sp = get_spouse( hh )
      sp.age = spouse_age
   end
   if earnings != -1
      head.income[wages] = earnings
   end
   if rent != -1
      hh.gross_rent = rent
   end
   make_eq_scales!( hh )
   return hh
end

#
# quickie for making a default-ish pid
#
function makePID( hid::Int, year = 2018, pno=1 )::BigInt
   get_pid( FRS, year, hid, pno )
end

#
# Messing around with periods to fit the Policy In Practice
# calculator.
#
const POLICY_IN_PRACTICE_WEEKS_PER_MONTH = 52/12
const PWPM = POLICY_IN_PRACTICE_WEEKS_PER_MONTH

"""
Compare PIP monthly calculations with our weekly, using their 52/12 weeks per month.
Throw an AssertionException if different, so we can use this directly
in a @test macro.
"""
function compare_w_2_m( uspw::Real, thempm::Real, ps :: Real = 1 ) :: Bool
   uspm = uspw*PWPM
   thempw = thempm/PWPM
   @assert to_nearest_p(uspw,thempw,ps) "us $(uspw)pw ($(uspm)pm) != $(thempm)pm"
   return true
end

sys21_22 = load_file( "../params/sys_2021_22.jl" )

	weeklyise!( sys21_22; wpy=52, wpm=PWPM  )
end

# ╔═╡ Cell order:
# ╠═bb28a1be-25cb-11ec-3a42-a3577c700396
# ╠═b9c7ccae-afbd-4e15-b508-44bc0803bdab
# ╠═9b44570c-bd0f-427a-bdbb-497737ede8cc

max_stars_repo_name: grahamstark/ScottishTaxBenefitModel

max_issues_repo_name: grahamstark/ScottishTaxBenefitModel

max_forks_repo_name: grahamstark/ScottishTaxBenefitModel


text:
function f_sel(t_o, aspl)

    # Number of time steps
    n_t = size(t_o)[1]

    # Interpolate time, aspl
    dt = 0.5
    n_ip = Int64(ceil((t_o[end]-t_o[1])/dt))

    T = eltype(aspl)
    t_ip = zeros(T, n_ip)
    for i in range(1, n_ip, step=1)
        t_ip[i] = t_o[1] + (i-1)*dt
    end

    # Interpolate the data
    f_aspl_interp = PCHIPInterpolation.Interpolator(t_o, reshape(aspl, (n_t, )))
    aspl_ip = f_aspl_interp.(t_ip)

    # Compute max. PNLT
    asplm = maximum(aspl_ip)

    # ICAO Annex 16 procedures (p132 - App. 2-18)
    f_int = 10 .^(aspl_ip / 10.)

    # Compute integration bounds
    D = 10 * log10(sum(f_int)) - asplm - 10 * log10(10. / dt)

    # Compute sel
    sel = asplm + D

    return sel
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
function multi_factorial(a)

    # function to compute multivariate factorials of the form
    # \bm{i!} = i_1!i_2!...i_N!

    return prod([factorial(b) for b in a])

end

function gamma_factorial(a, i)

    # function to compute gamma factors of the form
    # (α)_i = α(α+1)...(α+i-1)

    product = 1.0
    for j in 1:i
        product *= a+i-j
    end

    return product

end

function gamma_closure(sys::MomentEquations, binary_vars::Array{Int,1}=Int[])

    closure = OrderedDict() # leaving symbolic higher order terms
    closure_exp = OrderedDict() # expanding out all higher order terms fully

    N = sys.N
    if sys isa CentralMomentEquations
        M = copy(sys.M)
        μ = central_to_raw_moments(N, sys.m_order)
        μ_symbolic = define_μ(N, sys.q_order)
    else
        # consider covariances explicitly as M (without converting to raw moments)
        # and converting at the end leads to much more tractable expressions
        M = define_M(N, 2)
        μ = copy(sys.μ)
        μ_symbolic = copy(μ)
    end

    # obtain the shape and scale parameters
    # starting with symbolic α initially and substituting the respective
    # central/raw moments after the expressions are simplified

    #iter = Iterators.product(fill(1:N, 2)...)
    #α = define_temp_vars(:α, vcat(iter...))
    #x = define_temp_vars(:x, 1:N)

    #αs = Matrix{typeof(α[1])}(undef, N, N)
    #for (ind, (i,j)) in enumerate(iter)
    #    αs[i, j] = α[ind]
    #end

    @variables αs[1:N, 1:N]
    @variables x[1:N]

    αs = scalarize(αs)
    x = scalarize(x)
    for i in 1:N
        for j in i+1:N
            αs[j, i] = αs[i, j]
        end
    end

    for i in 1:N
        αs[i, i] = x[i]
        for k in 1:N
            if i != k
                αs[i, i] -= αs[i, k]
            end
        end
    end
    # the computations above only work with Num (maybe a bug?), hence `value` only now
    α = value.(αs)
    x = value.(x)

    symbolic_sub = Dict()
    β = Array{Any}(undef, N)
    for i in 1:N
        eᵢ = sys.iter_1[i]
        for j in i+1:N
            eⱼ = sys.iter_1[j]
            #symbolic_sub[αs[i, j]] = M[eᵢ .+ eⱼ] * μ[eᵢ] * μ[eⱼ] / M[2 .* eᵢ] / M[2 .* eⱼ]
            symbolic_sub[αs[i, j]] = M[eᵢ .+ eⱼ] * μ[eᵢ] * μ[eⱼ] * M[2 .* eᵢ]^-1 * M[2 .* eⱼ]^-1
        end
        #symbolic_sub[x[i]] = μ[eᵢ]^2 / M[2 .* eᵢ]
        symbolic_sub[x[i]] = μ[eᵢ]^2 * M[2 .* eᵢ]^-1
        #β[i] = M[2 .* eᵢ] / μ[eᵢ]
        β[i] = M[2 .* eᵢ] * μ[eᵢ]^-1

    end

    # construct the raw moments that follow the multivariate gamma distribution
    iters = Vector(undef, N)

    # TupleTools.sort comes in handy here
    unique_iter_q = unique(sort(i) for i in sys.iter_q)
    sub = Dict()

    iter_2nd_order = filter(x-> sum(x) == 2, sys.iter_m)
    # line below reproduces Lakatos et al. (2015) results - stronger assumptions
    #unique_iter_q = unique(sort(i) for i in vcat(sys.iter_m, sys.iter_q))

    for iter in unique_iter_q

        # construct iterator through the product of sums in the Eq. for raw moments
        for j in 1:N
            iter_j = Iterators.filter(x -> sum(x)==iter[j], sys.iter_all)
            iters[j] = iter_j
        end

        iter_k = Iterators.product(iters...)
        suma = 0.0

        multif_iter = multi_factorial(iter)
        for k in iter_k
            #factor1 = multif_iter / prod([multi_factorial(k_i) for k_i in k])
            factor1 = multif_iter * prod([multi_factorial(k_i) for k_i in k])^-1
            factor2 = 1.0
            for a in 1:N
                factor2 *= gamma_factorial(α[a, a], k[a][a])
                for b in a+1:N
                    factor2 *= gamma_factorial(α[a, b], k[a][b]+k[b][a])
                end
            end
            suma += factor1*factor2
        end

        suma = simplify(suma, expand=true)
        suma = substitute(suma, symbolic_sub)
        suma *= prod([β[j]^iter[j] for j in 1:N])
        μ[iter] = simplify(suma, expand=true)
        closure[μ_symbolic[iter]] = μ[iter]

        # Note that joint moments of multivariate distributions are symmetric in the sense
        # that, for example, M₀₃₀ is equivalent in form to M₃₀₀, only in their respective expressions
        # indices of all lower order moments have to be permuted correspondingly
        # let M₃₀₀ = μ₁₀₀+M₂₁₀, then M₀₃₀ = μ₀₁₀+M₁₂₀ (same form but indices were changed as 1->2 & 2->1)
        # Utilising the symmetry we can greatly reduce the computation time as no redundant
        # expression constructions/simplifications need to be performed anymore
        # TODO: implement for all closures?

        perms = collect(multiset_permutations(iter, length(iter)))[2:end]

        for iter_perm in perms

            iter_perm_ind = sortperm(sortperm(iter_perm))
            for i in sys.iter_1
                sub[μ_symbolic[i]] = μ_symbolic[i[iter_perm_ind]]
            end
            for i in iter_2nd_order
                sub[M[i]] = M[i[iter_perm_ind]]
            end
            #maybe faster to do substitution step by step without dict?
            iter_perm = Tuple(iter_perm)
            μ[iter_perm] = substitute(μ[iter], sub)
            closure[μ_symbolic[iter_perm]] = μ[iter_perm]
        end

    end

    # construct the corresponding truncated expressions of higher order
    # central moments from the obtained gamma raw moment expressions
    if sys isa CentralMomentEquations
        raw_to_central = raw_to_central_moments(N, sys.q_order, μ)
        central_to_raw = central_to_raw_moments(N, sys.q_order)
        closure_M = OrderedDict()
        for i in sys.iter_q
            closure_exp[M[i]] = raw_to_central[i]
            closure_M[M[i]] = simplify(closure[μ_symbolic[i]]-(central_to_raw[i]-M[i]))
        end
        closure = closure_M
    else
        raw_to_central = raw_to_central_moments(N, 2)
        M_to_μ = [M[i] => raw_to_central[i] for i in filter(x -> sum(x)==2, sys.iter_all)]

        # variance of bernoulli variable is expressed in terms of its mean
        #for ind in binary_vars
    #        eᵢ = sys.iter_1[ind]
    #        push!(M_to_μ, M[2 .* eᵢ] => μ[eᵢ] - μ[eᵢ]^2)
    #    end

        for i in sys.iter_q
            closure_exp[sys.μ[i]] = substitute(μ[i], M_to_μ)
            expr = substitute(closure[sys.μ[i]], M_to_μ)
            closure[sys.μ[i]] = simplify(expr)
        end
    end

    if isempty(binary_vars)
        return close_eqs(sys, closure_exp, closure, true)
    else
        # have to perform the bernoulli simplify in the end for gamma closure
        # as otherwise some symbolic terms do not cancel out and break it

        redundant_iter, redundant_eqs, iter_sub = bernoulli_reduce(sys, binary_vars)

        closed_eqs = Equation[]
        for (i, eq) in enumerate(get_eqs(sys.odes))
            if !(i in redundant_eqs)

                closed_rhs = substitute(eq.rhs, closure_exp)
                closed_rhs = expand(closed_rhs)
                closed_rhs = substitute(closed_rhs, iter_sub)
                closed_rhs = simplify(closed_rhs)

                push!(closed_eqs, Equation(eq.lhs, closed_rhs))

            end
        end

        iter_rm = intersect(sys.iter_q, redundant_iter)

        if sys isa CentralMomentEquations
            moments = M
        else
            moments = μ_symbolic
        end

        for i in iter_rm
            delete!(closure, moments[i])
        end

        vars = extract_variables(closed_eqs, N, sys.q_order)
        odename = Symbol(nameof(sys), "_gamma_closure")
        odes = ODESystem(closed_eqs, get_iv(sys.odes), vars, get_ps(sys.odes); name=odename)

        return ClosedMomentEquations(odes, closure, sys)

    end

end