text:
module Modems

export Modem, createmodem, modulate, demodulate,display,
       getbps,getM, destroy, count_biterrors, PSK, DPSK,
       ASK, QAM, phaseerror, evm

immutable Modem
	scheme::AbstractString
	M::Int
	modemObjPtr::Ptr{Void}
end

# Creates Modem object
# See http://liquidsdr.org/doc/modem/
# Uses liquid_getopt_str2mod(), modem_create(), modem_get_bps()
function createmodem(scheme::AbstractString)
	ms = ccall((:liquid_getopt_str2mod, "libliquid"), Int, (Ptr{UInt8},), scheme);
	if(ms==0)
		error("unknown modulation scheme : ", scheme );
	else
		modem = ccall((:modem_create, "libliquid"), Ptr{Void}, (Int,), ms);
		bps   = ccall((:modem_get_bps, "libliquid"), Int, (Ptr{Void}, ), modem);
		M = 1 << bps;
		return Modem(scheme, M, modem);
	end
end

function Base.display(modem::Modem)
	ccall((:modem_print, "libliquid"), Void, (Ptr{Void}, ), modem.modemObjPtr);
end

function getbps(modem::Modem)
	bps = ccall((:modem_get_bps, "libliquid"), Int, (Ptr{Void}, ), modem.modemObjPtr);
	return bps;
end

function getM(modem::Modem)
	return modem.M;
end

function destroy(modem::Modem)
	ccall((:modem_destroy, "libliquid"), Void, (Ptr{Void}, ), modem.modemObjPtr);
end

function modulate(modem::Modem , symbol::Int)
	y = Complex64[0];
	M = modem.M;
	if(symbol > (M-1))
		error("symbol exceeds constellation size (maximum = " , (M-1), ")");
	end
	ccall((:modem_modulate, "libliquid"), Void, (Ptr{Void}, Int, Ptr{Complex64}),
		modem.modemObjPtr , symbol, y);
	return y[1];
end

function modulate(modem::Modem, symbols::Array{Int})
	y = Array(Complex64,length(symbols));
	for i in range(1, length(symbols))
		y[i] = modulate(modem, symbols[i]);
	end
	return y;
end

function modulate(modem::Modem, symbols::Range{Int})
	y = Array(Complex64,length(symbols));
	index = 1
	for i in symbols
		y[index] = modulate(modem, i);
		index += 1;
	end
	return y;
end

function demodulate(modem::Modem, mod::Complex64)
	x = Array(Int, 1);
	ccall((:modem_demodulate, "libliquid"), Void, (Ptr{Void}, Complex64, Ptr{Int}),
		modem.modemObjPtr , mod, x );
	return x[1];
end

function demodulate(modem::Modem, mods::Array{Complex64})
	y = Array(Int,length(mods));
	for i in range(1, length(mods))
		y[i] = demodulate(modem.modemObjPtr, mods[i]);
	end
	return y;
end

function count_biterrors(x::Int, y::Int)
	err = ccall((:count_bit_errors, "libliquid"), Int, (Int, Int), x,y );
	return err;
end

function count_biterrors(x::Array{Int}, y::Array{Int})
	if(length(x) == length(y))
		errors = Array(Int, length(x))
		for i in range(1,length(x))
			errors[i] = count_biterrors(x[i], y[i]);
		end
		return errors;
	else
		error("x and y must be of same length")
	end
end

function phaseerror(modem::Modem)
	ccall((:modem_get_demodulator_phase_error, "libliquid"), Float64, (Ptr{Void},), modem.modemObjPtr);
end

function evm(modem::Modem)
	ccall((:modem_get_demodulator_evm, "libliquid"), Float64, (Ptr{Void}, ), modem.modemObjPtr);
end

end


text:
# Script to replicate table 2 (regression from equation 11)

using DataFrames, FixedEffectModels, RegressionTables, StatsBase, Statistics, LinearAlgebra, XLSX, PrettyTables

# -------------- Previous scripts -----------------------------------------------------------------------------------------------------------------------

include("src/data_wrangling.jl") # formats the raw input data obtained from WIOD

# -------------- Country-sector I/O-Table ---------------------------------------------------------------------------------------------------------------

function country_sector_level_GVC_indicators(df::DataFrame, year::Int64)

    ID, FD, GO, VA, IV, other = obtain_matrices(df, year) # N*S×N*S, N*S×N, N*S×1, N*S×1, N*S×N, N*S×1

    # net inventory correction
    iv_correction = [GO[i]/(GO[i]-IV[i, j]) for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys] # N*S×N, to match country wise inventory changes
    # notice difference in column size, need spread inventory correction of country column to country-sector columns
    # i.e. ID[1,1:35]/iv_correction[1,1], however problem with indexing => ceil(Int, 34/35)=1, ceil(Int, 36/35)=2
    ID_corrected = [ID[i, j]*iv_correction[i, ceil(Int, j/n_sectors)] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors] # N*S×N*S
    FD_corrected = [FD[i, j]*iv_correction[i, j] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys]
    FD_corrected = [sum(FD_corrected[i, :]) for i in 1:n_ctrys*n_sectors]

    replace!(ID_corrected, NaN => 0.0) # sometimes zeros are interpreted as NaN (particular for sector 35, private households with employed persons)
    replace!(FD_corrected, NaN => 0.0)

    # input/output coefficient matrices
    A = [ID_corrected[i, j] / GO[j] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors] # N*S×N*S, iterate output over destination country
    B = [ID_corrected[i, j] / GO[i] for i in 1:n_ctrys*n_sectors, j in 1:n_ctrys*n_sectors] # N*S×N*S, iterate output over reporting country

    replace!(A, NaN => 0.0) # sometimes zeros are interpreted as NaN (particular for sector 35, private households with employed persons)
    replace!(B, NaN => 0.0)

    # country-sector level statistics
    FD_GO = FD_corrected ./ GO # N*S×1
    VA_GO = VA ./ GO # N*S×1
    U_num = inv(I - A) * GO # N*S×1, numerator of equation (5)
    D_num = transpose(GO) * inv(I - B) # N*S×1, numerator of equation (7), notice difference in matrix algebra

    U = U_num ./ GO # N*S×1
    D = transpose(D_num) ./ GO # N*S×1, transpose to have column vector

    # collect data in a DataFrame
    summary_stats = hcat(fill(year, n_ctrys*n_sectors), repeat(ctrys, inner=n_sectors), repeat(1:35, outer=n_ctrys), FD_GO, VA_GO, U, D)
    t_country_sector_wide = DataFrame(summary_stats, :auto)
    rename!(t_country_sector_wide, [:year, :country, :sector, :FD_GO, :VA_GO, :U, :D])
    transform!(t_country_sector_wide, [:FD_GO, :VA_GO, :U, :D] .=> ByRow(x -> round(x, digits=3)) .=> [:FD_GO, :VA_GO, :U, :D])

    return t_country_sector_wide
end

# gather data for all years in a single data frame
t_tab2 = DataFrame(year=Int[], country=String[], sector=Int[], FD_GO=Float64[], VA_GO=Float64[], U=Float64[], D=Float64[])

for i in years
    append!(t_tab2, country_sector_level_GVC_indicators(df, i))
end

XLSX.writetable("clean/t_tab2.xlsx", t_tab2, overwrite=true) # export table

# -------------- Table A1 --------------------------------------------------------------------------------------------------------------------------------

# import data from local drive
t_tab2 = DataFrame(XLSX.readtable("clean/t_tab2.xlsx", "Sheet1")...) # still a lot of NaN
transform!(t_tab2, :country .=> ByRow(string) .=> :country, renamecols=false) # change type to String
transform!(t_tab2, [:year, :sector] .=> ByRow(Int) .=> [:year, :sector], renamecols=false) # change type to Int64
transform!(t_tab2, [:FD_GO, :VA_GO, :U, :D] .=> ByRow(Float64) .=> [:FD_GO, :VA_GO, :U, :D], renamecols=false) # change type to Float64

# missing: create table A.1 with summary statistics of country-industry GVC measures (how to treat NaN?)
# treat NaN = 0
t_tab2_long = stack(t_tab2, Not([:year, :country, :sector]))
transform!(t_tab2_long, :value => ByRow(x -> ifelse(isnan(x), 0.0, x)) => :value, renamecols=false) # set NaN = 0.0
transform!(t_tab2_long, :sector => ByRow(x -> ifelse(x in goods_sectors, "goods", "services")) => :sector_dummy, renamecols=false) # sector dummy for grouping
gdf = groupby(t_tab2_long, :variable)

t_tabA1_1 = combine(gdf, nrow, :value .=> [x -> percentile(x, 10), median, x -> percentile(x, 90), mean, std] .=> [:perc_10, :median, :perc_90, :mean, :std])
t_tabA1_1[:, :sector_dummy] .= "all industries"

gdf = groupby(t_tab2_long, [:variable, :sector_dummy])
t_tabA1_2 = combine(gdf, nrow, :value .=> [x -> percentile(x, 10), median, x -> percentile(x, 90), mean, std] .=> [:perc_10, :median, :perc_90, :mean, :std])

t_tabA1 = vcat(t_tabA1_2, t_tabA1_1) # combine tables
sort!(t_tabA1, [:variable, :sector_dummy])

open("images/summary_stats_country_industry_GVC.txt", "w") do f
    pretty_table(f,t_tabA1; backend = Val(:html), nosubheader=true, standalone=false) # only prints html code for table itself
end

# -------------- Table 2 --------------------------------------------------------------------------------------------------------------------------------

# use package: FixedEffectModels

# different ways of specifying the same regression, i.e. using either @formula, Symbols, Strings
# FixedEffectModels.reg(data_reg, @formula(FD_GO ~ year + fe(country_sector)), Vcov.cluster(:year, :country_sector))
# FixedEffectModels.reg(data_reg, term(:VA_GO) ~ term(:year) + fe(:country_sector), Vcov.cluster(:year, :country_sector))
# FixedEffectModels.reg(data_reg, term("VA_GO") ~ term("year") + fe(term("country_sector")), Vcov.cluster(:year, :country_sector))

#data_reg = subset(t_tab2, :VA_GO => ByRow(x -> !isnan(x))) # filter out rows with NaN
data_reg = transform(t_tab2, [:FD_GO, :VA_GO, :U, :D] .=> ByRow(x -> ifelse(isnan(x), 0.0, x)) .=> [:FD_GO, :VA_GO, :U, :D], renamecols=false) # set NaN = 0.0
# setting NaN = 0 considerably improves match with table 2
data_reg[:, :country_sector] .= data_reg[:, :country] .* "_" .* string.(data_reg[:, :sector]) # considerably improves R² compared to use columns individually

rr_FD_GO_1 = FixedEffectModels.reg(data_reg, @formula(FD_GO ~ year + fe(country_sector)), Vcov.cluster(:country_sector), save=true)
rr_FD_GO_2 = FixedEffectModels.reg(data_reg, @formula(FD_GO ~ year + fe(country_sector)), contrasts = Dict(:year => DummyCoding(base=1995)),
 Vcov.cluster(:country_sector), save=true)
rr_VA_GO_1 = FixedEffectModels.reg(data_reg, @formula(VA_GO ~ year + fe(country_sector)), Vcov.cluster(:country_sector), save=true)
rr_VA_GO_2 = FixedEffectModels.reg(data_reg, @formula(VA_GO ~ year + fe(country_sector)), contrasts = Dict(:year => DummyCoding(base=1995)),
 Vcov.cluster(:country_sector), save=true)
rr_U_1 = FixedEffectModels.reg(data_reg, @formula(U ~ year + fe(country_sector)), Vcov.cluster(:country_sector), save=true)
rr_U_2 = FixedEffectModels.reg(data_reg, @formula(U ~ year + fe(country_sector)), contrasts = Dict(:year => DummyCoding(base=1995)),
 Vcov.cluster(:country_sector), save=true)
rr_D_1 = FixedEffectModels.reg(data_reg, @formula(D ~ year + fe(country_sector)), Vcov.cluster(:country_sector), save=true)
rr_D_2 = FixedEffectModels.reg(data_reg, @formula(D ~ year + fe(country_sector)), contrasts = Dict(:year => DummyCoding(base=1995)),
 Vcov.cluster(:country_sector), save=true)

# use package: RegressionTables

# create .png out of ascii output to include in README.md as unformatted html output does not look good
RegressionTables.regtable(rr_FD_GO_1, rr_FD_GO_2, rr_VA_GO_1, rr_VA_GO_2, rr_U_1, rr_U_2, rr_D_1, rr_D_2;
 renderSettings = asciiOutput(), regression_statistics=[:nobs, :r2], print_fe_section=true, estimformat="%0.4f") # output on the console

# html output
RegressionTables.regtable(rr_FD_GO_1, rr_FD_GO_2, rr_VA_GO_1, rr_VA_GO_2, rr_U_1, rr_U_2, rr_D_1, rr_D_2;
 renderSettings = htmlOutput("images/table2_html.txt"), regression_statistics=[:nobs, :r2], print_fe_section=true, estimformat="%0.4f")

