text:
# ------ Smoothing functions often used to permit continuous differentiability -----

"""
    abs_smooth(x, delta_x)

Smooth out the absolute value function with a quadratic interval.
delta_x is the half width of the smoothing interval.
Typically usage is with gradient-based optimization.
"""
function abs_smooth(x, delta_x)

    if real(x) < delta_x && real(x) > -delta_x
        return x^2/(2*delta_x) + delta_x/2
    else
        return abs_cs_safe(x)
    end

end

"""
    ksmax(x, hardness=50)

Kreisselmeier–Steinhauser constraint aggregation function.  In the limit as `hardness`
goes to infinity the maximum function is returned. Is mathematically guaranteed to
overestimate the maximum function, i.e. `maximum(x) <= ksmax(x, hardness)`.
"""
function ksmax(x, hardness=50)
    k = maximum(x)
    return 1.0/hardness*log(sum(exp.(hardness*(x.-k)))) .+ k
end

"""
    ksmin(x, hardness=50)

Kreisselmeier–Steinhauser constraint aggregation function.  In the limit as `hardness`
goes to infinity the minimum function is returned. Is mathematically guaranteed to
underestimate the minimum function, i.e. `minimum(x) <= ksmin(x, hardness)`.
"""
ksmin(x, hardness=50) = -ksmax(-x, hardness)

"""
    ksmax_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1)

Kreisselmeier–Steinhauser constraint aggregation function using the adaptive
hardness proposed by Poon and Martins in "An adaptive approach to constraint
aggregation using adjoint sensitivity analysis".  This implementation uses
Newton's method rather than the secant method for increasing hardness values.
Some blending is also used to ensure that result is C1 continuous.
`smoothing_fraction` controls the smoothness of this blending.
"""
function ksmax_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1)
    # check if derivative of the KS function wrt hardness is within the tolerance
    deriv = ksmax_h(x, hardness)
    target_deriv = -tol
    if abs(deriv) > tol * (1-smoothing_fraction)
        # increase hardness by applying Newton's method once in log-log space
        dderiv = ksmax_hh(x, hardness)*hardness/deriv
        new_hardness = 10^(log10(hardness) + log10(target_deriv/deriv)/dderiv)
        hardness = quintic_blend(hardness, new_hardness, -deriv, -target_deriv, smoothing_fraction*tol)
    end
    # use the new hardness to compute ksmax
    return ksmax(x, hardness)
end

"""
    ksmin_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1)

Kreisselmeier–Steinhauser constraint aggregation function using the adaptive
hardness proposed by Poon and Martins in "An adaptive approach to constraint
aggregation using adjoint sensitivity analysis".  This implementation uses
Newton's method rather than the secant method for increasing hardness values.
Some blending is also used to ensure that result is C1 continuous.
`smoothing_fraction` controls the smoothness of this blending.
"""
ksmin_adaptive(x, hardness=50; tol=1e-6, smoothing_fraction=0.1) =
    -ksmax_adaptive(-x, hardness, tol=tol, smoothing_fraction=smoothing_fraction)

"""
    ksmax_h(x, hardness)

Computes the derivative of the Kreisselmeier–Steinhauser constraint aggregation
function with respect to `hardness`.
"""
function ksmax_h(x, hardness)
    k = maximum(x)
    tmp1 = exp.(hardness*(x.-k))
    tmp2 = sum((x.-k).*tmp1)
    tmp3 = sum(tmp1)
    tmp4 = 1.0/hardness*log(tmp3)
    return 1.0/hardness*(tmp2/tmp3 - tmp4)
end

"""
    ksmax_hh(x, hardness)

Computes the second derivative of the Kreisselmeier–Steinhauser constraint aggregation
function with respect to `hardness`.
"""
function ksmax_hh(x, hardness)
    k = maximum(x)
    tmp1 = exp.(hardness*(x.-k))
    tmp2 = sum((x.-k).*tmp1)
    tmp2_h = sum((x.-k).^2 .* tmp1)
    tmp3 = sum(tmp1)
    tmp3_h = sum((x.-k).*tmp1)
    tmp4 = 1.0/hardness*log(tmp3)
    tmp4_h = 1.0/hardness*(tmp2/tmp3 - tmp4)
    return -1.0/hardness^2*(tmp2/tmp3 - tmp4) +
        1.0/hardness*(tmp2_h/tmp3 - tmp2*tmp3_h/tmp3^2 - tmp4_h)
end

"""
    sigmoid(x)

Sigmoid function, implemented with branching to avoid NaNs
"""
function sigmoid(x)
    if x >= zero(x)
        z = exp(-x)
        return one(z) / (one(z) + z)
    else
        z = exp(x)
        return z / (one(z) + z)
    end
end

"""
    sigmoid_blend(f1x, f2x, x, xt, hardness=50)

Smoothly transitions the results of functions f1 and f2 using the sigmoid function,
with the transition between the functions located at `xt`. `hardness` controls the
sharpness of the transition between the two functions.
"""
function sigmoid_blend(f1x, f2x, x, xt, hardness=50)
    sx = sigmoid(hardness*(x-xt))
    return f1x + sx*(f2x-f1x)
end

"""
    cubic_blend(f1x, f2x, x, xt, delta_x)

Smoothly transitions the results of functions f1 and f2 using a cubic polynomial,
with the transition between the functions located at `xt`. delta_x is the half
width of the smoothing interval.  The resulting function is C1 continuous.
"""
function cubic_blend(f1x, f2x, x, xt, delta_x)
    if x <= xt - delta_x
        return f1x
    elseif x >= xt + delta_x
        return f2x
    else
        xp = (x-xt)/(2*delta_x) + 1/2
        sx = -2*xp^3 + 3*xp^2
        return f1x + sx*(f2x-f1x)
    end
end

"""
    quintic_blend(f1x, f2x, x, xt, delta_x)

Smoothly transitions the results of functions f1 and f2 using a quintic polynomial,
with the transition between the functions located at `xt`. delta_x is the half
width of the smoothing interval.  The resulting function is C2 continuous.
"""
function quintic_blend(f1x, f2x, x, xt, delta_x)
    if x <= xt - delta_x
        return f1x
    elseif x >= xt + delta_x
        return f2x
    else
        xp = (x-xt)/(2*delta_x) + 1/2
        sx = 6*xp^5 - 15*xp^4 + 10*xp^3
        return f1x + sx*(f2x-f1x)
    end
end

# TODO AN: add smooth max/min with cubic splines


text:
abstract type AbstractContinuousMarkovChain <: AbstractMarkovChain end

"""
    ContinuousMarkovChain(transition_matrix)
    ContinuousMarkovChain(state_space, transition_matrix)
    ContinuousMarkovChain(discrete_markov_chain)

Creates a new continuous Markov chain object.
This is also known as a Markov jump process.

Note that an irreducible finite continuous-time Markov chain is
always positive recurrent and its stationary distribution always
exists, is unique and is equal to the limiting distribution.

# Arguments
- `state_space`: The names of the states that make up the Markov chain.
- `transition_matrix`: The transition intensity matrix. Also known as the generator matrix.
- `discrete_markov_chain`: An instance of `DiscreteMarkovChain`.

# Examples
The following shows a basic Sunny-Cloudy-Rainy weather model.
```jldoctest ContinuousMarkovChain
using DiscreteMarkovChains
T = [
    -.1 0.1 0.0;
    0.5 -.8 0.3;
    0.1 0.4 -.5;
]
X = ContinuousMarkovChain(["Sunny", "Cloudy", "Rainy"], T)
println(state_space(X))

# output

["Sunny", "Cloudy", "Rainy"]
```

```jldoctest ContinuousMarkovChain
println(transition_matrix(X))

# output

[-0.1 0.1 0.0; 0.5 -0.8 0.3; 0.1 0.4 -0.5]
```

# References
1. [Wikipedia](https://en.wikipedia.org/wiki/Continuous-time_Markov_chain)
"""
struct ContinuousMarkovChain <: AbstractContinuousMarkovChain
    state_space
    transition_matrix
    function ContinuousMarkovChain(state_space, transition_matrix)
        check(state_space, transition_matrix, ContinuousMarkovChain)
        new(state_space, transition_matrix)
    end
end

required_row_sum(::Core.Type{<:AbstractContinuousMarkovChain}) = 0
required_row_sum(::Type{<:AbstractContinuousMarkovChain}, ::Type{T}) where T = zero(T)
function ContinuousMarkovChain(transition_matrix)
    return ContinuousMarkovChain(1:(size(transition_matrix)[1]), transition_matrix)
end

function characteristic_matrix(::AbstractContinuousMarkovChain)
    return LinearAlgebra.UniformScaling(0)
end

probability_matrix(x::AbstractContinuousMarkovChain) = exp(transition_matrix(x))

Base.convert(
    ::Type{ContinuousMarkovChain},
    x::AbstractDiscreteMarkovChain,
) = ContinuousMarkovChain(x)
function ContinuousMarkovChain(x::AbstractDiscreteMarkovChain)
    S = state_space(x)
    T = transition_matrix(x)

    if length(S) == 0
        return ContinuousMarkovChain(S, T)
    end

    Q = log(T)
    return ContinuousMarkovChain(S, Q)
end

Base.convert(
    ::Type{DiscreteMarkovChain},
    x::AbstractContinuousMarkovChain,
) = DiscreteMarkovChain(x)
function DiscreteMarkovChain(x::AbstractContinuousMarkovChain)
    S = state_space(x)
    Q = transition_matrix(x)

    if length(S) == 0
        return DiscreteMarkovChain(S, Q)
    end

    T = exp(Q)
    return DiscreteMarkovChain(S, T)
end

function embedded(x::AbstractContinuousMarkovChain)
    Q = transition_matrix(x)

    if size(Q)[1] == 0
        return DiscreteMarkovChain(state_space(x), Q)
    end

    D = LinearAlgebra.diag(Q)
    diag_Q_inv = LinearAlgebra.diagm(1 ./ D)
    T = -diag_Q_inv*Q
    T[LinearAlgebra.diagind(T)] .= 0  # Sometimes NaNs appear
    T[T .=== NaN] .= 0

    for row in 1:(size(Q)[1])
        if D[row] == 0
            T[row, :] .= 0
            T[row, row] = 1
        end
    end

    return DiscreteMarkovChain(state_space(x), T)
end

