text:
# elastic-net regularization

export ElasticNet

"""
**Elastic-net regularization**

    ElasticNet(μ=1, λ=1)

Returns the function
```math
f(x) = μ\\|x\\|_1 + (λ/2)\\|x\\|^2,
```
for nonnegative parameters `μ` and `λ`.
"""
struct ElasticNet{R <: Real, S <: Real} <: ProximableFunction
    mu::R
    lambda::S
    function ElasticNet{R, S}(mu::R, lambda::S) where {R <: Real, S <: Real}
        if lambda < 0 || mu < 0
            error("parameters `μ` and `λ` must be nonnegative")
        else
            new(mu, lambda)
        end
    end
end

is_separable(f::ElasticNet) = true
is_prox_accurate(f::ElasticNet) = true
is_convex(f::ElasticNet) = true

ElasticNet(mu::R=1, lambda::S=1) where {R <: Real, S <: Real} = ElasticNet{R, S}(mu, lambda)

function (f::ElasticNet)(x::AbstractArray{T}) where {
    R <: Real, T <: RealOrComplex{R}
}
    return f.mu * norm(x, 1) + f.lambda / R(2) * norm(x, 2)^2
end

function prox!(y::AbstractArray{R}, f::ElasticNet, x::AbstractArray{R}, gamma::R=R(1)) where {
    R <: Real,
}
    sqnorm2x = R(0)
    norm1x = R(0)
    gm = gamma * f.mu
    gl = gamma * f.lambda
    for i in eachindex(x)
        y[i] = (x[i] + (x[i] <= -gm ? gm : (x[i] >= gm ? -gm : -x[i])))/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y::AbstractArray{R}, f::ElasticNet, x::AbstractArray{R}, gamma::AbstractArray{R}) where {
    R <: Real
}
    sqnorm2x = R(0)
    norm1x = R(0)
    for i in eachindex(x)
        gm = gamma[i] * f.mu
        gl = gamma[i] * f.lambda
        y[i] = (x[i] + (x[i] <= -gm ? gm : (x[i] >= gm ? -gm : -x[i])))/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y::AbstractArray{Complex{R}}, f::ElasticNet, x::AbstractArray{Complex{R}}, gamma::R=R(1)) where {
    R <: Real
}
    sqnorm2x = R(0)
    norm1x = R(0)
    gm = gamma * f.mu
    gl = gamma * f.lambda
    for i in eachindex(x)
        y[i] = sign(x[i]) * max(0, abs(x[i]) - gm)/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function prox!(y::AbstractArray{Complex{R}}, f::ElasticNet, x::AbstractArray{Complex{R}}, gamma::AbstractArray{R}) where {
    R <: Real
}
    sqnorm2x = R(0)
    norm1x = R(0)
    for i in eachindex(x)
        gm = gamma[i] * f.mu
        gl = gamma[i] * f.lambda
        y[i] = sign(x[i]) * max(0, abs(x[i]) - gm)/(1 + gl)
        sqnorm2x += abs2(y[i])
        norm1x += abs(y[i])
    end
    return f.mu * norm1x + f.lambda / R(2) * sqnorm2x
end

function gradient!(y::AbstractArray{T}, f::ElasticNet, x::AbstractArray{T}) where {
    R <: Real, T <: RealOrComplex{R}
}
    # Gradient of 1 norm
    y .= f.mu .* sign.(x)
    # Gradient of 2 norm
    y .+= f.lambda .* x
    return f.mu * norm(x, 1) + f.lambda / R(2) * norm(x, 2)^2
end

fun_name(f::ElasticNet) = "elastic-net regularization"
fun_dom(f::ElasticNet) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::ElasticNet) = "x ↦ μ||x||_1 + (λ/2)||x||²"
fun_params(f::ElasticNet) = "μ = $(f.mu), λ = $(f.lambda)"

function prox_naive(f::ElasticNet, x::AbstractArray{T}, gamma::R=R(1)) where {
    R <: Real, T <: RealOrComplex{R}
}
    uz = max.(0, abs.(x) .- gamma * f.mu)/(1 + f.lambda * gamma)
    return sign.(x) .* uz, f.mu * norm(uz, 1) + f.lambda / R(2) * norm(uz)^2
end

function prox_naive(f::ElasticNet, x::AbstractArray{T}, gamma::AbstractArray{R}) where {
    R <: Real, T <: RealOrComplex{R}
}
    uz = max.(0, abs.(x) .- gamma .* f.mu)./(1 .+ f.lambda .* gamma)
    return sign.(x) .* uz, f.mu * norm(uz, 1) + f.lambda / R(2) * norm(uz)^2
end


text:
using Test

using PyPlot
using LaTeXStrings

using Oceananigans.Advection

# Define a few utilities for running tests and unpacking and plotting results
include("ConvergenceTests/ConvergenceTests.jl")

using .ConvergenceTests
using .ConvergenceTests.OneDimensionalGaussianAdvectionDiffusion: run_test
using .ConvergenceTests.OneDimensionalUtils: unpack_errors, defaultcolors, removespines