# regressions coincide very well, possible source of difference is the additional 300 observations of Antras and Chor (2018)
# have for VA, U, D (i.e. must have fewer NaN, but how?) => set NaN = 0

ext: jl


text:
mutable struct LinearApproximator <: AbstractApproximator
    basis
    param  # parameter
end

"""
n: input dimension
d: polynomial degree
m: output dimension
"""
function LinearApproximator(n::Int, d::Int, m::Int=1;
        with_bias=true, init_method=rand)
    basis = PolynomialBasis(n, d; with_bias=with_bias)
    w_length = basis(rand(n)) |> length
    param = m == 1 ? init_method(w_length) : init_method(w_length, m)
    LinearApproximator(basis, param)
end

function LinearApproximator(sym::Symbol)
    if sym == :tao_bian_nonlinear_VI_V̂
        basis = PolynomialBasis(2, 2; with_bias=false)
        param = basis(rand(2)) |> length |> zeros
        return LinearApproximator(basis, param)
    elseif sym == :tao_bian_nonlinear_VI_dV̂
        @polyvar x[1:2]
        @polyvar u
        Ψ = [
             x[1], x[2],  # deg: 1
             x[1]^2, x[1]*x[2], x[2]^2,  # deg: 2
             x[1]^3, x[1]^2 * x[2], x[1] * x[2]^2, x[2]^3,  # deg: 3
             x[1]^4, x[1]^3 * x[2], x[1]^2 * x[2]^2,
             x[1] * x[2]^3, x[2]^4,  # deg: 4
             x[1]*u^3, x[2]*u^3,
            ]
        basis = PolynomialBasis(Ψ, [x..., u])
        return LinearApproximator(basis, nothing)
    end
end

(a::LinearApproximator)(x) = a.param' * a.basis(x)


text:
using HorizonSideRobots

function mark_kross(r::Robot)
    for side1 in 0:3
        side2 = (side1+1)%4
        move_forward!(r, side1, side2)
        back!(r, inverse(side1),inverse(side2))
    end
    putmarker!(r)
end

function move_two!(r, side1, side2)
	move!(r, HorizonSide(side1))
	move!(r,  HorizonSide(side2))
end

function movements!(r, side1, side2)
    while !(isborder(r,HorizonSide(side1)) || isborder(r,HorizonSide(side2)))
        move_twice!(r, side1, side2)
        putmarker!(r)
    end
end

function back!(r, side1, side2)
    while ismarker(r)
        move_twice!(r, side1, side2)
    end
end

function inverse(side)
    return mod(side+2, 4)
end

mark_x_cross(r)

max_stars_repo_name: NikitaDestroyer/Basics

max_issues_repo_name: NikitaDestroyer/Basics

max_forks_repo_name: NikitaDestroyer/Basics


text:
"""
    struct GadgetCube{T} <: AbstractGadgetGeometry
        corner_lower_left::Vector{T}
        corner_upper_right::Vector{T}
    end

Defines a cube by `corner_lower_left` and `corner_upper_right`.
To be used for [`read_particles_in_geometry`](@ref)
"""
struct GadgetCube{T} <: AbstractGadgetGeometry
    corner_lower_left::Vector{T}
    corner_upper_right::Vector{T}
end

"""
    get_geometry_box_corners(cube::GadgetCube)

Returns a `Tuple` with the lower left and upper right corner of a box which contains the `cube`.
"""
function get_geometry_box_corners(cube::GadgetCube)
    cube.corner_lower_left, cube.corner_upper_right
end

"""
    function get_geometry_center(sphere::GadgetCube)

Returns the center of `cube`.
"""
get_geometry_center(cube::GadgetCube) =  1 // 2 .* cube.corner_lower_left .+ cube.corner_upper_right

"""
    get_geometry_mask(cube::GadgetCube, pos::Matrix{T}) where T

Returns the indices of all particles contained in the `cube`.
"""
function get_geometry_mask(cube::GadgetCube, pos::Matrix{T}) where T
    # by definition all particles are in the cube
    #return collect(1:size(pos,2))
    mask = @views @. ( ( cube.corner_lower_left[1] <= pos[1,:] <= cube.corner_upper_right[1] ) &
                       ( cube.corner_lower_left[2] <= pos[2,:] <= cube.corner_upper_right[2] ) &
                       ( cube.corner_lower_left[3] <= pos[3,:] <= cube.corner_upper_right[3] ) )

    return findall(mask)
end


text:
if R123_USE_AESNI
    key = rand(UInt128)
    r = ARS1x(key)
    r1 = ARS4x((unsafe_wrap(Array, Ptr{UInt32}(pointer_from_objref(key)), 4)...))
    @test r.x == rand(r1, UInt128)
    @test rand(r, UInt128) == rand(r1, UInt128)
end