function mean_recurrence_time(x::AbstractContinuousMarkovChain)
    Q = transition_matrix(x)
    T = transition_matrix(embedded(x))
    n = size(Q)[1]
    result = []

    for state in 1:n
        new_mat = copy(Q)
        new_mat[state, :] .= 0
        y = ContinuousMarkovChain(new_mat)
        states, _, _, _ = decompose(embedded(y))

        lag = -1/Q[state, state]
        values = mean_time_to_absorption(y)
        probabilities = T[state, setdiff(states, [state])]

        push!(result, lag+(values'*probabilities)[1, 1])
    end

    return result
end

function mean_first_passage_time(x::AbstractContinuousMarkovChain)
    Q = transition_matrix(x)
    n = size(Q)[1]
    result = Array{Any}(undef, n, n)

    for state in 1:n
        new_mat = copy(Q)
        new_mat[state, :] .= 0
        y = ContinuousMarkovChain(new_mat)
        states, _, _, _ = decompose(embedded(y))

        values = mean_time_to_absorption(y)

        result[setdiff(states, [state]), state] = values
    end

    result[LinearAlgebra.diagind(result)] .= 0

    return result
end

ext: jl


text:
T=9
tsname="joExtend"
@testset "$tsname" begin
for t=1:T # start test loop
    n=8*t
    l=rand(0:n)
    u=rand(0:n)
    m=n+l+u
    Az=joExtend(n,:zeros,pad_lower=l,pad_upper=u)
    Ab=joExtend(n,:border,pad_lower=l,pad_upper=u)
    Am=joExtend(n,:mirror,pad_lower=l,pad_upper=u)
    Ap=joExtend(n,:periodic,pad_lower=l,pad_upper=u)

    verbose && println("$tsname ($n,$l,$u)")
    @testset "$n x $m" begin
        @test isadjoint(Az)[1]
        @test isadjoint(Ab)[1]
        @test isadjoint(Am)[1]
        @test isadjoint(Ap)[1]
        @test islinear(Az)[1]
        @test islinear(Ab)[1]
        @test islinear(Am)[1]
        @test islinear(Ap)[1]
    end

end # end test loop
end


text:
opensmt_expr = Ptr{Void}
opensmt_context = Ptr{Void}

## Communication APIs
## ==================
opensmt_init() = ccall( (:opensmt_init, "libdreal"), Void, ())

opensmt_set_verbosity(ctx::opensmt_context, level::Cint) =
  ccall((:opensmt_set_verbosity, "libdreal"), Void, (Ptr{Void}, Cint), ctx, level)

opensmt_set_precision(ctx::opensmt_context, p::Float64) =
  ccall( (:opensmt_set_precision, "libdreal"), Ptr{Void}, (Ptr{Void}, Float64), ctx, p)

opensmt_get_precision(ctx::opensmt_context, p::Float64) =
  ccall( (:opensmt_get_precision, "libdreal"), Ptr{Void}, (Ptr{Void}, Float64), ctx, p)

opensmt_version() =
  ccall( (:opensmt_version, "libdreal"), Ptr{UInt8}, ())

opensmt_print_expr(e::opensmt_expr) =
  ccall( (:opensmt_print_expr, "libdreal"), Void, (Ptr{Void},), e)

opensmt_mk_context(l::Cuint) =
  ccall( (:opensmt_mk_context, "libdreal"), Ptr{Void}, (Cuint,), l)

opensmt_del_context(ctx::opensmt_context) =
  ccall( (:opensmt_del_context, "libdreal"), Void, (Ptr{Void},), ctx)

opensmt_reset(ctx::opensmt_context) =
  ccall( (:opensmt_reset, "libdreal"), Void, (Ptr{Void},), ctx)

opensmt_push(ctx::opensmt_context) =
  ccall( (:opensmt_push, "libdreal"), Void, (Ptr{Void},), ctx)

opensmt_pop(ctx::opensmt_context) =
  ccall( (:opensmt_pop, "libdreal"), Void, (Ptr{Void},), ctx)

opensmt_assert(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_assert, "libdreal"), Void,
    (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_check(ctx::opensmt_context) =
  ccall((:opensmt_check, "libdreal"), Cint, (Ptr{Void},), ctx)

opensmt_check_assump(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_check_assump, "libdreal"), Cuint, (Ptr{Void},Ptr{Void}), ctx, e)

opensmt_check_lim_assump(ctx::opensmt_context, e::opensmt_expr, i::Cuint) =
  ccall((:opensmt_check_assump, "libdreal"), Cuint, (Ptr{Void},Ptr{Void}), ctx, e, i)

opensmt_conflicts(ctx::opensmt_context) =
  ccall((:opensmt_conflicts, "libdreal"), Cuint, (Ptr{Void},), ctx)

opensmt_decisions(ctx::opensmt_context) =
  ccall((:opensmt_decisions, "libdreal"), Cuint, (Ptr{Void},), ctx)

opensmt_get_value(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_get_value, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_get_lb(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_get_lb, "libdreal"), Float64, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_get_ub(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_get_ub, "libdreal"), Float64, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_get_bool(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_get_bool, "libdreal"), Cint, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_prefer(e::opensmt_expr) =
  ccall((:opensmt_prefer, "libdreal"), Cuint, (Ptr{Void},), e)

opensmt_polarity(ctx::opensmt_context, e::opensmt_expr, pos::Cint) =
  ccall((:opensmt_polarity, "libdreal"), Cuint, (Ptr{Void}, Ptr{Void}), ctx, e, pos)

opensmt_print_model(ctx::opensmt_context, c::Ptr{UInt8}) =
  ccall((:opensmt_print_model, "libdreal"), Void, (Ptr{Void}, Ptr{UInt8}), ctx, c)

opensmt_print_proof(ctx::opensmt_context, c::Ptr{UInt8}) =
  ccall((:opensmt_print_proof, "libdreal"), Void, (Ptr{Void}, Ptr{UInt8}), ctx, c)

opensmt_print_interpolant(ctx::opensmt_context, c::Ptr{UInt8}) =
  ccall((:opensmt_print_interpolant, "libdreal"), Void, (Ptr{Void}, Ptr{UInt8}), ctx, c)

## Formula Construction APIs
## =========================

opensmt_mk_true(ctx::opensmt_context) =
  ccall((:opensmt_mk_true, "libdreal"), Ptr{Void}, (Ptr{Void},), ctx)

opensmt_mk_false(ctx::opensmt_context) =
  ccall((:opensmt_mk_false, "libdreal"), Ptr{Void}, (Ptr{Void},), ctx)

opensmt_mk_bool_var(ctx::opensmt_context, varname::Ptr{UInt8}) =
  ccall((:opensmt_mk_bool_var, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_bool_var(ctx,varname::String) = opensmt_mk_bool_var(ctx, pointer(varname))

opensmt_mk_real_var(ctx::opensmt_context, varname::String, lb::Float64, ub::Float64) =
  ccall((:opensmt_mk_real_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}, Float64, Float64), ctx, varname, lb, ub)

opensmt_mk_int_var(ctx::opensmt_context, varname::String, lb::Int64, ub::Int64) =
  ccall((:opensmt_mk_int_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}, Int64, Int64), ctx, varname, lb, ub)

opensmt_mk_unbounded_int_var(ctx::opensmt_context, varname::String) =
  ccall((:opensmt_mk_unbounded_int_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_forall_int_var(ctx::opensmt_context, varname::String, lb::Int64, ub::Int64) =
  ccall((:opensmt_mk_forall_int_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}, Int64, Int64), ctx, varname, lb, ub)

opensmt_mk_forall_unbounded_int_var(ctx::opensmt_context, varname::String) =
  ccall((:opensmt_mk_forall_unbounded_int_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_unbounded_real_var(ctx::opensmt_context, varname::String) =
  ccall((:opensmt_mk_unbounded_real_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_forall_real_var(ctx::opensmt_context, varname::String, lb::Float64, ub::Float64) =
  ccall((:opensmt_mk_forall_real_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}, Float64, Float64), ctx, varname, lb, ub)

opensmt_mk_forall_unbounded_real_var(ctx::opensmt_context, varname::String) =
  ccall((:opensmt_mk_forall_real_var, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_cost_var(ctx::opensmt_context, varname::Ptr{UInt8}) =
  ccall((:opensmt_mk_cost_var, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{UInt8}), ctx, varname)

opensmt_mk_cost_var(ctx,varname::String) = opensmt_mk_cost_var(ctx, pointer(varname))

opensmt_mk_or(ctx::opensmt_context, es::Vector{opensmt_expr}, i::Cuint) =
  ccall((:opensmt_mk_or, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{opensmt_expr}, Cuint), ctx, es, i)

opensmt_mk_or_2(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_or_2, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_or_3(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr, e3::opensmt_expr) =
  ccall((:opensmt_mk_or_3, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2, e3)

opensmt_mk_and(ctx::opensmt_context, es::Vector{opensmt_expr}, i::Cuint) =
  ccall((:opensmt_mk_and, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{opensmt_expr}, Cuint), ctx, es, i)

opensmt_mk_and_2(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_and_2, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_and_3(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr, e3::opensmt_expr) =
  ccall((:opensmt_mk_and_3, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2, e3)

opensmt_mk_eq(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_eq, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_ite(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr, e3::opensmt_expr) =
  ccall((:opensmt_mk_ite, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2, e3)

opensmt_mk_not(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_not, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_num(ctx::opensmt_context, c::Float64) =
  ccall((:opensmt_mk_num, "libdreal"), Ptr{Void},
        (Ptr{Void}, Float64), ctx, c)

opensmt_mk_num_from_string(ctx::opensmt_context, c::String) =
  ccall((:opensmt_mk_num_from_string, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{UInt8}), ctx, c)

opensmt_mk_plus(ctx::opensmt_context, es::Vector{opensmt_expr}, i::Cuint) =
  ccall((:opensmt_mk_plus, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{opensmt_expr}, Cuint), ctx, es, i)

opensmt_mk_plus_2(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_plus_2, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_plus_3(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr, e3::opensmt_expr) =
  ccall((:opensmt_mk_plus_3, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2, e3)

opensmt_mk_uminus(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_uminus, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_minus(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_minus, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_times(ctx::opensmt_context, es::Vector{opensmt_expr}, i::Cuint) =
  ccall((:opensmt_mk_times, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{opensmt_expr}, Cuint), ctx, es, i)

opensmt_mk_times_2(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_times_2, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_times_3(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr, e3::opensmt_expr) =
  ccall((:opensmt_mk_times_3, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2, e3)

opensmt_mk_div(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_div, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_lt(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_lt, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_leq(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_leq, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_gt(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_gt, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_geq(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_geq, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_abs(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_abs, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_exp(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_exp, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_log(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_log, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_sqrt(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_sqrt, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_pow(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_pow, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_sin(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_sin, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_cos(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_cos, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_tan(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_tan, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_asin(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_asin, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_acos(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_acos, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_atan(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_atan, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_sinh(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_sinh, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_cosh(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_cosh, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_tanh(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_tanh, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_atan2(ctx::opensmt_context, e::opensmt_expr) =
  ccall((:opensmt_mk_atan2, "libdreal"), Ptr{Void}, (Ptr{Void}, Ptr{Void}), ctx, e)

opensmt_mk_max(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_max, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

opensmt_mk_min(ctx::opensmt_context, e1::opensmt_expr, e2::opensmt_expr) =
  ccall((:opensmt_mk_min, "libdreal"), Ptr{Void},
        (Ptr{Void}, Ptr{Void}, Ptr{Void}), ctx, e1, e2)

## ODE
function opensmt_define_ode(
  ctx::opensmt_context, varname::String,
  es1::Vector{opensmt_expr}, es2::Vector{opensmt_expr}, n::Cuint)
  ccall((:opensmt_define_ode, "libdreal"),
        Void,
        (Ptr{Void},Ptr{UInt8}, Ptr{opensmt_expr}, Ptr{opensmt_expr}, Cuint),
        ctx, varname, es1, es2, n)
end

function opensmt_mk_integral(
  ctx::opensmt_context,
  vars_t::Vector{opensmt_expr},
  time_l::opensmt_expr,
  time_u::opensmt_expr,
  vars_0::Vector{opensmt_expr},
  n::Cuint,
  flowname::String)
  ccall((:opensmt_mk_integral, "libdreal"),
        Ptr{Void},
        (Ptr{Void}, Ptr{opensmt_expr}, Ptr{Void}, Ptr{Void}, Ptr{opensmt_expr}, Cuint,
        Ptr{UInt8}),
        ctx, vars_t, time_l, time_u, vars_0, n, flowname)
end


text:
abstract type AbstractNode{T, N, p} <: AbstractArray{T, N} end

@pure Base.size(::Type{Tnode}) where {Tnode <: AbstractNode} = convert.(Int, TreeSize(Tnode)[1])
@pure totalsize(::Type{Tnode}) where {Tnode <: AbstractNode} = totalsize(TreeSize(Tnode))

# null
@pure null(::Type{Tnode}) where {Tnode <: AbstractNode} = Tnode(undef)
null(x::AbstractNode) = null(typeof(x))
isnull(x::AbstractNode) = x === null(x)

# childtype
@pure childtype(::Type{<: AbstractNode}) = nothing
@pure childtype(::Type{<: AbstractNode{T}}) where {T <: AbstractNode} = T
childtype(x::AbstractNode) = childtype(typeof(x))

# leaftype
@pure leaftype(T::Type{<: AbstractNode}) = T
@pure leaftype(::Type{<: AbstractNode{T}}) where {T <: AbstractNode} = leaftype(T)
leaftype(x::AbstractNode) = leaftype(typeof(x))

# leafeltype
@pure leafeltype(T::Type{<: AbstractNode}) = eltype(leaftype(T))
leafeltype(x::AbstractNode) = leafeltype(typeof(x))

# leafblock
nleafblocks(node::AbstractNode) = Int.(totalsize(node) .÷ size(leaftype(node)))
leafblockunit(node::AbstractNode) = TreeSize(node)[end]

getmask(x::AbstractNode) = getmask(x.data)
isactive(x::AbstractNode, i...) = (@_propagate_inbounds_meta; isnull(x) ? false : getmask(x)[i...])

@inline function unsafe_getindex(x::AbstractNode, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds unsafe_getindex(x.data, i)
end

@inline function unsafe_setindex!(x::AbstractNode, v, i::Int)
    @boundscheck checkbounds(x, i)
    @inbounds unsafe_setindex!(x.data, v, i)
end

@inline function unsafe_getindex(x::AbstractNode, i::Int...)
    @boundscheck checkbounds(x, i...)
    @inbounds unsafe_getindex(x, Base._sub2ind(size(x), i...))
end

@inline function unsafe_setindex!(x::AbstractNode, v, i::Int...)
    @boundscheck checkbounds(x, i...)
    @inbounds unsafe_setindex!(x, v, Base._sub2ind(size(x), i...))
end

function allocate!(node::AbstractNode, mask::AbstractArray{Bool})
    checkbounds(CartesianIndices(totalsize(node)), CartesianIndices(mask))
    dims = totalsize(Base.tail(TreeSize(node)))
    for I in CartesianIndices(node)
        start = CartesianIndex(@. dims * ($Tuple(I) - 1) + 1)
        stop = min(CartesianIndex(@. $Tuple(start) + dims - 1), last(CartesianIndices(mask)))
        childmask = @view mask[start:stop]
        if any(childmask)
            child = allocate!(node, I)
            allocate!(child, childmask)
        end
    end
end

function cleanup!(x::AbstractNode)
    @inbounds for i in eachindex(x)
        if isactive(x, i)
            # if `i` is active, then child node is always not null
            childnode = unsafe_getindex(x, i)
            cleanup!(childnode)
            !any(getmask(childnode)) && delete!(x, i)
        else
            delete!(x, i)
        end
    end
    x
end

function nleaves(x::AbstractNode)
    isnull(x) && return 0
    count = 0
    @inbounds for i in eachindex(x)
        if isactive(x, i)
            child = unsafe_getindex(x, i)
            count += nleaves(child)
        end
    end
    count
end

struct Dynamic end
struct TreeSize{S}
    function TreeSize{S}() where {S}
        new{S::Tuple{Vararg{Union{Tuple{Vararg{Integer}}, Dynamic}}}}()
    end
end
@pure TreeSize(S::Tuple) = TreeSize{S}()
@pure Base.Tuple(::TreeSize{S}) where {S} = S

@pure nodesize(::Type{<: AbstractNode{<: Any, <: Any, Dynamic()}}) = Dynamic()
@pure nodesize(::Type{<: AbstractNode{<: Any, N, p}}) where {N, p} = nfill(Power2(p), Val(N))
@pure TreeSize(::Nothing) = TreeSize{()}()
@pure TreeSize(::Type{Tnode}) where {Tnode <: AbstractNode} = TreeSize{(nodesize(Tnode), Tuple(TreeSize(childtype(Tnode)))...)}()
@pure TreeSize(x::AbstractNode) = TreeSize(typeof(x))

@pure Base.length(::TreeSize{S}) where {S} = length(S)
@pure Base.getindex(::TreeSize{S}, i::Int) where {S} = S[i]
@pure Base.firstindex(::TreeSize) = 1
@pure Base.lastindex(x::TreeSize) = length(x)
@pure Base.tail(::TreeSize{S}) where {S} = length(S) == 1 ? TreeSize{(map(zero, S[1]),)}() : TreeSize{Base.tail(S)}()
@pure totalsize(::TreeSize{S}) where {S} = broadcast(*, S...)

Base.show(io::IO, ::TreeSize{S}) where {S} = print(io, "TreeSize", S)

struct Activated{P, I}
    parent::P
    index::I
end
@inline set!(x::Activated, v) = @inbounds unsafe_setindex!(x.parent, v, x.index)
@inline set!(x::Activated, name::Symbol, v) = @inbounds unsafe_setindex!(getproperty(x.parent, name), v, x.index)
@inline function allocate!(x::Activated, i...)
    @_propagate_inbounds_meta
    child = @inbounds unsafe_getindex(x.parent, x.index)
    allocate!(child, i...)
end
@inline function Base.setindex!(x::Activated, v, i...)
    @_propagate_inbounds_meta
    child = @inbounds unsafe_getindex(x.parent, x.index)
    setindex!(child, v, i...)
end


text:
using Printf

# This is a copy paste from the examples in CImGui

function TestFunc(clear_color=copy(CImGuiFrontEnd.default_clear_color))
    @cstatic show_demo_window=false show_another_window=false begin
        # show a simple window that we create ourselves.
        # we use a Begin/End pair to created a named window.
        @cstatic f=Cfloat(0.0) counter=Cint(0) begin
            CImGui.Begin("Hello, world!")  # create a window called "Hello, world!" and append into it.
            CImGui.Text("This is some useful text.")  # display some text
            @c CImGui.Checkbox("Demo Windows", &show_demo_window)  # edit bools storing our window open/close state
            @c CImGui.Checkbox("Another Window", &show_another_window)

            @c CImGui.SliderFloat("float", &f, 0, 1)  # edit 1 float using a slider from 0 to 1
            CImGui.ColorEdit3("clear color", clear_color)  # edit 3 floats representing a color
            CImGui.Button("Button") && (counter += 1)

            CImGui.SameLine()
            CImGui.Text("counter = $counter")
            CImGui.Text(@sprintf("Application average %.3f ms/frame (%.1f FPS)", 1000 / CImGui.GetIO().Framerate, CImGui.GetIO().Framerate))

            CImGui.End()
        end

        # show another simple window.
        if show_another_window
            @c CImGui.Begin("Another Window", &show_another_window)  # pass a pointer to our bool variable (the window will have a closing button that will clear the bool when clicked)
            CImGui.Text("Hello from another window!")
            CImGui.Button("Close Me") && (show_another_window = false;)
            CImGui.End()
        end

        if show_demo_window
            CImGui.ShowStyleSelector("asdf")
            CImGui.Begin("asdf2")
            CImGui.ShowStyleEditor()
            CImGui.End()
            CImGui.ShowFontSelector("asdf3")
        end
    end

end

DoGui(TestFunc)


text:
module DashAladinLite
using Dash

const resources_path = realpath(joinpath( @__DIR__, "..", "deps"))
const version = "0.1.1"

include("jl/dal_dashaladinlite.jl")

function __init__()
    DashBase.register_package(
        DashBase.ResourcePkg(
            "dash_aladin_lite",
            resources_path,
            version = version,
            [
                DashBase.Resource(
    relative_package_path = "third_party/aladin.min.js",
    external_url = "https://aladin.u-strasbg.fr/AladinLite/api/v2/latest/aladin.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_aladin_lite.min.js",
    external_url = "https://unpkg.com/dash_aladin_lite@0.1.1/dash_aladin_lite/dash_aladin_lite.min.js",
    dynamic = nothing,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "dash_aladin_lite.min.js.map",
    external_url = "https://unpkg.com/dash_aladin_lite@0.1.1/dash_aladin_lite/dash_aladin_lite.min.js.map",
    dynamic = true,
    async = nothing,
    type = :js
),
DashBase.Resource(
    relative_package_path = "third_party/aladin.min.css",
    external_url = "https://aladin.u-strasbg.fr/AladinLite/api/v2/latest/aladin.min.css",
    dynamic = nothing,
    async = nothing,
    type = :css
)
            ]
        )

    )
end
end


text:
# Basically loads the weights tests on the dev set after each epoch.
using CorpusLoaders
using MultiResolutionIterators
using TextAnalysis
using TextAnalysis: CRF, crf_loss
using WordTokenizers
using Embeddings
using Flux
using Flux: onehot, param, onehotbatch, Conv, LSTM, flip, Momentum, reset!, onecold
using Tracker
using BSON: @save, @load
using LinearAlgebra

println("Dependecies Loaded")

CHAR_EMBED_DIMS = 25
WORD_EMBED_DIMS = 50
DESIRED_TAG_SCHEME = "BIOES"

train_set = load(CoNLL(), "train") # training set
# test_set = load(CoNLL(), "test") # test set
dev_set = load(CoNLL(), "dev") # dev set

dataset = flatten_levels(train_set, lvls(CoNLL, :document)) |> full_consolidate
dev_dataset = flatten_levels(dev_set, lvls(CoNLL, :document)) |> full_consolidate

X_train = [CorpusLoaders.word.(sent) for sent in dataset]
Y_train = [CorpusLoaders.named_entity.(sent) for sent in dataset]
tag_scheme!.(Y_train, "BIO2", DESIRED_TAG_SCHEME)
@assert length.(X_train) == length.(Y_train)

X_dev = [CorpusLoaders.word.(sent) for sent in dev_dataset]
Y_dev = [CorpusLoaders.named_entity.(sent) for sent in dev_dataset]
tag_scheme!.(Y_dev, "BIO2", DESIRED_TAG_SCHEME)
@assert length.(X_dev) == length.(Y_dev)

words_vocab = unique(vcat(X_train...))
alphabets = unique(vcat(collect.(words_vocab)...))
labels = unique(vcat(Y_train...))
DENSE_OUT_SIZE = num_labels = length(labels)

#### Word Embeddings
# using these indices to generate the flux one hot input word embeddings vectors.

embtable = load_embeddings(GloVe)
get_word_index = Dict(word => ii for (ii, word) in enumerate(embtable.vocab))
get_word_from_index = Dict(value => key for (key, value) in get_word_index)
W_word_Embed = (embtable.embeddings)

# One Vec for unknown chars
UNK_Word = "<UNK>"
UNK_Word_Idx = length(get_word_index) + 1
embedding_vocab_length = length(get_word_index) + 1
@assert UNK_Word ∉ collect(keys(get_word_index))

get_word_index[UNK_Word] = UNK_Word_Idx
get_word_from_index[UNK_Word_Idx] = UNK_Word
W_word_Embed = hcat(W_word_Embed, rand(WORD_EMBED_DIMS))
W_word_Embed = Float32.(W_word_Embed)
@assert size(W_word_Embed, 2) == embedding_vocab_length

##### Char Embeddings
UNK_char = '¿'
@assert UNK_char ∉ alphabets
push!(alphabets, UNK_char)
sort!(unique!(alphabets))

# using these indices to generate the flux one hot input Char embeddings vectors.

W_Char_Embed = rand(CHAR_EMBED_DIMS, length(alphabets))
W_Char_Embed = Float32.(W_Char_Embed)
W_Char_Embed = param(W_Char_Embed)
get_char_index = Dict(char => ii for (ii, char) in enumerate(alphabets))
get_char_from_index = Dict(value => key for (key, value) in get_char_index)

println("Embeddings Made")

############# Creating onehoot vectors (useful for embeddings lookup), convenient
# TODO: Minibatches
onehotword(word) = onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length)
onehotchars(word) = onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets))
onehotlabel(label) = onehot(label, labels)
onehotinput(word) = (onehot(get(get_word_index, lowercase(word), get_word_index[UNK_Word]), 1:embedding_vocab_length),
                onehotbatch([get(get_char_index, c, get_char_index[UNK_char]) for c in word], 1:length(alphabets)))

oh_seq(arr, f) = [f(element) for element in arr]

# X_words_train = [oh_seq(sentence, onehotword) for sentence in X_train] # A Bunch of Sequences of words, i.e. sentences
# X_chars_train = [oh_seq(sentence, onehotchars) for sentence in X_train] # A Bunch Sequences of Array of Chars, done to prevent repeated computations.
X_input_train = [(oh_seq(sentence, onehotinput)) for sentence in X_train]
Y_oh_train = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_train]

Y_oh_train = (Y_oh_train)

X_input_dev = [oh_seq(sentence, onehotinput) for sentence in X_dev]
Y_oh_dev = [oh_seq(tags_sequence, onehotlabel) for tags_sequence in Y_dev]

@assert length.(X_train) == length.(X_input_train) == length.(Y_oh_train)

#################### MODEL ######################
# 1. Character Embeddings with CNNs
# 2. Obtain Word Embeddings, concatenate Char & Word Embeddings
# 3. Bi-LSTM
# 4. Softmax or CRF on each unit in sequence

CNN_OUTPUT_SIZE = 53
CONV_WINDOW_LENGTH = 3
LSTM_STATE_SIZE = 253
DROPOUT_RATE_CNN = DROPOUT_INPUT_EMBEDDING = DROPOUT_OUT_LAYER = 0.68

# 1. Dropout before conv, Max poll layer, PADDING on both sides, hyperparams = window size, output vector size.
conv1 = Conv((CHAR_EMBED_DIMS, CONV_WINDOW_LENGTH), 1=>CNN_OUTPUT_SIZE, pad=(0,2))
dropout1 = Dropout(DROPOUT_RATE_CNN)

char_features = Chain(x -> W_Char_Embed * x,
                      x -> reshape(x, size(x)..., 1,1),
                      dropout1,
                      conv1,
                      x -> maximum(x, dims=2),
                      x -> reshape(x, length(x),1))

# 2. Input embeddings:
# Maybe could use only the embeddings for the words needed

W_word_Embed = param(W_word_Embed)# For trainable
get_word_embedding(w) = W_word_Embed * w# works coz - onehot

# Dropout before LSTM
dropout_embed = Dropout(DROPOUT_INPUT_EMBEDDING)
input_embeddings((w, cs)) = dropout_embed(vcat(get_word_embedding(w), char_features(cs)))

# 3. Bi-LSTM

forward_lstm = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
backward = LSTM(CNN_OUTPUT_SIZE + WORD_EMBED_DIMS, LSTM_STATE_SIZE)
backward_lstm(x) = reverse(backward.(reverse(x)))

bilstm_layer(x) = vcat.(forward_lstm.(x), backward_lstm(x))

# TODO: Bias vectors are initialized to zero, except the bias bf for the forget gate in LSTM , which is initialized to 1.0

# 4. Softmax / CRF Layer
# Dropout after LSTM

dropout2 = Dropout(DROPOUT_OUT_LAYER)
d_out = Dense(LSTM_STATE_SIZE * 2, DENSE_OUT_SIZE + 2)

m = Chain(x -> input_embeddings.(x),
                bilstm_layer,
                x -> dropout2.(x),
                x -> d_out.(x))

using TextAnalysis: crf_loss, CRF

Flux.@treelike TextAnalysis.CRF
c = TextAnalysis.CRF(num_labels)

init_α = fill(-10000, (c.n + 2, 1))
init_α[c.n + 1] = 0

loss(x, y) =  crf_loss(c, m(x), y, init_α)

η = 0.005 # learning rate
β = 0.05 # rate decay
ρ = 0.9 # momentum

# TODO: rate decay
opt = Flux.Optimiser(ExpDecay(β), Momentum(η, ρ))
data = zip(X_input_train, Y_oh_train)

println("Model Built")

function load_weights(epoch, conv1, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c)

    conv_cpu = conv1 |> cpu
    W_word_cpu = W_word_Embed |> cpu
    W_char_cpu = W_Char_Embed |> cpu
    forward_lstm_cpu = forward_lstm |> cpu
    backward_lstm_cpu = backward |> cpu
    crf_cpu = c |> cpu
    d_cpu = d_out |> cpu

    @load "./weights/conv_cpu$(epoch).bson" conv_cpu
    @load "./weights/W_word_cpu$(epoch).bson" W_word_cpu
    @load "./weights/W_char_cpu$(epoch).bson" W_char_cpu
    @load "./weights/forward_lstm$(epoch).bson" forward_lstm_cpu
    @load "./weights/backward_lstm$(epoch).bson" backward_lstm_cpu
    @load "./weights/d_cpu.bson$(epoch)" d_cpu
    @load "./weights/crf.bson$(epoch)" crf_cpu

    return conv_cpu, W_word_cpu, W_char_cpu, forward_lstm_cpu, backward_lstm_cpu, d_cpu, crf_cpu
end

function try_outs_for_confusion_mat()
    Flux.testmode!(dropout1)
    Flux.testmode!(dropout2)
    Flux.testmode!(dropout_embed)

    confusion_matrix = zeros(Int, (19, 19))

    for (x, y) in zip(X_input_dev, Y_oh_dev)
        reset!(forward_lstm)
        reset!(backward)
        x_seq = m(x)
        ohs = TextAnalysis.viterbi_decode(c, x_seq, init_α)

        for d in zip(ohs, y)
            confusion_matrix[d[1].ix, d[2].ix] += 1
        end
    end

    reset!(forward_lstm)
    reset!(backward)
    Flux.testmode!(dropout1, false)
    Flux.testmode!(dropout2, false)
    Flux.testmode!(dropout_embed, false)

    return confusion_matrix
end

# org_idx = [1, 7, 8, 14]
# misc_idx = [3,11, 12, 13]
# per_idx = [4, 5, 9, 10]
# loc_idx = [6, 15, 16 ,17]
# o_idx = [2]

function find_precision(confusion_matrix)
    s = sum(confusion_matrix, dims=1)'
    dg = diag(confusion_matrix)
    labels = [[1, 7, 8, 14], [3,11, 12, 13], [4, 5, 9, 10], [6, 15, 16 ,17], [2]]
    label_wise(labels_indices) = sum([dg[i] for i in labels_indices]) / sum([s[i] for i in labels_indices])
    return sum([label_wise(labels_indices) for labels_indices in labels]) / 5
end

function find_recall(confusion_matrix)
    s = sum(confusion_matrix, dims=2)
    dg = diag(confusion_matrix)
    labels = [[1, 7, 8, 14], [3,11, 12, 13], [4, 5, 9, 10], [6, 15, 16 ,17], [2]]
    label_wise(labels_indices) = sum([dg[i] for i in labels_indices]) / sum([s[i] for i in labels_indices])
    return sum([label_wise(labels_indices) for labels_indices in labels]) / 5
end

function call_for_f1(epoch)
    conv1, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c =
            load_weights(epoch, conv1, W_word_Embed, W_Char_Embed, forward_lstm, backward, d_out, c)

    Flux.testmode!(dropout1)
    Flux.testmode!(dropout2)
    Flux.testmode!(dropout_embed)

    confusion_matrix = zeros(Int, (19, 19))
    for (x, y) in zip(X_input_dev, Y_oh_dev)
        reset!(forward_lstm)
        reset!(backward)
        x_seq = m(x)
        ohs = TextAnalysis.viterbi_decode(c, x_seq, init_α)

        for d in zip(ohs, y)
            confusion_matrix[d[1].ix, d[2].ix] += 1
        end
    end

    prec = find_precision(confusion_matrix)
    rc = find_recall(confusion_matrix)
    println("Precision and recall are:", prec, " ", rc)
    println("F1 is:", (2 * prec * rc) / (prec + rc))

    Flux.testmode!(dropout1, false)
    Flux.testmode!(dropout2, false)
    Flux.testmode!(dropout_embed, false)

end


text:
export Arity, Image, Signature, Closure
export Associative, Commutative

"""
    Arity <: Property

Representation of function `f` having the arity `n`, `Arity(f, n)`.
"""
struct Arity <: Property
    f
    n::Int
end

"""
    Image <: Property

Representation of expression `ex` having the mathematical image `S::AbstractSet`,
`Image(ex, S)`.
"""
struct Image{T} <: Property
    ex::T
    S::AbstractSet
end
isvalid(i::Image{<:Union{Symbolic,Variable,Expr}}, ctx::Context) =
    invoke(isvalid, Tuple{Property,Context}, i, ctx)
isvalid(i::Image, ctx::Context) =
    i.ex ∈ i.S || invoke(isvalid, Tuple{Property,Context}, i, ctx)
implies(i::Image, i′::Image, ::Context) = i.ex == i′.ex && i.S ⊆ i′.S

"""
    Signature <: Property

Representation of function `f` having signature `(x₁, x₂, …, xₙ) -> y` given `xᵢ` and `y` as
image sets, `Signature(f, [x₁, x₂, …, xₙ], y)`.
"""
struct Signature <: Property
    f
    args::Vector{AbstractSet}
    S::AbstractSet
end
implies(sig::Signature, a::Arity, ::Context) = a.f == sig.f && a.n == length(sig.args)
implies(sig::Signature, i::Image{Expr}, ctx::Context) =
    sig.S ⊆ i.S && i.ex.head === :call && length(i.ex.args) == 1 + length(sig.args) &&
        sig.f == i.ex.args[1] && all(isvalid(ctx), Image.(i.ex.args[2:end], sig.args))

"""
    Closure <: Property

Representation of function `f` being closed under set `S`, `Closure(f, S)`.
"""
struct Closure <: Property
    f
    S::AbstractSet
end
implies(p::Closure, q::Closure, ::Context) = p.f == q.f && q.S ⊆ p.S
implies(p::Closure, sig::Signature, ::Context) =
    p.f == sig.f && q.S ⊆ p.S && all(sig.args .⊆ Ref(p.S))
implies(p::Closure, i::Image{Expr}, ctx::Context) =
    p.S ⊆ i.S && i.ex.head === :call && !isempty(i.ex.args) &&
        p.f == i.ex.args[1] && all(isvalid(ctx), Image.(i.ex.args[2:end], Ref(p.S)))

"""
    Associative <: Property

Representation of function `f` having the associative property, `Associative(f)`.
"""
struct Associative <: Property
    f
end
implies(p::Associative, a::Arity, ::Context) = a.n == 2 && a.f == p.f

"""
    Commutative <: Property

Representation of function `f` having the commutative property, `Commutative(f)`.
"""
struct Commutative <: Property
    f
end


text:
__precompile__(true)

module Sorting

	function is_sorted(data::Array)
		sorted = true
		for i in 1 : length(data) - 1
			sorted &= data[i + 1] >= data[i]
		end

		return sorted
	end

	function is_sorted(data::Array, length::Int)
		sorted = true
		for i in 1 : length - 1
			sorted &= data[i + 1] >= data[i]
		end

		return sorted
	end

	function median!(a::Array)
		p, r = 1, length(a) + 1
		mid = fld(r + p, 2)
		if (r - p) % 2 == 0
			return (select_element!(a, mid) + select_element!(a, mid - 1)) / 2
		else
			return select_element!(a, mid)
		end
	end

	function median(a::Array)
		return median!(copy(a))
	end

	include("insertionsort.jl")
	include("bubblesort.jl")
	include("quicksort.jl")
	include("mergesort.jl")
	include("bisect.jl")
	include("select.jl")
	include("radix.jl")
	include("heapsort.jl")
	include("filesort.jl")
end


text:
"""
This function does nothing if the product is not of Capacity type.
"""
function update_capacity_revenues!(product::T,
                                scenario_name::String,
                                price_years::NamedTuple{(:start_year, :end_year),
                                    Tuple{Int64, Int64}},
                                capacity_revenues::Vector{Float64},
                                capacity_prices::Dict{String, AxisArrays.AxisArray{Float64, 1}},
                                max_cap::Float64) where T <: Product
    return capacity_revenues
end

"""
This function calculates expected capacity revenues for each project.
Returns capacity revenues.
"""
function update_capacity_revenues!(product::Capacity,
                                scenario_name::String,
                                price_years::NamedTuple{(:start_year, :end_year),
                                    Tuple{Int64, Int64}},
                                capacity_revenues::Vector{Float64},
                                capacity_prices::Dict{String, AxisArrays.AxisArray{Float64, 1}},
                                max_cap::Float64)
    capacity_revenues += capacity_prices[scenario_name][price_years[:start_year]:price_years[:end_year]] .*
                         get_accepted_perc(product)[scenario_name][price_years[:start_year]:price_years[:end_year]] *
                         max_cap *
                         get_derating(product)
    return capacity_revenues
end

"""
This function updates the expected profit array of the project
if the product is a capacity product.
Returns nothing.
"""
function update_profit!(project::P,
                        scenario_name::String,
                        product::Capacity,
                        operating_profit::AxisArrays.AxisArray{Float64, 2},
                        capacity_revenues::Vector{Float64},
                        capacity_forward_years::Int64,
                        rec_revenues::Vector{Float64},
                        price_years::NamedTuple{(:start_year, :end_year),
                                     Tuple{Int64, Int64}},
                        update_years::NamedTuple{(:start_year, :end_year),
                                      Tuple{Int64, Int64}},
                        iteration_year::Int64) where P <: Project{<: BuildPhase}

    capacity_revenue_start_year = max(1, capacity_forward_years - price_years[:start_year] + 1)

        for i in capacity_revenue_start_year:length(capacity_revenues)
                set_scenario_profit!(get_finance_data(project),
                            scenario_name,
                            get_name(product),
                            iteration_year,
                            update_years[:start_year] + i - 1,
                            capacity_revenues[i])
        end

    return
end


text:
"""
```
OneAssetHANK{T} <: AbstractCTModel{T}
```

The `OneAssetHANK` type defines the structure of the one-asset HANK model
originally written in MATLAB by Ben Moll and SeHyoun Ahn at
https://sehyoun.com/EXAMPLE_one_asset_HANK_web.html

### Fields

#### Parameters and Steady-States
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant model
  parameters.

* `steady_state::Vector{AbstractParameter}`: Model steady-state values, computed
  as a function of elements of `parameters`.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for all model
  parameters and steady-states to their indices in `parameters` and
  `steady_state`.

#### Inputs to Measurement and Equilibrium Condition Equations

The following fields are dictionaries that map human-readable names to row and
column indices in the matrix representations of of the measurement equation and
equilibrium conditions.

* `endogenous_states::OrderedDict{Symbol,Int}`: Maps each state to a column in
  the measurement and equilibrium condition matrices.

* `exogenous_shocks::OrderedDict{Symbol,Int}`: Maps each shock to a column in
  the measurement and equilibrium condition matrices.

* `expected_shocks::OrderedDict{Symbol,Int}`: Maps each expected shock to a
  column in the measurement and equilibrium condition matrices.

* `equilibrium_conditions::OrderedDict{Symbol,Int}`: Maps each equlibrium
  condition to a row in the model's equilibrium condition matrices.

* `endogenous_states_augmented::OrderedDict{Symbol,Int}`: Maps lagged states to
  their columns in the measurement and equilibrium condition equations. These
  are added after `gensys` solves the model.

* `observables::OrderedDict{Symbol,Int}`: Maps each observable to a row in the
  model's measurement equation matrices.

* `pseudo_observables::OrderedDict{Symbol,Int}`: Maps each pseudo-observable to
  a row in the model's pseudo-measurement equation matrices.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"one_asset_hank\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.

* `pseudo_observable_mappings::OrderedDict{Symbol,PseudoObservable}`: A
  dictionary that stores names and transformations to/from model units. See
  `PseudoObservable` for further details.
"""
mutable struct OneAssetHANK{T} <: AbstractCTModel{T}
    parameters::ParameterVector{T}                          # vector of all time-invariant model parameters
    steady_state::ParameterVector{T}                        # model steady-state values
    keys::OrderedDict{Symbol,Int}                           # human-readable names for all the model
                                                            # parameters and steady-states

    endogenous_states::OrderedDict{Symbol,Vector{Int}}      # these fields used to create matrices in the
    exogenous_shocks::OrderedDict{Symbol,Int}               # measurement & equilibrium condition equations
    expected_shocks::OrderedDict{Symbol,Vector{Int}}        #
    equilibrium_conditions::OrderedDict{Symbol,Vector{Int}} #
    endogenous_states_augmented::OrderedDict{Symbol,Int}    #
    observables::OrderedDict{Symbol,Int}                    #
    pseudo_observables::OrderedDict{Symbol,Int}             #

    spec::String                                            # Model specification number (eg "m990")
    subspec::String                                         # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                          # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                     # Settings/flags for testing mode
    rng::MersenneTwister                                    # Random number generator
    testing::Bool                                           # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
    pseudo_observable_mappings::OrderedDict{Symbol, PseudoObservable}
end

description(m::OneAssetHANK) = "Julia implementation of the one-asset HANK model
originally written in MATLAB by Ben Moll and SeHyoun Ahn at
https://sehyoun.com/EXAMPLE_one_asset_HANK_web.html: OneAssetHANK, $(m.subspec)"

"""
`init_model_indices!(m::OneAssetHANK)`

Arguments:
`m:: OneAssetHANK`: a model object

Description:
Initializes indices for all of `m`'s states, shocks, and equilibrium conditions.
"""
function init_model_indices!(m::OneAssetHANK)
    # Endogenous states
    endogenous_states = collect([:value_function, :inflation, :distribution,
                                 :monetary_policy, :w, :N, :C, :output, :B])

    # Exogenous shocks
    exogenous_shocks = collect([:mp_shock])

    # Expectations shocks
    expected_shocks = collect([:E_V, :E_π])

    # Equilibrium conditions
    equilibrium_conditions = collect([:value_function, :inflation, :distribution,
                                      :monetary_policy, :w, :N, :C, :Y, :B])

    # Additional states added after solving model
    # Lagged states and observables measurement error
    endogenous_states_augmented = []

    # Observables
    observables = keys(m.observable_mappings)

    # Pseudo-observables
    pseudo_observables = keys(m.pseudo_observable_mappings)

    # Assign dimensions
    m.endogenous_states[:value_function]  = collect(1:200)      # size of state space grid
    m.endogenous_states[:inflation]       = [201]               # inflation
    m.endogenous_states[:distribution]    = collect(202:400)    # size of grid minus 1 for ...
    m.endogenous_states[:monetary_policy] = [401]               # monetary policy shocks
    m.endogenous_states[:w]               = [402]               # static conditions from market-clearing
    m.endogenous_states[:N]               = [403]
    m.endogenous_states[:C]               = [404]
    m.endogenous_states[:output]          = [405]
    m.endogenous_states[:B]               = [406]
    m.expected_shocks[:E_V]               = m.endogenous_states[:value_function]
    m.expected_shocks[:E_π]               = m.endogenous_states[:inflation]

    for (i,k) in enumerate(endogenous_states_augmented)
        m.endogenous_states_augmented[k] = i + length(endogenous_states)
    end
    for (i,k) in enumerate(observables);          m.observables[k]          = i end
    for (i,k) in enumerate(exogenous_shocks);     m.exogenous_shocks[k]     = i end
    for (i,k) in enumerate(pseudo_observables);   m.pseudo_observables[k]   = i end
end

function OneAssetHANK(subspec::String="ss0";
                       custom_settings::Dict{Symbol, Setting} = Dict{Symbol, Setting}(),
                       testing = false)

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # Initialize empty model
    m = OneAssetHANK{Float64}(
            # Model parameters and steady state values # before 2nd arg was Vector{Float64}()
            Vector{AbstractParameter{Float64}}(), Vector{Float64}(),
            OrderedDict{Symbol,Int}(),

            # Model indices
            OrderedDict{Symbol,Vector{Int}}(), OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,Vector{Int}}(), OrderedDict{Symbol,Vector{Int}}(),
            OrderedDict{Symbol,Int}(),         OrderedDict{Symbol,Int}(),
            OrderedDict{Symbol,Int}(),

            spec,
            subspec,
            settings,
            test_settings,
            rng,
            testing,
            OrderedDict{Symbol,Observable}(),
            OrderedDict{Symbol,PseudoObservable}())

    # Initialize parameters
    init_parameters!(m)
    init_model_indices!(m)

    # Set settings
    model_settings!(m)
    default_test_settings!(m)

    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Test to see if necessary values are defined, throw error otherwise
    try
        @assert typeof(get_setting(m, :reduce_state_vars)) == Bool
        @assert typeof(get_setting(m, :reduce_v)) == Bool
    catch
        error("Need to specify in settings whether to perform distribution
               and/or value function reduction")
    end

    if get_setting(m, :reduce_state_vars) || get_setting(m, :reduce_v)
        try
            @assert get_setting(m, :n_jump_vars) > 0
            @assert get_setting(m, :n_state_vars) > 0
            @assert get_setting(m, :n_state_vars_unreduce) >= 0
        catch
            error("Need to enter values in settings for n_jump_vars, n_state_vars,
                   and/or n_state_vars_unreduce in Settings. These indicate the
                   number of jump and aggregate state variables.")
        end
    end

    if get_setting(m, :reduce_state_vars)
        try
            @assert get_setting(m, :krylov_dim) > 0
        catch
            error("Need to specify in settings desired dimension of
                   Krylov subspace, krylov_dim.")
        end
    end

    if get_setting(m, :reduce_v)
        try
            @assert get_setting(m, :n_prior) > 0 && typeof(get_setting(m, :n_prior)) == Int
            @assert get_setting(m, :n_post) > 0 && typeof(get_setting(m, :n_post)) == Int
            @assert in(:knots_dict, keys(m.settings))
            @assert in(:spline_grid, keys(m.settings))
        catch
            error("Need to specify in settings n_prior, n_post, knots_dict
                   and/or spline_grid for value function reduction.
                   See reduction.jl for their uses.")
        end
    end

    # Compute steady state
    steadystate!(m)
    return m
end

"""
```
init_parameters!(m::OneAssetHANK)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::OneAssetHANK)
    # Initialize parameters
    m <= parameter(:coefrra, 1.0, fixed=true,
                   description="Relative risk aversion",
                   tex_label="coefrra")

    m <= parameter(:frisch, 0.5, fixed=true,
                   description="Frisch elasticity",
                   tex_label="frisch")

    m <= parameter(:meanlabeff, 3.0, fixed=true,
                   description="meanlabeff: so that at h=1/3 output will be approximately = 1",
                   tex_label="meanlabeff")

    m <= parameter(:maxhours, 1.0, fixed=true,
                   description="maxhours:...",
                   tex_label="maxhours")

    # Production
    m <= parameter(:ceselast, 10., fixed=true,
                   description="CES elasticity",
                   tex_label="ceselast")

    m <= parameter(:priceadjust, 100., fixed=true,
                   description="priceadjust...",
                   tex_label="priceadjust")

    # Policy parameters
    m <= parameter(:taylor_inflation, 1.25, fixed=true,
                   description="Taylor rule coefficient on inflation",
                   tex_label="taylor_inflation")

    m <= parameter(:taylor_outputgap, 0., fixed=true,
                   description="Taylor rule coefficient on output",
                   tex_label="taylor_outputgap")

    m <= parameter(:labtax, 0.2, fixed=true,
                   description="Marginal tax rate on labor income",
                   tex_label="labtax")

    m <= parameter(:govbondtarget, 6., fixed=true,
                   description="govbondtarget: multiple of quarterly GDP",
                   tex_label="govbondtarget")

    m <= parameter(:lumptransferpc, 0.06, fixed=true,
                   description="lumptransferepc: 6% of quarterly GDP in steady state",
                   tex_label="lumptransferepc")

    m <= parameter(:govbcrule_fixnomB, 0., fixed=true,
                   description="govbcrule_fixnomB...",
                   tex_label="govbvrule_fixnomB")

    m <= parameter(:labdisutil, m[:meanlabeff] / ((0.75 ^(-m[:coefrra])) *
                                                  ((1. / 3.)^(1/m[:frisch]))), fixed=true,
                   description="Coefficient of labor disutility",
                   tex_label="labdisutil")

    # Aggregate shocks
    m <= parameter(:σ_MP, sqrt(0.05), (1e-20, 1e5), (1e-20, 1e5), ModelConstructors.Exponential(),
                   RootInverseGamma(4, .4), fixed=false,
                   description="Volatility of monetary policy shocks", tex_label="σ_MP")

    m <= parameter(:θ_MP, 0.25, fixed=true,
                   description="Rate of mean reversion in monetary policy shocks",
                   tex_label="θ_MP")

    m <= SteadyStateParameterArray(:V_ss, Vector{Float64}(),
                                   description = "Stacked steady-state value function")

    m <= SteadyStateParameter(:inflation_ss, NaN, description = "Steady-state rate of inflation")
    m <= SteadyStateParameterArray(:g_ss, Vector{Float64}(),
                                   description = "Stacked steady-state distribution")
    m <= SteadyStateParameter(:r_ss, NaN, description = "Steady-state real interest rate")
    m <= SteadyStateParameterArray(:u_ss, Vector{Float64}(), description = "Steady-state flow utility")
    m <= SteadyStateParameterArray(:c_ss, Vector{Float64}(),
                                   description = "Steady-state individual consumption rate")
    m <= SteadyStateParameterArray(:h_ss, Vector{Float64}(),
                                   description = "Steady-state individual labor hours")
    m <= SteadyStateParameterArray(:s_ss, Vector{Float64}(),
                                   description = "Steady-state individual savings rate")
    m <= SteadyStateParameter(:rnom_ss, NaN, description = "Steady-state nominal interest rate")
    m <= SteadyStateParameter(:B_ss, NaN, description = "Steady-state bond supply")
    m <= SteadyStateParameter(:N_ss, NaN, description = "Steady-state labor supply")
    m <= SteadyStateParameter(:Y_ss, NaN, description = "Steady-state output")
    m <= SteadyStateParameter(:labor_share_ss, NaN, description = "Frictionless labor share of income")
    m <= SteadyStateParameter(:w_ss, NaN, description = "Steady-state wage rate")
    m <= SteadyStateParameter(:profit_ss, NaN, description = "Steady-state profits")
    m <= SteadyStateParameter(:C_ss, NaN, description = "Steady-state aggregate consumption")
    m <= SteadyStateParameter(:T_ss, NaN, description = "Steady-state total taxes")
    m <= SteadyStateParameter(:G_ss, NaN, description = "Steady-state government spending")
    m <= SteadyStateParameter(:ρ_ss, NaN, description = "Steady-state discount rate")
end

"""
```
steadystate!(m::OneAssetHANK)
```

Calculates the model's steady-state values. `steadystate!(m)` must be called whenever
the parameters of `m` are updated.
"""
function steadystate!(m::OneAssetHANK)

    # Read in parameters
    coefrra        = m[:coefrra].value
    frisch         = m[:frisch].value
    meanlabeff     = m[:meanlabeff].value
    maxhours       = m[:maxhours].value
    ceselast       = m[:ceselast].value
    labtax         = m[:labtax].value
    govbondtarget  = m[:govbondtarget].value
    labdisutil     = m[:labdisutil].value
    lumptransferpc = m[:lumptransferpc].value

    # Read in grids
    I                = get_setting(m, :I)
    J                = get_setting(m, :J)
    a                = get_setting(m, :a)
    g_z              = get_setting(m, :g_z)
    zz               = get_setting(m, :zz)
    ymarkov_combined = get_setting(m, :ymarkov_combined)

    # Set necessary variables
    aa   = repeat(a, 1, J)
    amax = maximum(a)
    amin = minimum(a)

    # Read in initial rates
    iterate_r = get_setting(m, :iterate_r)
    r         = get_setting(m, :r0)
    rmin      = get_setting(m, :rmin)
    rmax      = get_setting(m, :rmax)
    iterate_ρ = get_setting(m, :iterate_ρ)
    ρ         = get_setting(m, :ρ0)
    ρmin      = get_setting(m, :ρmin)
    ρmax      = get_setting(m, :ρmax)

    # Read in approximation parameters
    Ir          = get_setting(m, :Ir)
    maxit_HJB   = get_setting(m, :maxit_HJB)
    tol_HJB     = get_setting(m, :tol_HJB)
    Δ_HJB       = get_setting(m, :Δ_HJB)
    maxit_kfe   = get_setting(m, :maxit_kfe)
    tol_kfe     = get_setting(m, :tol_kfe)
    Δ_kfe       = get_setting(m, :Δ_kfe)
    niter_hours = get_setting(m, :niter_hours)
    crit_S      = get_setting(m, :crit_S)

    # Initializing equilibrium objects
    labor_share_ss = (ceselast - 1) / ceselast
    w       = w_ss = labor_share_ss

    # compute initial guesses at steady state values given zz, labor_share_ss, etc.
    N_ss, Y_ss, B_ss, profit_ss, profshare, lumptransfer =
        calculate_ss_equil_vars(zz, labor_share_ss, meanlabeff, lumptransferpc, govbondtarget)

    # Initialize matrices for finite differences
    Vaf = Array{ComplexF64}(undef, I, J)
    Vab = Array{ComplexF64}(undef, I, J)

    cf  = Array{ComplexF64}(undef, I, J) # forward consumption difference
    hf  = Array{ComplexF64}(undef, I, J) # forward hours difference
    sf  = Array{ComplexF64}(undef, I, J) # forward saving difference
    cb  = Array{ComplexF64}(undef, I, J) # backward consumption difference
    hb  = Array{ComplexF64}(undef, I, J) # backward hours difference
    sb  = Array{ComplexF64}(undef, I, J) # backward saving difference
    c0  = Array{ComplexF64}(undef, I, J)
    A   = Array{ComplexF64}(undef, I*J, I*J)

    # Aswitch*v = \lambda_z(v(a,z') - v(a,z))
    # Captures expected change in value function due to jumps in z
    Aswitch = kron(ymarkov_combined, SparseMatrixCSC{ComplexF64}(LinearAlgebra.I, I, I))

    # Initialize steady state variables
    V  = Array{ComplexF64}(undef, I, J) # value function
    u  = Array{ComplexF64}(undef, I, J) # flow utility across state space
    s  = Array{ComplexF64}(undef, I, J) # savings across state space
    c  = Array{ComplexF64}(undef, I, J) # flow consumption
    h  = Array{ComplexF64}(undef, I, J) # flow hours of labor
    h0 = Array{ComplexF64}(undef, I, J) # guess of what h will be

    # Creates functions for computing flow utility, income earned, and labor done given
    # CRRA + frisch elasticity style labor disutility
    util, income, labor = construct_household_problem_functions(V, w, coefrra, frisch, labtax, labdisutil)

    # Setting up forward/backward difference grids for a. See helpers.jl
    daf, dab, azdelta = initialize_diff_grids(a, I, J)

    for ir = 1:Ir
        fill!(c,  complex(0.))
        fill!(h,  complex(1/3))
        fill!(h0, complex(1.))

        # Initial guess
        v = similar(h)

        for i in eachindex(h)
            inc  = income(h[i], zz[i], profshare[i], lumptransfer, r, aa[i]) # Get income
            v[i] = util(inc, h[i]) / ρ                                       # Value function guess
        end

        # Iterate HJB
        for ihjb = 1:maxit_HJB
            V = v
            Vaf, Vab, cf, hf, cb, hb = construct_initial_diff_matrices(V, Vaf, Vab,
                                                                       income, labor, h, h0, zz,
                                                                       profshare, lumptransfer,
                                                                       amax, amin,
                                                                       coefrra, r, daf,
                                                                       dab, maxhours)

            # Iterative method to find consistent forward/backward/neutral difference matrices for c and h
            cf, hf, cb, hb, c0, h0 = hours_iteration(income, labor, zz, profshare,
                                                     lumptransfer, aa,
                                                     coefrra, r, cf, hf, cb, hb, c0, h0,
                                                     maxhours, niter_hours)

            for i in eachindex(h0)
                c0[i] = income(h0[i], zz[i], profshare[i], lumptransfer, r, aa[i])
                sf[i] = income(hf[i], zz[i], profshare[i], lumptransfer, r, aa[i]) - cf[i]
                sb[i] = income(hb[i], zz[i], profshare[i], lumptransfer, r, aa[i]) - cb[i]
            end

            for j=1:J
                Vaf[I, j] = cf[I, j] ^ (-coefrra) # Forward difference for value function w.r.t. wealth
                Vab[1, j] = cb[1, j] ^ (-coefrra) # Backward difference for value function w.r.t. wealth
            end

            # See helpers.jl: this is a custom upwind function, not one from solve_hjb.jl
            V, A, u, h, c, s  = upwind(ρ, V, util, Aswitch, cf, cb, c0, hf, hb, h0, sf, sb, Vaf, Vab,
                                    daf, dab; Δ_HJB = Δ_HJB)

            # Check for convergence
            Vchange = V - v
            v = V

            err_HJB = maximum(abs.(Vchange))
            if err_HJB < tol_HJB
                break
            end
        end

        # Create initial guess for g0
        g0 = zeros(ComplexF64, I, J)
        # Assign stationary income distribution weight at a = 0, zero elsewhere
        g0[iszero.(a), :] = g_z
        # g_z is marginal distribution, so re-weight by some multiplier of Lebesgue measure
        g0 = g0 ./ reshape(azdelta, I, J)
        # Solve for distribution
        g = solve_kfe(A, g0, spdiagm(0 => azdelta),
                       maxit_kfe = maxit_kfe, tol_kfe = tol_kfe,
                       Δ_kfe = Δ_kfe)

        # Back out static conditions/steady state values given our value function and distribution
        N_ss, Y_ss, B_ss, profit_ss, profshare, lumptransfer, bond_err =
            calculate_ss_equil_vars(zz, h, g, azdelta, aa, labor_share_ss, meanlabeff,
                                lumptransferpc, govbondtarget)

        # Check bond market for market clearing
        r, rmin, rmax, ρ, ρmin, ρmax, clear_cond =
            check_bond_market_clearing(bond_err, crit_S, r, rmin, rmax, ρ, ρmin, ρmax,
                                       iterate_r, iterate_ρ)

        if clear_cond
            # Set steady state values
            m[:V_ss]           = real(vec(V))
            m[:inflation_ss]   = 0.
            m[:g_ss]           = real(g)
            m[:r_ss]           = r
            m[:u_ss]           = real(vec(u))
            m[:c_ss]           = real(vec(c))
            m[:h_ss]           = real(vec(h))
            m[:s_ss]           = real(vec(s))
            m[:rnom_ss]        = m[:r_ss].value + m[:inflation_ss].value
            m[:B_ss]           = sum(m[:g_ss].value .* vec(aa) .* azdelta)
            m[:N_ss]           = real(sum(vec(zz) .* m[:h_ss].value .* m[:g_ss].value .* azdelta))
            m[:Y_ss]           =  m[:N_ss]
            m[:labor_share_ss] = (ceselast - 1) / ceselast
            m[:w_ss]           =  m[:labor_share_ss].value
            m[:profit_ss]      = (1 - m[:labor_share_ss].value) .* m[:Y_ss].value
            m[:C_ss]           = sum(m[:c_ss].value .* m[:g_ss].value .* azdelta)
            m[:T_ss]           = real(lumptransfer)
            m[:ρ_ss]           = ρ
            m[:G_ss]           = labtax * m[:w_ss].value * m[:N_ss].value -
                                 m[:T_ss].value - m[:r_ss].value * m[:B_ss].value
            break
        end
    end
    return m
end

function model_settings!(m::OneAssetHANK)
    default_settings!(m)

    # Steady state r iteration
    m <= Setting(:r0, 0.005, "Initial guess for real interest rate")
    m <= Setting(:rmax, 0.08, "Maximum guess for real interest rate")
    m <= Setting(:rmin, 0.001, "Minimum guess for real interest rate")

    # Steady state ρ iteration
    m <= Setting(:ρ0, 0.02, "Initial guess for discount rate")
    m <= Setting(:ρmax, 0.05, "Maximum guess for discount rate")
    m <= Setting(:ρmin, 0.005, "Minimum guess for discount rate")

    # State space grid
    m <= Setting(:I, 100, "Number of grid points")
    m <= Setting(:J, 2, "Number of income states")
    m <= Setting(:amin, 0., "Minimum asset grid value")
    m <= Setting(:amax, 40., "Maximum asset grid value")
    m <= Setting(:agridparam, 1, "Bending coefficient: 1 for linear")
    m <= Setting(:a, construct_asset_grid(get_setting(m,:I), get_setting(m, :agridparam),
                        get_setting(m, :amin), get_setting(m,:amax)), "Asset grid")
    m <= Setting(:ygrid_combined, [0.2, 1])
    m <= Setting(:ymarkov_combined, [-0.5 0.5; 0.0376 -0.0376], "Markov transition parameters")
    m <= Setting(:g_z, compute_stationary_income_distribution(get_setting(m, :ymarkov_combined),
                       get_setting(m, :J)), "Stationary income distribution")
    m <= Setting(:zz, construct_labor_income_grid(get_setting(m, :ygrid_combined),
                       get_setting(m, :g_z), m[:meanlabeff].value, get_setting(m, :I)),
                       "Labor income grid repeated across asset dimension")
    m <= Setting(:z, get_setting(m, :zz)[1, :], "Labor income grid")

    # Number of variables
    m <= Setting(:n_jump_vars, get_setting(m, :I) * get_setting(m, :J) + 1, "Number of jump variables")
    m <= Setting(:n_state_vars, get_setting(m, :I) * get_setting(m, :J), "Number of state variables")
    m <= Setting(:n_state_vars_unreduce, 0, "Number of state variables not being reduced")
    # R: Inserted from old code; not certain necessary
    m <= Setting(:n_static_relations, 5, "Number of static relations: bond-market clearing, labor market
                                          clearing, consumption, output, total assets")
    m <= Setting(:n_vars, get_setting(m, :n_jump_vars) + get_setting(m, :n_state_vars)
                 + get_setting(m, :n_static_relations),
                 "Number of variables, total")

    # Steady state approximation
    m <= Setting(:Ir, 100, "Max number of iterations for steady state")
    m <= Setting(:maxit_HJB, 500, "Max number of iterations for HJB")
    m <= Setting(:tol_HJB, 1e-8, "Tolerance for HJB error")
    m <= Setting(:Δ_HJB, 1e6, "Multiplier for implicit upwind scheme of HJB")
    m <= Setting(:maxit_kfe, 1000, "Max number of iterations for kfe")
    m <= Setting(:tol_kfe, 1e-12, "Tolerance for kfe error")
    m <= Setting(:Δ_kfe, 1e6, "Multiplier for implicit upwind scheme of kfe")
    m <= Setting(:niter_hours, 10, "Max number of iterations for finding hours worked")
    m <= Setting(:iterate_r, false, "Iterate on real interest rate or not")
    m <= Setting(:iterate_ρ, true, "Iterate on discount rate or not")
    m <= Setting(:crit_S, 1e-5, "Tolerance for error in bond markets")

    # Reduction
    m <= Setting(:n_knots, 12, "Number of knot points")
    m <= Setting(:c_power, 1, "Amount of bending of knot point locations to make them nonuniform")
    m <= Setting(:n_post, length(get_setting(m, :zz)[1,:]),
                 "Number of dimensions that need to be approximated by spline basis")
    m <= Setting(:n_prior, 1,
                 "Number of dimensions approximated by spline basis that
                  were not used to compute the basis matrix")

    knots = collect(range(get_setting(m, :amin), stop=get_setting(m, :amax), length=get_setting(m, :n_knots)-1))
    knots = (get_setting(m, :amax) - get_setting(m, :amin)) /
            (2^get_setting(m, :c_power)-1) * ((knots .- get_setting(m, :amin))
                                              / (get_setting(m, :amax) -
                                                 get_setting(m, :amin)) .+ 1) .^ get_setting(m, :c_power) .+
                                                 get_setting(m, :amin) .- (get_setting(m, :amax) -
                                                                          get_setting(m, :amin)) /
                                                                          (2 ^ get_setting(m, :c_power) - 1)
    m <= Setting(:knots_dict, Dict(1 => knots),
                 "Location of knot points for each dimension for value function reduction")
    m <= Setting(:krylov_dim, 20, "Krylov reduction dimension")
    m <= Setting(:reduce_state_vars, true, "Reduce state variables or not") #turned off reduction
    m <= Setting(:reduce_v, true, "Reduce value function or not") #turned off reduction
    m <= Setting(:spline_grid, get_setting(m, :a), "Grid of knot points for spline basis")

    # Sampling method
    m <= Setting(:sampling_method, :SMC, "Set sampling method to SMC")

    # SMC Settings
    m <= Setting(:n_particles, 3000)
    m <= Setting(:n_Φ, 200)
    m <= Setting(:λ, 2.0)
    m <= Setting(:n_smc_blocks, 1)
    m <= Setting(:use_parallel_workers, true)
    m <= Setting(:step_size_smc, 0.5)
    m <= Setting(:n_mh_steps_smc, 1)
    m <= Setting(:resampler_smc, :polyalgo)
    m <= Setting(:target_accept, 0.25)
    m <= Setting(:mixture_proportion, .9)
    m <= Setting(:tempering_target, 0.95)
    m <= Setting(:resampling_threshold, .5)
    m <= Setting(:use_fixed_schedule, true)
    m <= Setting(:smc_iteration, 0)

    # Track lag
    m <= Setting(:track_lag, false, "Add first lag when constructing measurement equation")

    # Forecast
    m <= Setting(:use_population_forecast, true,
                 "Whether to use population forecasts as data")
    m <= Setting(:forecast_zlb_value, 0.13,
        "Value of the zero lower bound in forecast periods, if we choose to enforce it")

    # Simulating states
    m <= Setting(:state_simulation_freq, 2,
                 "How many states you want to simulate between states + 1")
end


text:
    application_mime_types = (

        "application/atom+xml", # Atom feeds
        "application/ecmascript",
        # ECMAScript/JavaScript; Defined in RFC 4329 (equivalent to application/javascript but with stricter processing rules)
        "application/EDI-X12", # EDI X12 data; Defined in RFC 1767
        "application/EDIFACT", # EDI EDIFACT data; Defined in RFC 1767
        "application/json", # JavaScript Object Notation JSON; Defined in RFC 4627
        "application/javascript", # ECMAScript/JavaScript; Defined in RFC 4329 (equivalent to application/ecmascript
        #   but with looser processing rules) It is not accepted in IE 8
        #   or earlier - text/javascript is accepted but it is defined as obsolete in RFC 4329.
        #   The "type" attribute of the <script> tag in HTML5 is optional and in practice
        #   omitting the media type of JavaScript programs is the most interoperable
        #   solution since all browsers have always assumed the correct
        #   default even before HTML5.                                                                                "application/octet-stream",  # Arbitrary binary data.[6] Generally speaking this type identifies files that are not associated with a specific application. Contrary to past assumptions by software packages such as Apache this is not a type that should be applied to unknown files. In such a case, a server or application should not indicate a content type, as it may be incorrect, but rather, should omit the type in order to allow the recipient to guess the type.[7]
        "application/ogg", # Ogg, a multimedia bitstream container format; Defined in RFC 5334
        "application/pdf", # Portable Document Format, PDF has been in use for document exchange
        #   on the Internet since 1993; Defined in RFC 3778
        "application/postscript", # PostScript; Defined in RFC 2046
        "application/rdf+xml", # Resource Description Framework; Defined by RFC 3870
        "application/rss+xml", # RSS feeds
        "application/soap+xml", # SOAP; Defined by RFC 3902
        "application/font-woff", # Web Open Font Format; (candidate recommendation; use application/x-font-woff
        #   until standard is official)
        "application/xhtml+xml", # XHTML; Defined by RFC 3236
        "application/xml-dtd", # DTD files; Defined by RFC 3023
        "application/xop+xml", # XOP
        "application/zip", # ZIP archive files; Registered[8]
        "application/gzip",         # Gzip, Defined in RFC 6713
    )

    audio_mime_types = (
        "audio/basic", # mulaw audio at 8 kHz, 1 channel; Defined in RFC 2046
        "audio/L24", # 24bit Linear PCM audio at 8-48 kHz, 1-N channels; Defined in RFC 3190
        "audio/mp4", # MP4 audio
        "audio/mpeg", # MP3 or other MPEG audio; Defined in RFC 3003
        "audio/ogg", # Ogg Vorbis, Speex, Flac and other audio; Defined in RFC 5334
        "audio/vorbis", # Vorbis encoded audio; Defined in RFC 5215
        "audio/vnd.rn-realaudio", # RealAudio; Documented in RealPlayer Help[9]
        "audio/vnd.wave", # WAV audio; Defined in RFC 2361
        "audio/webm",               # WebM open media format
    )

    image_mime_types = (
        "image/gif", # GIF image; Defined in RFC 2045 and RFC 2046
        "image/jpeg", # JPEG JFIF image; Defined in RFC 2045 and RFC 2046
        "image/pjpeg",
        # JPEG JFIF image; Associated with Internet Explorer; Listed in ms775147(v=vs.85) - Progressive JPEG, initiated before global browser support for progressive JPEGs (Microsoft and Firefox).
        "image/png", # Portable Network Graphics; Registered,[10] Defined in RFC 2083
        "image/svg+xml", # SVG vector image; Defined in SVG Tiny 1.2 Specification Appendix M
        "image/tiff", # Tag Image File Format (only for Baseline TIFF); Defined in RFC 3302
        "image/vnd.microsoft.icon",  # ICO image; Registered[11]
    )

    message_mime_types = (
        "message/http", # Defined in RFC 2616
        "message/imdn+xml", # IMDN Instant Message Disposition Notification; Defined in RFC 5438
        "message/partial", # Email; Defined in RFC 2045 and RFC 2046
        "message/rfc822",  # Email; EML files, MIME files, MHT files, MHTML files; Defined in RFC 2045 and RFC 2046
    )

    model_mime_types = (
        "model/example", # Defined in RFC 4735
        "model/iges", # IGS files, IGES files; Defined in RFC 2077
        "model/mesh", # MSH files, MESH files; Defined in RFC 2077, SILO files
        "model/vrml", # WRL files, VRML files; Defined in RFC 2077
        "model/x3d+binary", # X3D ISO standard for representing 3D computer graphics, X3DB binary files
        "model/x3d+vrml", # X3D ISO standard for representing 3D computer graphics, X3DV VRML files
        "model/x3d+xml",  # X3D ISO standard for representing 3D computer graphics, X3D XML files
    )

    multipart_mime_types = (
        "multipart/mixed", # MIME Email; Defined in RFC 2045 and RFC 2046
        "multipart/alternative", # MIME Email; Defined in RFC 2045 and RFC 2046
        "multipart/related", # MIME Email; Defined in RFC 2387 and used by MHTML (HTML mail)
        "multipart/form-data", # MIME Webform; Defined in RFC 2388
        "multipart/signed", # Defined in RFC 1847
        "multipart/encrypted",  # Defined in RFC 1847
    )

    text_mime_types = (
        "text/cmd", # commands; subtype resident in Gecko browsers like Firefox 3.5
        "text/css", # Cascading Style Sheets; Defined in RFC 2318
        "text/csv", # Comma-separated values; Defined in RFC 4180
        "text/html", # HTML; Defined in RFC 2854
        "text/javascript",
        # (Obsolete): JavaScript; Defined in and obsoleted by RFC 4329 in order to discourage its usage in favor of application/javascript. However, text/javascript is allowed in HTML 4 and 5 and, unlike application/javascript, has cross-browser support. The "type" attribute of the <script> tag in HTML5 is optional and there is no need to use it at all since all browsers have always assumed the correct default (even in HTML 4 where it was required by the specification).
        "text/plain", # Textual data; Defined in RFC 2046 and RFC 3676
        "text/vcard", # vCard (contact information); Defined in RFC 6350
        "text/xml",  # Extensible Markup Language; Defined in RFC 3023
    )

    video_mime_types = (
        "video/mpeg", # MPEG-1 video with multiplexed audio; Defined in RFC 2045 and RFC 2046
        "video/mp4", # MP4 video; Defined in RFC 4337
        "video/ogg", # Ogg Theora or other video (with audio); Defined in RFC 5334
        "video/quicktime", # QuickTime video; Registered[12]
        "video/webm", # WebM Matroska-based open media format
        "video/x-matroska", # Matroska open media format
        "video/x-ms-wmv", # Windows Media Video; Documented in Microsoft KB 288102
        "video/x-flv",  # Flash video (FLV files)
    )

    audio_file_extenstions = (
        "flac",
        "mp3",
        "wav",
    )

    image_file_extenstions = (
        "bmp",
        "gif",
        "jpeg",
        "jpg",
        "png",
        "tiff",
    )

    text_file_extensions = (
        "css",
        "csv",
        "html",
        "js",
        "json",
        "txt",
    )

    video_file_extensions = (
        "mp4",
        "avi",
        "mov",
        "webm",
    )

    function mime_type(category="None")

        mime_types = Dict(
        "application"=>application_mime_types[rand(1:length(application_mime_types))] ,
        "audio"=>audio_mime_types[rand(1:length(audio_mime_types))],
        "image"=>image_mime_types[rand(1:length(image_mime_types))],
        "message"=>message_mime_types[rand(1:length(message_mime_types))],
        "model"=>model_mime_types[rand(1:length(model_mime_types))],
        "multipart"=>multipart_mime_types[rand(1:length(multipart_mime_types))],
        "text"=>text_mime_types[rand(1:length(text_mime_types))],
        "video"=>video_mime_types[rand(1:length(video_mime_types))] )

        if category == "None"
        category= [key for key in  keys(mime_types)][rand(1:8)]
        end
        mime_types[category]
    end

        #param category: application|audio|image|message|model|multipart|text|video

    function file_extension( category="None")

       file_extensions = Dict(
        "audio"=>  audio_file_extenstions[rand(1:length(audio_file_extenstions))],
        "image"=>  image_file_extenstions[rand(1:length(image_file_extenstions))],
        "text"=>  text_file_extensions[rand(1:length(text_file_extensions))],
        "video"=>  video_file_extensions[rand(1:length(video_file_extensions))]
    )
        #:param category: audio|image|text|video
       if category == "None"
        category= [key for key in  keys(file_extensions)][rand(1:4)]
        end

     file_extensions[category]

    end

  function file_name( category="None", extencion="None")

        #:param category: audio|image|text|video
        #:param extension: file extension

    if extencion == "None"
     extencion = file_extension(category)
    end

    filename = lexify("????")

    "$(filename).$(extencion) "

    end


text:
ENV["PYCALL_JL_RUNTIME_PYTHON"] = Sys.which("python")

using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))

using ArgParse
using ProgressMeter
using Random
using RiskSensitiveSAC
import RobotOS.Time

################ COMMAND LINE ARGUMENTS ################
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "params_file_name"
            help = "parameter file name in default_params (without .jl)"
            arg_type = String
            required = true
        "test_case_file_name_1"
            help = "a test case file name in test_cases (without .jl)"
            arg_type = String
            required = true
        "test_case_file_name_2"
            help = "a test case file name in test_cases (without .jl)"
            arg_type = String
        "test_case_file_name_3"
            help = "a test case file name in test_cases (without .jl)"
            arg_type = String
        "--sigma_risk"
            help = "σ_risk value (risk sensitivity parameter)"
            arg_type = Float64
        "--alpha_col"
            help = "α_col value (collision cost peak parameter)"
            arg_type = Float64
        "--lambda_col"
            help = "λ_col value (collision cost bandwidth parameter)"
            arg_type = Float64
        "--num_samples"
            help = "number of prediction samples (not used in bic controller or oracle predictor)"
            arg_type = Int64
        "--nominal_search_only"
            help = "use nominal search results only instead of rssac"
            action = :store_true
        "--nominal_base_only"
            help = "use nominal base only instead of rssac"
            action = :store_true
        "--deterministic_prediction"
            help = "perform deterministic prediction"
            action = :store_true
        "--future_conditional"
            help = "condition prediction on future nominal control"
            action = :store_true
        "--use_bic"
            help = "use buffered input cell controller instead of rssac"
            action = :store_true
        "--use_crowd_nav"
            help = "use crowd nav controller instead of rssac"
            action = :store_true
        "--rng_seed"
            help = "rng seed for prediction (and simulation if in synthetic scene mode)"
            arg_type = Int64
        "--verbose"
            help = "verbose (true or false)"
            arg_type = Bool
            default = true
        "--show_log"
            help = "show log after evaluation is done"
            action = :store_true
        "--num_runs"
            help = "number of random simulation runs. default = 10"
            arg_type = Int64
            default = 10
        "--global_seed"
            help = "use global rng seed"
            action = :store_true
        "--randomize_goal"
            help = "perturb goal location with unit Gaussian"
            action = :store_true
    end
    return parse_args(s)
end

function main()
    parsed_args = parse_commandline()
    params_file_name = parsed_args["params_file_name"];
    test_case_file_name_1 = parsed_args["test_case_file_name_1"]
    test_case_file_name_2 = parsed_args["test_case_file_name_2"]
    test_case_file_name_3 = parsed_args["test_case_file_name_3"]
    sigma_risk = parsed_args["sigma_risk"]
    alpha_col = parsed_args["alpha_col"]
    lambda_col = parsed_args["lambda_col"]
    num_samples = parsed_args["num_samples"]
    nominal_search_only = parsed_args["nominal_search_only"]
    nominal_base_only = parsed_args["nominal_base_only"]
    deterministic_prediction = parsed_args["deterministic_prediction"]
    future_conditional = parsed_args["future_conditional"]
    use_bic = parsed_args["use_bic"]
    use_crowd_nav = parsed_args["use_crowd_nav"]
    verbose = parsed_args["verbose"]
    show_log = parsed_args["show_log"]
    randomize_goal = parsed_args["randomize_goal"]

    if parsed_args["global_seed"]
        rng = Random.GLOBAL_RNG
    else
        rng = MersenneTwister(593803);
    end
    seed_pool = Int[];

    @showprogress for run_id = 1:parsed_args["num_runs"]
        rng_seed = nothing;
        while true
            rng_seed_candidate = abs(rand(rng, Int32));
            if !in(rng_seed_candidate, seed_pool)
                rng_seed = rng_seed_candidate;
                push!(seed_pool, rng_seed);
                break;
            end
        end
        global newARGS = [string(params_file_name),
                          string(test_case_file_name_1)];
        if !isnothing(test_case_file_name_2)
            push!(newARGS, string(test_case_file_name_2))
        end
        if !isnothing(test_case_file_name_3)
            push!(newARGS, string(test_case_file_name_3))
        end
        push!(newARGS, "--run_id")
        push!(newARGS, string(run_id))
        push!(newARGS, "--rng_seed")
        push!(newARGS, string(rng_seed))
        if !isnothing(sigma_risk)
            push!(newARGS, "--sigma_risk")
            push!(newARGS, string(sigma_risk))
        end
        if !isnothing(alpha_col)
            push!(newARGS, "--alpha_col")
            push!(newARGS, string(alpha_col))
        end
        if !isnothing(lambda_col)
            push!(newARGS, "--lambda_col")
            push!(newARGS, string(lambda_col))
        end
        if !isnothing(num_samples)
            push!(newARGS, "--num_samples")
            push!(newARGS, string(num_samples))
        end
        if nominal_search_only
            push!(newARGS, "--nominal_search_only")
        end
        if nominal_base_only
            push!(newARGS, "--nominal_base_only")
        end
        if deterministic_prediction
            push!(newARGS, "--deterministic_prediction")
        end
        if future_conditional
            push!(newARGS, "--future_conditional")
        end
        if use_bic
            push!(newARGS, "--use_bic")
        end
        if use_crowd_nav
            push!(newARGS, "--use_crowd_nav")
        end
        if !verbose
            push!(newARGS, "--verbose")
            push!(newARGS, string(false))
        end
        if show_log
            push!(newARGS, "--show_log")
        end
        if randomize_goal
            push!(newARGS, "--randomize_goal")
        end
        try
            include("run_single_evaluation.jl")
        catch
            @warn "Run $(run_id) has errored. Skipping..."
            continue;
        end
        GC.gc();
    end
    # println(seed_pool)
end

main()


text:
# simulate MUSIC with a random array
using DirectionFinding

begin
    fc = 1e9
    λ = 299.79e6/fc
    a1 = 0.95; a2 = 6;
    N₀ = 0.01
    ra = randomarray(5, 1)
    s1 = SignalSource(a1, halfsine_train(1e6))
    s2 = SignalSource(a2, halfsine_train(1e6))
    sim = Simulation(ra, fc, N₀, s1, s2)
    mps = MUSIC(sim, 2.00001e6, 1000)
    ϕ = 0:0.01:2π-0.01
    peaks = findpeaks(mps, ϕ, n=2)
    println("Angles: $a1 -- $a2")
    println("Found: $(peaks[1]) -- $(peaks[2])")
    plot(ϕ, mps.(ϕ))
end


text:
function petk03()
#
# M. Petkovic testing polynomials, p134
#
    p1 = [1.0;-1-99*im/70];
    p2 = [1;2;3];
    p2 = conv(p2,p2);
    p = conv(p1,p2);
    p = conv(p,[1;1]);
    z = [-1.00000000000000 + 1.41421356237309im 2;
         -1.00000000000000 - 1.41421356237309im 2;
         -1-99*im/70 1;
         -1 1];
    p, PolyZeros(z)
end


text:
### Statistics

rank_inverse_transform(dist::Distribution, x) = quantile(dist,  (tiedrank(x) - 0.5) / length(x))

inverse_normal_transform(x::Array) = rank_inverse_transform(Normal(), x)

function inverse_normal_transform(x::DataArray)
	y  = copy(x)
	y[!isna(y)] = rank_inverse_transform(Normal(), dropna(x))
	return y
end

## FIXME: These functions needs cleanup!

function mean_imputation!(data::DataFrame,key::Symbol)
	mean_value = mean(map(Float32, x[find(!isna(x[:,key])),key]))
	data[key] = map(1:nrow(data)) do i
		if isna(data[i,key])
			value = group_mean_dict[data[i,groupby]]
			isnan(value) ? NA : value
		else
			data[i,key]
		end
	end
end

function group_mean_imputation!(data::DataFrame,key::Symbol, groupby::Symbol)
	group_mean = by(data, groupby, x->mean(map(Float32, x[find(!isna(x[:,key])),key])))
	group_mean_dict = Dict(zip(group_mean[groupby], group_mean[:x1]))
	data[key] = map(1:nrow(data)) do i
		if isna(data[i,key])
			value = group_mean_dict[data[i,groupby]]
			isnan(value) ? NA : value
		else
			data[i,key]
		end
	end
end


text:
sp = Specie(ρ=2200.,r=30.e-6,c=8433.,volfrac = 0.1)
sp = Specie(ρ=2200,r=30e-6,c=8433.,volfrac = 0.1)

ω = 30.0e6
# 2pi*2250.0/(ω*30e-6)

medium = Medium(ρ=2329.0,c = 2250.0)
kT = sqrt(wavenumber_low_volfrac(ω, medium, sp; verbose=true, tol=0.00002))

ωs = LinRange(0.01,1.0,60)*30.0e6
kTs = sqrt(wavenumber_low_volfrac(ωs, medium, sp))

speed_arr = ωs./real(kTs)
atten_arr = imag(kTs)

using Plots
pyplot()

p1 = plot(ωs./real(medium.c),speed_arr,  ylabel="wave speed", xlabel="k");
p2 = plot(ωs./real(medium.c), atten_arr, ylabel="attenuation", xlabel="k");
plot(p1,p2)

## An example where we vary the species

height=450
pyplot(linewidth=3, size=(2*height,height))
# pyplot(linewidth=3)
Plots.scalefontsizes(1.5)

# for fixed total volfrac fraction
medium = Medium(ρ=1.0,c = 1.0)
ωs = LinRange(0.01,1.0,60)
volfrac = 0.25
r1 = 0.5
r2 = 1.5
N=5
vols = LinRange(0.0,1.0,N)*volfrac

kTs_arr = [
  begin
    sp1 = Specie(0.0, r1; volfrac=vols[i])
    sp2 = Specie(Inf, r1; volfrac=volfrac-vols[i])
    [ sqrt(wavenumber_low_volfrac(ω, medium, [sp1,sp2])) for ω in ωs]
  end
for i = 1:N];

speed_arr = [ ωs./real(kTs) for kTs in kTs_arr]
atten_arr = imag(kTs_arr)

labs = reshape( map(v -> "void vol = $(Int(round(100*v)))%",vols),1, length(vols));
p1 = plot(ωs./real(medium.c), speed_arr, labels=labs, ylabel="wave speed" ,xlabel="k");
p2 = plot(ωs./real(medium.c), atten_arr, labels=labs, xlabel="k", ylabel="attenuation");

plot(p1,p2)
savefig("varying_voids.png")

Plots.scalefontsizes(1./1.5)


text:
using Plots
pyplot()
x = range(0, stop=2*pi, length=50)
println(x)
p = plot(x = x, y = sin(x), size =(900,300))

plot(x = x, y = sin(x))

plot(rand(100) / 3,reg=true,fill=(0,:green))
scatter!(rand(100),markersize=6,c=:orange)

plot(rand(100) / 3,reg=true,fill=(0,:green))
scatter!(rand(100),markersize=6,c=:orange)

plot(y = cumsum(randn(1000, 1)))


text:
function is_locally_isometric_kirschmer(L::QuadLat, M::QuadLat, p::NfOrdIdl)
  R = base_ring(L)
  base_ring(L) != base_ring(M) && throw(error("Lattices must have the same base ring"))
  order(p) != R && throw(error("Ideal must be in the base ring of the lattices"))
  d = rank(L)
  if d != rank(M)
    return false
  elseif d == 0
    return true
  end

  if minimum(p) != 2
    SL = _genus_symbol_kirschmer(L, p)
    SM = _genus_symbol_kirschmer(M, p, uniformizer = uniformizer(SL))
    return (data(SL) == data(SM))::Bool
  end

  if !is_rationally_isometric(L, M, p)
    return false
  end

  dimL1, sL1, wL1, aL1, fL1, G1 = data(_genus_symbol_kirschmer(L, p))
  dimL2, sL2, wL2, aL2, fL2, G2 = data(_genus_symbol_kirschmer(M, p))

  if (dimL1 != dimL2) || (sL1 != sL2) || (wL1 != wL2)
    return false
  end

  uL1 = [valuation(aL1[k], p) for k in 1:length(aL1)]
  uL2 = [valuation(aL2[k], p) for k in 1:length(aL2)]
  if uL1 != uL2
    return false
  end

  bL = [aL1[k]//aL2[k] for k in 1:length(aL1)];
  qL = [quadratic_defect(bL[k], p) for k in 1:length(bL)]

  if reduce((x, y) -> (x || y), qL[k] < wL1[k] - uL2[k] for k in 1:length(qL))
    return false
  end

  e = ramification_index(p)
  # GL, GM: Gram matrix
  d1 = [ diagonal_matrix([G1[i] for i in 1:t]) for t in 1:length(G1)]
  d2 = [ diagonal_matrix([G2[i] for i in 1:t]) for t in 1:length(G2)]

  for i in 1:length(fL1)
    detquot = det(d1[i])//det(d2[i])
    if valuation(detquot, p) != 0
      return false
    end
    if quadratic_defect(detquot, p) < fL1[i]
      return false
    end
    if (fL1[i] > 2*e + uL1[i+1] - wL2[i+1]) && !(_comp_hasse(d1[i], diagonal_matrix(d2[i], matrix(base_ring(d2[i]), 1, 1, [aL2[i+1]])), p))
      return false
    end
    if (fL1[i] > 2*e + uL1[i  ] - wL2[i  ]) && !(_comp_hasse(d1[i], diagonal_matrix(d2[i], matrix(base_ring(d2[i]), 1, 1, [aL2[i]])), p))
      return false
    end
  end

  return true
end

# nrows(G2) = nrows(G1) + 1
function _comp_hasse(G1, G2, p)
  G1o = diagonal(_gram_schmidt(G1, identity)[1])
  G2o = diagonal(_gram_schmidt(G2, identity)[1])
  push!(G1o, prod(G1o) * prod(G2o))
  return _hasse_invariant(G1o, p) == _hasse_invariant(G2o, p)
end