""" Run advection test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions, advection_scheme)

    # Determine safe time-step
           Lx = 2.5
    stop_time = 1e-3
            h = Lx / maximum(resolutions)
           Δt = min(1e-2 * h / U, 1e-3 * h^2 / κ)

    # Run the tests
    results = [run_test(Nx=Nx, Δt=Δt, advection=advection_scheme, stop_iteration=1,
                        U=U, κ=κ) for Nx in resolutions]

    return results
end

#####
##### Run test
#####

advection_schemes = (CenteredSecondOrder(), CenteredFourthOrder(), UpwindBiasedThirdOrder(), WENO5())

U = -3
κ = 1e-8
Nx = 2 .^ (7:10)

results = Dict()
for scheme in advection_schemes
    t_scheme = typeof(scheme)
    results[t_scheme] = run_convergence_test(κ, U, Nx, scheme)
end

rate_of_convergence(::CenteredSecondOrder) = 2
rate_of_convergence(::CenteredFourthOrder) = 4
rate_of_convergence(::UpwindBiasedThirdOrder) = 3
rate_of_convergence(::WENO5) = 5
rate_of_convergence(::WENO{K}) where K = 2K-1

fig, ax = subplots()

@testset "tmp" begin
for (j, scheme) in enumerate(advection_schemes)
    t_scheme = typeof(scheme)
    name = string(t_scheme)
    roc = rate_of_convergence(scheme)

    u_L₁, v_L₁, cx_L₁, cy_L₁, u_L∞, v_L∞, cx_L∞, cy_L∞ = unpack_errors(results[typeof(scheme)])

    atol = t_scheme == CenteredSecondOrder ? 0.05 :
           t_scheme == CenteredFourthOrder ? 0.12 :
           t_scheme == WENO5               ? 0.40 : Inf

    test_rate_of_convergence(u_L₁,  Nx, expected=-roc, atol=atol, name=name*" u_L₁")
    test_rate_of_convergence(v_L₁,  Nx, expected=-roc, atol=atol, name=name*" v_L₁")
    test_rate_of_convergence(cx_L₁, Nx, expected=-roc, atol=atol, name=name*" cx_L₁")
    test_rate_of_convergence(cy_L₁, Nx, expected=-roc, atol=atol, name=name*" cy_L₁")
    test_rate_of_convergence(u_L∞,  Nx, expected=-roc, atol=atol, name=name*" u_L∞")
    test_rate_of_convergence(v_L∞,  Nx, expected=-roc, atol=atol, name=name*" v_L∞")
    test_rate_of_convergence(cx_L∞, Nx, expected=-roc, atol=atol, name=name*" cx_L∞")
    test_rate_of_convergence(cy_L∞, Nx, expected=-roc, atol=atol, name=name*" cy_L∞")

    @test  u_L₁ ≈  v_L₁
    @test cx_L₁ ≈ cy_L₁
    @test  u_L∞ ≈  v_L∞
    @test cx_L∞ ≈ cy_L∞

    common_kwargs = (linestyle="None", color=defaultcolors[j], mfc="None", alpha=0.8)
    loglog(Nx,  u_L₁; basex=2, marker="o", label="\$L_1\$-norm, \$u\$ $name", common_kwargs...)
    loglog(Nx,  v_L₁; basex=2, marker="2", label="\$L_1\$-norm, \$v\$ $name", common_kwargs...)
    loglog(Nx, cx_L₁; basex=2, marker="*", label="\$L_1\$-norm, \$x\$ tracer $name", common_kwargs...)
    loglog(Nx, cy_L₁; basex=2, marker="+", label="\$L_1\$-norm, \$y\$ tracer $name", common_kwargs...)

    loglog(Nx,  u_L∞; basex=2, marker="1", label="\$L_\\infty\$-norm, \$u\$ $name", common_kwargs...)
    loglog(Nx,  v_L∞; basex=2, marker="_", label="\$L_\\infty\$-norm, \$v\$ $name", common_kwargs...)
    loglog(Nx, cx_L∞; basex=2, marker="^", label="\$L_\\infty\$-norm, \$x\$ tracer $name", common_kwargs...)
    loglog(Nx, cy_L∞; basex=2, marker="s", label="\$L_\\infty\$-norm, \$y\$ tracer $name", common_kwargs...)

    label = raw"\sim N_x^{-" * "$roc" * raw"}" |> latexstring
    loglog(Nx, cx_L₁[1] .* (Nx[1] ./ Nx) .^ roc, color=defaultcolors[j], basex=2, alpha=0.8, label=label)

    xlabel(L"N_x")
    ylabel("\$L\$-norms of \$ | c_\\mathrm{sim} - c_\\mathrm{analytical} |\$")
    removespines("top", "right")
    lgd = legend(loc="upper right", bbox_to_anchor=(1.4, 1.0), prop=Dict(:size=>6))

    if j == length(advection_schemes)
        filepath = joinpath(@__DIR__, "figs", "one_dimensional_convergence.png")
        savefig(filepath, dpi=480, bbox_extra_artists=(lgd,), bbox_inches="tight")
        close(fig)
    end
end
end

ext: jl


text:
abstract type ImprovementAquisition{S} <: SingleObjAquisition{S} end

function update_parameters!(a::ImprovementAquisition{Min}, krg::Kriging{N,1}) where {N}
    yo = krg.yscaler \ [minimum(krg.gps[1].y)]
    return a.τ = yo[1]
end

function update_parameters!(a::ImprovementAquisition{Max}, krg::Kriging{N,1}) where {N}
    yo = krg.yscaler \ [maximum(krg.gps[1].y)]
    return a.τ = yo[1]
end

improvement(a::ImprovementAquisition{Min}, y) = a.τ - y

improvement(a::ImprovementAquisition{Max}, y) = y - a.τ

mutable struct ProbabilityOfImprovement{S} <: ImprovementAquisition{S}
    τ::Float64
end

ProbabilityOfImprovement(s::Symbol) = ProbabilityOfImprovement(sym2sense(s))
ProbabilityOfImprovement(s::EGOSense) = ProbabilityOfImprovement{s}(0.0)

function objective(a::ProbabilityOfImprovement{Min}, model::Kriging{N,1}, x) where {N}
    m, v = predict_full(model, x)
    μ = m[1]
    σ² = v[1]
    if σ² == 0
        float(μ < a.τ)
    else
        Φ((a.τ - μ) / σ²)
    end
end

function objective(a::ProbabilityOfImprovement{Max}, model::Kriging{N,1}, x) where {N}
    m, v = predict_full(model, x)
    μ = m[1]
    σ² = v[1]
    if σ² == 0
        float(μ > a.τ)
    else
        Φ((μ - a.τ) / σ²)
    end
end
mutable struct ExpectedImprovement{S} <: ImprovementAquisition{S}
    τ::Float64
end

ExpectedImprovement(s::Symbol) = ExpectedImprovement(sym2sense(s))
ExpectedImprovement(s::EGOSense) = ExpectedImprovement{s}(0.0)

function objective(a::ExpectedImprovement{Min}, model::Kriging{N,1}, x) where {N}
    m, v = predict_full(model, x)
    μ = m[1]
    σ² = v[1]
    δ = a.τ - μ
    if σ² == 0
        max(δ, 0.0)
    else
        σ = √σ²
        n = δ / σ
        δ * Φ(n) + σ * φ(n)
    end
end

function objective(a::ExpectedImprovement{Max}, model::Kriging{N,1}, x) where {N}
    m, v = predict_full(model, x)
    μ = m[1]
    σ² = v[1]
    δ = μ - a.τ
    if σ² == 0
        max(δ, 0.0)
    else
        σ = √σ²
        n = δ / σ
        δ * Φ(n) + σ * φ(n)
    end
end


text:
using DSGE, Test, JLD2

# What do you want to do?
check_steady_state = true
check_jacobian     = true
check_solution     = true
check_irfs         = true

path = dirname(@__FILE__)

het = HetDSGEGovDebt(testing_gamma = true, ref_dir = HETDSGEGOVDEBT)
het <= Setting(:steady_state_only, true)
steadystate!(het)

m = RepDSGEGovDebt(het)

# Steady-state computation
if check_steady_state
    steadystate!(m, het)

    file = jldopen("$path/../../heterogeneous/het_dsge_gov_debt/reference/steady_state.jld2", "r")
    saved_Rk   = read(file, "Rk")
    saved_ω    = read(file, "omega")
    saved_kl   = read(file, "kl")
    saved_k    = read(file, "k")
    saved_y    = read(file, "y")
    saved_T    = read(file, "T")

    close(file)

    file = jldopen("$path/reference/steady_state.jld2", "r")
    saved_c    = read(file, "c")
    saved_ell  = read(file, "l")
    saved_β    = read(file, "beta")
    close(file)

    @testset "Check steady state outputs" begin
        @test saved_Rk ≈ m[:Rkstar].value
        @test saved_ω  ≈ m[:ωstar].value
        @test saved_kl ≈ m[:klstar].value
        @test saved_k  ≈ m[:kstar].value
        @test saved_y  ≈ m[:ystar].value
        @test saved_T  ≈ m[:Tstar].value
        @test saved_β  ≈ m[:βstar].value
        @test isapprox(saved_c,   m[:cstar].value; atol = 0.02)
        @test isapprox(saved_ell, m[:lstar].value; atol = 0.01)
    end
end

# Jacobian computation
if check_jacobian
    steadystate!(m, het)
    JJ = DSGE.jacobian(m)

    file = jldopen("$path/reference/jacobian.jld2", "r")
    saved_JJ = read(file, "JJ")

    nvars = read(file, "nvars")

    KKP  = read(file, "KKP")
    LRRP = read(file, "LRRP")
    LIIP = read(file, "LIIP")
    LYP = read(file, "LYP")
    LWP = read(file, "LWP")
    LXP = read(file, "LXP")
    BGP = read(file, "BGP")
    BP = read(file, "BP")
    GP = read(file, "GP")
    ZP = read(file, "ZP")
    MUP = read(file, "MUP")
    LAMWP = read(file, "LAMWP")
    LAMFP = read(file, "LAMFP")
    MONP = read(file, "MONP")
    RRP = read(file, "RRP")
    IIP = read(file, "IIP")
    TTP = read(file, "TTP")
    WP = read(file, "WP")
    HHP = read(file, "HHP")
    PIP = read(file, "PIP")
    PIWP = read(file, "PIWP")
    LAMP = read(file, "LAMP")
    YP = read(file, "YP")
    XP = read(file, "XP")
    MCP = read(file, "MCP")
    QP = read(file, "QP")
    RKP = read(file, "RKP")
    ELLP = read(file, "ELLP")
    TGP = read(file, "TGP")

    KK = KKP + nvars
    LRR = LRRP + nvars
    LII = LIIP + nvars
    LY = LYP + nvars
    LW = LWP + nvars
    LX = LXP + nvars
    BG = BGP + nvars
    B = BP + nvars
    G = GP + nvars
    Z = ZP + nvars
    MU = MUP + nvars
    LAMW = LAMWP + nvars
    LAMF = LAMFP + nvars
    MON = MONP + nvars
    RR = RRP + nvars
    II = IIP + nvars
    TT = TTP + nvars
    W = WP + nvars
    HH = HHP + nvars
    PI = PIP + nvars
    PIW = PIWP + nvars
    LAM = LAMP + nvars
    Y = YP + nvars
    X = XP + nvars
    MC = MCP + nvars
    Q = QP + nvars
    RK = RKP + nvars
    ELL = ELLP + nvars
    TG = TGP + nvars

    F1 = read(file, "F1")
    F5 = read(file, "F5")
    F6 = read(file, "F6")
    F7 = read(file, "F7")
    F8 = read(file, "F8")
    F9 = read(file, "F9")
    F10 = read(file, "F10")
    F11 = read(file, "F11")
    F12 = read(file, "F12")
    F13 = read(file, "F13")
    F14 = read(file, "F14")
    F15 = read(file, "F15")
    F16 = read(file, "F16")
    F17 = read(file, "F17")
    F18 = read(file, "F18")
    F40 = read(file, "F40")
    F41 = read(file, "F41")
    F19 = read(file, "F19")
    F20 = read(file, "F20")
    F24 = read(file, "F24")
    F31 = read(file, "F31")
    F32 = read(file, "F32")
    F33 = read(file, "F33")
    F34 = read(file, "F34")
    F35 = read(file, "F35")
    F36 = read(file, "F36")
    F37 = read(file, "F37")
    F38 = read(file, "F38")
    F39 = read(file, "F39")
    close(file)

    @testset "Euler equation" begin
        @test saved_JJ[F1,ELLP] ≈ JJ[F1, ELLP]
        @test saved_JJ[F1,ZP]   ≈ JJ[F1, ZP]
        @test isapprox(saved_JJ[F1,B], JJ[F1, B]; atol=0.1)
        @test saved_JJ[F1,ELL]  ≈ JJ[F1, ELL]
        @test isapprox(saved_JJ[F1,RR], JJ[F1, RR]; atol=0.1)
    end

    @testset "mkt clearing" begin
        @test saved_JJ[F5,Y]   ≈ JJ[F5, Y]
        @test saved_JJ[F5,G]   ≈ JJ[F5, G]
        @test saved_JJ[F5,X]   ≈ JJ[F5, X]
        @test isapprox(saved_JJ[F5,ELL], JJ[F5, ELL]; atol = 0.1)
    end

    @testset "lambda = average marginal utility" begin
        @test isapprox(saved_JJ[F6,LAM], JJ[F6, LAM]; atol = 0.1)
        @test isapprox(saved_JJ[F6,ELL], JJ[F6, ELL]; atol = 0.1)
    end

    @testset "transfer" begin
        @test saved_JJ[F7,TT] ≈ JJ[F7,TT]
        @test saved_JJ[F7,RK] ≈ JJ[F7,RK]
        @test saved_JJ[F7,KK] ≈ JJ[F7,KK]
        @test saved_JJ[F7,Z]  ≈ JJ[F7,Z]
        @test saved_JJ[F7,X]  ≈ JJ[F7,X]
        @test saved_JJ[F7,MC] ≈ JJ[F7,MC]
        @test saved_JJ[F7,TG] ≈ JJ[F7,TG]
    end

    @testset "investment" begin
        @test saved_JJ[F8,Q]  ≈ JJ[F8,Q]
        @test saved_JJ[F8,MU] ≈ JJ[F8,MU]
        @test saved_JJ[F8,XP] ≈ JJ[F8,XP]
        @test saved_JJ[F8,ZP] ≈ JJ[F8,ZP]
        @test saved_JJ[F8,X]  ≈ JJ[F8,X]
        @test saved_JJ[F8,LX] ≈ JJ[F8,LX]
        @test saved_JJ[F8,Z]  ≈ JJ[F8,Z]
    end

    @testset "tobin's q" begin
        @test saved_JJ[F9,LAM]  ≈ JJ[F9,LAM]
        @test saved_JJ[F9,LAMP] ≈ JJ[F9,LAMP]
        @test saved_JJ[F9,Q]    ≈ JJ[F9,Q]
        @test saved_JJ[F9,ZP]   ≈ JJ[F9,ZP]
        @test saved_JJ[F9,RKP]  ≈ JJ[F9,RKP]
        @test saved_JJ[F9,QP]   ≈ JJ[F9,QP]
    end

    @testset "capital accumulation" begin
        @test saved_JJ[F10,KKP] ≈ JJ[F10,KKP]
        @test saved_JJ[F10,KK]  ≈ JJ[F10,KK]
        @test saved_JJ[F10,Z]   ≈ JJ[F10,Z]
        @test saved_JJ[F10,MU]  ≈ JJ[F10,MU]
        @test saved_JJ[F10,X]   ≈ JJ[F10,X]
    end

    @testset "wage phillips curve" begin
        @test saved_JJ[F11,PIW]  ≈ JJ[F11,PIW]
        @test saved_JJ[F11,LAMW] ≈ JJ[F11,LAMW]
        @test saved_JJ[F11,HH]   ≈ JJ[F11,HH]
        @test saved_JJ[F11,W]    ≈ JJ[F11,W]
        @test saved_JJ[F11,PIWP] ≈ JJ[F11,PIWP]
        @test saved_JJ[F11,LAM]  ≈ JJ[F11,LAM]
    end

    @testset "price phillips curve" begin
        @test saved_JJ[F12,PI]   ≈ JJ[F12,PI]
        @test saved_JJ[F12,MC]   ≈ JJ[F12,MC]
        @test saved_JJ[F12,LAMF] ≈ JJ[F12,LAMF]
        @test saved_JJ[F12,PIP]  ≈ JJ[F12,PIP]
    end

    @testset "marginal cost" begin
        @test saved_JJ[F13,MC] ≈ JJ[F13,MC]
        @test saved_JJ[F13,W]  ≈ JJ[F13,W]
        @test saved_JJ[F13,RK] ≈ JJ[F13,RK]
    end

    @testset "gdp" begin
        @test saved_JJ[F14,Y]  ≈ JJ[F14,Y]
        @test saved_JJ[F14,Z]  ≈ JJ[F14,Z]
        @test saved_JJ[F14,KK] ≈ JJ[F14,KK]
        @test saved_JJ[F14,HH] ≈ JJ[F14,HH]
    end

    @testset "optimal k/l ratio" begin
        @test saved_JJ[F15,RK] ≈ JJ[F15,RK]
        @test saved_JJ[F15,W]  ≈ JJ[F15,W]
        @test saved_JJ[F15,HH] ≈ JJ[F15,HH]
        @test saved_JJ[F15,KK] ≈ JJ[F15,KK]
        @test saved_JJ[F15,Z]  ≈ JJ[F15,Z]
    end

    @testset "taylor rule" begin
        @test saved_JJ[F16,II]  ≈ JJ[F16,II]
        @test saved_JJ[F16,LII] ≈ JJ[F16,LII]
        @test saved_JJ[F16,PI]  ≈ JJ[F16,PI]
        @test saved_JJ[F16,Y]   ≈ JJ[F16,Y]
        @test saved_JJ[F16,LY]  ≈ JJ[F16,LY]
        @test saved_JJ[F16,Z]   ≈ JJ[F16,Z]
        @test saved_JJ[F16,MON] ≈ JJ[F16,MON]
    end

    @testset "fisher eqn" begin
        @test saved_JJ[F17,RR]  ≈ JJ[F17,RR]
        @test saved_JJ[F17,PIP] ≈ JJ[F17,PIP]
        @test saved_JJ[F17,II]  ≈ JJ[F17,II]
    end

    @testset "wage inflation" begin
        @test saved_JJ[F18,PIW] ≈ JJ[F18,PIW]
        @test saved_JJ[F18,PI]  ≈ JJ[F18,PI]
        @test saved_JJ[F18,Z]   ≈ JJ[F18,Z]
        @test saved_JJ[F18,W]   ≈ JJ[F18,W]
        @test saved_JJ[F18,LW]  ≈ JJ[F18,LW]
    end

    @testset "fiscal rule" begin
        @test saved_JJ[F40,TG]  ≈ JJ[F40,TG]
        @test saved_JJ[F40,RR]  ≈ JJ[F40,RR]
        @test saved_JJ[F40,BG]  ≈ JJ[F40,BG]
        @test saved_JJ[F40,Z]   ≈ JJ[F40,Z]
        @test saved_JJ[F40,Y]   ≈ JJ[F40,Y]
        @test saved_JJ[F40,G]   ≈ JJ[F40,G]
    end

    @testset "govt budget constraint" begin
        @test saved_JJ[F41,BGP] ≈ JJ[F41,BGP]
        @test saved_JJ[F41,RR]  ≈ JJ[F41,RR]
        @test saved_JJ[F41,BG]  ≈ JJ[F41,BG]
        @test saved_JJ[F41,Z]   ≈ JJ[F41,Z]
        @test saved_JJ[F41,Y]   ≈ JJ[F41,Y]
        @test saved_JJ[F41,G]   ≈ JJ[F41,G]
        @test saved_JJ[F41,TG]  ≈ JJ[F41,TG]
    end

    @testset "update lagged variables" begin
        @test saved_JJ[F19,LRRP] ≈ JJ[F19,LRRP]
        @test saved_JJ[F19,RR]   ≈ JJ[F19,RR]

        @test saved_JJ[F20,LIIP] ≈ JJ[F20,LIIP]
        @test saved_JJ[F20,II]   ≈ JJ[F20,II]

        @test saved_JJ[F24,LYP]  ≈ JJ[F24,LYP]
        @test saved_JJ[F24,Y]    ≈ JJ[F24,Y]

        @test saved_JJ[F31,LWP] ≈ JJ[F31,LWP]
        @test saved_JJ[F31,W]   ≈ JJ[F31,W]

        @test saved_JJ[F32,LXP] ≈ JJ[F32,LXP]
        @test saved_JJ[F32,X]   ≈ JJ[F32,X]
    end

    @testset "discount factor shock" begin
        @test saved_JJ[F33,BP] ≈ JJ[F33,BP]
        @test saved_JJ[F33,B]  ≈ JJ[F33,B]
    end

    @testset "g/y shock" begin
        @test saved_JJ[F34,GP] ≈ JJ[F34,GP]
        @test saved_JJ[F34,G]  ≈ JJ[F34,G]
    end

    @testset "tfp growth shock" begin
        @test saved_JJ[F35,ZP] ≈ JJ[F35,ZP]
        @test saved_JJ[F35,Z]  ≈ JJ[F35,Z]
    end

    @testset "investment shock" begin
        @test saved_JJ[F36,MUP] ≈ JJ[F36,MUP]
        @test saved_JJ[F36,MU]  ≈ JJ[F36,MU]
    end

    @testset "wage mkup shock" begin
        @test saved_JJ[F37,LAMWP] ≈ JJ[F37,LAMWP]
        @test saved_JJ[F37,LAMW]  ≈ JJ[F37,LAMW]
    end

    @testset "price mkup shock" begin
        @test saved_JJ[F38,LAMFP] ≈ JJ[F38,LAMFP]
        @test saved_JJ[F38,LAMF]  ≈ JJ[F38,LAMF]
    end

    @testset "monetary policy shock" begin
        @test saved_JJ[F39,MONP] ≈ JJ[F39,MONP]
        @test saved_JJ[F39,MON]  ≈ JJ[F39,MON]
    end
end

# Solution
if check_solution
    steadystate!(m, het)
    gx, hx, _ = klein(m)

    file = jldopen("$path/reference/solve.jld2", "r")
    saved_gx   = read(file, "gx")
    saved_hx   = read(file, "hx")
    close(file)

    @show sum(abs.(saved_gx - gx))
    @show sum(abs.(saved_hx - hx))

    @testset "Check solve outputs" begin
        @test saved_gx  ≈ gx
        @test saved_hx  ≈ hx
    end
end

# IRFs
if check_irfs
    endo = m.endogenous_states
    endo_aug = m.endogenous_states_augmented
    exo  = m.exogenous_shocks
    system = compute_system(m)
    states, obs, pseudo = impulse_responses(m, system, flip_shocks = true)

    file = jldopen("$path/reference/irfs.jld2", "r")
    IRFxp = read(file, "IRFxp")
    IRFC  = read(file, "IRFC")
    IRFyp = read(file, "IRFyp")

    nkp  = read(file, "nkp")
    ZP   = read(file, "ZP")
    ELLP = read(file, "ELLP")
    RRP  = read(file, "RRP")
    IIP  = read(file, "IIP")
    WP   = read(file, "WP")
    PIP  = read(file, "PIP")
    TTP  = read(file, "TTP")
    YP   = read(file, "YP")
    PIWP = read(file, "PIWP")
    XP   = read(file, "XP")
    KKP  = read(file, "KKP")
    QP   = read(file, "QP")
    RKP  = read(file, "RKP")
    HHP  = read(file, "HHP")
    TGP  = read(file, "TGP")
    BGP  = read(file, "BGP")
    MCP  = read(file, "MCP")
    LAMP = read(file, "LAMP")
    close(file)

    @testset "Check IRFs" begin
        @test states[endo[:z′_t], :, exo[:z_sh]] ≈ IRFxp[ZP, :]
        @test states[endo[:l′_t], :, exo[:z_sh]] ≈ IRFyp[ELLP - nkp, :]
        @test states[endo[:R′_t], :, exo[:z_sh]] ≈ IRFyp[RRP  - nkp, :]
        @test states[endo[:i′_t], :, exo[:z_sh]] ≈ IRFyp[IIP  - nkp, :]
        @test states[endo[:w′_t], :, exo[:z_sh]] ≈ IRFyp[WP   - nkp, :]
        @test states[endo[:π′_t], :, exo[:z_sh]] ≈ IRFyp[PIP  - nkp, :]
        @test states[endo[:t′_t], :, exo[:z_sh]] ≈ IRFyp[TTP  - nkp, :]
        @test states[endo[:y′_t], :, exo[:z_sh]] ≈ IRFyp[YP   - nkp, :]
        @test states[endo[:π_w′_t], :, exo[:z_sh]] ≈ IRFyp[PIWP - nkp, :]
        @test states[endo[:I′_t], :, exo[:z_sh]] ≈ IRFyp[XP   - nkp, :]
        @test states[endo[:k′_t], :, exo[:z_sh]] ≈ IRFxp[KKP, :]
        @test states[endo[:Q′_t], :, exo[:z_sh]] ≈ IRFyp[QP   - nkp, :]
        @test states[endo[:capreturn′_t], :, exo[:z_sh]] ≈ IRFyp[RKP - nkp, :]
        @test states[endo[:L′_t], :, exo[:z_sh]] ≈ IRFyp[HHP  - nkp, :]
        @test states[endo[:tg′_t], :, exo[:z_sh]] ≈ IRFyp[TGP - nkp, :]
        @test isapprox(states[endo[:bg′_t], :, exo[:z_sh]], IRFxp[BGP, :], atol = 1e-6)

        # Consumption IRF tests
        # contemporaneous consumption is the same as -l′_t
        @test states[endo[:l′_t], :, exo[:z_sh]] ≈ -IRFC
        @test states[endo_aug[:c_t1], :, exo[:z_sh]][2:end] ≈ IRFC[1:end-1]

        @test states[endo[:mc′_t], :, exo[:z_sh]] ≈ IRFyp[MCP - nkp, :]
        @test states[endo[:margutil′_t], :, exo[:z_sh]] ≈ IRFyp[LAMP - nkp, :]
    end
end


text:
# Simple kernel for matrix multiplication
@kernel function matmul_kernel!(a, b, c)
    i, j = @index(Global, NTuple)

    # creating a temporary sum variable for matrix multiplication
    tmp_sum = zero(eltype(c))
    for k = 1:size(a)[2]
        @inbounds tmp_sum += a[i,k] * b[k, j]
    end

    c[i,j] = tmp_sum
end

function matmul_testsuite(backend, ArrayT)

    matmul = matmul_kernel!(backend(), (32, 32))
    a = ArrayT(rand(128, 256))
    b = ArrayT(rand(256, 128))
    c = ArrayT(zeros(128, 128))
    wait(matmul(a, b, c, ndrange=size(c)))

    @test c ≈ a*b

    dmatmul = Enzyme.autodiff(matmul)
    da = similar(a)
    da .= 0
    db = similar(b)
    db .= 0
    dc = similar(c)
    dc .= 1
    c .= 0

    compare_dc = copy(dc)
    wait(dmatmul(
        Duplicated(a, da),
        Duplicated(b, db),
        Duplicated(c, dc), ndrange=size(c)))

    @test da ≈ compare_dc * b'
    @test db ≈ a' * compare_dc
end


text:
module HolyTraitsDispatch

using BinaryTraits
using BinaryTraits.Prefix: Is, Not
using Test

# Traits used in this test
import Base: getindex
@trait Indexable
@implement Is{Indexable} by getindex(_, v::Integer)
@assign AbstractArray with Is{Indexable}
@assign AbstractString with Is{Indexable}

import Base: collect
@trait Collectable
@implement Is{Collectable} by collect(_)
@assign Base.Generator with Is{Collectable}

function test()

    @testset "General dispatch" begin
        @traitfn head(v::Is{Indexable}) = :index     # v[1]
        @traitfn head(v::Not{Indexable}) = :notindex # first(v[1])
        @test head([1,2,3]) == :index
        @test head("abcde") == :index
        @test head(i for i in 4:6) == :notindex
    end

    @testset "Keyword args handling" begin
        @traitfn function increment(x::Is{Indexable}, i::Int; by = 1)
            x[i] += by
        end
        @test increment([1,2,3], 2; by = 2) == 4
    end

    @testset "Where-clause handling" begin
        @traitfn function addfirst(y::Vector{T}, x::Is{Indexable}) where {T <: AbstractFloat}
            return y .+ x[1]
        end
        @test addfirst([1.0, 2.0, 3.0], [1, 2, 3]) == [2.0, 3.0, 4.0]
    end

    @testset "Unnamed args" begin
        @traitfn unamed(v::Is{Indexable}, ::Int) = :int       # unamed regular position arg
        @traitfn unamed(::Is{Indexable}, v::String) = :string # unamed trait arg
        @traitfn unamed(::Is{Indexable}, ::Bool) = :bool      # unamed all!
        @test unamed([1,2,3], 1) == :int
        @test unamed([1,2,3], "abc") == :string
        @test unamed([1,2,3], true) == :bool
    end

    @testset "Duck typed arg" begin
        @traitfn duckarg(::Is{Indexable}, v) = :quack
        @test duckarg([1,2,3], true) == :quack
    end

    @testset "Multi-trait dispatch" begin
        # How do we use more than one trait for the same argument?
        # There are 2^n cases (n = number of traits).
        # However, it can be simplified using `BinaryTrait{T}` when cases overlap.
        @traitfn seek(v::Is{Indexable},  w::BinaryTrait{Collectable}, i::Integer)  = :index
        @traitfn seek(v::Not{Indexable}, w::Is{Collectable}, i::Integer) = :collect
        @traitfn seek(v::Not{Indexable}, w::Not{Collectable}, i::Integer) = :error
        seek(v, i::Integer) = seek(v, v, i)   # write a custom dispatcher that duplicates the arg

        @test seek([1,2,3], 2) == :index
        @test seek((i for i in 4:6), 2) == :collect
        @test seek(123, 1) == :error
    end
end

end # module

using Test
@testset "Holy Traits dispatch" begin
    import .HolyTraitsDispatch
    HolyTraitsDispatch.test()
end


text:
using Combinatorics: combinations

test = """
    [[[0,[5,8]],[[1,7],[9,6]]],[[4,[1,2]],[[1,4],2]]]
    [[[5,[2,8]],4],[5,[[9,9],0]]]
    [6,[[[6,2],[5,6]],[[7,6],[4,7]]]]
    [[[6,[0,7]],[0,9]],[4,[9,[9,0]]]]
    [[[7,[6,4]],[3,[1,3]]],[[[5,5],1],9]]
    [[6,[[7,3],[3,2]]],[[[3,8],[5,7]],4]]
    [[[[5,4],[7,7]],8],[[8,3],8]]
    [[9,3],[[9,9],[6,[4,9]]]]
    [[2,[[7,7],7]],[[5,8],[[9,3],[0,2]]]]
    [[[[5,2],5],[8,[3,7]]],[[5,[7,5]],[4,4]]]
    """

function parse_input(input::AbstractString)
    lines = split(input, '\n'; keepempty = false)
    return eval(Expr(:vect, Meta.parse.(lines)...))
end

abstract type LinkedTree{T} end

LinkedTree(v::AbstractVector) = LinkedNode(LinkedTree.(v)...)
LinkedTree(v) = LinkedLeaf(v)

mutable struct LinkedLeaf{T} <: LinkedTree{T}
    v::T
    prev::Ref
    next::Ref

    LinkedLeaf(v::T) where {T} = new{T}(v, Ref(nothing), Ref(nothing))
end

mutable struct LinkedNode{T} <: LinkedTree{T}
    v::T                        # grabage can
    left::LinkedTree{T}
    right::LinkedTree{T}
    first::Ref
    last::Ref

    function LinkedNode(left::LinkedTree{T}, right::LinkedTree{T}) where {T}
        t = new{T}(zero(T), left, right, first(left), last(right))
        first(left)[].prev = Ref(t)
        last(left)[].next = first(right)
        first(right)[].prev = last(left)
        last(right)[].next = Ref(t)
        return t
    end
end

Base.first(leaf::LinkedLeaf) = Ref(leaf)
Base.first(node::LinkedNode) = node.first
Base.last(leaf::LinkedLeaf) = Ref(leaf)
Base.last(node::LinkedNode) = node.last

Base.show(io::IO, leaf::LinkedLeaf) = print(io, leaf.v)
function Base.show(io::IO, node::LinkedNode)
    print(io, '[')
    print(io, node.left)
    print(io, ',')
    print(io, node.right)
    print(io, ']')
end

as_list(leaf::LinkedLeaf) = leaf.v
as_list(node::LinkedNode) = [as_list(node.left), as_list(node.right)]

magnitude(v::AbstractVector) = magnitude(LinkedTree(v))
magnitude(leaf::LinkedLeaf) = leaf.v
magnitude(node::LinkedNode) = 3magnitude(node.left) + 2magnitude(node.right)

function explode!(v::AbstractVector; ok::Bool = true)
    (t, ok_) = explode!(LinkedTree(v))
    @assert ok_ == ok
    return as_list(t)
end

explode!(tree::LinkedTree) = explode!(tree, 0)
explode!(leaf::LinkedLeaf, depth) = (leaf, false)
function explode!(node::LinkedNode, depth)
    if depth == 4
        node.left.prev[].v += node.left.v
        node.right.next[].v += node.right.v
        return (LinkedLeaf(0), true)
    end
    (left, ok) = explode!(node.left, depth + 1)
    ok && return (LinkedNode(left, node.right), true)
    (right, ok) = explode!(node.right, depth + 1)
    ok && return (LinkedNode(node.left, right), true)
    return (node, false)
end

function split!(leaf::LinkedLeaf)
    i = leaf.v
    return i >= 10 ? (LinkedTree([i ÷ 2, i - i ÷ 2]), true) : (leaf, false)
end

function split!(node::LinkedNode)
    (left, ok) = split!(node.left)
    ok && return (LinkedNode(left, node.right), true)
    (right, ok) = split!(node.right)
    ok && return (LinkedNode(node.left, right), true)
    return (node, false)
end

function reduce!(tree::LinkedTree)
    while true
        (tree, ok) = explode!(tree)
        ok && continue
        (tree, ok) = split!(tree)
        !ok && break
    end
    return tree
end

add(left, right) = as_list(add!(LinkedTree(left), LinkedTree(right)))
add!(left::LinkedTree, right::LinkedTree) = reduce!(LinkedNode(left, right))

add_all(v) = as_list(foldl(add!, LinkedTree.(v)))

q1(v) = magnitude(foldl(add!, LinkedTree.(v)))

function q2(v)
    pairs = combinations(v, 2)
    return map(pairs) do (a, b)
        m1 = add!(LinkedTree(a), LinkedTree(b)) |> magnitude
        m2 = add!(LinkedTree(b), LinkedTree(a)) |> magnitude
        return max(m1, m2)
    end |> findmax |> first
end

@assert magnitude([[9,1], [1,9]]) == 129

@assert explode!([[[[[9,8],1],2],3],4]) == [[[[0,9],2],3],4]
@assert explode!([7,[6,[5,[4,[3,2]]]]]) == [7,[6,[5,[7,0]]]]
@assert explode!([[6,[5,[4,[3,2]]]],1]) == [[6,[5,[7,0]]],3]
@assert explode!([[3,[2,[1,[7,3]]]],[6,[5,[4,[3,2]]]]]) ==
    [[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]
@assert explode!([[3,[2,[8,0]]],[9,[5,[4,[3,2]]]]]) ==
    [[3,[2,[8,0]]],[9,[5,[7,0]]]]

@assert add([[[[4,3],4],4],[7,[[8,4],9]]], [1,1]) ==
    [[[[0,7],4],[[7,8],[6,0]]],[8,1]]
@assert add([[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]],
            [7,[[[3,7],[4,3]],[[6,3],[8,8]]]]) ==
                [[[[4,0],[5,4]],[[7,7],[6,0]]],[[8,[7,7]],[[7,9],[5,0]]]]

@assert add_all([[i,i] for i in 1:4]) == [[[[1,1],[2,2]],[3,3]],[4,4]]
@assert add_all([[i,i] for i in 1:5]) == [[[[3,0],[5,3]],[4,4]],[5,5]]
@assert add_all([[i,i] for i in 1:6]) == [[[[5,0],[7,4]],[5,5]],[6,6]]
@assert add_all([[[[0,[4,5]],[0,0]],[[[4,5],[2,6]],[9,5]]],
                    [7,[[[3,7],[4,3]],[[6,3],[8,8]]]],
                    [[2,[[0,8],[3,4]]],[[[6,7],1],[7,[1,6]]]],
                    [[[[2,4],7],[6,[0,5]]],[[[6,8],[2,8]],[[2,1],[4,5]]]],
                    [7,[5,[[3,8],[1,4]]]],
                    [[2,[2,2]],[8,[8,1]]],
                    [2,9],
                    [1,[[[9,3],9],[[9,0],[0,7]]]],
                    [[[5,[7,4]],7],1],
                    [[[[4,2],2],6],[8,7]]]) ==
                        [[[[8,7],[7,7]],[[8,6],[7,7]]],[[[0,7],[6,6]],[8,7]]]

let v = parse_input(test)
    @assert add_all(v) ==
        [[[[6,6],[7,6]],[[7,7],[7,0]]],[[[7,7],[7,7]],[[7,8],[9,9]]]]
    @assert q1(v) == 4140
    @assert q2(v) == 3993
end

v = parse_input(read("day18.in", String))
println("Q1: ", q1(v))
println("Q2: ", q2(v))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
"""
    CoordinateSystem(; length=1.0)