text:
function load_pbmc68k(local_path)
    @load local_path data
    expr_mat = data[:all_data][Symbol("17820")][:hg19]

    obs = DataFrame(barcode=expr_mat[:barcodes])
    var = DataFrame(gene_id=expr_mat[:genes], gene_symbol=expr_mat[:gene_symbols])
    prof = Profile(SparseMatrixCSC(expr_mat[:mat]'), obs, var)
    prof
end


text:
import EarthOrientation
using ReferenceFrameRotations: angleaxis_to_dcm, ddcm

import ..sec2rad
using AstroTime: value, julian, julian_twopart, TTEpoch, UT1Epoch, TDBEpoch
using ..EarthAttitude:
    celestial_to_intermediate,
    earth_rotation_angle,
    polar_motion,
    s06,
    tio_locator,
    xy06

export
    CIRF,
    ITRF,
    TIRF,
    cirf,
    itrf,
    tirf

struct CIRF <: InertialFrame end
const cirf = CIRF()
from_sym(::Val{:CIRF}) = cirf

struct TIRF <: RotatingFrame end
const tirf = TIRF()
from_sym(::Val{:TIRF}) = tirf

struct ITRF <: RotatingFrame end
const itrf = ITRF()
from_sym(::Val{:ITRF}) = itrf

function Rotation(::ICRF, ::CIRF, ep::Epoch)
    m = precession_nutation(TTEpoch(ep))
    Rotation{icrf, cirf}(m')
end

function Rotation(::CIRF, ::ICRF, ep::Epoch)
    m = precession_nutation(TTEpoch(ep))
    Rotation{cirf, icrf}(m)
end

function Rotation(::CIRF, ::TIRF, ep::Epoch)
    ut1 = UT1Epoch(ep)
    era = earth_rotation_angle(value.(julian_twopart(ut1))...)
    rate = rotation_rate(earth, TDBEpoch(ep))
    m = angleaxis_to_dcm(era, [0, 0, 1])
    m′ = ddcm(m, [0, 0, rate])
    Rotation{cirf, tirf}(m, m′)
end

function Rotation(::TIRF, ::CIRF, ep::Epoch)
    ut1 = UT1Epoch(ep)
    era = earth_rotation_angle(value.(julian_twopart(ut1))...)
    rate = rotation_rate(earth, TDBEpoch(ep))
    m = angleaxis_to_dcm(-era, [0, 0, 1])
    m′ = ddcm(m, [0, 0, -rate])
    Rotation{tirf, cirf}(m, m′)
end

function polarmotion(ep::TTEpoch)
    xp, yp = EarthOrientation.polarmotion(value(julian(ep)))
    xp = sec2rad(xp)
    yp = sec2rad(yp)
    sp00 = tio_locator(value.(julian_twopart(ep))...)
    reshape(polar_motion(xp, yp, sp00), (3,3))
end

function precession_nutation(ep::TTEpoch)
    dx, dy = EarthOrientation.precession_nutation00(value(julian(ep)))
    jd1, jd2 = value.(julian_twopart(ep))
    x, y = xy06(jd1, jd2)
    s = s06(jd1, jd2, x, y)
    x += sec2rad(dx/1000.0)
    y += sec2rad(dy/1000.0)
    reshape(celestial_to_intermediate(x, y, s), (3,3))
end

function Rotation(::TIRF, ::ITRF, ep::Epoch)
    m = polarmotion(TTEpoch(ep))
    Rotation{tirf, itrf}(m')
end

function Rotation(::ITRF, ::TIRF, ep::Epoch)
    m = polarmotion(TTEpoch(ep))
    Rotation{itrf, tirf}(m)
end


text:
for (f, FT, arg) in ((:-, typeof(-), Number),
                     (:+, typeof(+), Real),
                     (:*, typeof(*), Real))
    @eval begin
        function Base.broadcasted(b::DefaultArrayStyle{1}, ::$FT, x::$arg, r::AbstractAxis)
            return maybe_unsafe_reconstruct(r, Base.broadcasted(b, $f, x, parent(r)); keys=keys(r))
        end
        function Base.broadcasted(b::DefaultArrayStyle{1}, ::$FT, r::AbstractAxis, x::$arg)
            return maybe_unsafe_reconstruct(r, Base.broadcasted(b, $f, parent(r), x); keys=keys(r))
        end
    end
end

#=
for (f, FT, arg) in ((:+, typeof(+), Real),)
    @eval begin
        function Base.broadcasted(::DefaultArrayStyle{1}, ::$FT, x::$arg, r::AbstractAxis)
            return unsafe_reconstruct(r, broadcast($f, x, eachindex(r)))
        end
        function Base.broadcasted(::DefaultArrayStyle{1}, ::$FT, r::AbstractAxis, x::$arg)
            return unsafe_reconstruct(r, broadcast($f, eachindex(r), x))
        end
    end
end
=#

function Broadcast.broadcast_shape(
    shape1::Tuple,
    shape2::Tuple{Vararg{<:AbstractAxis}},
    shapes::Tuple...
)

    return Broadcast.broadcast_shape(_bcs(shape1, shape2), shapes...)
end

function Broadcast.broadcast_shape(
    shape1::Tuple{Vararg{<:AbstractAxis}},
    shape2::Tuple,
    shapes::Tuple...
)

    return Broadcast.broadcast_shape(_bcs(shape1, shape2), shapes...)
end
function Broadcast.broadcast_shape(
    shape1::Tuple{Vararg{<:AbstractAxis}},
    shape2::Tuple{Vararg{<:AbstractAxis}},
    shapes::Tuple...
)

    return Broadcast.broadcast_shape(_bcs(shape1, shape2), shapes...)
end

# _bcs consolidates two shapes into a single output shape
_bcs(::Tuple{}, ::Tuple{}) = ()
_bcs(::Tuple{}, newshape::Tuple) = (newshape[1], _bcs((), tail(newshape))...)
_bcs(shape::Tuple, ::Tuple{}) = (shape[1], _bcs(tail(shape), ())...)
function _bcs(shape::Tuple, newshape::Tuple)
    return (combine_axis(first(shape), first(newshape)), _bcs(tail(shape), tail(newshape))...)
end

_parent_if_axis(x::AbstractAxis) = parent(x)

###
### cat
###
cat_indices(x, y) = set_length(indices(x), length(x) + length(y))

cat_keys(x::AbstractVector, y::AbstractRange) = StaticRanges.grow_last(y, length(x))

cat_keys(x::AbstractRange, y::AbstractVector) = StaticRanges.grow_last(x, length(y))

cat_keys(x::AbstractRange, y::AbstractRange) = StaticRanges.grow_last(x, length(y))

cat_keys(x::AbstractVector, y::AbstractVector) = cat_keys!(promote(x, y)...,)

cat_keys(x::AbstractVector{T}, y::AbstractVector{T}) where {T} = cat_keys!(copy(x), y)

function cat_keys!(x::AbstractVector{T}, y::AbstractVector{T}) where {T}
    for x_i in x
        if x_i in y
            error("Element $x_i appears in both collections in call to cat_axis!(collection1, collection2). All elements must be unique.")
        end
    end
    return vcat(x, y)
end

function cat_axis(x::Axis, y::AbstractUnitRange, inds=cat_indices(x, y))
    return maybe_unsafe_reconstruct(x, inds; keys=cat_keys(keys(x), keys(y)))
end

function cat_axis(x::Axis, y::Axis, inds=cat_indices(x, y))
    return maybe_unsafe_reconstruct(x, inds; keys=cat_keys(keys(x), keys(y)))
end

function cat_axis(x::AbstractUnitRange, y::Axis, inds=cat_indices(x, y))
    return maybe_unsafe_reconstruct(y, inds; keys=cat_keys(keys(x), keys(y)))
end

function cat_axis(x::AbstractUnitRange, y::AbstractUnitRange, inds=cat_indices(x, y))
    if x isa AbstractAxis
        return unsafe_reconstruct(x, inds)
    else
        return unsafe_reconstruct(y, inds)
    end
end

#=
    cat_axes(x::AbstractArray, y::AbstractArray, xy::AbstractArray, dims)

Produces the appropriate set of axes where `x` and `y` are the arrays that were
concatenated over `dims` to produce `xy`. The appropriate indices of each axis
are derived from from `xy`.
=#
@inline function cat_axes(x::AbstractArray, y::AbstractArray, xy::AbstractArray{T,N}, dims) where {T,N}
    ntuple(Val(N)) do i
        if i in dims
            cat_axis(axes(x, i), axes(y, i), axes(xy, i))
        else
            combine_axis(axes(x, i), axes(y, i), axes(xy, i))
        end
    end
end
# TODO do these work?
vcat_axes(x::AbstractArray, y::AbstractArray, xy::AbstractArray) = cat_axes(x, y, xy, 1)

hcat_axes(x::AbstractArray, y::AbstractArray, xy::AbstractArray) = cat_axes(x, y, xy, 2)

###
### combine
###
# LinearIndices indicates that keys are not formally defined so the collection
# that isn't LinearIndices is used. If both are LinearIndices then take the underlying
# OneTo as the new keys.
combine_keys(x, y) = _combine_keys(keys(x), keys(y))
_combine_keys(x, y) = promote_axis_collections(x, y)
_combine_keys(x,                y::LinearIndices) = x
_combine_keys(x::LinearIndices, y               ) = y
_combine_keys(x::LinearIndices, y::LinearIndices) = first(y.indices)

# TODO I still really don't like this solution but the result seems better than any
# alternative I've seen out there
#=
Problem: String, Symbol, Second, etc. don't promote neatly with Int (default key value)
Solution: Symbol(element), Second(element) for promotion
Exception: String(element) doesn't work so string(element) has to be used
=#

promote_axis_collections(x::X, y::X) where {X} = x

function promote_axis_collections(x::LinearIndices{1}, y::Y) where {Y}
    return promote_axis_collections(x.indices[1], y)
end
function promote_axis_collections(x::X, y::LinearIndices{1}) where {X}
    return promote_axis_collections(x, y.indices[1])
end
function promote_axis_collections(x::LinearIndices{1}, y::LinearIndices{1})
    return promote_axis_collections(x.indices[1], y.indices[1])
end
function promote_axis_collections(x::X, y::Y) where {X,Y}
    if promote_rule(X, Y) <: Union{}
        Z = promote_rule(Y, X)
    else
        Z = promote_rule(X, Y)
    end
    if Z <: Union{}
        Tx = eltype(X)
        Ty = eltype(Y)
        Tnew = promote_type(Tx, Ty)
        if Tnew == Any
            if is_key(Tx)
                if Tx <: AbstractString
                    return promote_axis_collections(x, string.(y))
                else
                    return promote_axis_collections(x, Tx.(y))
                end
            else
                if is_key(Ty)
                    if Ty <: AbstractString
                        return promote_axis_collections(string.(x), y)
                    else
                        return promote_axis_collections(Ty.(x), y)
                    end
                else
                    error("No method available for promoting keys of type $Tx and $Ty.")
                end
            end
        else
            return promote_axis_collections(Tnew.(x), y)
        end
    else
        return Z(x)
    end
end

###
### combine_axis
###
combine_axis(x, y) = combine_axis(x, y, _combine_indices(x, y))
combine_axis(x, y, inds) = SimpleAxis(inds)
combine_axis(x, y::AbstractAxis, inds) = combine_axis(y, x, inds)
combine_axis(x::AbstractAxis, y, inds) = unsafe_reconstruct(x, inds)
function combine_axis(x::AbstractAxis, y::AbstractAxis, inds)
    if y isa AbstractOffsetAxis
        return unsafe_reconstruct(x, unsafe_reconstruct(parent(y), inds))
    else
        return unsafe_reconstruct(x, unsafe_reconstruct(y, inds))
    end
end

combine_axis(x::SimpleAxis, y::AbstractAxis, inds) = combine_axis(parent(x), y, inds)
combine_axis(x::AbstractAxis, y::SimpleAxis, inds) = combine_axis(x, parent(y), inds)
combine_axis(x::SimpleAxis, y::SimpleAxis, inds) = combine_axis(parent(x), parent(y), inds)

# Axis
combine_axis(x::Axis, y, inds) = resize_last(x, inds)
combine_axis(x::Axis, y::SimpleAxis, inds) = combine_axis(x, parent(y), inds)
function combine_axis(x::Axis, y::AbstractAxis, inds)
    return resize_last(x, combine_axis(parent(x), y, inds))
end
# TODO check this stuff
function combine_axis(x::Axis, y::Axis, inds)
    return Axis(
        combine_keys(x, y),
        combine_axis(parent(x), parent(y), inds)
    )
end

_combine_indices(x, y) = One():_combine_length(static_length(x), static_length(y))
function _combine_length(::StaticInt{X}, ::StaticInt{Y}) where {X,Y}
    return StaticInt(_combine_length(X, Y))
end
_combine_length(::StaticInt{X}, y) where {X} = _combine_length(X, y)
_combine_length(x, ::StaticInt{Y}) where {Y} = _combine_length(x, Y)
function _combine_length(x, y)
    if x == y || y == 1
        return x
    elseif x == 1
        return y
    else
        throw(DimensionMismatch("arrays could not be broadcast to a common size;" *
                                " got a dimension with lengths $(x) and $(y)"))
    end
end


text:
function istriustrict(A::Matrix{T}) where {T}
    m, n = size(A)
    @inbounds for j = 1:min(n,m-1), i = j:m
        if A[i,j] ≠ zero(T)
            return false
        end
    end
    return true
end

function istrilstrict(A::Matrix{T}) where {T}
    m, n = size(A)
    @inbounds for j = 2:n, i = 1:min(j,m)
        if A[i,j] ≠ zero(T)
            return false
        end
    end
    return true
end

"Copy the first dimension of a 2D array y into a 1D array x."
function simd_copy_xy_first!(x::Array{T,1}, y::Array{T,2}, j) where {T}
    @assert length(x) == size(y, 1)
    @inbounds for i in eachindex(x)
        x[i] = y[i,j]
    end
    nothing
end

"Copy the first dimension of a 3D array y into a 1D array x."
function simd_copy_xy_first!(x::Array{T,1}, y::Array{T,3}, j, k) where {T}
    @assert length(x) == size(y, 1)
    @inbounds for i in eachindex(x)
        x[i] = y[i,j,k]
    end
    nothing
end

"Copy the first dimension of a 3D array y into a 1D array x."
function simd_copy_xy_first!(x::Array{T,2}, y::Array{T,3}, k) where {T}
    @assert size(x,1) == size(y, 1)
    @assert size(x,2) == size(y, 2)
    @inbounds for i in 1:size(x,1)
        for j=1:size(x,2)
            x[i,j] = y[i,j,k]
        end
    end
    nothing
end

"Copy a 1D array x into the first dimension of a 2D array y."
function simd_copy_yx_first!(x::Array{T,1}, y::Array{T,2}, j) where {T}
    @assert length(x) == size(y, 1)
    @inbounds for i in 1:size(y, 1)
        y[i,j] = x[i]
    end
    nothing
end

"Copy a 1D array x into the second dimension of a 2D array y."
function simd_copy_yx_second!(x::Array{T,1}, y::Array{T,2}, j) where {T}
    @assert length(x) == size(y, 2)
    @inbounds for i in 1:size(y, 2)
        y[j,i] = x[i]
    end
    nothing
end

"Copy a 1D array x into the first dimension of a 3D array y."
function simd_copy_yx_first!(x::Array{T,1}, y::Array{T,3}, j, k) where {T}
    @assert length(x) == size(y, 1)
    @inbounds for i in 1:size(y, 1)
        y[i,j,k] = x[i]
    end
    nothing
end

"Copy a 1D array x into the first dimension of a 3D array y."
function simd_copy_yx_first!(x::Array{T,2}, y::Array{T,3}, k) where {T}
    @assert size(x, 1) == size(y, 1)
    @assert size(x, 2) == size(y, 2)
    @inbounds for j in 1:size(y, 2)
        for i in 1:size(y, 1)
            y[i,j,k] = x[i,j]
        end
    end
    nothing
end

function simd_copy_yx_first_last!(x::Array{T,2}, y::Array{T,3}, j) where {T}
    @assert size(x, 1) == size(y, 1)
    @assert size(x, 2) == size(y, 3)
    @inbounds for k=1:size(y, 3)
        for i in 1:size(y, 1)
            y[i,j,k] = x[i,k]
        end
    end
    nothing
end

function simd_axpy!(a::T, x::DenseMatrix{T}, y::DenseMatrix{T}, e::DenseMatrix{T}) where {T}
    @assert size(x) == size(y) == size(e)

    local err::T
    local ty::T

    @inbounds for i in eachindex(x,y,e)
        err = e[i] + a*x[i]
        ty  = y[i]
        y[i] = err + ty
        e[i] = err + (ty - y[i])
    end
    nothing
end

function simd_aXbpy!(a, b::Vector{T}, X::Matrix{T}, y::Vector{T}) where {T}
    @assert length(y) == size(X, 1)
    @assert length(b) == size(X, 2)
    local ty::T
    @inbounds for i in eachindex(y)
        ty = zero(T)
        for j=eachindex(b)
            ty += X[i,j] * b[j]
        end
        y[i] = a*ty
    end
    nothing
end

function simd_abXpy!(a::T, b::Vector{T}, X::Matrix{T}, y::Vector{T}) where {T}
    @assert length(y) == size(X, 2)
    @assert length(b) == size(X, 1)
    local ty::T
    @inbounds for i in eachindex(y)
        ty = zero(T)
        for j=eachindex(b)
            ty += b[j] * X[j,i]
        end
        y[i] += a*ty
    end
    nothing
end

function simd_mult!(w::Vector{T}, X::Matrix, y::Vector) where {T}
    @assert length(w) == size(X, 1)
    @assert length(y) == size(X, 2)
    local tw::T
    @inbounds for i in eachindex(w)
        tw = zero(T)
        for j=eachindex(y)
            tw += X[i,j] * y[j]
        end
        w[i] = tw
    end
    nothing
end

function simd_mult!(w::Vector{T}, y::Vector, X::Matrix) where {T}
    @assert length(w) == size(X, 2)
    @assert length(y) == size(X, 1)
    local tw::T
    @inbounds for i in eachindex(w)
        tw = 0
        for j=eachindex(y)
            tw += y[j] * X[j,i]
        end
        w[i] = tw
    end
    nothing
end

function simd_mult!(W::Matrix{T}, X::Matrix, Y::Matrix) where {T}
    @assert size(W,1) == size(X, 1)
    @assert size(W,2) == size(Y, 2)
    @assert size(X,2) == size(Y, 1)
    local tw::T
    @inbounds for j in 1:size(W,2)
        for i in 1:size(W,1)
            tw = 0
            for k in 1:size(X,2)
                tw += X[i,k] * Y[k,j]
            end
            W[i,j] = tw
        end
    end
    nothing
end


text:
module ShaleDrillingLikelihood_DynamicDrillModel_ComputeFullPayoff

DOBTIME = false
PRINTSTUFF = false
DOPROFILE = false

using ShaleDrillingLikelihood
using Test
using SparseArrays
using BenchmarkTools
using Calculus
using Random
using Profile
# using ProfileView
using InteractiveUtils

using Base.Iterators: product, OneTo

using ShaleDrillingLikelihood: DCDPEmax,
    DCDPTmpVars,
    anticipate_t1ev,
    _dmax,
    solve_vf_all!,
    update_static_payoffs!,
    NoValueFunction,
    ValueFunction,
    ValueFunctionArrayOnly,
    value_function,
    EV,
    dEV,
    ObservationDrill,
    check_flow_grad,
    SimulationDraw,
    idx_ρ,
    actionspace,
    discounted_dynamic_payoff!,
    full_payoff!,
    update_interpolation!

println("print to keep from blowing up")

@testset "Dynamic Drilling Model" begin

    Random.seed!(1234)
    f = DrillReward(DrillingRevenue(Constrained(),NoTrend(),NoTaxes()), DrillingCost_constant(), ExtensionCost_Constant())
    # @show theta

    nψ, dmx, nz =  13, 3, 5

    # psi
    ψs = range(-4.5; stop=4.5, length=nψ)

    # z transition
    zs = ( range(-4.5; stop=4.5, length=nz), )
    nz = prod(length.(zs))
    ztrans = spdiagm(-1 => fill(1.0/3, nz-1), 0 => fill(1.0/3, nz), 1 => fill(1.0/3, nz-1) )
    ztrans[1,1] = ztrans[1,2] = ztrans[end,end-1] = ztrans[end,end] = 0.5

    # wp = LeasedProblemContsDrill(dmx,4,5,3,2)
    wp = LeasedProblem(dmx,4,5,3,2)

    # ichars
    ichar = (4.0, 0.25,)

    # ddm object
    ddm_no_t1ev   = DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, false)
    ddm_with_t1ev = DynamicDrillModel(f, 0.9, wp, zs, ztrans, ψs, true)

    @testset "VF Structs" begin
        vfao = ValueFunctionArrayOnly(ddm_no_t1ev)
        vf = ValueFunction(ddm_no_t1ev)
        @test EV(vfao) == EV(vf)
        @test dEV(vfao) == dEV(vf)
    end

    tmpv = DCDPTmpVars(ddm_no_t1ev)

    @testset "VF interpolation" begin

        num_parms = _nparm(f)
        gradtmp = zeros(num_parms)
        fdgrad = similar(gradtmp)

        function solve_vf_and_fullpayoff!(gradtmp, tmpv, ddm, d, obs, theta, uv, ichars, dograd=false)
            θρ = theta[idx_ρ(reward(ddm))]
            sim = SimulationDraw(uv..., θρ)
            # fill!(value_function(ddm), 0)
            solve_vf_all!(tmpv, ddm, theta, ichars, dograd)
            update_interpolation!(value_function(ddm), dograd)
            return full_payoff!(gradtmp, d, obs, theta, sim, dograd)
        end

        for ddm in (ddm_with_t1ev, ddm_no_t1ev)
            println("doing ddm with anticipate: $(ddm.anticipate_t1ev)")
            wp = statespace(ddm)
            for i in OneTo(length(wp))
                for d in actionspace(wp, i)

                    theta = randn(num_parms)
                    z = clamp(randn(), minimum(zs[1]), maximum(zs[1]))
                    obs = ObservationDrill(ddm, ichar, (z,), d, i)
                    uv = randn(2)

                    check_flow_grad(reward(ddm), d, obs, theta)

                    fill!(gradtmp, 0)
                    f(θ) = solve_vf_and_fullpayoff!(gradtmp, tmpv, ddm, d, obs, θ, uv, ichar, false)
                    fd = Calculus.gradient(f, theta)
                    @test all(gradtmp .== 0)

                    solve_vf_and_fullpayoff!(gradtmp, tmpv, ddm, d, obs, theta, uv, ichar, true)
                    @test all(isfinite.(gradtmp))
                    d > 0 && @test any(gradtmp .!= 0)
                    # fd ≈ gradtmp || @show fd .- gradtmp, i, d, z  # prefilter???
                    @test fd ≈ gradtmp
                end
            end
        end
    end

end

end # module


text:
VERSION > v"0.4-" && __precompile__()

module Lazy

############
# Utilities
############

include("macros.jl")
include("dynamic.jl")
include("tail.jl")

export @listable

import Base: *, ==, +, -

using Compat

macro listable(f)
  if typeof(f) == Expr && f.head == :tuple
    return Expr(:block, [:(@listable $f) for f in f.args]...)
  end
  f = esc(f)
  :($f(ls::List...) = map($f, ls...))
end

########
# Types
########

import Base: isempty, first, colon

export List, list, @lazy,
       prepend, rest

abstract List

type EmptyList <: List
end

type LinkedList <: List
  first
  rest::List
end

type LazyList <: List
  list::List
  realised::Bool
  f::Function

  LazyList(f) = new(EmptyList(), false, f)
end

# Empty List

first(::EmptyList) = nothing
rest(l::EmptyList) = l
isempty(::EmptyList) = true

# Lists

prepend(x, l::List) = LinkedList(x, l)
colon(x, xs::List) = prepend(x, xs)
colon(x,y,xs::List) = x:prepend(y,xs)

list() = EmptyList()
list(x, xs...) = x:list(xs...)

first(l::LinkedList) = l.first
rest(l::LinkedList) = l.rest

isempty(::LinkedList) = false

# Lazy Lists

realise(xs::LazyList) =
  xs.realised? xs.list : (xs.realised = true; xs.list = xs.f())

for f in [:first :rest :isempty]
  @eval $f(l::LazyList) = $f(realise(l))
end

macro lazy(code)
  :(LazyList(() -> seq($(esc(code)))))
end

########
# Usage
########

include("liblazy.jl")
include("collections.jl")

# -----
# Eager
# -----

export dorun, doall, foreach

@rec dorun(xs::List) = isempty(xs) ? nothing : dorun(rest(xs))
doall(xs::List) = (dorun(xs); xs)

foreach(f, ls::List...) = map(f, ls...) |> dorun

# -------
# Interop
# -------

import Base: getindex, setindex!, start, next, done

getindex(l::List, i::Int) = i <= 1 ? first(l) : rest(l)[i-1]
getindex(l::List, r::UnitRange) = take(r.len, drop(r.start - 1, l))

setindex!(xs::LinkedList, v, i::Integer) = i <= 1 ? xs.first = v : (rest(xs)[i-1] = v)
setindex!(xs::LazyList, v, i::Integer) = i <= 1 ? realise(xs)[1] = v : (rest(xs)[i-1] = v)

# Iteration over a list holds on to the head
start(xs::List) = xs
done(::List, xs::List) = isempty(xs)
next(::List, xs::List) = first(xs), rest(xs)

###########
# Printing
###########

Base.show(io::IO, xs::List) =
  foreach(x->print(io,x), ["("] * interpose(xs, " ") * [")"])

function Base.writemime(io::IO, ::MIME"text/plain", xs::List)
  isempty(xs) && return println(io, "List()")

  print(io, "List:")
  for x in xs
    print(io, "\n  ")
    show(io, x)
  end
end

# Some example lists

@listable +, -, Base.exp, Base.log, Base.sin, Base.cos, Base.tan

const fibs = @lazy 0:big(1):(fibs + drop(1, fibs))

isprime(n) =
  @>> primes begin
    takewhile(x -> x<=sqrt(n))
    map(x -> n % x == 0)
    any; !
  end

const primes = filter(isprime, range(2))

end


text:
module NEOSServer

@warn(
    "Make sure you comply with the NEOS terms of use: " *
    "http://www.neos-server.org/neos/termofuse.html",
)

import AmplNLWriter
import Base64
import HTTP
import LightXML

"""
	Server(email)

Construct a `Server` object. The `email` argument should take the users
valid email address (required for solvers like CPLEX).
"""
struct Server
    user_agent::String
    host::String
    content_type::String
    email::String
    function Server(email::String)
        return new("NEOS.jl", "https://neos-server.org:3333", "text/xml", email)
    end
end

"""
	struct Job
		number::Int
		password::String
	end

Struct representing a NEOS job.
"""
struct Job
    number::Int
    password::String
end

# ==============================================================================
# 	NEOS API Methods
#
# Taken from https://neos-server.org/neos/xml-rpc.html
# ==============================================================================

_add_text(value, arg::String) = LightXML.add_text(value, arg)
function _add_text(value, arg)
    return LightXML.add_text(LightXML.new_child(value, "int"), string(arg))
end

function _build_xml(name::String, args...)
    xml = LightXML.XMLDocument()
    mthd = LightXML.create_root(xml, "methodCall")
    mname = LightXML.new_child(mthd, "methodName")
    LightXML.add_text(mname, name)
    if length(args) > 0
        params = LightXML.new_child(mthd, "params")
        for a in args
            param = LightXML.new_child(params, "param")
            value = LightXML.new_child(param, "value")
            _add_text(value, a)
        end
    end
    return string(xml)
end

function _get_values!(values, node)
    types = ["int", "i4", "string", "double", "base64", "dateTime.iso8601"]
    for child in LightXML.child_nodes(node)
        if LightXML.name(child) in types
            push!(values, LightXML.content(child))
        else
            _get_values!(values, child)
        end
    end
end

function _api_method(s::Server, name::String, args...)
    xml = _build_xml(name, args...)
    headers = [
        "user-agent" => s.user_agent,
        "host" => s.host,
        "content-type" => s.content_type,
        "content-length" => string(length(xml)),
    ]
    res = HTTP.request("POST", s.host, headers, xml)
    if res.status != 200
        error("XML-RPC failed with code: $(res.status)")
    end
    parameters = Any[]
    xml = LightXML.parse_string(String(res.body))
    xroot = LightXML.root(xml)
    _get_values!(parameters, xroot)
    return parameters
end

function _decode_to_string(s)
    return String(Base64.base64decode(replace(s, "\n" => "")))
end

neos_help(s::Server) = _api_method(s, "help")[1]

neos_welcome(s::Server) = _api_method(s, "welcome")[1]

neos_version(s::Server) = _api_method(s, "version")[1]

neos_ping(s::Server) = _api_method(s, "ping")[1]

neos_printQueue(s::Server) = _api_method(s, "printQueue")[1]

function neos_getSolverTemplate(
    s::Server,
    category::String,
    solvername::String,
    inputMethod::String,
)
    ret = _api_method(s, "getSolverTemplate", category, solvername, inputMethod)
    return ret[1]
end

neos_listAllSolvers(s::Server) = _api_method(s, "listAllSolvers")

neos_listCategories(s::Server) = _api_method(s, "listCategories")

function neos_listSolversInCategory(s::Server, category::String)
    return _api_method(s, "listSolversInCategory", category)
end

function neos_submitJob(s::Server, xmlstring::String)
    res = _api_method(s, "submitJob", xmlstring)
    return Job(parse(Int, res[1]), res[2])
end

function neos_getJobStatus(s::Server, j::Job)
    return _api_method(s, "getJobStatus", j.number, j.password)[1]
end

function neos_getCompletionCode(s::Server, j::Job)
    return _api_method(s, "getCompletionCode", j.number, j.password)[1]
end

function neos_getJobInfo(s::Server, j::Job)
    return _api_method(s, "getJobInfo", j.number, j.password)
end

function neos_killJob(s::Server, j::Job)
    return _api_method(s, "killJob", j.number, j.password)
end

function neos_getFinalResults(s::Server, j::Job)
    ret = _api_method(s, "getFinalResults", j.number, j.password)
    return _decode_to_string(ret[1])
end

function neos_emailFinalResults(s::Server, j::Job)
    return _api_method(s, "emailFinalresults", j.number, j.password)
end

function neos_getIntermediateResults(s::Server, j::Job, offset::Int)
    ret = _api_method(s, "getIntermediateResults", j.number, j.password, offset)
    return _decode_to_string(ret[1]), parse(Int, ret[2])
end

function neos_getFinalResultsNonBlocking(s::Server, j::Job)
    ret = _api_method(s, "getFinalResultsNonBlocking", j.number, j.password)
    return _decode_to_string(ret[1])
end

function neos_getIntermediateResultsNonBlocking(s::Server, j::Job, offset::Int)
    ret = _api_method(
        s,
        "getIntermediateResultsNonBlocking",
        j.number,
        j.password,
        offset,
    )
    return _decode_to_string(ret[1]), parse(Int, ret[2])
end

function neos_getOutputFile(s::Server, j::Job, fileName::String)
    ret = _api_method(s, "getOutputFile", j.number, j.password, fileName)
    return _decode_to_string(ret[1])
end

# ==============================================================================
# 	NEOS.Optimizer interface
# ==============================================================================

const _SUPPORTED_SOLVERS = Dict(
    "Ipopt" => "nco",
    "SNOPT" => "nco",
    "Gurobi" => "milp",
    "FICO-Xpress" => "milp",
    "CPLEX" => "milp",
    "MOSEK" => "milp",
    "Knitro" => "nco",
    "COPT" => "milp",
    "OCTERACT" => "minco",
)

function Optimizer(; email::String, solver::String, kwargs...)
    category = get(_SUPPORTED_SOLVERS, solver, nothing)
    if category === nothing
        error(
            "NEOS.Optimizer only supports the following solvers: " *
            join(collect(keys(_SUPPORTED_SOLVERS)), ", "),
        )
    end
    cmd = _SolverCommand(solver, category, Server(email))
    return AmplNLWriter.Optimizer(cmd; kwargs...)
end

struct _SolverCommand <: AmplNLWriter.AbstractSolverCommand
    solver::String
    category::String
    server::Server
end

function AmplNLWriter.call_solver(
    solver::_SolverCommand,
    nl_filename::String,
    options::Vector{String},
    ::IO,
    stdout::IO,
)::String
    if length(options) > 0
        @warn("Options are not supported by NEOS.jl")
    end
    xml = """
    <document>
    <client>NEOS.jl</client>
    <solver>$(solver.solver)</solver>
    <category>$(solver.category)</category>
    <inputMethod>NL</inputMethod>
    <email>$(solver.server.email)</email>
    <model>$(read(nl_filename, String))</model>
    </document>
    """
    job = neos_submitJob(solver.server, xml)
    polling_period = 1.0
    offset = 0
    while true
        ret, new_offset =
            neos_getIntermediateResultsNonBlocking(solver.server, job, offset)
        if new_offset > offset
            println(stdout, ret)
        end
        offset = new_offset
        if neos_getJobStatus(solver.server, job) == "Done"
            break
        end
        sleep(polling_period)
        polling_period = min(60.0, 2 * polling_period)
    end
    sol = neos_getOutputFile(solver.server, job, "ampl.sol")
    sol_file = replace(nl_filename, "model.nl" => "model.sol")
    write(sol_file, sol)
    return sol_file
end

# ==============================================================================
# 	Export NEOS API
# ==============================================================================

for sym in names(@__MODULE__; all = true)
    if startswith("$(sym)", "neos_")
        @eval export $sym
    end
end

end


text:
using Performance
using BenchmarkTools
using Test

@testset "arrays" begin
    ##
    ## Arrays are good examples of containers for illustrating:
    ## https://docs.julialang.org/en/v1/manual/performance-tips/#man-performance-abstract-container
    ##
    seq3 = sequence_of_three(17)
    @test seq3 == [17, 18, 19]
    @test @inferred(seq3[2]) == 18
    seq3 = sequence_of_three(17.5)
    @test seq3 == [17.5, 18.5, 19.5]
    @test @inferred(seq3[2]) == 18.5

    # Improve on this "bad" implementation. The starter code in `src`
    # is identical to the bad one, but you should improve it.
    # If you need a hint, see https://julialang.org/blog/2018/08/union-splitting/
    function alternate_items_bad(nrepeats::Integer, item1, item2)
        out = []
        for i = 1:nrepeats
            push!(out, item1, item2)
        end
        return out
    end
    vbad = alternate_items_bad(10^4, 3, 4)
    bbad = @benchmark(sum($vbad), samples=100, evals=10)
    vgood = alternate_items(10^4, 3, 4)
    @test isconcretetype(eltype(vgood))
    bgood = @benchmark(sum($vgood), samples=100, evals=10)
    @test 10*minimum(bgood).time < minimum(bbad).time
    # Let's see what the actual performance ratio is (probably a *lot* better than 10x):
    println("arrays.jl: actual performance ratio: ", minimum(bbad).time / minimum(bgood).time)
    # But make sure your implementation can handle this
    v = alternate_items(10^4, missing, 4)
    @test @inferred(sum(skipmissing(v))) == 4*10^4
    # Note: while this example might seem artificial, there turn out
    # to be many cases where you want to create a struct that can hold
    # one of two types in a given field, or an array that can store
    # `nothing` or `missing` in addition to a concrete type.

end


text:
######################################################
#
# Binomial Cascade Plot
#
# Date:   2014-09-29
# Author: Hessel Tuinhof
#
# Plot several conservative and canonical binomial
# cascades.
#
# Required packages: 	PyPlot
# Required functions: 	BinomialConservative.jl
#					  	BinomialCanonical.jl
#
######################################################

# load packages & functions
using PyPlot
require("../src/cascades/BinomialConservative.jl")
require("../src/cascades/BinomialCanonical.jl")

# define colors
rouge = (0.79, 0.0, 0.1)
turquoise = (0, 0.41, 0.41)

# set seed
srand(314159265358979323846)

###
### Plot conservative binomial cascade (k = 5)
###

# configure parameters
m0 = 0.6
num_cascades = 5

# plot cascade
figure()
mu = BinomialConservative(m0, num_cascades)
y_coordinates = 2 * push!(mu, mu[end])
x_coordinates = linspace(0, 1, 2^num_cascades + 1)
step(x_coordinates, y_coordinates, where = "post", linewidth = 1, color = rouge)
ylim(ymin = 0, ymax = maximum(y_coordinates) + (maximum(y_coordinates)/100)*5) # add 5% to the maximum of the y-values

###
### Plot conservative binomial cascade (k = 10)
###

# configure parameters
m0 = 0.6
num_cascades = 10

# plot cascade
figure()
mu = BinomialConservative(m0, num_cascades)
y_coordinates = 2 * push!(mu, mu[end])
x_coordinates = linspace(0, 1, 2^num_cascades + 1)
step(x_coordinates, y_coordinates, where = "post", linewidth = 0.6, color = rouge)
ylim(ymin = 0, ymax = maximum(y_coordinates) + (maximum(y_coordinates)/100)*5) # add 5% to the maximum of the y-values

###
### Plot canonical binomial cascade (k = 5)
###

# configure parameters
m0 = 0.6
num_cascades = 5

# plot cascade
figure()
mu = BinomialCanonical(m0, num_cascades)
y_coordinates = 2 * push!(mu, mu[end])
x_coordinates = linspace(0, 1, 2^num_cascades + 1)
step(x_coordinates, y_coordinates, where = "post", linewidth = 1, color = turquoise)
ylim(ymin = 0, ymax = maximum(y_coordinates) + (maximum(y_coordinates)/100)*5) # add 5% to the maximum of the y-values

###
### Plot canonical binomial cascade (k = 10)
###

# configure parameters
m0 = 0.6
num_cascades = 10

# plot cascade
figure()
mu = BinomialCanonical(m0, num_cascades)
y_coordinates = 2 * push!(mu, mu[end])
x_coordinates = linspace(0, 1, 2^num_cascades + 1)
step(x_coordinates, y_coordinates, where = "post", linewidth = 0.6, color = turquoise)
ylim(ymin = 0, ymax = maximum(y_coordinates) + (maximum(y_coordinates)/100)*5) # add 5% to the maximum of the y-values

max_stars_repo_name: hesseltuinhof/masterthesis

max_issues_repo_name: hesseltuinhof/masterthesis

max_forks_repo_name: hesseltuinhof/masterthesis


text:
using HorizonSideRobots

function chessall(r::Robot)
    s = 0
    w = 0

    while !isborder(r,Sud)
        move!(r,Sud)
        s+=1
    end
    while !isborder(r,West)
        move!(r,West)
        w+=1
    end

    if w%2==0
        shax(r)
    else
        move!(r,Ost)
        shax(r)
    end

    while !isborder(r,Sud)
        move!(r,Sud)
    end
    while !isborder(r,West)
        move!(r,West)
    end

    for a in 1:w
        move!(r,Ost)
    end
    for a in 1:s
        move!(r,Nord)
    end
end

function shax(r::Robot)
    putmarker!(r)
    i = 0
    while !isborder(r,Nord)
        if !isborder(r,Sud)
            move!(r,Nord)
        end
        i = obr(r,i)
        i = stripe(r,Ost,i)
        if !isborder(r,Nord)
            move!(r,Nord)
            i = obr(r,i)
            i = stripe(r,West,i)
        end
    end
end

function obr(r,i)
    if i == 1
        i = 0
    else
        i = 1
        putmarker!(r)
    end
    return i
end

function stripe(r,side,i)
    while !isborder(r,side)
        if i == 1
            move!(r,side)
            i = 0
        else
            move!(r,side)
            putmarker!(r)
            i = 1
        end
    end
    return i
end

chessall(r)

max_stars_repo_name: Cheburkova/mirea-progs

max_issues_repo_name: Cheburkova/mirea-progs

max_forks_repo_name: Cheburkova/mirea-progs


text:
# Copyright 2016, 2017, 2018 Martin Holters
# See accompanying license file.

using ACME

# Model of "Der Birdie", see
# http://diy.musikding.de/wp-content/uploads/2013/06/birdieschalt.pdf
# for schematics.
# Note: Contrary to the schematics, in the real circuit as well as in this
# model, C5 and D1 are connected directly to the power supply; R6 is in series
# with LED1 only, which has been omitted from the model, though.

function birdie(::Type{Circuit}; vol=nothing)
    @circuit begin
        j3 = voltagesource(9), [-] ⟷ gnd, [+] ⟷ vcc # 9V power supply
        c5 = capacitor(100e-6), [1] ⟷ gnd, [2] ⟷ vcc
        d1 = diode(is=350e-12, η=1.6), [-] ⟷ vcc, [+] ⟷ gnd
        j1 = voltagesource(), [-] ⟷ gnd # input
        r1 = resistor(1e6), [1] ⟷ j1[+], [2] ⟷ gnd
        c1 = capacitor(2.2e-9), [1] ⟷ j1[+]
        r2 = resistor(43e3), [1] ⟷ c1[2], [2] ⟷ gnd
        r3 = resistor(430e3), [1] ⟷ c1[2], [2] ⟷ vcc
        t1 = bjt(:npn, isc=154.1e-15, ise=64.53e-15, ηc=1.10, ηe=1.06, βf=500, βr=12),
             [base] ⟷ c1[2]
        r4 = resistor(390), [1] ⟷ t1[emitter], [2] ⟷ gnd
        r5 = resistor(10e3), [1] ⟷ t1[collector], [2] ⟷ vcc
        c3 = capacitor(2.2e-9), [1] ⟷ t1[collector]
        p1 = potentiometer(100e3, (vol == nothing ? () : (vol,))...), [1] ⟷ gnd, [3] ⟷ c3[2]
        j2 = voltageprobe(), [-] ⟷ gnd, [+] ⟷ p1[2]# output
    end
end

birdie(::Type{DiscreteModel}=DiscreteModel; vol=nothing, fs=44100, solver=nothing) =
    solver === nothing ?
        DiscreteModel(birdie(Circuit, vol=vol), 1//fs) :
        DiscreteModel(birdie(Circuit, vol=vol), 1//fs, solver)


text:
function cleanup()::Void
  pToClean=procs()
  for i in 2:length(pToClean)
    rmprocs(pToClean[i])
  end

  return nothing
end

#adds processes as desired
function initializePar(workers::Int = WORKERS; reset::Bool=false, fixRandom::Bool = true)::Void
  #if the reset flag is active we kill all active processes besides the main
  #note this takes a while
  if reset == true
    cleanup()
  end

  #create the processes
  if nworkers() < workers || (nworkers()==1 && workers == 1)
    addprocs(nworkers()==1?workers-nworkers()+1:workers-nworkers())
  end

  #make sure all processes have access to everything udrf
  #@everywhere eval( :(using FLMethods, FLOptimizer, FLData))
  #@everywhere eval( :(using FLOptimizerPar))
  #@everywhere eval( :(using FLScripts))
  @everywhere eval( :(using Endowment))

  return nothing
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
#=
---------------------------------------------------------------------
This function computes the integer array ISIDE which is needed to
create the FACE arrays.

Written by F.X. Giraldo on July 12, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

function create_side(intma,bsido,Npoin,Ne,Nboun,Nface,Np,DFloat)

    #global arrays
    iside = zeros(Int64,4,Nface)
    jeside= zeros(Int64,4,Ne)

    #local arrays
    lwher = zeros(Int64,Npoin,1)
    lhowm = zeros(Int64,Npoin,1)
    icone = zeros(Int64,5*Npoin,1)
    inode = zeros(Int64,4,1)
    jnode = zeros(Int64,4,1)

    #Fix lnode
    inode[1]=1
    inode[2]=Np
    inode[3]=Np
    inode[4]=1
    jnode[1]=1
    jnode[2]=1
    jnode[3]=Np
    jnode[4]=Np

    #count how many elements own each node
    for in0=1:4
        for e=1:Ne
            I=intma[inode[in0],jnode[in0],e]
            lhowm[I]=lhowm[I] + 1
        end #ie
    end #in0

    #track elements owning each node
    lwher[1]=0
    for I=2:Npoin
        lwher[I]=lwher[I-1] + lhowm[I-1]
    end #ip

    #another tracker array
    lhowm = zeros(Int64,Npoin,1)
    for in0=1:4
        for e=1:Ne
            I=intma[inode[in0],jnode[in0],e]
            lhowm[I]=lhowm[I] + 1
            jloca=lwher[I] + lhowm[I]
            icone[jloca]=e
        end #ie
    end #in

    #LOOP OVER THE NODES
    iloca=0
    for I=1:Npoin
        iloc1=iloca
        iele=lhowm[I]

        if (iele != 0 )
            iwher=lwher[I]

            #LOOP OVER THOSE ELEMENTS SURROUNDING NODE I
            ip1=I
            for el=1:iele
                e=icone[iwher+el]

                #find out position of ip in intma
                i0=0
                for in0=1:4
                    i0=in0
                    ipt=intma[inode[in0],jnode[in0],e]
                    if (ipt == I)
                        break
                    end
                end #in0
                in0=i0

                #Check Edge of Element E which claims I
                j=0
                for jnod = 1:2:3
                    iold=0
                    j=j+1
                    in2=in0 + jnod
                    if (in2 > 4)
                        in2=in2-4
                    end
                    ip2=intma[inode[in2],jnode[in2],e]
                    if (ip2 >= ip1)

                        #check whether side is old | new
                        jloca=0
                        if (iloca != iloc1)
                            for s=iloc1+1:iloca
                                jloca=s
                                if (iside[2,s] == ip2)
                                    iold=1
                                    break
                                end
                            end #s
                        end #iloca

                        if (iold == 0)
                            #NEW SIDE
                            iloca=iloca + 1
                            iside[1,iloca]=ip1
                            iside[2,iloca]=ip2
                            iside[2+j,iloca]=e
                        elseif (iold == 1)
                            #OLD SIDE
                            iside[2+j,jloca]=e
                        end #iold
                    end #ip2
                end #jnod
            end #iel

            #Perform some Shifting to order the nodes of a side in CCW direction
            for s=iloc1+1:iloca
                if (iside[3,s] == 0)
                    iside[3,s]=iside[4,s]
                    iside[4,s]=0
                    iside[1,s]=iside[2,s]
                    iside[2,s]=ip1
                end #iside
            end #s
        end #if iele
    end #ip

    if (iloca != Nface)
        println( "Error in CREATE_SIDE. iloca = ",iloca," Nface = ",Nface)
        exit
    end

    #RESET THE BOUNDARY MARKERS
    for s=1:Nface
        if (iside[4,s] == 0)
        il=iside[1,s]
        ir=iside[2,s]
        e=iside[3,s]
        for ib=1:Nboun
            ibe=bsido[3,ib]
            ibc=bsido[4,ib]
            if (ibe == e)
                ilb=bsido[1,ib]
                irb=bsido[2,ib]
                if  (ilb == il && irb == ir)
                    iside[4,s]=-ibc
                    break
                end #ilb
            end #ibe
        end #ib
        end #iside
    end #s

    #FORM ELEMENT/SIDE CONNECTIVITY ARRAY
    for s=1:Nface
        el=iside[3,s]
        er=iside[4,s]
        is1=iside[1,s]
        is2=iside[2,s]

        #LEFT SIDE
        for in0=1:4
            i1=intma[inode[in0],jnode[in0],el]
            in1=in0 + 1
            if (in1 > 4)
                in1=1
            end #in1
            i2=intma[inode[in1],jnode[in1],el]
            if ((is1 == i1) && (is2 == i2))
                jeside[in0,el]=s
            end #is1
        end #in0

        #RIGHT SIDE
        if (er > 0)
            for in0=1:4
                i1=intma[inode[in0],jnode[in0],er]
                in1=in0 + 1
                if (in1 > 4)
                    in1=1
                end #in1
                i2=intma[inode[in1],jnode[in1],er]
                if ((is1 == i2) && (is2 == i1))
                    jeside[in0,er]=s
                end #is1
            end #in
        end #ier
    end #s

    return(iside,jeside)
end


text:
module TestMessages

using ZulipReminderBot
using ZulipReminderBot: Adapter, ZulipRequest
using ZulipReminderBot: process_list
import ZulipReminderBot
using Test
using Dates
using Base64

import DBInterface

struct MockDB end
struct MockDBStmt
    s::String
end

function ZulipReminderBot.prepare(db::Adapter{MockDB, MockDB}, query)
    MockDBStmt(query)
end

function DBInterface.execute(stmt::MockDBStmt, x)
    if stmt.s == "SELECT * FROM senders WHERE id = ?" && x == (1, )
        return [(id = 1, v1 = "Europe/London")]
    elseif stmt.s == "SELECT * FROM senders WHERE id = ?" && x == (2, )
        return [(id = 1, v1 = "UTC+4")]
    elseif stmt.s == "SELECT * FROM messages WHERE msg_sender_id = ?" && x == (1, )
        return [(id = 123, cts = 0, ets = 1622571921000, s1 = "Stream1", t1 = "Topic1", tp = "private", sid = 1, cn = base64encode("Hello"))]
    elseif stmt.s == "SELECT * FROM messages WHERE msg_sender_id = ?" && x == (2, )
        return [(id = 123, cts = 0, ets = 1622571921000, s1 = "Stream1", t1 = "Topic1", tp = "private", sid = 1, cn = base64encode("Hello")),
               (id = 124, cts = 0, ets = 1622571931000, s1 = "Stream1", t1 = "Topic1", tp = "private", sid = 1, cn = base64encode("Hello2"))]
    elseif stmt.s == "SELECT * FROM messages WHERE msg_sender_id = ?" && x == (3, )
        return []
    end
end

function DBInterface.close!(stmt::MockDBStmt) end

@testset "test list command" begin
    zreq = ZulipRequest()
    zreq.message.sender_id = 1
    db = Adapter(MockDB(), MockDB)

    res = process_list(zreq, db, nothing)
    @test res == "**id:** 123\n**scheduled:** 2021-06-01 19:25:21 +01:00\nHello"

    zreq.message.sender_id = 2
    res = process_list(zreq, db, nothing)
    @test res == "**id:** 123\n**scheduled:** 2021-06-01 22:25:21 +04:00\nHello\n---\n**id:** 124\n**scheduled:** 2021-06-01 22:25:31 +04:00\nHello2"
    zreq.message.sender_id = 3
    res = process_list(zreq, db, nothing)
    @test res == "No messages scheduled"
end

end # module


text:
# Generate 3D data for the overthrust model
# Author: Philipp Witte, pwitte@eoas.ubc.ca
# Date: May 2017
#

Pkg.checkout("opesciSLIM","development")
using PyCall, PyPlot, JLD, HDF5, opesciSLIM.TimeModeling

# Load velocity model
n,d,o,m = h5open("../data/overthrust_3D_true_model.h5","r") do file
	read(file, "n", "d", "o", "m")
end

# Set up model structure
model = Model((n[1],n[2],n[3]), (d[1],d[2],d[3]), (o[1],o[2],o[3]), m)	# need n,d,o as tuples and m0 as array

# Bound constraints
vmin = ones(Float32,model.n) + 0.3
vmax = ones(Float32,model.n) + 5.5

# Slowness squared [s^2/km^2]
mmin = vec((1./vmax).^2.)
mmax = vec((1./vmin).^2.)

# Set up source grid
nsrc = 97^2
xsrc = 400.f0:200.f0:19600.f0	# 500 m space from boundary
ysrc = 400.f0:200.f0:19600.f0
zsrc = 50.f0
(xsrc,ysrc,zsrc) = setup_3D_grid(xsrc,ysrc,zsrc)
xsrc = convertToCell(xsrc)
ysrc = convertToCell(ysrc)
zsrc = convertToCell(zsrc)
timeS = 3000.f0
dtS = 4.f0
srcGeometry = Geometry(xsrc,ysrc,zsrc;dt=dtS,t=timeS)

## Set up receivers
xmin = 100.f0; xmax = 19900.f0
ymin = 100.f0; ymax = 19900.f0
offset = 6000.f0
nxrec = 240
nyrec = 240
xrec = Array{Any,1}(nsrc)
yrec = Array{Any,1}(nsrc)
zrec = Array{Any}(nsrc)
for j=1:nsrc
	xsrc[j] - offset < xmin ? xlocal1 = xmin : xlocal1 = xsrc[j] - offset
	xsrc[j] + offset > xmax ? xlocal2 = xmax : xlocal2 = xsrc[j] + offset
	ysrc[j] - offset < ymin ? ylocal1 = ymin : ylocal1 = ysrc[j] - offset
	ysrc[j] + offset > ymax ? ylocal2 = ymax : ylocal2 = ysrc[j] + offset
	xlocal = xlocal1:50.f0:xlocal2
	ylocal = ylocal1:50.f0:ylocal2
	zlocal = 500.f0
	(xlocal,ylocal,zlocal) = setup_3D_grid(xlocal,ylocal,zlocal)
	xrec[j] = xlocal
	yrec[j] = ylocal
	zrec[j] = zlocal
end
timeR = 3000.f0	# receiver recording time [ms]
dtR = 4.f0	# receiver sampling interval
recGeometry = Geometry(xrec,yrec,zrec;dt=dtR,t=timeR)

# Set up source
f0 = 0.008
wavelet = source(timeS,dtS,f0)

# Info for modeling operators
ntComp = get_computational_nt(srcGeometry,recGeometry,model)
info = Info(prod(model.n),nsrc,ntComp)

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

# Set up modeling options
opt = Options(limit_m = true,
			  save_data_to_disk = true,
			  file_path = "../data",
 			  file_name = "overthrust_3D_shot_",
			  save_rate = 14.
			  )

# Set up operators
Pr = joProjection(info,recGeometry)
Ps = joProjection(info,srcGeometry)
F = joModeling(info,model; options=opt)
q = joData(srcGeometry,wavelet)

# Generate data
dobs = Pr*F*Ps'*q

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
struct Foo
    x
    y
end

mutable struct Bar
    x::GPlot.Option{Int}
    y::GPlot.Option{Int}
    z::Bool
end

@testset "▶ /utils                      " begin
    g = GPlot.GLE()
    gp = GPlot.Gnuplot()
    @test typeof(g) <: GPlot.Backend
    @test typeof(gp) <: GPlot.Backend
    @test g.io isa IOBuffer
    @test gp.io isa IOBuffer

    "blah" |> g
    @test String(take!(g)) == "blah " # NOTE extra space
    "blih" |> (g, gp)
    @test String(take!(g)) == String(take!(gp)) == "blih "

    @test isnothing(nothing)
    @test GPlot.isdef(1)
    @test !GPlot.isdef(nothing)

    @test GPlot.isanydef(Foo(1, 2))
    @test GPlot.isanydef(Foo(1, nothing))
    @test !GPlot.isanydef(Foo(nothing, nothing))

    @test GPlot.round3d(pi) == 3.142
    @test GPlot.col2str(colorant"red") == "rgba(1.0,0.0,0.0,1.0)"

    @test GPlot.vec2str([1,2,3]) == "1 2 3 "
    @test GPlot.vec2str(["a","b","c"]) == "\"a\" \"b\" \"c\" "
    @test GPlot.svec2str(["a","b","c"]) == "a,b,c"
    @test GPlot.svec2str(["a"]) == "a"
    @test GPlot.svec2str(("a$i" for i ∈ 1:3)) == "a1,a2,a3"

    @test t"\sin(x)+\cos(x)" == "\\sin(x)+\\cos(x)"
    _α = 5
    @test t"\sin(##_α x)" == "\\sin(5 x)"

    # FL
    @test G.fl(nothing) === nothing
    @test G.fl(missing) === missing
    @test G.fl(2) == 2.0

    # CSV_writer
    f, _ = mktemp()
    z = zip([missing, 1, Inf, 3, missing, NaN])
    G.csv_writer(f, z, true)
    ff = read(f, String)
    ff == "?\n1\n?\n3\n?\n?\n"
end

@testset "▶ /utils2                     " begin
    f = Figure("blah", reset=true)
    @test gcf().id == f.id
    @test gcf().size == f.size
    @test gca() === nothing

    m = G.MarkerStyle("circle", 0.5, colorant"blue")
    @test G.str(m) == "mk_circle_rgba_0_0_0_0_1_0_1_0_"

    set_palette([colorant"blue", colorant"red"])
    @test G.GP_ENV["PALETTE"] == [colorant"blue", colorant"red"]
    @test G.GP_ENV["SIZE_PALETTE"] == 2
end


text:
const HDF_FILENAME = "simulation_store.h5"
const HDF_SIMULATION_ROOT_PATH = "simulation"
const EMULATION_MODEL_PATH = "$HDF_SIMULATION_ROOT_PATH/emulation_model"
const OPTIMIZER_STATS_PATH = "optimizer_stats"
const SERIALIZED_KEYS_PATH = "serialized_keys"

# This only applies if chunks are enabled, and that will only likely happen if we enable
# compression.
# The optimal number of chunks to store in memory will vary widely.
# The HDF docs recommend keeping chunk byte sizes between 10 KiB - 1 MiB.
# We want to make it big enough to compress duplicate values.
# The downside to making this larger is that any read causes the
# entire chunk to be read.
# If one variable has 10,000 components and each value is a Float64 then one row would
# consume 10,000 * 8 = 78 KiB
DEFAULT_MAX_CHUNK_BYTES = 128 * KiB

"""
Stores simulation data in an HDF file.
"""
mutable struct HdfSimulationStore <: SimulationStore
    file::HDF5.File
    params::SimulationStoreParams
    # The key order is the problem execution order.
    dm_data::OrderedDict{Symbol, DatasetContainer{HDF5Dataset}}
    em_data::DatasetContainer{HDF5Dataset}
    # The key is the problem name.
    optimizer_stats_datasets::Dict{Symbol, HDF5.Dataset}
    optimizer_stats_write_index::Dict{Symbol, Int}
    cache::OptimizationOutputCaches
end

get_initial_time(store::HdfSimulationStore) = get_initial_time(store.params)

function HdfSimulationStore(file_path::AbstractString, mode::AbstractString)
    if !(mode == "w" || mode == "r")
        throw(IS.ConflictingInputsError("mode can only be 'w' or 'r'"))
    end

    if isfile(file_path) && mode == "w"
        throw(IS.ConflictingInputsError("$file_path already exists"))
    end

    file = HDF5.h5open(file_path, mode)
    if mode == "w"
        HDF5.create_group(file, HDF_SIMULATION_ROOT_PATH)
        @debug "Created store" file_path
    end

    store = HdfSimulationStore(
        file,
        SimulationStoreParams(),
        OrderedDict{Symbol, DatasetContainer{HDF5Dataset}}(),
        DatasetContainer{HDF5Dataset}(),
        Dict{Symbol, HDF5.Dataset}(),
        Dict{Symbol, Int}(),
        OptimizationOutputCaches(),
    )
    mode == "r" && _deserialize_attributes!(store)

    finalizer(_check_state, store)
    return store
end

"""
Construct and open an HdfSimulationStore.

When reading or writing results in a program you should use the method that accepts a
function in order to guarantee that the file handle gets closed.

# Arguments

  - `directory::AbstractString`: Directory containing the store file
  - `mode::AbstractString`: Mode to use to open the store file
  - `filename::AbstractString`: Base name of the store file

# Examples

```julia
# Assumes a simulation has been executed in the './rts' directory with these parameters.
path = "./rts"
problem = :ED
var_name = :P__ThermalStandard
timestamp = DateTime("2020-01-01T05:00:00")
store = open_store(HdfSimulationStore, path)
df = PowerSimulations.read_result(DataFrame, store, model, :variables, var_name, timestamp)
```
"""
function open_store(
    ::Type{HdfSimulationStore},
    directory::AbstractString,
    mode="r";
    filename=HDF_FILENAME,
)
    return HdfSimulationStore(joinpath(directory, filename), mode)
end

function open_store(
    func::Function,
    ::Type{HdfSimulationStore},
    directory::AbstractString,
    mode="r";
    filename=HDF_FILENAME,
)
    store = nothing
    try
        store = HdfSimulationStore(joinpath(directory, filename), mode)
        return func(store)
    finally
        if store !== nothing
            close(store)
        end
    end
end

function Base.close(store::HdfSimulationStore)
    flush(store)
    HDF5.close(store.file)
    empty!(store.cache)
    @debug "Close store file handle" store.file
end

function Base.isopen(store::HdfSimulationStore)
    return store.file === nothing ? false : HDF5.isopen(store.file)
end

function Base.flush(store::HdfSimulationStore)
    for (key, output_cache) in store.cache.data
        _flush_data!(output_cache, store, key, false)
        @assert !has_dirty(output_cache) "$key has dirty cache after flushing"
    end

    flush(store.file)
    @debug "Flush store"
    return
end

get_params(store::HdfSimulationStore) = store.params

function get_decision_model_params(store::HdfSimulationStore, model_name::Symbol)
    return get_decision_model_params(get_params(store), model_name)
end

function get_emulation_model_params(store::HdfSimulationStore)
    return get_emulation_model_params(get_params(store))
end

function get_container_key_lookup(store::HdfSimulationStore)
    function _get_lookup()
        root = _get_root(store)
        buf = IOBuffer(root[SERIALIZED_KEYS_PATH][:])
        return Serialization.deserialize(buf)
    end
    isopen(store) && return _get_lookup()

    store.file = HDF5.h5open(store.file.filename, "r")
    try
        return _get_lookup()
    finally
        HDF5.close(store.file)
    end
end

"""
Return the problem names in order of execution.
"""
list_decision_models(store::HdfSimulationStore) = keys(get_dm_data(store))

"""
Return the fields stored for the `problem` and `container_type` (duals/parameters/variables).
"""
function list_decision_model_keys(
    store::HdfSimulationStore,
    model::Symbol,
    container_type::Symbol,
)
    container = getfield(get_dm_data(store)[model], container_type)
    return keys(container)
end

function list_emulation_model_keys(store::HdfSimulationStore, container_type::Symbol)
    container = getfield(get_em_data(store), container_type)
    return keys(container)
end

function write_optimizer_stats!(
    store::HdfSimulationStore,
    model::OperationModel,
    ::DecisionModelIndexType,
)
    stats = get_optimizer_stats(model)
    model_name = get_name(model)
    dataset = _get_dataset(OptimizerStats, store, model_name)

    # Uncomment for performance measures of HDF Store
    #TimerOutputs.@timeit RUN_SIMULATION_TIMER "Write optimizer stats" begin
    dataset[:, store.optimizer_stats_write_index[model_name]] = to_matrix(stats)
    #end

    store.optimizer_stats_write_index[model_name] += 1
    return
end

function write_optimizer_stats!(
    store::HdfSimulationStore,
    model::OperationModel,
    ::EmulationModelIndexType,
)
    return
end

"""
Read the optimizer stats for a problem execution.
"""
function read_optimizer_stats(
    store::HdfSimulationStore,
    simulation_step::Int,
    model_name::Symbol,
    execution_index::Int,
)
    optimizer_stats_write_index =
        (simulation_step - 1) *
        store.params.decision_models_params[model_name].num_executions + execution_index
    dataset = _get_dataset(OptimizerStats, store, model_name)
    return OptimizerStats(dataset[:, optimizer_stats_write_index])
end

"""
Return the optimizer stats for a problem as a DataFrame.
"""
function read_optimizer_stats(store::HdfSimulationStore, model_name)
    dataset = _get_dataset(OptimizerStats, store, model_name)
    data = permutedims(dataset[:, :])
    stats = [to_namedtuple(OptimizerStats(data[i, :])) for i in axes(data)[1]]
    return DataFrames.DataFrame(stats)
end

function initialize_problem_storage!(
    store::HdfSimulationStore,
    params::SimulationStoreParams,
    dm_problem_reqs::Dict{Symbol, SimulationModelStoreRequirements},
    em_problem_reqs::SimulationModelStoreRequirements,
    flush_rules::CacheFlushRules,
)
    store.params = params
    root = store.file[HDF_SIMULATION_ROOT_PATH]
    problems_group = _get_group_or_create(root, "decision_models")
    store.cache = OptimizationOutputCaches(flush_rules)
    @info "Initialize store cache" get_min_flush_size(store.cache) get_max_size(store.cache)
    initial_time = get_initial_time(store)
    container_key_lookup = Dict{String, OptimizationContainerKey}()
    for (problem, problem_params) in store.params.decision_models_params
        get_dm_data(store)[problem] = DatasetContainer{HDF5Dataset}()
        problem_group = _get_group_or_create(problems_group, string(problem))
        for type in STORE_CONTAINERS
            group = _get_group_or_create(problem_group, string(type))
            for (key, reqs) in getfield(dm_problem_reqs[problem], type)
                !should_write_resulting_value(key) && continue
                name = encode_key_as_string(key)
                dataset = _create_dataset(group, name, reqs)
                # Columns can't be stored in attributes because they might be larger than
                # the max size of 64 KiB.
                col = _make_column_name(name)
                HDF5.write_dataset(group, col, string.(reqs["columns"]))
                column_dataset = group[col]
                datasets = getfield(get_dm_data(store)[problem], type)
                datasets[key] = HDF5Dataset(
                    dataset,
                    column_dataset,
                    get_resolution(problem_params),
                    initial_time,
                )
                add_output_cache!(
                    store.cache,
                    problem,
                    key,
                    get_rule(flush_rules, problem, key),
                )
                container_key_lookup[encode_key_as_string(key)] = key
            end
        end

        num_stats = params.num_steps * params.decision_models_params[problem].num_executions
        columns = fieldnames(OptimizerStats)
        num_columns = length(columns)
        dataset = HDF5.create_dataset(
            problem_group,
            OPTIMIZER_STATS_PATH,
            HDF5.datatype(Float64),
            HDF5.dataspace((num_columns, num_stats)),
        )
        HDF5.attributes(dataset)["columns"] = [string(x) for x in columns]
        store.optimizer_stats_datasets[problem] = dataset
        store.optimizer_stats_write_index[problem] = 1
        @debug "Initialized optimizer_stats_datasets $problem ($num_columns, $num_stats)"
    end

    emulation_group = _get_group_or_create(root, "emulation_model")
    for emulation_params in values(store.params.emulation_model_params)
        for type in STORE_CONTAINERS
            group = _get_group_or_create(emulation_group, string(type))
            for (key, reqs) in getfield(em_problem_reqs, type)
                name = encode_key_as_string(key)
                dataset = _create_dataset(group, name, reqs)
                # Columns can't be stored in attributes because they might be larger than
                # the max size of 64 KiB.
                col = _make_column_name(name)
                HDF5.write_dataset(group, col, string.(reqs["columns"]))
                column_dataset = group[col]
                datasets = getfield(store.em_data, type)
                datasets[key] = HDF5Dataset(
                    dataset,
                    column_dataset,
                    get_resolution(emulation_params),
                    initial_time,
                )
                container_key_lookup[encode_key_as_string(key)] = key
            end
        end
    end
    buf = IOBuffer()
    Serialization.serialize(buf, container_key_lookup)
    seek(buf, 0)
    root[SERIALIZED_KEYS_PATH] = buf.data

    # This has to run after problem groups are created.
    _serialize_attributes(store)
    return
end

log_cache_hit_percentages(x::HdfSimulationStore) = log_cache_hit_percentages(x.cache)

"""
Return DataFrame, DenseAxisArray, or Array for a model result at a timestamp.
"""
function read_result(
    ::Type{DataFrames.DataFrame},
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
)
    data, columns = _read_data_columns(store, model_name, key, index)

    if (ndims(data) < 2 || size(data)[1] == 1) && size(data)[2] != size(columns)[1]
        data = reshape(data, length(data), 1)
    end
    return DataFrames.DataFrame(data, columns)
end

function read_result(
    ::Type{DenseAxisArray},
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
)
    data, columns = _read_data_columns(store, model_name, key, index)
    return DenseAxisArray(permutedims(data), columns, 1:size(data)[1])
end

function read_result(
    ::Type{<:Array},
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::Union{DecisionModelIndexType, EmulationModelIndexType},
)
    if is_cached(store.cache, model_name, key, index)
        data = _read_result(store.cache, model_name, key, index)
    else
        # PERF: If this will be commonly used then we need to remove reading of columns.
        data, _ = _read_result(store, model_name, key, index)
    end
end

function read_results(
    store::HdfSimulationStore,
    key::OptimizationContainerKey;
    index::Union{Nothing, EmulationModelIndexType}=nothing,
    len::Union{Nothing, Int}=nothing,
)
    dataset = _get_em_dataset(store, key)
    @assert_op ndims(dataset.values) == 2
    if isnothing(index)
        @assert_op(isnothing(len))
        data = dataset.values[:, :]
    elseif isnothing(len)
        data = dataset.values[index:end, :]
    else
        data = dataset.values[index:(index + len - 1), :]
    end
    columns = get_column_names(key, dataset)
    @assert_op size(data)[2] == length(columns)
    return DataFrames.DataFrame(data, columns)
end

function get_emulation_model_dataset_size(
    store::HdfSimulationStore,
    key::OptimizationContainerKey,
)
    dataset = _get_em_dataset(store, key)
    return size(dataset.values)
end

function _read_result(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
)
    !isopen(store) && throw(ArgumentError("store must be opened prior to reading"))
    model_params = get_emulation_model_params(store)

    if index > model_params.num_executions
        throw(
            ArgumentError(
                "index = $index cannot be larger than $(model_params.num_executions)",
            ),
        )
    end

    dataset = _get_em_dataset(store, key)
    dset = dataset.values
    # Uncomment for performance checking
    #TimerOutputs.@timeit RUN_SIMULATION_TIMER "Read dataset" begin
    @assert_op ndims(dset) == 2
    data = dset[index, :]
    #end
    columns = get_column_names(key, dataset)
    data = permutedims(data)
    @assert_op size(data)[2] == length(columns)
    @assert_op size(data)[1] == 1
    return data, columns
end

function _read_result(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
)
    simulation_step, execution_index = _get_indices(store, model_name, index)
    return _read_result(store, model_name, key, simulation_step, execution_index)
end

function _read_result(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    simulation_step::Int,
    execution_index::Int,
)
    !isopen(store) && throw(ArgumentError("store must be opened prior to reading"))

    model_params = get_decision_model_params(store, model_name)
    num_executions = model_params.num_executions
    if execution_index > num_executions
        throw(
            ArgumentError(
                "execution_index = $execution_index cannot be larger than $num_executions",
            ),
        )
    end

    dataset = _get_dm_dataset(store, model_name, key)
    dset = dataset.values
    row_index = (simulation_step - 1) * num_executions + execution_index
    columns = get_column_names(key, dataset)

    # Uncomment for performance checking
    #TimerOutputs.@timeit RUN_SIMULATION_TIMER "Read dataset" begin
    num_dims = ndims(dset)
    if num_dims == 3
        data = dset[:, :, row_index]
        #elseif num_dims == 4
        #    data = dset[:, :, :, row_index]
    else
        error("unsupported dims: $num_dims")
    end
    #end

    return data, columns
end

"""
Write a decision model result for a timestamp to the store.
"""
function write_result!(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
    ::Dates.DateTime,
    data,
)
    output_cache = get_output_cache(store.cache, model_name, key)

    cur_size = get_size(store.cache)
    add_result!(output_cache, index, to_matrix(data), is_full(store.cache, cur_size))

    if get_dirty_size(output_cache) >= get_min_flush_size(store.cache)
        discard = !should_keep_in_cache(output_cache)

        # PERF: A potentially significant performance improvement would be to queue several
        # flushes and submit them in parallel.
        size_flushed = _flush_data!(output_cache, store, model_name, key, discard)

        @debug "flushed data" LOG_GROUP_SIMULATION_STORE key size_flushed discard cur_size
    end

    # Disabled because this is currently a noop.
    #if is_full(store.cache)
    #    _flush_data!(store.cache, store)
    #end

    @debug "write_result" get_size(store.cache) encode_key_as_string(key)
    return
end

"""
Write an emulation model result for an execution index value and the timestamp of the update
"""
function write_result!(
    store::HdfSimulationStore,
    ::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    simulation_time::Dates.DateTime,
    data::Matrix{Float64},
)
    dataset = _get_em_dataset(store, key)
    _write_dataset!(dataset.values, data, index:index)
    set_last_recorded_row!(dataset, index)
    set_update_timestamp!(dataset, simulation_time)
    return
end

function write_result!(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    simulation_time::Dates.DateTime,
    data,
)
    write_result!(store, model_name, key, index, simulation_time, to_matrix(data))
    return
end

function _check_state(store::HdfSimulationStore)
    if has_dirty(store.cache)
        error("BUG!!! dirty cache is present at shutdown: $(store.file)")
    end
end

function _compute_chunk_count(dims, dtype; max_chunk_bytes=DEFAULT_MAX_CHUNK_BYTES)
    bytes_per_element = sizeof(dtype)

    if length(dims) == 2
        size_row = dims[1] * bytes_per_element
    elseif length(dims) == 3
        size_row = dims[1] * dims[2] * bytes_per_element
    elseif length(dims) == 4
        size_row = dims[1] * dims[2] * dims[3] * bytes_per_element
    else
        error("unsupported dims = $dims")
    end

    chunk_count = minimum((trunc(max_chunk_bytes / size_row), dims[end]))
    if chunk_count == 0
        error(
            "HDF Max Chunk Bytes is smaller than the size of a row. Please increase it. " *
            "max_chunk_bytes=$max_chunk_bytes dims=$dims " *
            "size_row=$size_row",
        )
    end

    chunk_dims = [x for x in dims]
    chunk_dims[end] = chunk_count
    return chunk_dims
end

function _create_dataset(group, name, reqs)
    dataset = HDF5.create_dataset(
        group,
        name,
        HDF5.datatype(Float64),
        HDF5.dataspace(reqs["dims"]),
        # We are choosing to optimize read performance in the first implementation.
        # Compression would slow that down.
        #chunk = _compute_chunk_count(reqs["dims"], Float64),
        #shuffle = (),
        #deflate = 3,
    )
    @debug "Created dataset for" group name size(dataset)
    return dataset
end

function _deserialize_attributes!(store::HdfSimulationStore)
    container_key_lookup = get_container_key_lookup(store)
    group = store.file["simulation"]
    initial_time = Dates.DateTime(HDF5.read(HDF5.attributes(group)["initial_time"]))
    step_resolution =
        Dates.Millisecond(HDF5.read(HDF5.attributes(group)["step_resolution_ms"]))
    num_steps = HDF5.read(HDF5.attributes(group)["num_steps"])
    store.params = SimulationStoreParams(initial_time, step_resolution, num_steps)
    empty!(get_dm_data(store))
    for model in HDF5.read(HDF5.attributes(group)["problem_order"])
        problem_group = store.file["simulation/decision_models/$model"]
        model_name = Symbol(model)
        store.params.decision_models_params[model_name] = ModelStoreParams(
            HDF5.read(HDF5.attributes(problem_group)["num_executions"]),
            HDF5.read(HDF5.attributes(problem_group)["horizon"]),
            Dates.Millisecond(HDF5.read(HDF5.attributes(problem_group)["interval_ms"])),
            Dates.Millisecond(HDF5.read(HDF5.attributes(problem_group)["resolution_ms"])),
            HDF5.read(HDF5.attributes(problem_group)["base_power"]),
            Base.UUID(HDF5.read(HDF5.attributes(problem_group)["system_uuid"])),
        )
        get_dm_data(store)[model_name] = DatasetContainer{HDF5Dataset}()
        for type in STORE_CONTAINERS
            group = problem_group[string(type)]
            for name in keys(group)
                if !endswith(name, "columns")
                    dataset = group[name]
                    column_dataset = group[_make_column_name(name)]
                    resolution =
                        get_resolution(get_decision_model_params(store, model_name))
                    item = HDF5Dataset(dataset, column_dataset, resolution, initial_time)
                    container_key = container_key_lookup[name]
                    getfield(get_dm_data(store)[model_name], type)[container_key] = item
                    add_output_cache!(
                        store.cache,
                        model_name,
                        container_key,
                        CacheFlushRule(),
                    )
                end
            end
        end

        store.optimizer_stats_datasets[model_name] = problem_group[OPTIMIZER_STATS_PATH]
        store.optimizer_stats_write_index[model_name] = 0
    end

    em_group = _get_emulation_model_path(store)
    model_name = Symbol(HDF5.read(HDF5.attributes(em_group)["name"]))
    resolution = Dates.Millisecond(HDF5.read(HDF5.attributes(em_group)["resolution_ms"]))
    store.params.emulation_model_params[model_name] = ModelStoreParams(
        HDF5.read(HDF5.attributes(em_group)["num_executions"]),
        HDF5.read(HDF5.attributes(em_group)["horizon"]),
        Dates.Millisecond(HDF5.read(HDF5.attributes(em_group)["interval_ms"])),
        resolution,
        HDF5.read(HDF5.attributes(em_group)["base_power"]),
        Base.UUID(HDF5.read(HDF5.attributes(em_group)["system_uuid"])),
    )
    for type in STORE_CONTAINERS
        group = em_group[string(type)]
        for name in keys(group)
            if !endswith(name, "columns")
                dataset = group[name]
                column_dataset = group[_make_column_name(name)]
                item = HDF5Dataset(dataset, column_dataset, resolution, initial_time)
                container_key = container_key_lookup[name]
                getfield(store.em_data, type)[container_key] = item
                add_output_cache!(store.cache, model_name, container_key, CacheFlushRule())
            end
        end
    end
    # TODO: optimizer stats are not being written for EM.

    @debug "deserialized store params and datasets" store.params
end

function _serialize_attributes(store::HdfSimulationStore)
    params = store.params
    group = store.file["simulation"]
    HDF5.attributes(group)["problem_order"] =
        [string(k) for k in keys(params.decision_models_params)]
    HDF5.attributes(group)["initial_time"] = string(params.initial_time)
    HDF5.attributes(group)["step_resolution_ms"] =
        Dates.Millisecond(params.step_resolution).value
    HDF5.attributes(group)["num_steps"] = params.num_steps

    for problem in keys(params.decision_models_params)
        problem_group = store.file["simulation/decision_models/$problem"]
        HDF5.attributes(problem_group)["num_executions"] =
            params.decision_models_params[problem].num_executions
        HDF5.attributes(problem_group)["horizon"] =
            params.decision_models_params[problem].horizon
        HDF5.attributes(problem_group)["resolution_ms"] =
            Dates.Millisecond(params.decision_models_params[problem].resolution).value
        HDF5.attributes(problem_group)["interval_ms"] =
            Dates.Millisecond(params.decision_models_params[problem].interval).value
        HDF5.attributes(problem_group)["base_power"] =
            params.decision_models_params[problem].base_power
        HDF5.attributes(problem_group)["system_uuid"] =
            string(params.decision_models_params[problem].system_uuid)
    end

    if !isempty(params.emulation_model_params)
        em_params = first(values(params.emulation_model_params))
        emulation_group = store.file["simulation/emulation_model"]
        HDF5.attributes(emulation_group)["name"] =
            string(first(keys(params.emulation_model_params)))
        HDF5.attributes(emulation_group)["num_executions"] = em_params.num_executions
        HDF5.attributes(emulation_group)["horizon"] = em_params.horizon
        HDF5.attributes(emulation_group)["resolution_ms"] =
            Dates.Millisecond(em_params.resolution).value
        HDF5.attributes(emulation_group)["interval_ms"] =
            Dates.Millisecond(em_params.interval).value
        HDF5.attributes(emulation_group)["base_power"] = em_params.base_power
        HDF5.attributes(emulation_group)["system_uuid"] = string(em_params.system_uuid)
    end
    return
end

function _flush_data!(
    cache::OptimizationOutputCache,
    store::HdfSimulationStore,
    model_name,
    key::OptimizationContainerKey,
    discard,
)
    return _flush_data!(cache, store, OptimizationResultCacheKey(model_name, key), discard)
end

function _flush_data!(
    cache::OptimizationOutputCache,
    store::HdfSimulationStore,
    cache_key::OptimizationResultCacheKey,
    discard::Bool,
)
    !has_dirty(cache) && return 0
    dataset = _get_dm_dataset(store, cache_key)
    timestamps, data = get_dirty_data_to_flush!(cache)
    num_results = length(timestamps)
    @assert_op num_results == size(data)[end]
    end_index = dataset.write_index + length(timestamps) - 1
    write_range = (dataset.write_index):end_index
    # Enable only for development and benchmarking
    #TimerOutputs.@timeit RUN_SIMULATION_TIMER "Write $(key.key) array to HDF" begin
    _write_dataset!(dataset.values, data, write_range)
    #end

    discard && discard_results!(cache, timestamps)

    dataset.write_index += num_results
    size_flushed = cache.size_per_entry * num_results

    @debug "Flushed cache results to HDF5" LOG_GROUP_SIMULATION_STORE cache_key size_flushed num_results get_size(
        store.cache,
    )
    return size_flushed
end

function _get_dataset(::Type{OptimizerStats}, store::HdfSimulationStore, model_name)
    return store.optimizer_stats_datasets[model_name]
end

function _get_dataset(::Type{OptimizerStats}, store::HdfSimulationStore)
    return store.optimizer_stats_datasets
end

function _get_em_dataset(store::HdfSimulationStore, key::OptimizationContainerKey)
    return getfield(get_em_data(store), get_store_container_type(key))[key]
end

function _get_dm_dataset(store::HdfSimulationStore, model_name::Symbol)
    return get_dm_data(store)[model_name]
end

function _get_dm_dataset(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
)
    return getfield(get_dm_data(store)[model_name], get_store_container_type(key))[key]
end

function _get_dm_dataset(store::HdfSimulationStore, key::OptimizationResultCacheKey)
    return _get_dm_dataset(store, key.model, key.key)
end

function _get_group_or_create(parent, group_name)
    if haskey(parent, group_name)
        group = parent[group_name]
    else
        group = HDF5.create_group(parent, group_name)
        @debug "Created group" group
    end

    return group
end

_make_column_name(name) = string(name) * "__columns"

function _get_indices(store::HdfSimulationStore, model_name::Symbol, timestamp)
    time_diff = Dates.Millisecond(timestamp - store.params.initial_time)
    step = time_diff ÷ store.params.step_resolution + 1
    if step > store.params.num_steps
        throw(
            ArgumentError("timestamp = $timestamp is beyond the simulation: step = $step"),
        )
    end
    problem_params = store.params.decision_models_params[model_name]
    initial_time = store.params.initial_time + (step - 1) * store.params.step_resolution
    time_diff = timestamp - initial_time
    if time_diff % problem_params.interval != Dates.Millisecond(0)
        throw(ArgumentError("timestamp = $timestamp is not a valid problem timestamp"))
    end
    execution_index = time_diff ÷ problem_params.interval + 1
    return step, execution_index
end

_get_root(store::HdfSimulationStore) = store.file[HDF_SIMULATION_ROOT_PATH]
_get_emulation_model_path(store::HdfSimulationStore) = store.file[EMULATION_MODEL_PATH]

function _read_column_names(::Type{OptimizerStats}, store::HdfSimulationStore)
    dataset = _get_dataset(OptimizerStats, store)
    return HDF5.read(HDF5.attributes(dataset), "columns")
end

function _read_data_columns(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
)
    if is_cached(store.cache, model_name, key, index)
        data = read_result(store.cache, model_name, key, index)
        column_dataset = _get_dm_dataset(store, model_name, key).column_dataset
        columns = column_dataset[:]
    else
        data, columns = _read_result(store, model_name, key, index)
    end

    return data, columns
end

function _read_data_columns(
    store::HdfSimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
)
    # TODO: Enable once the cache is in use for em_data
    #if is_cached(store.cache, model_name, key, index)
    # data = read_result(store.cache, model_name, key, index)
    #columns = _get_em_dataset(store, model_name, key).column_dataset[:]
    #else
    #    data, columns = _read_result(store, model_name, key, index)
    #end

    return _read_result(store, model_name, key, index)
end

function _read_length(::Type{OptimizerStats}, store::HdfSimulationStore)
    dataset = _get_dataset(OptimizerStats, store)
    return HDF5.read(HDF5.attributes(dataset), "columns")
end

function _write_dataset!(
    dataset,
    array::Matrix{Float64},
    row_range::UnitRange{Int64},
    ::Val{3},
)
    dataset[:, 1, row_range] = array
    @debug "wrote dataset" dataset row_range
    return
end

function _write_dataset!(
    dataset,
    array::Matrix{Float64},
    row_range::UnitRange{Int64},
    ::Val{2},
)
    dataset[row_range, :] = array
    @debug "wrote dataset" dataset row_range
    return
end

function _write_dataset!(dataset, array::Matrix{Float64}, row_range::UnitRange{Int64})
    _write_dataset!(dataset, array, row_range, Val{ndims(dataset)}())
    return
end

function _write_dataset!(dataset, array::Array{Float64, 3}, row_range::UnitRange{Int64})
    dataset[:, :, row_range] = array
    @debug "wrote dataset" dataset row_range
    return
end