Generate a new visual shape representing a coordinate system.
The reference frame = Object3D frame coincides with the coordinate system.

# Arguments
- `length` defines the length of the coordinate system axes.

# Notes
- CoordinateSystem is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct CoordinateSystem <: Modia3D.AbstractShape
    length::Float64

    function CoordinateSystem(; length=1.0)
        @assert(length >= 0.0)
        new(length)
    end
end

"""
    Grid(; axis=3, length=1.0, width=0.6, distance=0.1, lineWidth=1.0)

Generate a new visual shape representing a planar grid.
The reference frame = Object3D frame is located at the center of the grid.

# Arguments
- `axis` defines the normal axis of the grid: 1=x; 2=y; 3=z.
- `length` defines the length of the grid.
- `width` defines the width of the grid.
- `distance` defines the distance of the grid lines.
- `lineWidth` defines the line width of the grid.

The dimension directions depend on `axis` by circularly shift:

| `axis` | `length` | `width` |
|:------:|:--------:|:-------:|
|   1    |   y      |   z     |
|   2    |   z      |   x     |
|   3    |   x      |   y     |

# Notes
- Grid is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct Grid <: Modia3D.AbstractShape
    axis::Int
    length::Float64
    width::Float64
    distance::Float64
    lineWidth::Float64

    function Grid(; axis=3, length=1.0, width=0.6, distance=0.1, lineWidth=1.0)
        @assert(1 <= axis <= 3)
        @assert(length >= 0.0)
        @assert(width >= 0.0)
        @assert(0.0 < distance <= min(length, width))
        @assert(lineWidth > 0.0)
        new(axis, length, width, distance, lineWidth)
    end
end

"""
    Spring(; axis=3, length=1.0, diameter=1.0, wireDiameter=0.05, windings=5)

Generate a new visual shape representing a helical spiral spring.
The reference frame = Object3D frame is located at the center of the base circle of the spring.

# Arguments
- `axis` defines the axis of the spring: 1=x; 2=y; 3=z.
- `length` defines the length of the spring.
- `diameter` defines the diameter of the spring.
- `wireDiameter` defines the diameter of the spring wire.
- `windings` defines the number of windings of the spring.

# Notes
- SimVis visualizes whole-numbered values of `windings`.
- Spring is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct Spring <: Modia3D.AbstractShape
    axis::Int
    length::Float64
    diameter::Float64
    wireDiameter::Float64
    windings::Float64

    function Spring(; axis=3, length=1.0, diameter=1.0, wireDiameter=0.05, windings=5)
        @assert(1 <= axis <= 3)
        @assert(length >= 0.0)
        @assert(diameter >= 0.0)
        @assert(0.0 < wireDiameter < min(diameter, length))
        new(axis, length, diameter, wireDiameter, windings)
    end
end

"""
    GearWheel(; axis=3, diameter=1.0, length=1.0, innerDiameter=0.0, angle=0.0, teeth=20)

Generate a new visual shape representing a gearwheel.
The reference frame = Object3D frame is located at the center of the gearwheel.

# Arguments
- `axis` defines the rotation axis of the gearwheel: 1=x; 2=y; 3=z.
- `diameter` defines the diameter of the gearwheel.
- `length` defines the length of the gearwheel.
- `innerDiameter` defines the inner diameter of the gearwheel (where `innerDiameter=0` defines a full gearwheel).
- `angle` defines the bevel angle of the gearwheel (where `angle=0` defines a cylindrical gearwheel).
- `teeth` defines the number of teeth of the gearwheel.

# Notes
- GearWheel is [not supported by animation export](https://github.com/ModiaSim/PrivateModia3D.jl/issues/77).
"""
mutable struct GearWheel <: Modia3D.AbstractShape
    axis::Int
    diameter::Float64
    length::Float64
    innerDiameter::Float64
    angle::Float64
    teeth::Int

    function GearWheel(; axis=3, diameter=1.0, length=1.0, innerDiameter=0.0, angle=0.0, teeth=20)
        @assert(1 <= axis <= 3)
        @assert(diameter >= 0.0)
        @assert(length >= 0.0)
        @assert(0.0 <= innerDiameter < diameter)
        @assert(isinteger(teeth))
        @assert(teeth > 0)
        new(axis, diameter, length, innerDiameter, angle, teeth)
    end
end

getShapeKind(shape::Sphere)    = SphereKind
getShapeKind(shape::Ellipsoid) = EllipsoidKind
getShapeKind(shape::Box)       = BoxKind
getShapeKind(shape::Cylinder)  = CylinderKind
getShapeKind(shape::Cone)      = ConeKind
getShapeKind(shape::Capsule)   = CapsuleKind
getShapeKind(shape::Beam)      = BeamKind
getShapeKind(shape::FileMesh)  = FileMeshKind
getShapeKind(shape::Spring)    = SpringKind
getShapeKind(shape::CoordinateSystem) = CoordinateSystemKind
getShapeKind(shape::Grid)      = GridKind
getShapeKind(shape::GearWheel) = GearWheelKind
getShapeKind(shape::TextShape) = TextKind
getShapeKind(shape)            = UndefinedShapeKind


text:
using Base.Meta: isexpr

# A parser for a simple human-readable of an MOI model.
# This should be thought of as a compact way to write small models
# for tests, and not an exchange format.
#
# variables: x, y, z
# minobjective: 2x + 3y
# con1: x + y <= 1
# con2: [x,y] in Set
#
# special labels: variables, minobjective, maxobjective
# everything else denotes a constraint with a name
# all constraints must be named
# "x - y" does NOT currently parse, needs to be written as "x + -1.0*y"
# "x^2" does NOT currently parse, needs to be written as "x*x"

struct ParsedScalarAffineTerm
    coefficient::Float64
    variable_index::Symbol
end

struct ParsedScalarAffineFunction
    terms::Vector{ParsedScalarAffineTerm}
    constant::Float64
end

struct ParsedVectorAffineTerm
    output_index::Int64
    scalar_term::ParsedScalarAffineTerm
end

struct ParsedVectorAffineFunction
    terms::Vector{ParsedVectorAffineTerm}
    constant::Vector{Float64}
end

struct ParsedScalarQuadraticTerm
    coefficient::Float64
    variable_index_1::Symbol
    variable_index_2::Symbol
end

struct ParsedScalarQuadraticFunction
    affine_terms::Vector{ParsedScalarAffineTerm}
    quadratic_terms::Vector{ParsedScalarQuadraticTerm}
    constant::Float64
end

struct ParsedVectorQuadraticTerm
    output_index::Int64
    scalar_term::ParsedScalarQuadraticTerm
end

struct ParsedVectorQuadraticFunction
    affine_terms::Vector{ParsedVectorAffineTerm}
    quadratic_terms::Vector{ParsedVectorQuadraticTerm}
    constant::Vector{Float64}
end

struct ParsedSingleVariable
    variable::Symbol
end

struct ParsedVectorOfVariables
    variables::Vector{Symbol}
end

# Not written with any considerations for performance
function parsefunction(ex)
    if isa(ex, Symbol)
        return ParsedSingleVariable(ex)
    elseif isexpr(ex, :vect)
        if all(s -> isa(s, Symbol), ex.args)
            return ParsedVectorOfVariables(copy(ex.args))
        else
            singlefunctions = parsefunction.(ex.args)
            affine_terms = ParsedVectorAffineTerm[]
            quadratic_terms = ParsedVectorQuadraticTerm[]
            constant = Float64[]
            for (outindex, f) in enumerate(singlefunctions)
                if isa(f, ParsedSingleVariable)
                    push!(
                        affine_terms,
                        ParsedVectorAffineTerm(
                            outindex,
                            ParsedScalarAffineTerm(1.0, f.variable),
                        ),
                    )
                    push!(constant, 0.0)
                elseif isa(f, ParsedScalarAffineFunction)
                    append!(
                        affine_terms,
                        ParsedVectorAffineTerm.(outindex, f.terms),
                    )
                    push!(constant, f.constant)
                else
                    @assert isa(f, ParsedScalarQuadraticFunction)
                    append!(
                        affine_terms,
                        ParsedVectorAffineTerm.(outindex, f.affine_terms),
                    )
                    append!(
                        quadratic_terms,
                        ParsedVectorQuadraticTerm.(outindex, f.quadratic_terms),
                    )
                    push!(constant, f.constant)
                end
            end
            if length(quadratic_terms) == 0
                return ParsedVectorAffineFunction(affine_terms, constant)
            else
                return ParsedVectorQuadraticFunction(
                    affine_terms,
                    quadratic_terms,
                    constant,
                )
            end
        end
    else
        # only accept Expr(:call, :+, ...), no recursive expressions
        # TODO: generalize. x - y + z would be useful
        if isexpr(ex, :call) && ex.args[1] == :*
            # handle 2x as (+)(2x)
            ex = Expr(:call, :+, ex)
        end
        if ex isa Number
            ex = Expr(:call, :+, ex)
        end
        @assert isexpr(ex, :call)
        if ex.args[1] != :+
            error("Expected `+`, got `$(ex.args[1])`.")
        end
        affine_terms = ParsedScalarAffineTerm[]
        quadratic_terms = ParsedScalarQuadraticTerm[]
        constant = 0.0
        for subex in ex.args[2:end]
            if isexpr(subex, :call) && subex.args[1] == :*
                if length(subex.args) == 3
                    # constant * variable
                    @assert isa(subex.args[2], Number)
                    @assert isa(subex.args[3], Symbol)
                    push!(
                        affine_terms,
                        ParsedScalarAffineTerm(subex.args[2], subex.args[3]),
                    )
                else
                    # constant * variable * variable for quadratic
                    @assert length(subex.args) == 4 "Multiplication with more than three terms not supported"
                    @assert isa(subex.args[2], Number)
                    @assert isa(subex.args[3], Symbol)
                    @assert isa(subex.args[4], Symbol)
                    if subex.args[3] == subex.args[4]
                        coefficient = 2 * subex.args[2]
                    else
                        coefficient = subex.args[2]
                    end
                    push!(
                        quadratic_terms,
                        ParsedScalarQuadraticTerm(
                            coefficient,
                            subex.args[3],
                            subex.args[4],
                        ),
                    )
                end
            elseif isa(subex, Symbol)
                push!(affine_terms, ParsedScalarAffineTerm(1.0, subex))
            else
                @assert isa(subex, Number)
                constant += subex
            end
        end
        if length(quadratic_terms) == 0
            return ParsedScalarAffineFunction(affine_terms, constant)
        else
            return ParsedScalarQuadraticFunction(
                affine_terms,
                quadratic_terms,
                constant,
            )
        end
    end
end

# see tests for examples
function separatelabel(ex)
    if isexpr(ex, :call) && ex.args[1] == :(:)
        return ex.args[2], ex.args[3]
    elseif isexpr(ex, :tuple)
        ex = copy(ex)
        @assert isexpr(ex.args[1], :call) && ex.args[1].args[1] == :(:)
        label = ex.args[1].args[2]
        ex.args[1] = ex.args[1].args[3]
        return label, ex
    elseif isexpr(ex, :call)
        ex = copy(ex)
        @assert isexpr(ex.args[2], :call) && ex.args[2].args[1] == :(:)
        label = ex.args[2].args[2]
        ex.args[2] = ex.args[2].args[3]
        return label, ex
    else
        error("Unrecognized expression $ex")
    end
end

function parsedtoMOI(model, s::Symbol)
    index = MOI.get(model, MOI.VariableIndex, String(s))
    if index === nothing
        error("Invalid variable name $s.")
    end
    return index
end

# Used for Vector{Symbol}, Vector{ParsedScalarAffineTerm}, Vector{ParsedVectorAffineTerm},
# Vector{ParsedScalarQuadraticTerm} and Vector{ParsedVectorQuadraticTerm}
parsedtoMOI(model, s::Vector) = parsedtoMOI.(model, s)

parsedtoMOI(model, s::Union{Float64,Int64}) = s

for typename in [
    :ParsedScalarAffineTerm,
    :ParsedScalarAffineFunction,
    :ParsedVectorAffineTerm,
    :ParsedVectorAffineFunction,
    :ParsedScalarQuadraticTerm,
    :ParsedScalarQuadraticFunction,
    :ParsedVectorQuadraticTerm,
    :ParsedVectorQuadraticFunction,
    :ParsedSingleVariable,
    :ParsedVectorOfVariables,
]
    moiname = Meta.parse(replace(string(typename), "Parsed" => "MOI."))
    fields = fieldnames(eval(typename))
    constructor = Expr(
        :call,
        moiname,
        [
            Expr(
                :call,
                :parsedtoMOI,
                :model,
                Expr(:., :f, Base.Meta.quot(field)),
            ) for field in fields
        ]...,
    )
    @eval parsedtoMOI(model, f::$typename) = $constructor
end

function loadfromstring!(model, s)
    parsedlines = filter(ex -> ex != nothing, Meta.parse.(split(s, "\n")))

    for line in parsedlines
        label, ex = separatelabel(line)
        if label == :variables
            if isexpr(ex, :tuple)
                for v in ex.args
                    vindex = MOI.add_variable(model)
                    MOI.set(model, MOI.VariableName(), vindex, String(v))
                end
            else
                @assert isa(ex, Symbol)
                vindex = MOI.add_variable(model)
                MOI.set(model, MOI.VariableName(), vindex, String(ex))
            end
        elseif label == :maxobjective
            f = parsedtoMOI(model, parsefunction(ex))
            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
            MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        elseif label == :minobjective
            f = parsedtoMOI(model, parsefunction(ex))
            MOI.set(model, MOI.ObjectiveFunction{typeof(f)}(), f)
            MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
        else
            # constraint
            @assert isexpr(ex, :call)
            f = parsedtoMOI(model, parsefunction(ex.args[2]))
            if ex.args[1] == :in
                # Could be safer here
                set = Core.eval(MOI, ex.args[3])
            elseif ex.args[1] == :<=
                set = MOI.LessThan(ex.args[3])
            elseif ex.args[1] == :>=
                set = MOI.GreaterThan(ex.args[3])
            elseif ex.args[1] == :(==)
                set = MOI.EqualTo(ex.args[3])
            else
                error("Unrecognized expression $ex")
            end
            F, S = typeof(f), typeof(set)
            if !MOI.supports_constraint(model, F, S)
                throw(MOI.UnsupportedConstraint{F,S}())
            end
            cindex = MOI.add_constraint(model, f, set)
            MOI.set(model, MOI.ConstraintName(), cindex, String(label))
        end
    end
end


text:
import Base.rand

export Ellipsoid,
       shape_matrix

"""
    Ellipsoid{N<:AbstractFloat, VN<:AbstractVector{N},
              MN<:AbstractMatrix{N}} <: AbstractCentrallySymmetric{N}

Type that represents an ellipsoid.

It is defined as the set

```math
E = \\left\\{ x ∈ \\mathbb{R}^n : (x-c)^T Q^{-1} (x-c) ≤ 1 \\right\\},
```
where ``c \\in \\mathbb{R}^n`` is its *center* and ``Q \\in \\mathbb{R}^{n×n}``
its *shape matrix*, which should be a positive definite matrix.
An ellipsoid can also be characterized as the image of a Euclidean ball by an
invertible linear transformation. It is the higher-dimensional generalization
of an ellipse.

### Fields

- `center`       -- center of the ellipsoid
- `shape_matrix` -- real positive definite matrix, i.e. it is equal to its transpose
                    and ``x^\\mathrm{T}Qx > 0`` for all nonzero ``x``

## Notes

By default, the inner constructor checks that the given shape matrix is positive
definite. Use the flag `check_posdef=false` to disable this check.

### Examples

An ellipsoid can be created passing its center and shape matrix (which should be
positive definite).

For instance, we create a two-dimensional ellipsoid with center `[1, 1]`:

```jldoctest ellipsoid_constructor
julia> using LinearAlgebra

julia> E = Ellipsoid(ones(2), Diagonal([2.0, 0.5]))
Ellipsoid{Float64,Array{Float64,1},Diagonal{Float64,Array{Float64,1}}}([1.0, 1.0], [2.0 0.0; 0.0 0.5])
```

If the center is not specified, it is assumed that it is the origin. For instance,
a three-dimensional ellipsoid centered at the origin and the shape matrix being
the identity can be created with:

```jldoctest ellipsoid_constructor
julia> E = Ellipsoid(Matrix(1.0I, 3, 3))
Ellipsoid{Float64,Array{Float64,1},Array{Float64,2}}([0.0, 0.0, 0.0], [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])

julia> dim(E)
3
```
The center and shape matrix of the ellipsoid can be retrieved with the functions
`center` and `shape_matrix` respectively:

```jldoctest ellipsoid_constructor
julia> center(E)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> shape_matrix(E)
3×3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0
```

The function `an_element` returns some element of the ellipsoid:

```jldoctest ellipsoid_constructor
julia> an_element(E)
3-element Array{Float64,1}:
 0.0
 0.0
 0.0

julia> an_element(E) ∈ E
true
```

We can evaluate its support vector in a given direction, say `[1, 1, 1]`:

```jldoctest ellipsoid_constructor
julia> σ(ones(3), E)
3-element Array{Float64,1}:
 0.5773502691896258
 0.5773502691896258
 0.5773502691896258
```
"""
struct Ellipsoid{N<:AbstractFloat, VN<:AbstractVector{N},
                 MN<:AbstractMatrix{N}} <: AbstractCentrallySymmetric{N}
    center::VN
    shape_matrix::MN

    # default constructor with dimension check
    function Ellipsoid(c::VN, Q::MN;
                       check_posdef::Bool=true) where {N<:AbstractFloat,
                                                       VN<:AbstractVector{N},
                                                       MN<:AbstractMatrix{N}}

        @assert length(c) == checksquare(Q) "the length of the center and the size "
            "of the shape matrix do not match; they are $(length(c)) and $(size(Q)) respectively"

        if check_posdef
            isposdef(Q) || throw(ArgumentError("an ellipsoid's shape matrix " *
                                               "must be positive definite"))
        end
        return new{N, VN, MN}(c, Q)
    end
end

isoperationtype(::Type{<:Ellipsoid}) = false
isconvextype(::Type{<:Ellipsoid}) = true

# convenience constructor for ellipsoid centered in the origin
function Ellipsoid(Q::AbstractMatrix{N}; check_posdef::Bool=true) where {N<:AbstractFloat}
    # TODO: use similar vector type for the center, see #2032
    return Ellipsoid(zeros(N, size(Q, 1)), Q; check_posdef=check_posdef)
end

# --- AbstractCentrallySymmetric interface functions ---

"""
    center(E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the center of the ellipsoid.

### Input

- `E` -- ellipsoid

### Output

The center of the ellipsoid.
"""
function center(E::Ellipsoid{N}) where {N<:AbstractFloat}
    return E.center
end

"""
    shape_matrix(E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the shape matrix of the ellipsoid.

### Input

- `E` -- ellipsoid

### Output

The shape matrix of the ellipsoid.
"""
function shape_matrix(E::Ellipsoid{N}) where {N<:AbstractFloat}
    return E.shape_matrix
end

# --- LazySet interface functions ---

"""
    σ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the support vector of an ellipsoid in a given direction.

### Input

- `d` -- direction
- `E` -- ellipsoid

### Output

Support vector in the given direction.

### Algorithm

Let ``E`` be an ellipsoid of center ``c`` and shape matrix ``Q = BB^\\mathrm{T}``.
Its support vector along direction ``d`` can be deduced from that of the unit
Euclidean ball ``\\mathcal{B}_2`` using the algebraic relations for the support
vector,

```math
σ_{B\\mathcal{B}_2 ⊕ c}(d) = c + Bσ_{\\mathcal{B}_2} (B^\\mathrm{T} d)
= c + \\dfrac{Qd}{\\sqrt{d^\\mathrm{T}Q d}}.
```
"""
function σ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}
    if norm(d, 2) == zero(N)
        return an_element(E)
    end
    α = sqrt(dot(d, E.shape_matrix * d))
    return E.center .+ E.shape_matrix * d ./ α
end

"""
    ρ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}

Return the support function of an ellipsoid in a given direction.

### Input

- `d` -- direction
- `E` -- ellipsoid

### Output

The support function of the ellipsoid in the given direction.

### Algorithm

The support value is ``cᵀ d + ‖Bᵀ d‖₂`` where ``c`` is the center and
``Q = B Bᵀ`` is the shape matrix of `E`.
"""
function ρ(d::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}
    return dot(center(E), d) + sqrt(inner(d, E.shape_matrix, d))
end

"""
    ∈(x::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}

Check whether a given point is contained in an ellipsoid.

### Input

- `x` -- point/vector
- `E` -- ellipsoid

### Output

`true` iff `x ∈ E`.

### Algorithm

The point ``x`` belongs to the ellipsoid of center ``c`` and shape matrix ``Q``
if and only if

```math
(x-c)^\\mathrm{T} Q^{-1} (x-c) ≤ 1.
```
"""
function ∈(x::AbstractVector{N}, E::Ellipsoid{N}) where {N<:AbstractFloat}
    @assert length(x) == dim(E)
    w, Q = x-E.center, E.shape_matrix
    return dot(w, Q \ w) ≤ 1
end

"""
    rand(::Type{Ellipsoid}; [N]::Type{<:AbstractFloat}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random ellipsoid.

### Input

- `Ellipsoid` -- type for dispatch
- `N`         -- (optional, default: `Float64`) numeric type
- `dim`       -- (optional, default: 2) dimension
- `rng`       -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed`      -- (optional, default: `nothing`) seed for reseeding

### Output

A random ellipsoid.

### Algorithm

The center is a normally distributed vector with entries of mean 0 and standard
deviation 1.

The idea for the shape matrix comes from
[here](https://math.stackexchange.com/a/358092).
The matrix is symmetric positive definite, but also diagonally dominant.

```math
Q =  \\frac{1}{2}(S + S^T) + nI,
```
where ``n`` = `dim` (defaults to 2), and ``S`` is a ``n \\times n`` random
matrix whose coefficients are uniformly distributed in the interval ``[-1, 1]``.
"""
function rand(::Type{Ellipsoid};
              N::Type{<:AbstractFloat}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    rng = reseed(rng, seed)
    center = randn(rng, N, dim)
    # random entries in [-1, 1]
    # this needs a bit of code because 'rand' only samples from [0, 1]
    shape_matrix = Matrix{N}(undef, dim, dim)
    for j in 1:dim
        for i in 1:dim
            entry = rand(N)
            if rand(Bool)
                entry = -entry
            end
            shape_matrix[i, j] = entry
        end
    end
    # make diagonally dominant
    shape_matrix = N(0.5) * (shape_matrix + shape_matrix') +
                   Matrix{N}(dim*I, dim, dim)
    return Ellipsoid(center, shape_matrix)
end

"""
    translate(E::Ellipsoid{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:AbstractFloat}

Translate (i.e., shift) an ellipsoid by a given vector.

### Input

- `E`     -- ellipsoid
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated ellipsoid.

### Notes

The shape matrix is shared with the original ellipsoid if `share == true`.

### Algorithm

We add the vector to the center of the ellipsoid.
"""
function translate(E::Ellipsoid{N}, v::AbstractVector{N}; share::Bool=false
                  ) where {N<:AbstractFloat}
    @assert length(v) == dim(E) "cannot translate a $(dim(E))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    c = center(E) + v
    shape_matrix = share ? E.shape_matrix : copy(E.shape_matrix)
    return Ellipsoid(c, shape_matrix)
end


text:
using Cubature
# using Cuba
# using MonteCarloIntegration

include("vegas.jl")

function koopman(g,prob,u0,p,args...;kwargs...)
    g(solve(remake(prob,u0=u0,p=p),args...;kwargs...))
end

function koopman_cost(u0s,ps,g,prob,args...;maxevals=0,
                      ireltol = 1e-2, iabstol=1e-2, use_vegas=false, kwargs...)
    n = length(u0s)
    function _f(x)
    u0 = x[1:n]
    p  = x[n+1:end]
    k = koopman(g,prob,u0,p,args...;kwargs...)
    w = prod(pdf(a,b) for (a,b) in zip(u0s,u0))*
        prod(pdf(a,b) for (a,b) in zip(ps,p))
    k*w
    end
    xs = [u0s;ps]
    st = minimum.(xs)
    en = maximum.(xs)
    @show st
    @show en

    if use_vegas
        # vegas((x,f) -> f[1] = _f(st .+ x ./ (en .- st)) * sum(en .- st), rtol = ireltol)
        vegas(_f, st, en)
    else
        hcubature(_f, minimum.(xs), maximum.(xs);
            reltol=ireltol, abstol=iabstol, maxevals = maxevals)
    end
end

function montecarlo_cost(u0s,ps,g,prob,args...;num_monte,kwargs...)
  prob_func = function (prob,i,repeat)
    remake(prob,u0=rand.(u0s),p=rand.(ps))
  end
  output_func = (sol,i) -> (g(sol),false)
  monte_prob = MonteCarloProblem(prob;
                                 output_func = output_func,
                                 prob_func = prob_func)
  mean(solve(monte_prob,args...;num_monte=num_monte,kwargs...).u)
end


text:
using SentinelArrays, Test, Random

@testset "SentinelArrays" begin

# this is a very rare corner case
# we put this test first because it relies on
# the seeded random number generator for thread 1
# the issue is if we had an underlying array of non-sentinel values
# then tried to `setindex!` to the _NEXT_ chosen sentinel value
# then the element getting set would end up `missing` (because it == the sentinel value)
# instead of the sentinel value getting cycled to something else
Random.seed!(SentinelArrays.RNG[1], 0)
x = SentinelVector{Int64}(undef, 1)
x[1] = 1
x.sentinel = 1369352191816061504
x[1] = 1369352191816061504
@test x[1] == 1369352191816061504

x = SentinelVector{Int}(undef, 10)
fill!(x, missing)
@test length(x) == 10
@test all(isequal.(x.data, x.sentinel))

@test x[1] === missing
# force recoding
sent = x.sentinel
x[1] = x.sentinel
@test x[1] === sent
r = reverse(x)
@test x.sentinel == r.sentinel

@test size(x) == (10,)
x[1] = 3
@test x[1] === 3

resize!(x, length(x) + 1)
@test x[end] === missing

x = SentinelVector{Int64}(undef, 1)
x[1] = missing
@test x[1] === missing

x = SentinelArray{Float64, 2}(undef, 10, 10)
@test size(x) == (10, 10)
x = SentinelArray{Float64}(undef, 10, 10)
@test size(x) == (10, 10)

x = SentinelArray(fill(3.0, 10, 10))
y = convert(SentinelArray{Int64}, x)
@test size(y) == (10, 10)
@test y isa SentinelArray{Int64}
@test all(y .=== Int64(3))

y = convert(SentinelVector{Int64}, x)
@test size(y) == (10, 10)
@test y isa SentinelArray{Int64}
@test all(y .=== Int64(3))

y = convert(SentinelArray, fill(Int64(3), 10, 10))
@test size(y) == (10, 10)
@test size(y, 1) == 10
@test axes(y, 1) == Base.OneTo(10)
@test stride(y, 1) == 1
@test strides(y) == (1, 10)
@test y isa SentinelArray{Int64}
@test all(y .=== Int64(3))

x = SentinelVector{Union{Bool, Missing}}(undef, 1, missing, missing)
@test x[1] === missing
x[1] = true
@test x[1] === true
x[1] = missing
@test x[1] === missing

x = SentinelArray(["$i" for i = 1:10])
deleteat!(x, [9, 10])
@test length(x) == 8
@test x == ["$i" for i = 1:8]
deleteat!(x, [true, false, true, false, true, false, true, false])
@test length(x) == 4
@test x == ["2", "4", "6", "8"]

x = SentinelVector{String}(undef, 10)
@test x[1] === missing
x[1] = "hey"
@test x[1] == "hey"
x[1] = missing
@test x[1] === missing

@test SentinelArrays.newsentinel!(x) === nothing

@test all(x .=== copy(x))
@test length(empty!(x)) == 0
@test length(resize!(x, 10)) == 10
@test x[1] === missing

push!(x, missing)
@test x[end] === missing
push!(x, "ho")
@test x[end] == "ho"

deleteat!(x, length(x))
@test x[end] === missing
@test length(x) == 11

empty!(x)
append!(x, ["hey", "ho", missing])
@test length(x) == 3
@test isequal(x, ["hey", "ho", missing])

pushfirst!(x, missing)
@test x[1] === missing
prepend!(x, [missing, "first"])
@test x[1] === missing
@test x[2] == "first"

@test popfirst!(x) === missing
@test x[1] == "first"

insert!(x, 1, missing)
@test x[1] === missing
insert!(x, length(x) + 1, "pirate")
@test x[end] == "pirate"

@test splice!(x, length(x), ["pirate2"]) == "pirate"
@test splice!(x, length(x)) == "pirate2"
@test splice!(x, length(x), ["pirate3", "pirate4"]) === missing

@test splice!(x, (length(x)-1):length(x), ["pirate5", "pirate6"]) == ["pirate3", "pirate4"]
@test splice!(x, (length(x)-1):length(x), ["pirate7"]) == ["pirate5", "pirate6"]
@test splice!(x, length(x), ["pirate8", "pirate9"]) == "pirate7"
@test splice!(x, (length(x)-1):length(x), ["pirate10", "pirate11", "pirate12"]) == ["pirate8", "pirate9"]
@test splice!(x, (length(x)-2):length(x)) == ["pirate10", "pirate11", "pirate12"]

t = [SentinelVector{Int64}(undef, 10), SentinelVector{Int64}(undef, 10), SentinelVector{Int64}(undef, 5)]
sent = t[1].sentinel
@test all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...; force=false)
# make sure nothing got recoded
@test all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...; force=true)
# make sure all got recoded
@test all(x->x.sentinel != sent, t)

# force recode of just one vector
sent = t[1].sentinel
t[2][1] = t[2].sentinel
@test !all(x->x.sentinel == sent, t)
SentinelArrays.newsentinel!(t...)
# make sure all got recoded to same
@test all(x->x.sentinel == t[1].sentinel, t)

A = SentinelArray(Int64[i for i = 1:10])
B = SentinelVector{Int64}(undef, 10)

C = vcat(A, B)
@test C[1:10] == collect(1:10)
@test all(C[11:20] .=== missing)

append!(A, B)
@test A[1:10] == collect(1:10)
@test all(A[11:20] .=== missing)

A = SentinelArray(Int64[i for i = 1:10])
B = SentinelVector{Int64}(undef, 10)
# force B to recode
B[1] = B.sentinel
B[1] = missing

C = vcat(A, B)
@test C[1:10] == collect(1:10)
@test all(C[11:20] .=== missing)

append!(A, B)
@test A[1:10] == collect(1:10)
@test all(A[11:20] .=== missing)

# make sure we bail when can't find a new automatic sentinel
A = SentinelArray([i for i = 0x00:0xff])
@test_throws SentinelCollisionError setindex!(A, A.sentinel, 1)

@test_throws ErrorException SentinelVector{Bool}(undef, 1)

t = SentinelVector{Tuple{Int32, Int32}}(undef, 1)
@test t.data[1] === t.sentinel

t = SentinelMatrix{Float64}(undef, (10, 10))
@test size(t) == (10, 10)

t = SentinelArray(collect(1:10))
pushfirst!(t, 3, 2, 1)
@test t[1:3] == [3, 2, 1]

prepend!(t, (i for i = 1:3 if i > 0))
@test t[1:3] == [1, 2, 3]

@test pop!(t) == 10
empty!(t)
@test_throws ArgumentError pop!(t)
@test_throws ArgumentError popfirst!(t)

x = SentinelArray(ones(3))

# broadcasting
@test x .+ x == 2 * ones(3)
@test x .+ x .+ x == 3 * ones(3)
@test (x .= 0 .* x .+ 7) == 7 * ones(3)

v = SentinelArray(zeros(3,))
m = SentinelArray(ones(3, 3))
s = 0

@test v .+ m == ones(3, 3) == m .+ v
@test s .+ m == ones(3, 3) == m .+ s
@test s .+ v .+ m == ones(3, 3) == m .+ s .+ v

casts = (
    SentinelArray,  # Named Matrix
    x->SentinelArray(x[:, 1]),  # Named Vector
    x->SentinelArray(x[:, 1:1]),  # Named Single Column Matrix
    identity, # Matrix
    x->x[:, 1], # Vector
    x->x[:, 1:1], # Single Column Matrix
    first, # Scalar
    )
for (T1, T2, T3) in Iterators.product(casts, casts, casts)
    all(isequal(identity), (T1, T2, T3)) && continue
    !any(isequal(SentinelArray), (T1, T2, T3)) && continue

    total = T1(ones(3, 6)) .+ T2(2ones(3, 6)) .+ T3(3ones(3, 6))
    @test total == 6ones(3, 6)
end

s = SentinelVector{Int}(undef, 2)
@test isequal(max.(s, 0), [missing, missing])
@test (s .=== missing) isa BitArray

s = SentinelArray(collect(1:5))
c = ChainedVector([s, s, s])
c2 = copy(c)
@test length(c) == length(c2)
@test c2 isa Vector

# deleteat! of SentinelArray w/ underlying ChainedVector w/ UndefInitializer
x = SentinelArray(ChainedVector([Vector{String}(undef, 5)]))
deleteat!(x, 1)
@test length(x) == 4

end # @testset

@testset "MissingVector" begin

x = MissingVector(10)
@test all(x .=== missing)
@test length(x) == 10
@test length(x) isa Int
@test Base.IndexStyle(x) == Base.IndexLinear()

y = similar(x, Missing, 5)
@test length(y) == 5

y = empty(x)
@test length(y) == 0

x[1] = missing
x[end] = missing
@test x[1] === missing
@test x[end] === missing

y = similar(x)
@test typeof(y) == MissingVector
@test length(y) == 10
y = similar(x, 20)
@test typeof(y) == MissingVector
@test length(y) == 20

@test isequal(copy(x), x)
empty!(x)
@test length(x) == 0
@test isequal(copy(x), x)

@test_throws ArgumentError resize!(x, -1)
resize!(x, 10)
@test length(x) == 10

push!(x, missing)
@test x[end] === missing
empty!(x)
push!(x, missing)
@test x[1] === x[end] === missing

pushfirst!(x, missing)
@test x[1] === missing
empty!(x)
pushfirst!(x, missing)
@test x[1] === x[end] === missing
pushfirst!(x, missing)

@test pop!(x) === missing
@test popfirst!(x) === missing
@test isempty(x)

@test_throws BoundsError insert!(x, 0, missing)
@test_throws BoundsError insert!(x, 2, missing)
insert!(x, 1, missing)
@test x[1] === missing
insert!(x, 1, missing)
@test x[1] === missing

x = MissingVector(10)
y = MissingVector(10)

z = vcat(x, y)
@test length(z) == 20

empty!(x)
z = vcat(x, y)
@test isequal(z, y)

x = MissingVector(10)
append!(x, y)
@test length(x) == 20

x = MissingVector(10)
append!(x, (missing for _ = 1:10))
@test length(x) == 20

empty!(x)
append!(x, y)
@test isequal(x, y)

x = MissingVector(10)
y = MissingVector(10)

prepend!(y, x)
@test length(y) == 20

y = MissingVector(10)
prepend!(y, (missing for _ = 1:10))
@test length(y) == 20

empty!(y)
prepend!(y, x)
@test isequal(y, x)

x = MissingVector(10)
deleteat!(x, 1)
@test x[1] === missing
deleteat!(x, 1:4)
@test length(x) == 5
deleteat!(x, [2, 4])
@test length(x) == 3

m = similar(x)
@test length(m) == length(x)
m = similar(x, Missing)
@test length(m) == length(x)
@test typeof(m[1:3]) == typeof(m)

deleteat!(x, [true, true, false])
@test length(x) == 1
empty!(x)
@test_throws ArgumentError pop!(x)
@test_throws ArgumentError popfirst!(x)

m = MissingVector(5)
c = ChainedVector([m, m, m])
c2 = copy(c)
@test length(c) == length(c2)
@test c2 isa MissingVector

deleteat!(c2, Int[])
@test length(c2) == 15

end

include("chainedvector.jl")


text:
#------------------------------------------------------------------------------#
#                                   Package                                    #
#------------------------------------------------------------------------------#

"""
# 1. Description

    module HermitePolyTensor

Symbolic Hermite Polynomial Tensors — or N-dimensional Hermite polynomials [1] —
in one to three Euclidean space dimensions.

# 2. Usage

```julia-repl
julia> import HermitePolyTensor

```

# 3. References

[1]: H. Grad, “Note on N-dimensional  Hermite  polynomials,”  Communications  on
Pure and Applied Mathematics, vol. 2, no. 4, pp. 325–330, 1949.
"""
module HermitePolyTensor

#------------------------------------------------------------------------------#
#                                   Imports                                    #
#------------------------------------------------------------------------------#

using Reexport

#------------------------------------------------------------------------------#
#                                   Includes                                   #
#------------------------------------------------------------------------------#

include("version.jl")

include("HPT.jl")
@reexport using .HPT

end # module


text:
module TestExpressionsLocals

import Test: @testset, @test
import Volcanito: unescape_locals

input_output = (
    (
        input=:(a + b + $(Expr(:$, :c))),
        output=:(a + b + c),
    ),
)

@testset "unescape_locals" begin
    for pair in input_output
        @test unescape_locals(pair.input) == pair.output
    end
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 1e957cea-7f1e-11eb-065e-6b9e6f9f0976
md"# Linear Regression"

# ╔═╡ 79c2e252-7f22-11eb-21d4-e1836ee79caa

# ╔═╡ a1e6948c-7f21-11eb-2485-654e59594d7c

# ╔═╡ 9570fb28-7f21-11eb-2c7a-ffa16a382522

# ╔═╡ Cell order:
# ╠═1e957cea-7f1e-11eb-065e-6b9e6f9f0976
# ╠═79c2e252-7f22-11eb-21d4-e1836ee79caa
# ╠═a1e6948c-7f21-11eb-2485-654e59594d7c
# ╠═9570fb28-7f21-11eb-2c7a-ffa16a382522

max_stars_repo_name: AI-Pree/ML-and-algorithms

max_issues_repo_name: AI-Pree/ML-and-algorithms

max_forks_repo_name: AI-Pree/ML-and-algorithms


text:
# This file was generated, do not modify it. # hide
mμσ = time_normal((μ,σ) -> Normal(μ,σ))
dμσ = time_normal((μ,σ) -> Dists.Normal(μ,σ; check_args=false))

mμ1 = time_normal((μ,σ) -> Normal(μ=μ))
dμ1 = time_normal((μ,σ) -> Dists.Normal(μ))

m0σ = time_normal((μ,σ) -> Normal(σ=σ))
d0σ = time_normal((μ,σ) -> Dists.Normal(0.0,σ; check_args=false))

m01 = time_normal((μ,σ) -> Normal())
d01 = time_normal((μ,σ) -> Dists.Normal())

max_stars_repo_name: cscherrer/InformativePrior

max_issues_repo_name: cscherrer/InformativePrior

max_forks_repo_name: cscherrer/InformativePrior


text:
const Region3_vt =(
    n = [0.106_580_700_285_13E1,
    -0.157_328_452_902_39E2,
     0.209_443_969_743_07E2,
    -0.768_677_078_787_16E1,
     0.261_859_477_879_54E1,
    -0.280_807_811_486_20E1,
     0.120_533_696_965_17E1,
    -0.845_668_128_125_02E-2,
    -0.126_543_154_777_14E1,
    -0.115_244_078_066_81E1,
     0.885_210_439_843_18,
    -0.642_077_651_816_07,
     0.384_934_601_866_71,
    -0.852_147_088_242_06,
     0.489_722_815_418_77E1,
    -0.305_026_172_569_65E1,
     0.394_205_368_791_54E-1,
     0.125_584_084_243_08,
    -0.279_993_296_987_10,
     0.138_997_995_694_60E1,
    -0.201_899_150_235_70E1,
    -0.821_476_371_739_63E-2,
    -0.475_960_357_349_23,
     0.439_840_744_735_00E-1,
    -0.444_764_354_287_39,
     0.905_720_707_197_33,
     0.705_224_500_879_67,
     0.107_705_126_263_32,
    -0.329_136_232_589_54,
    -0.508_710_620_411_58,
    -0.221_754_008_730_96E-1,
     0.942_607_516_650_92E-1,
     0.164_362_784_479_61,
    -0.135_033_722_413_48E-1,
    -0.148_343_453_524_72E-1,
     0.579_229_536_280_84E-3,
     0.323_089_047_037_11E-2,
     0.809_648_029_962_15E-4,
    -0.165_576_797_950_37E-3,
    -0.449_238_990_618_15E-4],

I = [0,
     0,
     0,
     0,
     0,
     0,
     0,
     0,
     1,
     1,
     1,
     1,
     2,
     2,
     2,
     2,
     2,
     2,
     3,
     3,
     3,
     3,
     3,
     4,
     4,
     4,
     4,
     5,
     5,
     5,
     6,
     6,
     6,
     7,
     8,
     9,
     9,
     10,
     10,
     11],

J = [0,
     0,
     1,
     2,
     7,
     10,
     12,
     23,
     2,
     6,
     15,
     17,
     0,
     2,
     6,
     7,
     22,
     26,
     0,
     2,
     4,
     16,
     26,
     0,
     2,
     4,
     26,
     1,
     3,
     26,
     0,
     2,
     26,
     2,
     26,
     2,
     26,
     0,
     1,
     26])

     #==
function Region3_ρ(Output::Symbol, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J

    δ = ρ / ρstar
    τ = Tstar / T

    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    ϕ_δδ = -n[1]/(δ^2) + sum(n[i]*I[i]*(I[i]-1)*(δ^(I[i]-2))*(τ^J[i]) for i=2:40)
    ϕ_τ  =  sum(n[i]*(δ^I[i])*J[i]*(τ^(J[i]-1)) for i=2:40)
    ϕ_ττ =  sum(n[i]*(δ^I[i])*J[i]*(J[i]-1)*(τ^(J[i]-2)) for i=2:40)
    ϕ_δτ =  sum(n[i]*I[i]*(δ^(I[i]-1))*J[i]*(τ^(J[i]-1)) for i=2:40)

    if Output == :SpecificG            #kJ/kg
        return IF97_R*T*(δ*ϕ_δ + ϕ)
    elseif Output == :SpecificF        #kJ/kg
        return IF97_R*T*ϕ
    elseif Output == :Pressure         #MPa
        return ρ*IF97_R*T*δ*ϕ_δ/1000
    elseif Output == :SpecificU        #kJ/kg
        return IF97_R*T*τ*ϕ_τ
    elseif Output == :SpecificS        #kJ/kgK
        return IF97_R*(τ*ϕ_τ - ϕ)
    elseif Output == :SpecificH        #kJ/kg
        return IF97_R*T*(τ*ϕ_τ + δ*ϕ_δ)
    elseif Output == :SpecificCP       #kJ/kgK
        return IF97_R*(-τ^2*ϕ_ττ+((δ*ϕ_δ - δ*τ*ϕ_δτ)^2)/(2*δ*ϕ_δ + δ^2*ϕ_δδ))
    elseif Output == :SpecificCV       #kJ/kgK
        return IF97_R*(-τ^2*ϕ_ττ)
    elseif Output == :SpeedOfSound     #m/s
        return sqrt(1000*IF97_R*T*(2*δ*ϕ_δ + δ^2*ϕ_δδ - ((δ*ϕ_δ - δ*τ*ϕ_δτ)^2)/(τ^2*ϕ_ττ)))
    else
        throw(DomainError(Output, "Unknown value requested."))
    end
end
==#
function mass_gibbs_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    res = IF97_R*T*(δ*ϕ_δ + ϕ)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_helmholtz_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)

    res = IF97_R*T*ϕ
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function pressure_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    res = ρ*IF97_R*T*δ*ϕ_δ/1000
    return convert_unit(u"MPa",u"Pa",res)
end

function mass_internal_energy_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)
    ϕ_τ  =  sum(n[i]*(δ^I[i])*J[i]*(τ^(J[i]-1)) for i=2:40)
    res = IF97_R*T*τ*ϕ_τ
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_entropy_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)
    ϕ_τ  =  sum(n[i]*(δ^I[i])*J[i]*(τ^(J[i]-1)) for i=2:40)
    res = IF97_R*(τ*ϕ_τ - ϕ)
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function mass_enthalpy_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ    =  n[1]*log(δ) + sum(n[i]*(δ^I[i])*(τ^J[i]) for i=2:40)
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    ϕ_τ  =  sum(n[i]*(δ^I[i])*J[i]*(τ^(J[i]-1)) for i=2:40)
    res = IF97_R*T*(τ*ϕ_τ + δ*ϕ_δ)
    return convert_unit(u"kJ/kg",u"J/kg",res)
end

function mass_cp_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ_ττ =  sum(n[i]*(δ^I[i])*J[i]*(J[i]-1)*(τ^(J[i]-2)) for i=2:40)
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    ϕ_δτ =  sum(n[i]*I[i]*(δ^(I[i]-1))*J[i]*(τ^(J[i]-1)) for i=2:40)
    ϕ_δδ = -n[1]/(δ^2) + sum(n[i]*I[i]*(I[i]-1)*(δ^(I[i]-2))*(τ^J[i]) for i=2:40)

    res = IF97_R*(-τ^2*ϕ_ττ+((δ*ϕ_δ - δ*τ*ϕ_δτ)^2)/(2*δ*ϕ_δ + δ^2*ϕ_δδ))
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function mass_cv_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ_ττ =  sum(n[i]*(δ^I[i])*J[i]*(J[i]-1)*(τ^(J[i]-2)) for i=2:40)

    res = IF97_R*(-τ^2*ϕ_ττ)
    return convert_unit(u"kJ/(kg*K)",u"J/(kg*K)",res)
end

function sound_speed_impl(mt::SingleVT,model::IF97Region{:r3}, _ρ, _T)
    ρ = normalize_units(_ρ)
    T = normalize_units(_T)
    ρstar = IF97_ρc_inv      #kg/m3
    Tstar = IF97_Tc   #K
    n = Region3_vt.n
    I = Region3_vt.I
    J = Region3_vt.J
    δ = ρ * ρstar
    τ = Tstar / T
    ϕ_ττ =  sum(n[i]*(δ^I[i])*J[i]*(J[i]-1)*(τ^(J[i]-2)) for i=2:40)
    ϕ_δ  =  n[1]/δ + sum(n[i]*I[i]*(δ^(I[i]-1))*(τ^J[i]) for i=2:40)
    ϕ_δτ =  sum(n[i]*I[i]*(δ^(I[i]-1))*J[i]*(τ^(J[i]-1)) for i=2:40)
    ϕ_δδ = -n[1]/(δ^2) + sum(n[i]*I[i]*(I[i]-1)*(δ^(I[i]-2))*(τ^J[i]) for i=2:40)

    res = sqrt(1000*IF97_R*T*(2*δ*ϕ_δ + δ^2*ϕ_δδ - ((δ*ϕ_δ - δ*τ*ϕ_δτ)^2)/(τ^2*ϕ_ττ)))
    return res
end

#==
function Region3(Output::Symbol, P, T)
    ρ0 = 1000*P/(R*T) #Starting value from ideal gas

    f(ρ) = Region3_ρ(:Pressure, ρ, T) - P
    ρ = Roots.find_zero(f, ρ0)

    if Output == :SpecificV
        return 1.0 / ρ
    else
        return Region3_ρ(Output, ρ, T)
    end
end
==#

for op in [:mass_gibbs_impl,:mass_helmholtz_impl,
    :mass_enthalpy_impl,:mass_internal_energy_impl,:mass_entropy_impl,
    :mass_cv_impl,:mass_cp_impl,:sound_speed_impl]

    @eval begin
        function $op(mt::SinglePT,model::IF97Region{:r3},p,t)
            P =convert_unit(u"Pa",u"MPa",p)
            T = normalize_units(t)
            ρ0 = 1000*P/(IF97_R*T) #Starting value from ideal gas
            f(ρ) = pressure_impl(QuickStates.ρt(), model,ρ, T) - normalize_units(p)
            ρ = Roots.find_zero(f, ρ0)
            return $op(QuickStates.ρt(), model,ρ, T)
        end
    end

end

function mass_volume_impl(mt::SinglePT,model::IF97Region{:r3},p,t)
    P =convert_unit(u"Pa",u"MPa",p)
    T = normalize_units(t)
    ρ0 = 1000*P/(IF97_R*T) #Starting value from ideal gas
    f(ρ) = pressure_impl(QuickStates.ρt(), model,ρ, T) - normalize_units(p)
    ρ = Roots.find_zero(f, ρ0)
    return inv(ρ)
end


text:
using DataFrames
using Missings
using EcoBase

"""
    Metacommunity{FP, ARaw, AProcessed, Part, Sim}

Metacommunity type, representing a whole metacommunity containing a
single community or a collection of subcommunities. The metacommunity
of individuals *may* be further partitioned into smaller groups. For
instance this may be an ecosystem, which consists of a series of
subcommunities. The AbstractPartition subtype within it stores
relative abundances of different types, e.g. species, and also allows
for similarity between individuals.

# Constructor:

Metacommunity(abundances::AbstractArray,
              part::AbstractPartition,
              types::AbstractTypes)

# Members:

- `abundances` the abundance matrix for the metacommunity.

- `partition` the instance of the AbstractPartition subtype, containing the
  subcommunities.

- `types` The instance of the AbstractTypes subtype, from which
  similarities between individuals can be calculated.

- `ordinariness` A cache of the ordinariness of the individuals in the
  Partition. Should only be accessed through
  getordinariness!(::Metacommunity), which will populate the cache if
  it has not yet been calculated.

"""
mutable struct Metacommunity{FP, ARaw, AProcessed, Sim, Part} <:
    Diversity.API.AbstractMetacommunity{FP, ARaw, AProcessed, Sim, Part}
    rawabundances::ARaw
    processedabundances::AProcessed
    scale::FP
    types::Sim
    partition::Part
    ordinariness::Union{AProcessed, Missing}

    function Metacommunity{FP, ARaw, AProcessed,
                           Sim, Part}(abundances::ARaw,
                                      matrix::AProcessed,
                                      types::Sim,
                                      part::Part) where
        {FP <: AbstractFloat, ARaw <: AbstractArray, AProcessed <: AbstractArray{FP},
         Sim <: AbstractTypes, Part <: AbstractPartition}
        mcmatch(matrix, types, part) ||
            error("Type or size mismatch between abundance array, " *
                  "partition and type list")
        processedabundances, scale = _calcabundance(types, matrix)
        new{FP, ARaw, AProcessed, Sim, Part}(abundances, processedabundances, scale,
                                             types, part, missing)
    end
end

function Metacommunity(abundances::ARaw, meta::Meta) where
    {ARaw <: AbstractArray, Meta <: AbstractMetacommunity}
    types = gettypes(meta)
    part = getpartition(meta)
    mat = reshape(abundances, counttypes(types), countsubcommunities(part))
    if sum(mat) ≉ one(eltype(mat))
        @warn "Abundances not normalised to 1, correcting..."
        mat = mat / sum(mat)
    end
    return Metacommunity{eltype(mat), ARaw, typeof(mat),
                         typeof(types), typeof(part)}(abundances, mat, types, part)
end

function Metacommunity(abundances::V,
                       types::Sim = UniqueTypes(size(abundances, 1)),
                       part::Part = Onecommunity()) where
    {V <: AbstractVector, Sim <: AbstractTypes, Part <: AbstractPartition}
    mat = reshape(abundances / sum(abundances), length(abundances), 1)
    return Metacommunity{eltype(mat), V, typeof(mat),
                         typeof(types), typeof(part)}(abundances, mat, types, part)
end

function Metacommunity(abundances::V,
                       types::Sim = UniqueTypes(size(abundances, 1)),
                       part::Part = Onecommunity()) where
    {FP <: AbstractFloat, V <: AbstractVector{FP},
     Sim <: AbstractTypes, Part <: AbstractPartition}
    mat = reshape(abundances, length(abundances), 1)
    if sum(mat) ≉ one(eltype(mat))
        @warn "Abundances not normalised to 1, correcting..."
        mat = mat / sum(mat)
    end
    return Metacommunity{eltype(mat), typeof(abundances), typeof(mat),
                         typeof(types), typeof(part)}(abundances, mat, types, part)
end

function Metacommunity(abundances::M,
                       types::Sim = UniqueTypes(size(abundances, 1)),
                       part::Part = Subcommunities(size(abundances, 2))) where
    {M <: AbstractMatrix, Sim <: AbstractTypes, Part <: AbstractPartition}
    mat = abundances / sum(abundances)
    return Metacommunity{eltype(mat), M, typeof(mat),
                         typeof(types), typeof(part)}(abundances, mat, types, part)
end

function Metacommunity(abundances::M,
                       types::Sim = UniqueTypes(size(abundances, 1)),
                       part::Part = Subcommunities(size(abundances, 2))) where
    {FP <: AbstractFloat, M <: AbstractMatrix{FP},
     Sim <: AbstractTypes, Part <: AbstractPartition}
    mat = abundances
    if sum(mat) ≉ one(eltype(mat))
        @warn "Abundances not normalised to 1, correcting..."
        mat = mat / sum(mat)
    end
    return Metacommunity{eltype(mat), typeof(abundances), typeof(mat),
                         typeof(types), typeof(part)}(abundances, mat, types, part)
end

function Metacommunity(abundances::V, zmatrix::M) where
    {FP <: AbstractFloat, V <: AbstractVector, M <: AbstractMatrix{FP}}
    return Metacommunity(abundances, GeneralTypes(zmatrix), Onecommunity())
end

function Metacommunity(abundances::MU, zmatrix::M) where
    {FP <: AbstractFloat, MU <: AbstractMatrix{FP}, M <: AbstractMatrix{FP}}
    return Metacommunity(abundances, GeneralTypes(zmatrix),
                         Subcommunities(size(abundances, 2)))
end

Metacommunity(asm::EcoBase.AbstractAssemblage) =
    hassimilarity(asm) ?
        Metacommunity(occurrences(asm), _calcsimilarity(asm)) :
        Metacommunity(occurrences(asm))

import Diversity.API._gettypes
_gettypes(meta::Metacommunity) = meta.types

import Diversity.API._getpartition
_getpartition(meta::Metacommunity) = meta.partition

import Diversity.API._getabundance
function _getabundance(meta::Metacommunity, raw::Bool)
    return raw ? meta.rawabundances : meta.processedabundances
end

import Diversity.API._getordinariness!
function _getordinariness!(meta::Metacommunity)
    if ismissing(meta.ordinariness)
        meta.ordinariness = _calcordinariness(meta.types,
                                              meta.processedabundances,
                                              meta.scale)
    end
    return meta.ordinariness
end

import Diversity.API._getscale
_getscale(m::Metacommunity) = m.scale


text:
function SplineInterpolator()
    return SplineInterpolator(())
end

function interpolate(obj::SplineInterpolator, arg0::Vector{RealFieldElement}, arg1::Vector{RealFieldElement})
    return jcall(obj, "interpolate", FieldPolynomialSplineFunction, (Vector{RealFieldElement}, Vector{RealFieldElement}), arg0, arg1)
end

function interpolate(obj::SplineInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble})
    return jcall(obj, "interpolate", PolynomialSplineFunction, (Vector{jdouble}, Vector{jdouble}), arg0, arg1)
end


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

function build_default_parameter_dictionary()::Dict{String,Any}

    # initailize -
    parameter_dictionary = Dict{String,Any}()
    system_type_flag = 	:HUMAN

    try

        # get the biophysical parameters for this system type -
        biophysical_parameters_dictionary = get_biophysical_parameters_dictionary_for_system(system_type_flag)

        # initial condition array -
        initial_condition_array = [

        ];

        # ======= DO NOT EDIT BELOW THIS LINE ========================================================== #
        parameter_dictionary["biophysical_parameters_dictionary"] = biophysical_parameters_dictionary
        parameter_dictionary["initial_condition_array"] = initial_condition_array
        return parameter_dictionary
        # ======= DO NOT EDIT ABOVE THIS LINE ========================================================== #
    catch error
        rethrow(error)
    end
end

max_stars_repo_name: varnerlab/VLCoVGenomeReplicationModel

max_issues_repo_name: varnerlab/VLCoVGenomeReplicationModel

max_forks_repo_name: varnerlab/VLCoVGenomeReplicationModel


text:
using InteractiveDynamics
using DynamicalSystems, GLMakie
using OrdinaryDiffEq

# Lorenz
ds = Systems.lorenz()

u1 = [10,20,40.0]
u2 = [10,20,40.0 + 1e-3]
u3 = [20,10,40.0]
u0s = [u1, u2, u3]

idxs = (1, 2, 3)
diffeq = (alg = Tsit5(), dtmax = 0.01)

figure, obs = interactive_evolution(
    ds, u0s; idxs, tail = 1000, diffeq
)

# %% Lorenz with timeseries
ds = Systems.lorenz()

u1 = [10,20,40.0]
u2 = [10,20,40.0 + 1e-3]
u3 = [20,10,40.0]
u0s = [u1, u2, u3]

idxs = (1, 2, 3)
diffeq = (alg = Tsit5(), dtmax = 0.01)

figure, obs = interactive_evolution_timeseries(
    ds, u0s; idxs, tail = 1000, diffeq
)

# %% Lorenz with timeseries and parameters
ps = Dict(
    1 => 1:0.1:30,
    2 => 10:0.1:100,
    3 => 1:0.01:10.0,
)
pnames = Dict(1 => "σ", 2 => "ρ", 3 => "β")

lims = (
    (-30, 30),
    (-30, 30),
    (0, 100),
)

ds = Systems.lorenz()

u1 = [10,20,40.0]
u3 = [20,10,40.0]
u0s = [u1, u3]

idxs = (1, 2, 3)
diffeq = (alg = Tsit5(), dtmax = 0.01)

figure, obs = interactive_evolution_timeseries(
    ds, u0s, ps; idxs, tail = 1000, diffeq, pnames, lims
)

# %% towel
ds = Systems.towel()
u0s = [0.1ones(3) .+ 1e-3rand(3) for _ in 1:3]
idxs = (1, 2, 3)

figure, obs = interactive_evolution_timeseries(
    ds, u0s; idxs, tail = 10000,
)

# %% SM
ds = Systems.standardmap()
u0s = [[0.1, 0.1], [2.5, 0.4], [1.88, 3.25]]
lims = ((0, 2π), (0, 2π))

figure, obs = interactive_evolution(
    ds, u0s; tail = 100000, lims
)

# %% Henon helies

ds = Systems.henonheiles()

u0s = [[0.0, -0.25, 0.42081, 0.0],
[0.0, 0.1, 0.5, 0.0],
[0.0, -0.31596, 0.354461, 0.0591255]]

diffeq = (alg = Vern9(), dtmax = 0.01)
idxs = (1, 2, 4)
colors = Makie.to_color.(["#233B43", "#499cbf", "#E84646"])

figure, obs = interactive_evolution(
    ds, u0s; idxs, tail = 10000, colors, diffeq
)
# main.scene[Axis][:names, :axisnames] = ("q₁", "q₂", "p₂")