text:
function post(𝒜::ASB07,
              𝑃::InitialValueProblem{<:AbstractContinuousSystem},
              𝑂::Options)
    # =================================
    # Initialization and discretization
    # =================================

    𝑂 = merge(𝒜.options.defaults, 𝑂, 𝒜.options.specified)
    δ, T = 𝑂[:δ], 𝑂[:T]
    N = round(Int, T / δ)

    # compute and unrwap discretized system
    𝑃_discrete = discretize(𝑃, δ; algorithm="interval_matrix",
                            order=𝑂[:order_discretization],
                            set_operations=𝑂[:set_operations_discretization])
    Ω0, Φ = 𝑃_discrete.x0, 𝑃_discrete.s.A
    if ! (genmat(Ω0) isa Matrix)
        Ω0 = Zonotope(center(Ω0), Matrix(genmat(Ω0)))
    end

    # ====================
    # Flowpipe computation
    # ====================

    # preallocate output
    T = 𝑂[:set_operations_discretization] == "zonotope" ?
        Zonotope{Float64, Vector{Float64}, Matrix{Float64}} :
        LazySet{Float64}
    Rsets = Vector{ReachSet{T}}(undef, N)

    max_order = 𝑂[:max_order]

    info("Reachable States Computation...")
    @timing begin
    if inputdim(𝑃_discrete) == 0
        U = nothing
    else
        U = inputset(𝑃_discrete)
    end
    reach_ASB07!(Rsets, Ω0, U, Φ, N, δ, max_order)
    end # timing

    Rsol = ReachSolution(Flowpipe(Rsets), 𝑂)

    # ==========
    # Projection
    # ==========

    if 𝑂[:project_reachset]
        info("Projection...")
        Rsol = @timing project(Rsol)
    end

    return Rsol
end


text:
using FastTransforms, LinearAlgebra, Test

include("specialfunctionstests.jl")
include("chebyshevtests.jl")
include("quadraturetests.jl")
include("libfasttransformstests.jl")
include("nuffttests.jl")
include("fftBigFloattests.jl")
include("paduatests.jl")
include("gaunttests.jl")
include("hermitetests.jl")
include("toeplitztests.jl")
include("clenshawtests.jl")

ext: jl


text:
using SialiaCLI
using Test

@testset "SialiaCLI.jl" begin
    # Write your own tests here.
end


text:
using FEniCS
using PyPlot

mu = 1
rho = 1
lambda = 1
g = 9.81

# Create mesh and define function space
mesh = UnitSquareMesh(20, 20) # TODO: use RectangleMesh(0.0, 0.0, 10.0, 4.0, 10, 10)
V = VectorFunctionSpace(mesh, "P", 1)
c = Constant((0, 0))
bc = DirichletBC(V, c, "on_boundary && x[1]<1E-14") # What BCs are bing set???
tol = 1E-14

# Define strain
function epsilon(u)
     0.5 * (nabla_grad(u) + Transpose(nabla_grad(u)))
end

# Define stress
function sigma(u)
     lambda * nabla_div(u) * Identity(d) + 2 * mu * epsilon(u)
end

u = TrialFunction(V)
d = geometric_dimension(u) # space dimension
v = TestFunction(V)
f = Constant((0, -rho*g)) # Vector of uniform body force
T = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v) * dx + dot(T, v) * ds
u = FeFunction(V)
lvsolve(a, L, u, bc)

# s = sigma(u) - (1.0/3.0) * tr(sigma(u)) * Identity(d)
V = FunctionSpace(mesh, "P", 1)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plothandle = FEniCS.Plot(u_magnitude)
colorbar(plothandle)


text:
# Tests for the integrate_timestep! function

# Used in finite difference derivatives
function compute_timestep(params, t0, a, b, points, h, ia, ib, it)
    rstar = params[end]
    xc = components(points[ib,it,:,1]./rstar, h)
    yc = components(points[ib,it,:,2]./rstar, h)
    trans = transit_init(params[ib - 1], zero(typeof(t0)), params[end - 2:end-1], false)
    integrate_timestep!(t0, a, b, xc, yc, trans, ia)
    return ia.I_of_f[1]
end

# Used in finite difference derivatives
function compute_timestep_nbody(coords, s_copy, tt, trans, ia, ic, intr, a, b, it, ib)
    s = deepcopy(s_copy)
    coords = reshape(coords, 7, s.n)
    for (i, col) in enumerate(eachcol(coords))
        s.x[:,i] .= col[1:3]
        s.v[:,i] .= col[4:6]
        s.m[i]    = col[end]
    end

    pd = compute_pd(s, ic, tt, intr, grad=false);
    rstar = big(0.00465047 * 0.1192) # Trappist-1 (Rstar/AU)
    normalize_points!(pd.points, rstar);

    t0 = pd.times[it]
    xc = components(pd.points[ib,it,:,1], pd.h)
    yc = components(pd.points[ib,it,:,2], pd.h)

    integrate_timestep!(t0, a, b, xc, yc, trans, ia)
    return ia.I_of_f[1]
end

# Test the derivatives of the integrated exposure with respect to the initial
# Cartesian coordinates, radius ratios, and limbdarkening coefficents
function test_integrate_timestep_derivatives(n)
    # Setup initial conditions and compute points
    BJD = 7250.0; t0_ic = 7.0
    tmax = 2.0
    k = get_radius_ratios_trappist(n);
    u_n = get_limdark_coeffs_trappist();
    ic = setup_ICs(n, BJD, t0_ic);
    intr = setup_integrator(ic, tmax);
    tt = compute_transit_times(ic, intr);
    pd = compute_pd(ic, tt, intr, grad=true);

    # Normalize points to stellar radius
    rstar = 0.00465047 * 0.1192 # Trappist-1 (Rstar/AU)
    normalize_points!(pd.points, rstar);
    normalize_points!(pd.dpoints, rstar);

    for it in 1:length(pd.times)
        # Make a fake exposure time
        dt = 2 / 60 / 24 # 2 minute cadence in days
        ib = pd.bodies[it]
        t0 = pd.times[it]
        texp = t0 + dt / 2

        # Get a transit structure
        trans = transit_init(k[ib - 1], 0.0, u_n, false)

        # Compute x and y components for impact parameter calculation
        xc = components(pd.points[ib,it,:,1], pd.h)
        yc = components(pd.points[ib,it,:,2], pd.h)

        # Integrate a timestep over a -> b
        a = texp - dt / 2
        b = texp + dt / 2
        tol = 1e-15
        maxdepth = 40
        ia = IntegralArrays(1, maxdepth, tol)
        integrate_timestep!(t0, a, b, xc, yc, trans, ia)
        flux_int = ia.I_of_f[1]

        # Make bigfloat variants
        ic_big = setup_ICs(n, big(BJD), big(t0_ic));
        intr_big = setup_integrator(ic_big, big(tmax));
        tt_big = compute_transit_times(ic_big, intr_big);
        pd_big = compute_pd(ic_big, tt_big, intr_big, grad=false);
        ia_big = IntegralArrays(1, maxdepth, big(tol))

        # Make sure we get the same answer
        params = [k...,u_n...,rstar]
        flux_big = compute_timestep(big.(params), big(t0), big(a), big(b), copy(pd_big.points), pd_big.h, ia_big, ib, it)
        @test flux_big ≈ flux_int

        # Compute finite diff derivatives with respect to the transit parameters
        grad_num = grad(central_fdm(5, 1), p -> compute_timestep(p, big(t0), big(a), big(b), copy(pd_big.points), pd_big.h, ia_big, ib, it), big.(params))[1]

        # Compute analytic derivatives
        inds = [1,2,3,4,5,6,7]
        dxc = [components(pd.dpoints[ib,it,:,1,k,i], pd.h) for i in inds, k in 1:ic.nbody][:]
        dyc = [components(pd.dpoints[ib,it,:,2,k,i], pd.h) for i in inds, k in 1:ic.nbody][:]
        dbdq0 = zeros(7 * ic.nbody)
        ki = ib - 1
        trans_grad = transit_init(k[ki], 0.0, u_n, true)
        n_params = length(dbdq0) + length(k) + length(u_n) + 1 # flux
        ia_grad = IntegralArrays(n_params + 1, maxdepth, tol)
        integrate_timestep!(t0, a, b, xc, yc, dxc, dyc, trans_grad, ia_grad, dbdq0, ki)

        # Check flux again
        @test ia_grad.I_of_f[1] ≈ flux_int

        # Check derivative wrt radius ratios
        @test isapprox(ia_grad.I_of_f[2 + length(dbdq0):end-4], Float64.(grad_num[1:n-1]), norm=x -> maximum(abs.(x)))

        # Now the derivatives wrt limbdark coefficients
        dfdu = trans_grad.dgdu' * ia_grad.I_of_f[end-trans_grad.n-1:end-1]
        @test isapprox(dfdu, Float64.(grad_num[n:end-1]), norm=x -> maximum(abs.(x)))

        # Now the derivatives wrt the stellar radius
        dfdr = ia_grad.I_of_f[end] / rstar
        @test isapprox(dfdr, Float64.(grad_num[end]), norm=x -> maximum(abs(x)), atol=1e-8)

        # Finally, derivatives wrt the Nbody initial conditions
        s_copy = deepcopy(State(ic_big))
        params_nbody = vcat([[x...,v...,m] for (x, v, m) in zip(eachcol(s_copy.x), eachcol(s_copy.v), s_copy.m)]...)

        # Again, check flux
        trans = transit_init(k[ib - 1], 0.0, u_n, false)
        flux_nbody = compute_timestep_nbody(params_nbody, State(ic), tt, trans, ia, deepcopy(ic), intr, a, b, it, ib)
        @test flux_nbody ≈ flux_int

        # Compute finite FiniteDifferences
        trans_big = transit_init(big(k[ki]), big(0.0), big.(u_n), false)
        grad_nbody_num = grad(central_fdm(12, 1), p -> compute_timestep_nbody(p, s_copy, tt_big, trans_big, ia_big, ic_big, intr_big, big(a), big(b), it, ib), big.(params_nbody))[1]

        # Compare with analytic
        @test isapprox(ia_grad.I_of_f[2:1 + length(dbdq0)], Float64.(grad_nbody_num), norm=x -> maximum(abs.(x)))
    end
end

@testset "Integrated Flux" begin
    test_integrate_timestep_derivatives(3) # This runs for a while...
end


text:
import ThermoState
import ThermoState: @spec_str, FromState
import Unitful: @u_str

Unitful.@derived_dimension __MassDensity Unitful.𝐌/Unitful.𝐋^3
Unitful.@derived_dimension __MolDensity Unitful.𝐍/Unitful.𝐋^3
Unitful.@derived_dimension __MassVolume Unitful.𝐋^3/Unitful.𝐌
Unitful.@derived_dimension __MolVolume Unitful.𝐋^3/Unitful.𝐍

const __VolumeKind = Union{__MassDensity,__MolDensity,__MassVolume,__MolVolume,Unitful.Volume}

#real, use st argument
function autospec(λ::Real,st)
    return ThermoState.spec(st,λ)
end

#with units, use units
function autospec(λ::Unitful.Temperature,st)
    return ThermoState.spec(spec"T",λ)
end

function autospec(λ::Unitful.Pressure,st)
    return ThermoState.spec(spec"p",λ)
end

function autospec(λ::Unitful.Volume,st)
    return ThermoState.spec(spec"total_v",λ)
end

function autospec(λ::__MolDensity,st)
    return ThermoState.spec(spec"mol_rho",λ)
end

function autospec(λ::__MassDensity,st)
    return ThermoState.spec(spec"mass_rho",λ)
end

function autospec(λ::__MolVolume,st)
    return ThermoState.spec(spec"mol_v",λ)
end

function autospec(λ::__MassVolume,st)
    return ThermoState.spec(spec"mass_v",λ)
end

function autospec(t::AbstractVector{<:Real},st)
    if length(t) == 1
        return ThermoState.spec(spec"moles",only(t))
    else
        return ThermoState.spec(spec"n",t)
    end
end

function autospec(t::AbstractVector{<:Unitful.Amount},st)
    if length(t) == 1
        return ThermoState.spec(spec"moles",only(t))
    else
        return ThermoState.spec(spec"n",t)
    end
end

function autospec(t::AbstractVector{<:Unitful.Mass},st)
    if length(t) == 1
        return ThermoState.spec(spec"mass",only(t))
    else
        return ThermoState.spec(spec"m",t)
    end
end

ThermoState.molecular_weight(model::EoSModel) = mw(model)
#transforms vt or pt model to ThermoState.state
function standarize(model,λ,t,z)
    spec_t = autospec(t,spec"T")
    spec_λ = autospec(λ,spec"p") #using pressure by default
    spec_z = autospec(z,spec"n")
    st = ThermoState.state(spec_λ,spec_t,spec_z)
end

function state_to_vt(model,st::ThermoState.ThermodynamicState)
    return state_to_vt(ThermoState.state_type(st),model,st)
end

function state_to_vt(::ThermoState.QuickStates.SingleVT,model,st::ThermoState.ThermodynamicState)
    v = ThermoState.total_volume(FromState(),st,u"m^3",mw(model)) |> only
    T = ThermoState.temperature(FromState(),st,u"K") |> only
    _z = ThermoState.moles(FromState(),st,u"mol",mw(model))
    z = SA[_z]
    return v,T,z
end

function state_to_vt(::ThermoState.QuickStates.MultiVT,model,st::ThermoState.ThermodynamicState)
    v = ThermoState.total_volume(FromState(),st,u"m^3",mw(model)) |> only
    T = ThermoState.temperature(FromState(),st,u"K") |> only
    z = ThermoState.mol_number(FromState(),st,u"mol",mw(model))
    return v,T,z
end
function state_to_pt(model,st::ThermoState.ThermodynamicState)
    return state_to_pt(ThermoState.state_type(st),model,st)
end

function state_to_pt(::ThermoState.QuickStates.SinglePT,model,st::ThermoState.ThermodynamicState)
    p = ThermoState.pressure(FromState(),st,u"Pa") |> only
    T = ThermoState.temperature(FromState(),st,u"K") |> only
    _z = ThermoState.moles(FromState(),st,u"mol",mw(model))
    z = SA[_z]
    return p,T,z
end
function state_to_pt(::ThermoState.QuickStates.MultiPT,model,st::ThermoState.ThermodynamicState)
    p = ThermoState.pressure(FromState(),st,u"Pa") |> only
    T = ThermoState.temperature(FromState(),st,u"K") |> only
    z = ThermoState.mol_number(FromState(),st,u"mol",mw(model))
    return p,T,z
end


text:
### TEST FOR THE TRANSMISSION COEFFICIENT FOR DIFFERENT SOLITON VELOCITIES AND BARRIER STRENGTHS
## PLOTTING RESULTS IN NORMALIZED COORDINATES AS GARDNER et al
using SolitonBEC
using Printf
using Plots
using Elliptic
using Distributed
#using FileIO
hbar = 6.62607015e-34 / (2 * pi)
width = 1.42e-6
mass = 6.941 * 1.660539e-27
omega_perp = 2 * pi * 710
l_perp = sqrt(hbar / (mass * omega_perp))
N = 4e3
as = -0.21e-9

interaction_g = abs(2*hbar^2 * as / mass / l_perp^2)
ggg = 2*hbar*omega_perp*abs(as)

print("\ninteraction_g: ", interaction_g)

# mass as total mass or mass per atom??
time_unit = hbar^2/(mass * ggg * N)
space_unit = hbar^3/(mass * ggg^2 * N^2)
velocity_unit = ggg * N / hbar
print("\n\ttime_unit: ", time_unit)
print("\n\tspace_unit: ", space_unit)
print("\n\tvelocity unit: ", velocity_unit, " m/s")

# --------- Simulation ---------
khaykovich_gpe = Simulation(
  "NPSE",
  "barrier",
)

γ = 0 * 1.77e-11
# --------- Apparata ---------
std_apparatus = Apparatus(
  mass, #m (conversion AMU -> kg)
  as, # as
  omega_perp, # ω_perp
  N, #N
  γ, #γ
  2 * pi * 4,# ω_z
)

# Energy = 1 / 2 * (std_apparatus.as^2 / (width^2) + std_apparatus.as^2 / (2 * l_perp^2) + width^2 / std_apparatus.as^2 - interaction_g * std_apparatus.as^3 / (l_perp * width^2)) * hbar * omega_perp
L = 2 * space_unit

v0 = 1*velocity_unit
x0 =-L/4
# --------- InitialStates ---------
function initial_state_velocity(velocity::Float64)
  init = InitialState1 = InitialState(
    "sech", #type
    width*1.8, # width
    velocity, #v0
    x0,
  )
  return init
end

# --------- Numerics ---------
function adaptive_numerics(velocity::Float64, L, x0, velocity_unit)
  if (velocity==0)
    T = abs(x0)/ velocity_unit * 2
  else
    T = abs(x0)/ velocity * 2
  end
  num = Numerics(
    T, #T
    T * 2e-4, #dt
    L, #S
    L * 2e-4, #ds
  )
  return num
end

#--------- Potential ---------

function barrier_height(energy::Float64)
  r = Potential(
    "barrier", # type
    L * 1e-3, #width
    0, #position
    energy, #energy
    0, #ϵ
    10e-6 #a
  )
  return r
end

energy_unit = mass * N^2 * ggg^2 / hbar^2
print("\n\tenergy unit: ", energy_unit, " J")

GSEnergy = energy_unit * (-N/24)
#print("\nground state energy: ", GSEnergy/hbar, " hbar\n")
phys_vel = 4.918e-3
normd_vel = phys_vel*hbar/ggg/N
Energy = GSEnergy + energy_unit * normd_vel^2*N/2
#print("\ntraveling state energy: ", Energy/hbar, " hbar\n")

## ==================== Transmission grid configuration
configs = []
num_barr = 200
barrier_list = LinRange(0, 15036, num_barr)
velocity_list = LinRange(0, 1, num_barr)

for vel in velocity_list
  for barrier_energy in barrier_list
    potential = barrier_height(barrier_energy * hbar)
    state = initial_state_velocity(vel * velocity_unit)
    numerics = adaptive_numerics(vel * velocity_unit, L, x0, velocity_unit)
    push!(configs, (numerics, khaykovich_gpe, std_apparatus, potential, state))
  end
end

mem_limit = 15000000000 #byte
pyplot()
plt_width = 800
plt_height = 600

## Soliton - barrier collision --------------------------------------------
#run_dynamics(configs[20*20 - 5]...)

T = zeros(Float64, length(velocity_list), length(barrier_list))

nth = Threads.nthreads() #print number of threads
print("number of threads: ", nth)

Threads.@threads for iv in axes(velocity_list, 1)

    for (ib, barrier_energy) in enumerate(barrier_list)

    (numerics, sim, app, pot, state) = configs[(iv-1) * num_barr + ib]

    # potential space index
    coeffs = get_coefficients(sim, app, pot, state)
    time, space, ψ, ψ_spect = ssfm_solve(numerics, coeffs)
    #potential_idx = Int64(floor((pot.position+numerics.S/2) / numerics.ds))

    T[iv, ib] = sum(abs.(ψ[Int(floor(length(space)/2)):end, end]) .^ 2 * numerics.ds)
    print("\nT[", iv, ", ", ib ,"] = ", T[iv, ib])

    #counter = sum(abs.(ψ[1:Int(floor(length(space)/2)), end]) .^ 2 * numerics.ds)
    #print("\n\t Total integral: ", T[iv, ib]+counter)
    # fig = plot(space, abs.(ψ[:, end]).^2, title="Transmission heatmap",
    # xlabel="barrier",
    # ylabel="velocity",
    # reuse=false,
    # size=(plt_width, plt_height))
    # #plot!(abs.(ψ[Int(floor(length(space)/2)):end, end]))
    # #plot!(space, abs.(coeffs.β.(space)))
    # display(fig)
  end
end

fig1 = plot(title="Transmission heatmap",
  xlabel="barrier",
  ylabel="velocity",
  reuse=false,
  size=(plt_width, plt_height))
heatmap!(T)
display(fig1)
savefig("T200.pdf")
savefig("T200-raster.png")

#  @save "T40.jld2" T
#save("T100.jld2", T)

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
@testset "Single stage sequential tests" begin
    template_ed = get_template_nomin_ed_simulation()
    c_sys = PSB.build_system(PSITestSystems, "c_sys5_uc")
    problems = SimulationProblems(
        ED = OperationsProblem(template_ed, c_sys, optimizer = ipopt_optimizer),
    )
    test_sequence = SimulationSequence(
        problems = problems,
        intervals = Dict("ED" => (Hour(24), Consecutive())),
        ini_cond_chronology = InterProblemChronology(),
    )
    sim_single = Simulation(
        name = "consecutive",
        steps = 2,
        problems = problems,
        sequence = test_sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    build_out = build!(sim_single)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim_single)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "All stages executed - No Cache" begin
    duals = [:CopperPlateBalance]
    template_uc = get_template_basic_uc_simulation()
    template_ed = get_template_nomin_ed_simulation()
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    problems = SimulationProblems(
        UC = OperationsProblem(
            template_uc,
            c_sys5_hy_uc;
            optimizer = GLPK_optimizer,
            balance_slack_variables = true,
        ),
        ED = OperationsProblem(
            template_ed,
            c_sys5_hy_ed;
            optimizer = ipopt_optimizer,
            constraint_duals = duals,
            # Needed do to inconsistency in the test data
            balance_slack_variables = true,
        ),
    )

    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
            ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                variable_source_problem = PSI.ACTIVE_POWER,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    sim = Simulation(
        name = "no_cache",
        steps = 2,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )

    build_out = build!(sim; console_level = Logging.Info)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Simulation Single Stage with Cache" begin
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    template = get_template_hydro_st_ed()
    problems = SimulationProblems(
        ED = OperationsProblem(template, c_sys5_hy_ed; optimizer = ipopt_optimizer),
    )

    single_sequence = SimulationSequence(
        problems = problems,
        intervals = Dict("ED" => (Hour(1), Consecutive())),
        cache = Dict(("ED",) => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY)),
        ini_cond_chronology = IntraProblemChronology(),
    )

    sim_single_wcache = Simulation(
        name = "cache_st",
        steps = 2,
        problems = problems,
        sequence = single_sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    # Disabled due to https://github.com/NREL-SIIP/PowerSystemCaseBuilder.jl/issues/11
    #build_out = build!(sim_single_wcache)
    #@test build_out == PSI.BuildStatus.BUILT
    #execute_out = execute!(sim_single_wcache)
    #@test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Simulation with 2-Stages and Cache" begin
    template_uc = get_template_hydro_st_uc()
    template_ed = get_template_hydro_st_ed()
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    problems = SimulationProblems(
        UC = OperationsProblem(template_uc, c_sys5_hy_uc; optimizer = GLPK_optimizer),
        ED = OperationsProblem(template_ed, c_sys5_hy_ed; optimizer = GLPK_optimizer),
    )

    sequence_cache = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
            ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                variable_source_problem = PSI.ACTIVE_POWER,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(
            ("UC",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON),
            ("UC", "ED") => StoredEnergy(PSY.HydroEnergyReservoir, PSI.ENERGY),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    sim_cache = Simulation(
        name = "cache",
        steps = 2,
        problems = problems,
        sequence = sequence_cache,
        simulation_folder = mktempdir(cleanup = true),
    )
    build_out = build!(sim_cache)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim_cache)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Test Recedin Horizon Chronology" begin
    template_uc = get_template_basic_uc_simulation()
    template_ed = get_template_nomin_ed_simulation()
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    problems = SimulationProblems(
        UC = OperationsProblem(template_uc, c_sys5_hy_uc; optimizer = GLPK_optimizer),
        ED = OperationsProblem(
            template_ed,
            c_sys5_hy_ed;
            optimizer = ipopt_optimizer,
            # Added because of data issues
            balance_slack_variables = true,
        ),
    )

    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => RecedingHorizon()),
        intervals = Dict(
            "UC" => (Hour(24), RecedingHorizon()),
            "ED" => (Minute(60), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )

    sim = Simulation(
        name = "receding_horizon",
        steps = 2,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )
    build_out = build!(sim)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim)
    @test execute_out == PSI.RunStatus.SUCCESSFUL
end

@testset "Test Simulation Utils" begin
    template_uc = get_template_basic_uc_simulation()
    template_ed = get_template_nomin_ed_simulation()
    set_device_model!(template_ed, HydroEnergyReservoir, HydroDispatchReservoirBudget)
    c_sys5_hy_uc = PSB.build_system(PSITestSystems, "c_sys5_hy_uc")
    c_sys5_hy_ed = PSB.build_system(PSITestSystems, "c_sys5_hy_ed")
    problems = SimulationProblems(
        UC = OperationsProblem(
            template_uc,
            c_sys5_hy_uc;
            optimizer = GLPK_optimizer,
            constraint_duals = [:CopperPlateBalance],
        ),
        ED = OperationsProblem(
            template_ed,
            c_sys5_hy_ed;
            optimizer = ipopt_optimizer,
            # Added because of data issues
            balance_slack_variables = true,
            constraint_duals = [:CopperPlateBalance],
        ),
    )

    sequence = SimulationSequence(
        problems = problems,
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :ThermalStandard) => SemiContinuousFF(
                binary_source_problem = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
            ("ED", :devices, :HydroEnergyReservoir) => IntegralLimitFF(
                variable_source_problem = PSI.ACTIVE_POWER,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterProblemChronology(),
    )
    sim = Simulation(
        name = "aggregation",
        steps = 2,
        problems = problems,
        sequence = sequence,
        simulation_folder = mktempdir(cleanup = true),
    )

    build_out = build!(sim; console_level = Logging.Info)
    @test build_out == PSI.BuildStatus.BUILT
    execute_out = execute!(sim)
    @test execute_out == PSI.RunStatus.SUCCESSFUL

    @testset "Verify simulation events" begin
        file = joinpath(PSI.get_simulation_dir(sim), "recorder", "simulation.log")
        @test isfile(file)
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step = 1,
        )
        @test length(events) == 0
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step = 2,
        )
        @test length(events) == 10
        PSI.show_simulation_events(
            devnull,
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step = 2,
        )
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim);
            step = 1,
            problem = 1,
        )
        @test length(events) == 0
        events = PSI.list_simulation_events(
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step = 2,
            problem = 1,
        )
        @test length(events) == 10
        PSI.show_simulation_events(
            devnull,
            PSI.InitialConditionUpdateEvent,
            PSI.get_simulation_dir(sim),
            ;
            step = 2,
            problem = 1,
        )
    end

    # @testset "Check Serialization - Deserialization of Sim" begin
    #     path = mktempdir()
    #     files_path = PSI.serialize_simulation(sim; path = path)
    #     deserialized_sim = Simulation(files_path, stage_info)
    #     build_out = build!(deserialized_sim)
    #     @test build_out == PSI.BuildStatus.BUILT
    #     for stage in values(PSI.get_stages(deserialized_sim))
    #         @test PSI.is_stage_built(stage)
    #     end
    # end

    # @testset "Test print methods" begin
    #     list = [sim, sim.sequence, sim.stages["UC"]]
    #     _test_plain_print_methods(list)
    # end
end


text:
@testset "malicious html" begin
    # https://github.com/rgrove/sanitize/blob/master/test/test_malicious_html.rb

  @testset "comments" begin
    @testset "should not allow script injection via conditional comments" begin
      @test "" == HTMLSanitizer.sanitize("""""")
    end
  end

  @testset "" begin
    @testset "should not be possible to inject JS via a malformed event attribute" begin
      whitelist = deepcopy(HTMLSanitizer.WHITELIST)
      append!(whitelist[:elements], ["body", "head"])
      @test """""" == HTMLSanitizer.sanitize("""""", isfragment = false, whitelist = whitelist)
    end
  end

  @testset "" begin
    @testset "should not be possible to inject an iframe using an improperly closed tag" begin
      @test "" == HTMLSanitizer.sanitize("""<iframe src=http://ha.ckers.org/scriptlet.html <""")
    end
  end

  @testset "" begin
    @testset "should not be possible to inject JS via an unquoted  src attribute" begin
      @test "" == HTMLSanitizer.sanitize("""""")
    end

    @testset "should not be possible to inject JS using grave accents as  src delimiters" begin
      @test "" == HTMLSanitizer.sanitize("""<img src=`javascript:alert('XSS')`>""")
    end

    @testset "should not be possible to inject ">""")
    end

    @testset "should not be possible to inject protocol-based JS" begin
      @test "" == HTMLSanitizer.sanitize("""""")

      @test "" == HTMLSanitizer.sanitize("""""")

      @test "" == HTMLSanitizer.sanitize("""""")

      # Encoded tab character.
      @test "" == HTMLSanitizer.sanitize("""""")

      # Encoded newline.
      @test "" == HTMLSanitizer.sanitize("""""")

      # Encoded carriage return.
      @test "" == HTMLSanitizer.sanitize("""""")

      # Spaces plus meta char.
      @test "" == HTMLSanitizer.sanitize("""""")

      # Mixed spaces and tabs.
      @test "" == HTMLSanitizer.sanitize("""""")
    end

    @testset "should not be possible to inject protocol-based JS via whitespace" begin
      @test "" == HTMLSanitizer.sanitize("""""")
    end

    @testset "should not be possible to inject JS using a half-open  tag" begin
      @test """""" == HTMLSanitizer.sanitize("""<img src="javascript:alert('XSS')" """)
    end
  end
end


text:
#!/usr/bin/env julia

#Import Statements for MPC
import Pkg; Pkg.activate(@__DIR__);
Pkg.instantiate()
include("Arthur.jl")
using RigidBodyDynamics
using StaticArrays
using RobotDynamics
using Rotations
using LinearAlgebra
using ForwardDiff, FiniteDiff
using Altro
using TrajectoryOptimization

# Black magic to make Julia and Python happy
using Distributed

using Pkg
Distributed.@everywhere using Pkg

Distributed.@everywhere begin
  ENV["PYTHON"] = "/usr/bin/python3"
  Pkg.build("PyCall")
end

using PyCall
Distributed.@everywhere using PyCall

#Import Statements for RobotOS
using RobotOS

#For Testing Purposes
@rosimport geometry_msgs.msg: Point
rostypegen()
using .geometry_msgs.msg

# function callback(msg::Point,pub_obj::Publisher{Point})
#     pt_msg = Point(10,10,10)
#     publish(pub_obj,pt_msg)
# end

# function loop(pub_obj)
#     loop_rate = Rate(5.0)
#     while ! is_shutdown()
#         pt_msg = Point(rand(),rand(),rand())
#         publish(pub_obj,pt_msg)
#         rossleep(loop_rate)
#     end
# end

# function main()
#     init_node("julia2ros")
#     pub = Publisher{Point}("julia_out",queue_size=10)
#     sub = Subscriber{Point}("julia_in",callback,(pub,),queue_size=10)
#     loop(pub)
# end

# if !isinteractive()
#     main()
# end

@rosimport trajectory_msgs.msg: JointTrajectoryPoint
rostypegen()
using .trajectory_msgs.msg

JointTrajectoryOutput = JointTrajectoryPoint([1,1,1],[1,1,1],[1,1,1],[1,1,1],Duration(1))

function callback(msg::Point,pub_obj::Publisher{JointTrajectoryPoint})
    # JointTrajectory = msg.positions
    #Do magic here

    # JointTrajectoryOutput .= JointTrajectoryPoint([0,0,0],[0,0,0],[0,0,0],[0,0,0],Duration(0))

    publish(pub_obj,JointTrajectoryOutput)
end

function loop(sub_obj, pub_obj)
    loop_rate = Rate(5.0)
    while ! is_shutdown()
        # publish(pub_obj,JointTrajectoryOutput)
        RobotOS._run_callbacks(sub_obj)
        rossleep(loop_rate)
    end
end

function main()
    init_node("julia2ros")
    pub = Publisher{JointTrajectoryPoint}("julia_out",queue_size=10)
    sub = Subscriber{Point}("julia_in",callback,(pub,),queue_size=10)
    loop(sub, pub)
end

if !isinteractive()
    main()
end


text:
function dump_perf_buf(buf::IOBuffer)
    seekstart(buf)
    while !eof(buf)
        header = unpack(buf, perf_event_header)
        @show header.event_type
        @show header.misc
        println(PERF_EVENT_TYPE[header.event_type])
        pos = position(buf)
        if header.event_type == PERF_RECORD_ITRACE_START
            println(unpack(buf, itrace_start))
        elseif header.event_type == PERF_RECORD_AUX
            println(unpack(buf, aux_record))
        elseif header.event_type == PERF_RECORD_COMM
            println(unpack(buf, itrace_start))
            println(bytestring(read(buf, UInt8, header.size - sizeof(perf_event_header) - sizeof(itrace_start))))
        end
        seek(buf, pos + header.size - sizeof(perf_event_header))
    end
end


text:
@testset "Integral Data Generation (Univariate)" begin
    m = InfiniteModel();
    @infinite_parameter(m, t in [-Inf, Inf], supports = [0, .5, 1])
    @infinite_parameter(m, x[1:2] in [-Inf, Inf])
    # test _trapezoid_coeff
    @testset "_trapezoid_coeff" begin
        @test InfiniteOpt.MeasureToolbox._trapezoid_coeff([1,2,3]) == [0.5, 1., 0.5]
        @test_throws ErrorException InfiniteOpt.MeasureToolbox._trapezoid_coeff([1])
    end
    # test generate_integral_data (trapezoid)
    @testset "generate_integral_data (trapezoid)" begin
        @test generate_integral_data(t, 0, 1, UniTrapezoid()) ==
                    FunctionalDiscreteMeasureData(t, IOMT._trapezoid_coeff, 0, All,
                               NoGenerativeSupports(),default_weight, 0, 1, false)
    end
    # test _ensure_independent_param
    @testset "_ensure_independent_param" begin
        @test IOMT._ensure_independent_param(t, Automatic()) isa Nothing
        @test_throws ErrorException IOMT._ensure_independent_param(x[1], Automatic())
    end
    # test generate_integral_data (Gauss-Legendre)
    @testset "generate_integral_data (Gauss-Legendre)" begin
        @test generate_integral_data(t, 0, 1, GaussLegendre()) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, 1, GaussLegendre()).coefficients == FastGaussQuadrature.gausslegendre(3)[2] * 1/2
    end
    # test generate_integral_data (Gauss-Lobatto)
    @testset "generate_integral_data (Gauss-Lobato)" begin
        @test generate_integral_data(t, 0, 1, GaussLobatto()) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, 1, GaussLobatto()).coefficients == FastGaussQuadrature.gausslobatto(3)[2] * 1/2
    end
    # test _lobatto_coeff
    @testset "_lobatto_coeff" begin
        @test InfiniteOpt.MeasureToolbox._lobatto_coeff([1,2], 2) == [0.5, 0.5]
        @test_throws ErrorException InfiniteOpt.MeasureToolbox._lobatto_coeff([1], 3)
    end
    # test generate_integral_data (FEGaussLobatto)
    @testset "generate_integral_data (FEGaussLobatto)" begin
        @test_throws ErrorException generate_integral_data(t, 0, Inf, FEGaussLobatto())
        @test_throws ErrorException generate_integral_data(t, 0, 1, FEGaussLobatto(), num_nodes = 1)
        @test_throws ErrorException generate_integral_data(x[1], 0, 1, FEGaussLobatto(), num_nodes = 3)
        @test generate_integral_data(t, 0, 1, FEGaussLobatto()) isa FunctionalDiscreteMeasureData
        @test generate_integral_data(t, 0, 1, FEGaussLobatto()).label == All
        expected = [0.08333333333333333, 0.3333333333333333, 0.16666666666666666,
                    0.3333333333333333, 0.08333333333333333]
        @test all(abs.(generate_integral_data(t, 0, 1, FEGaussLobatto()).coeff_function([0, 0.5, 1]) .- expected) .< 0.00005)
        @variable(m, flob, Infinite(t))
        integral(flob, t, 0, 1, eval_method = FEGaussLobatto(), num_nodes = 3)
        add_generative_supports(t)
        lobb_supps = supports(t, label = All)
        @test all(abs.(lobb_supps .- [0, .25, .5, .75, 1] .< .000005))
    end
    # test generate_integral_data (Gauss-Radau)
    @testset "generate_integral_data (Gauss-Radau)" begin
        @test generate_integral_data(t, 0, 1, GaussRadau()) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, 1, GaussRadau()).coefficients == FastGaussQuadrature.gaussradau(3)[2] * 1/2
    end
    # test generate_integral_data (Gauss-Jacobi)
    @testset "generate_integral_data (Gauss-Jacobi)" begin
        @test generate_integral_data(t, 0, 1, GaussJacobi(5.0, 4.0)) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, 1, GaussJacobi(5.0, 4.0)).coefficients == FastGaussQuadrature.gaussjacobi(3, 5.0, 4.0)[2] * 1/2
        @test_throws ErrorException GaussJacobi(-5, -5)
    end
     # test generate_integral_data (Gauss-Chebyshev)
    @testset "generate_integral_data (Gauss-Chebyshev)" begin
        @test generate_integral_data(t, 0, 1, GaussChebyshev(3)) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, 1, GaussChebyshev(3)).coefficients == FastGaussQuadrature.gausschebyshev(3, 3)[2] * 1/2
        @test_throws ErrorException GaussChebyshev(5)
    end
    # test generate_integral_data (Gauss-Laguerre)
    @testset "generate_integral_data (Gauss-Laguerre)" begin
        @test generate_integral_data(t, -Inf, 0, GaussLaguerre()) isa DiscreteMeasureData
        @test generate_integral_data(t, 0, Inf, GaussLaguerre()) isa DiscreteMeasureData
    end
    # test generate_integral_data (Gauss-Hermite)
    @testset "generate_integral_data (Gauss-Hermite)" begin
        @test generate_integral_data(t, -Inf, Inf, GaussHermite()) isa DiscreteMeasureData
    end
    # test generate_integral_data (Quadrature)
    @testset "generate_integral_data (Quadrature)" begin
        @test generate_integral_data(t, 0, 1, Quadrature()) isa FunctionalDiscreteMeasureData
        @test generate_integral_data(t, -Inf, 0, Quadrature()) isa DiscreteMeasureData
        @test generate_integral_data(t, -Inf, Inf, Quadrature()) isa DiscreteMeasureData
    end
    # test handling of inappropriate quadrature method use
    @testset "generate_integral_data (improper quadrature method usage)" begin
        warn = "The `GaussLegendre` method can only be applied on finite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussLegendre())) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussLegendre())) isa DiscreteMeasureData
        warn = "The `GaussLobatto` method can only be applied on finite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussLobatto())) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussLobatto())) isa DiscreteMeasureData
        warn = "The `GaussRadau` method can only be applied on finite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussRadau())) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussRadau())) isa DiscreteMeasureData
        warn = "The `GaussJacobi` method can only be applied on finite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussJacobi(5.0, 4.0))) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussJacobi(5.0, 4.0))) isa DiscreteMeasureData
        warn = "The `GaussChebyshev` method can only be applied on finite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussChebyshev(3))) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussChebyshev(3))) isa DiscreteMeasureData
        warn = "`GaussLaguerre` quadrature can only be applied on semi-infinite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, 0, 1, GaussLaguerre())) isa DiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, Inf, GaussLaguerre())) isa DiscreteMeasureData
        warn = "`GaussHermite` quadrature can only be applied on infinite intervals, " *
               "switching to an appropriate method."
        @test (@test_logs (:warn, warn) generate_integral_data(t, 0, 1, GaussHermite())) isa FunctionalDiscreteMeasureData
        @test (@test_logs (:warn, warn) generate_integral_data(t, -Inf, 0, GaussHermite())) isa DiscreteMeasureData
    end
    # test generate_integral_data (uniform Monte Carlo)
    @testset "generate_integral_data (UniMCSampling())" begin
        @test generate_integral_data(t, 0, 1, UniMCSampling()) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException generate_integral_data(t, -Inf, 0, UniMCSampling())
        warn = "Cannot specify a nonzero `num_supports` when making an " *
               "integral with repsect to an individual dependent parameter. " *
               "Setting `num_supports = 0`."
        @test_logs (:warn, warn) generate_integral_data(x[1], 0, 1, UniMCSampling(), num_supports = 1)
    end
    # test generate_integral_data (independent Monte Carlo)
    @testset "generate_integral_data (UniIndepMCSampling())" begin
        @test generate_integral_data(t, 0, 1, UniIndepMCSampling()) isa DiscreteMeasureData
        @test_throws ErrorException generate_integral_data(t, -Inf, 0, UniIndepMCSampling())
        @test_throws ErrorException generate_integral_data(t, -Inf, Inf, UniIndepMCSampling())
    end
    # test generate_integral_data (Automatic)
    @testset "generate_integral_data (Automatic)" begin
        @test generate_integral_data(t, 0, 1, Automatic()) isa FunctionalDiscreteMeasureData
        @test generate_integral_data(t, -Inf, 1, Automatic()) isa DiscreteMeasureData
        @test generate_integral_data(t, -Inf, Inf, Automatic()) isa DiscreteMeasureData
        @test generate_integral_data(x[1], 0, 1, Automatic()) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException generate_integral_data(x[1], -Inf, 0, Automatic())
        @test_throws ErrorException generate_integral_data(x[1], -Inf, Inf, Automatic())
    end
    # test fallback
    @testset "generate_integral_data (fallback)" begin
        @test_throws ErrorException generate_integral_data(x, 1, 2, Quadrature())
    end
end

@testset "Integral Data Generation (Multivariate)" begin
    m = InfiniteModel();
    @infinite_parameter(m, x[1:2] in [-Inf, Inf])
    # test generate_integral_data (MultiMCSampling)
    @testset "generate_integral_data (MultiMCSampling)" begin
        @test generate_integral_data(x, [0, 0], [1, 1], MultiMCSampling()) isa FunctionalDiscreteMeasureData
        @test_throws ErrorException generate_integral_data(x, [0, 0], [Inf, 1], MultiMCSampling())
        @test_throws ErrorException("`MultiMCSampling` is not applicable to (semi-)infinite intervals.") generate_integral_data(x, [0, -Inf], [1, 1], MultiMCSampling())
    end
    # test _make_multi_mc_supports
    # @testset "_make_multi_mc_supports" begin
    #     dx = dispatch_variable_ref.(x)
    #     @test IOMT._make_multi_mc_supports(dx, [0, 0], [1, 1], 5) isa Matrix{Float64}
    #     @test size(IOMT._make_multi_mc_supports(dx, [0, 0], [1, 1], 5)) == (2, 5)
    # end
    # test generate_integral_data (MultiIndepMCSampling)
    @testset "generate_integral_data (MultiIndepMCSampling)" begin
        @test generate_integral_data(x, [0, 0], [1, 1], MultiIndepMCSampling()) isa DiscreteMeasureData
        @test_throws ErrorException generate_integral_data(x, [0, 0], [Inf, 1], MultiIndepMCSampling())
        @test_throws ErrorException generate_integral_data(x, [0, -Inf], [1, 1], MultiIndepMCSampling())
    end
    # test generate_integral_data (Automatic)
    @testset "generate_integral_data (Automatic)" begin
        @test generate_integral_data(x, [0, 0], [1, 1], Automatic()) isa FunctionalDiscreteMeasureData
    end
end

@testset "Univariate Integrals" begin
    m = InfiniteModel();
    @infinite_parameter(m, t in [0, 10])
    @variable(m, inf1, Infinite(t))
    @testset "uni_integral_defaults" begin
        @test uni_integral_defaults() == IOMT.UniIntegralDefaults
    end
    @testset "[set/clear]_uni_integral_defaults" begin
        @test uni_integral_defaults() == Dict{Symbol, Any}(:eval_method => Automatic())
        @test set_uni_integral_defaults(num_nodes = 5, new_kwarg = true) isa Nothing
        @test uni_integral_defaults()[:num_nodes] == 5
        @test uni_integral_defaults()[:new_kwarg]
        @test clear_uni_integral_defaults() isa Nothing
        @test uni_integral_defaults() == Dict{Symbol, Any}(:eval_method => Automatic())
    end
    @testset "integral" begin
        @test InfiniteOpt._index_type(integral(inf1, t)) == MeasureIndex
        @test_throws ErrorException integral(inf1, t, -1, 1)
        @test_throws ErrorException integral(inf1, t, 9, 11)
        @test_throws ErrorException integral(inf1, t, 3, 1)
        @test name(integral(inf1, t)) == "integral"
    end
    @testset "∫" begin
        @test InfiniteOpt._index_type(∫(inf1, t)) == MeasureIndex
        @test_throws ErrorException ∫(inf1, t, -1, 1)
        @test_throws ErrorException ∫(inf1, t, 9, 11)
        @test_throws ErrorException ∫(inf1, t, 3, 1)
        @test name(∫(inf1, t)) == "integral"
    end
end

@testset "Multivariate Integrals" begin
    m = InfiniteModel();
    @infinite_parameter(m, x[1:2] in [0, 10])
    @infinite_parameter(m, xi[1:2] in [0, 5])
    @infinite_parameter(m, y[2:3] in [0, 1])
    @variable(m, inf1, Infinite(x))
    @variable(m, inf2, Infinite(xi))
    @variable(m, inf3, Infinite(y))
    @testset "multi_integral_defaults" begin
        @test multi_integral_defaults() == IOMT.MultiIntegralDefaults
    end
    @testset "[set/clear]_multi_integral_defaults" begin
        @test multi_integral_defaults() == Dict{Symbol, Any}(:eval_method => Automatic())
        @test set_multi_integral_defaults(num_supports = 5, new_kwarg = true) isa Nothing
        @test multi_integral_defaults()[:num_supports] == 5
        @test multi_integral_defaults()[:new_kwarg]
        @test clear_multi_integral_defaults() isa Nothing
        @test multi_integral_defaults() == Dict{Symbol, Any}(:eval_method => Automatic())
    end
    @testset "integral" begin
        @test InfiniteOpt._index_type(integral(inf1, x)) == MeasureIndex
        @test InfiniteOpt._index_type(integral(inf2, xi)) == MeasureIndex
        @test InfiniteOpt._index_type(integral(inf1, x, 1, 2)) == MeasureIndex
        @test InfiniteOpt._index_type(integral(inf2, xi, 1, 2)) == MeasureIndex
        @test_throws ErrorException integral(inf3, y, [0, 0], [1, 1])
        @test_throws ErrorException integral(inf1, x, 3, 1)
        @test_throws ErrorException integral(inf1, x, -1, 1)
        @test_throws ErrorException integral(inf1, x, 9, 11)
        @test name(integral(inf1, x)) == "integral"
    end
    @testset "∫" begin
        @test InfiniteOpt._index_type(∫(inf1, x)) == MeasureIndex
        @test InfiniteOpt._index_type(∫(inf2, xi)) == MeasureIndex
        @test InfiniteOpt._index_type(∫(inf1, x, 1, 2)) == MeasureIndex
        @test InfiniteOpt._index_type(∫(inf2, xi, 1, 2)) == MeasureIndex
        @test_throws ErrorException ∫(inf3, y, [0, 0], [1, 1])
        @test_throws ErrorException ∫(inf1, x, 3, 1)
        @test_throws ErrorException ∫(inf1, x, -1, 1)
        @test_throws ErrorException ∫(inf1, x, 9, 11)
        @test name(∫(inf1, x)) == "integral"
    end
end

@testset "Macro" begin
    m = InfiniteModel();
    @infinite_parameter(m, t in [0, 10])
    @infinite_parameter(m, x[1:2] in [0, 10])
    @variable(m, inf1, Infinite(t))
    @variable(m, inf2, Infinite(x))
    @testset "@integral" begin
        @test InfiniteOpt._index_type(@integral(inf1, t)) == MeasureIndex
        @test InfiniteOpt._index_type(@integral(inf2, x)) == MeasureIndex
        @test InfiniteOpt._index_type(@integral(inf2, t)) == MeasureIndex
        @test InfiniteOpt._index_type(@integral(inf1, t, 3, 7)) == MeasureIndex
        @test InfiniteOpt._index_type(@integral(inf2, x, [0, 2], [3, 5])) == MeasureIndex
        @test_macro_throws ErrorException @integral(inf1, t, 1)
    end
    @testset "@∫" begin
        @test InfiniteOpt._index_type(@∫(inf1, t)) == MeasureIndex
        @test InfiniteOpt._index_type(@∫(inf2, x)) == MeasureIndex
        @test InfiniteOpt._index_type(@∫(inf2, t)) == MeasureIndex
        @test InfiniteOpt._index_type(@∫(inf1, t, 3, 7)) == MeasureIndex
        @test InfiniteOpt._index_type(@∫(inf2, x, [0, 2], [3, 5])) == MeasureIndex
        @test_macro_throws ErrorException @∫(inf1, t, 1)
    end
end


text:
# ps3000_set_trigger.jl

export ps3000_set_trigger, ps3000_set_trigger2, PS3000Source, PS3000Direction

function ps3000_set_trigger(handle::Int16, source, threshold, direction,
                             delay::Int, auto_trigger_ms)
  ps_status = ccall((:ps3000_set_trigger, ps3000driver), Int16,
    (Int16, Int16, Int16, Int16, Int16, Int16),
    handle, source, threshold, direction, delay, auto_trigger_ms)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end

function ps3000_set_trigger(handle::Int16, source, threshold, direction,
                             delay, auto_trigger_ms)
  ps_status = ccall((:ps3000_set_trigger2, ps3000driver), Int16,
    (Int16, Int16, Int16, Int16, Cfloat, Int16),
    handle, source, threshold, direction, delay, auto_trigger_ms)
  if ps_status == 0
    error("ps3000 Error: one of the parameters is out of range")
  end
  return nothing
end

function ps3000_set_trigger2(handle::Int16, source, threshold, direction,
                             delay, auto_trigger_ms)
  ps3000_set_trigger(handle::Int16, source, threshold, direction,
                     delay, auto_trigger_ms)
end

baremodule PS3000Source
  const CHANNEL_A = 0
  const CHANNEL_B = 1
  const CHANNEL_C = 2
  const CHANNEL_D = 3
  const EXTERNAL  = 4
  const NONE      = 5
end # module

baremodule PS3000Direction
  const RISING = 0
  const FALLING = 1
end


text:
struct GQL end

Base.convert(::Type{String}, t::RBNF.Token) = t.str
Base.convert(::Type{Symbol}, t::RBNF.Token{:name}) = Symbol(t.str)
Base.convert(::Type{Symbol}, t::RBNF.Token{:reserved}) = Symbol(t.str)
Base.convert(::Type{Int}, t::RBNF.Token{:int_value}) = Base.parse(Int, t.str)
Base.convert(::Type{Float64}, t::RBNF.Token{:float_value}) = Base.parse(Float64, t.str)
function Base.convert(::Type{String}, t::RBNF.Token{:str})
    startswith(t.str, "\"\"\"") ? String(t.str[4:end-3]) : String(t.str[2:end-1])
end
function Base.convert(::Type{String}, t::RBNF.Token{:single_quote_string_value})
    string(t.str[2:end-1])
end
function Base.convert(::Type{String}, t::RBNF.Token{:triple_quote_string_value})
    string(t.str[4:end-3])
end

RBNF.@parser GQL begin
    ignore{whitespace,comment,comma}

    @grammar
    atom = (name | float_value | int_value | string_value | punctuator)

    document::Document := [definitions = definition{*}]
    definition = (executable_definition | type_system_definition | type_system_extension)

    executable_definition = (operation_definition | fragment_definition)

    operation_definition = (
        selection_set |
        [operation_type, name.?, variable_definitions.?, directives.?, selection_set]
    )

    operation_type = ("query" | "mutation" | "subscription")

    selection_set::SelectionSet := ['{', selections = selection{*}, '}']

    selection = (field | fragment_spread | inline_fragment)

    field::Field := [
        alias = alias.?,
        name = name,
        arguments = arguments.?,
        directives = directives.?,
        selection_set = selection_set.?,
    ]

    alias::Alias := [name = name, ':']

    arguments::Arguments := ['(', arguments = argument{*}, ')']
    argument::Argument := [name = name, ':', value = value]

    fragment_spread::FragmentSpread :=
        ["...", fragment_name = fragment_name, directives = directives.?]
    inline_fragment::InlineFragment := [
        "...",
        type_condition = type_condition.?,
        directives = directives.?,
        selection_set = selection_set,
    ]
    fragment_definition::FragmentDefinition := [
        "fragment",
        fragment_name = fragment_name,
        type_condition = type_condition,
        directives = directives.?,
        selection_set = selection_set,
    ]

    # but not on
    fragment_name = name
    type_condition::TypeCondition := ["on", named_type = named_type]

    value = (
        variable |
        int_value |
        float_value |
        string_value |
        bool_value |
        null_value |
        enum_value |
        list_value |
        object_value
    )

    bool_value = ("true" | "false")
    null_value = "null"

    # but not true, false, or null
    enum_value = name

    list_value::ListValue := ['[', values = value{*}.?, ']']
    object_value::ObjectValue := ['{', object_fields = object_field{*}.?, '}']

    object_field::ObjectField := [name = name, ':', value = value]

    variable_definitions::VariableDefinitions :=
        ['(', variable_definitions = variable_definition{*}, ')']
    variable_definition::VariableDefinition :=
        [name = variable, ':', type = typ, value = default_value.?]
    variable::Variable := ['$', name = name]
    default_value::DefaultValue := ['=', value = value]

    typ = (list_type | named_type)
    named_type::NamedType := [type = name, non_null = '!'.?]
    list_type::ListType := ['[', type = typ, ']', non_null = '!'.?]

    directives::Directives := [directives = directive{*}]
    directive::Directive := ['@', name = name, arguments = arguments.?]

    type_system_definition = (schema_definition | type_definition | directive_definition)

    type_system_extension = (schema_extension | type_extension)

    schema_definition::SchemaDefinition := [
        "schema",
        directives = directives.?,
        '{',
        operation_type_definitions = operation_type_definition{*},
        '}',
    ]
    schema_extension = (
        ["extend", "schema", directives.?, '{', operation_type_definition{*}, '}'] |
        ["extend", "schema", directives]
    )

    operation_type_definition::OperationTypeDefinition :=
        [operation_type = operation_type, ':', named_type = named_type]

    type_definition::TypeDefinition := [
        description = string_value.?,
        type = (
            scalar_type_definition |
            object_type_definition |
            interface_type_definition |
            union_type_definition |
            enum_type_definition |
            input_object_type_definition
        ),
    ]

    # TODO: extensions might not work
    type_extension = (
        scalar_type_extension |
        object_type_extension |
        interface_type_extension |
        union_type_extension |
        enum_type_extension |
        input_object_type_extension
    )

    scalar_type_definition::ScalarTypeDefinition :=
        ["scalar", name = name, directives = directives.?]
    scalar_type_extension::ScalarTypeExtension :=
        ["extend", "scalar", name = name, directives = directives]

    object_type_definition::ObjectTypeDefinition := [
        "type",
        name = name,
        implements_interfaces = implements_interfaces.?,
        directives = directives.?,
        fields_definition = fields_definition.?,
    ]
    object_type_extension = (
        ["extend", "type", name, implements_interfaces.?, directives.?, fields_definition] |
        ["extend", "type", name, implements_interfaces.?, directives] |
        ["extend", "type", name, implements_interfaces]
    )

    # TODO: recursion
    implements_interfaces = ["implements", '&'.?, type = named_type]

    # implements_interfaces = (
    #                          ["implements", '&'.?, named_type] |
    #                          [implements_interfaces, '&', named_type]
    #                         )

    fields_definition::FieldsDefinition :=
        ['{', field_definitions = field_definition{*}, '}']

    field_definition::FieldDefinition := [
        description = string_value.?,
        name = name,
        arguments_definition = arguments_definition.?,
        ':',
        type = typ,
        directives = directives.?,
    ]

    arguments_definition::ArgumentsDefinition :=
        ['(', input_value_definitions = input_value_definition{*}, ')']

    input_value_definition::InputValueDefinition := [
        description = string_value.?,
        name = name,
        ':',
        type = typ,
        default_value = default_value.?,
        directives = directives.?,
    ]

    interface_type_definition::InterfaceTypeDefinition := [
        "interface",
        name = name,
        directives = directives.?,
        fields_definition = fields_definition.?,
    ]

    interface_type_extension = (
        ["extend", "interface", name, directives.?, fields_definition] |
        ["extend", "interface", name, directives]
    )

    union_type_definition::UnionTypeDefinition :=
        ["union", name = name, directives = directives.?, types = union_member_types.?]

    union_member_types = @direct_recur begin
        init = [['=', '|'.?, named_type] % (x -> x[3])]
        prefix = [recur..., ('|', named_type) % (x -> x[2])]
    end

    union_type_extension = (
        ["extend", "union", name, directives.?, union_member_types] |
        ["extend", "union", name, directives]
    )

    enum_type_definition::EnumTypeDefinition :=
        ["enum", name = name, directives = directives.?, enums = enum_values_definition.?]

    enum_values_definition::EnumValuesDefinition :=
        ['{', values = enum_value_definition{*}, '}']
    enum_value_definition::EnumValueDefinition :=
        [description = string_value.?, value = enum_value, directives = directives.?]
    enum_type_extension = (
        ["extend", "enum", name, directives.?, enum_values_definition] |
        ["extend", "enum", name, directives]
    )

    input_object_type_definition::InputObjectTypeDefinition := [
        "input",
        name = name,
        directives = directives.?,
        fields = input_fields_definition.?,
    ]

    input_fields_definition::InputFieldsDefinition :=
        ['{', input_value_definitions = input_value_definition{*}, '}']

    input_object_type_extension = (
        ["extend", "input", name, directives.?, input_fields_definition] |
        ["extend", "input", name, directives]
    )

    directive_definition::DirectiveDefinition := [
        description = string_value.?,
        "directive",
        '@',
        name = name,
        arguments_definition = arguments_definition.?,
        "on",
        directive_locations = directive_locations,
    ]

    # TODO: recursion
    directive_locations =
        (['|'.?, directive_location] | [directive_locations, '|', directive_location])

    directive_location = (executable_directive_location | type_system_directive_location)

    executable_directive_location = (
        "QUERY" |
        "MUTATION" |
        "SUBSCRIPTION" |
        "FIELD" |
        "FRAGMENT_DEFINITION" |
        "FRAGMENT_SPREAD" |
        "INLINE_FRAGMENT"
    )
    type_system_directive_location = (
        "SCHEMA" |
        "SCALAR" |
        "OBJECT" |
        "FIELD_DESCRIPTION" |
        "ARGUMENT_DEFINITION" |
        "INTERFACE" |
        "UNION" |
        "ENUM" |
        "ENUM_VALUE" |
        "INPUT_OBJECT" |
        "INPUT_FIELD_DEFINITION"
    )

    string_value = triple_quote_string_value | single_quote_string_value

    @token
    comment := r"\G#.*(\r|\n)?"
    whitespace := r"\G\s+"
    comma := r"\G,"

    name := r"\G[A-Za-z][_0-9A-Za-z]*"
    float_value := r"\G([0-9]+\.[0-9]*|[0-9]*\.[0.9]+)([eE][-+]?[0-9]+)?"
    int_value := r"\G([1-9]+[0-9]*|0)"
    triple_quote_string_value := @quote ("\"\"\"", "\\\"\"\"", "\"\"\"")
    single_quote_string_value := @quote ("\"", "\\\"", "\"")

    punctuator := r"\G(!|\$|\(|\)|\.\.\.|:|=|@|\[|\]|{|}|\|)"
end


text:
# include model data -
include("SolveBalances.jl")

using PyCall
@pyimport numpy as np

using Debug

# some global parameters -
BIG = 1e10
SMALL = 1e-6

# Globally load data so we don't have load on each iteration -
MEASURED_ARRAY_C3a_O1 = readdlm("./data/Shaw2015_Fig2a_C3a.txt",',')
MEASURED_ARRAY_C5a_O1 = readdlm("./data/Shaw2015_Fig3ai_C5a.txt",',')

MEASURED_ARRAY_C3a_O2 = readdlm("./data/Shaw2015_Fig2e_C3a.txt",',')
MEASURED_ARRAY_C5a_O2 = readdlm("./data/Shaw2015_Fig3c_C5a_original.txt",',')

MEASURED_ARRAY_C3a_O3 = readdlm("./data/Shaw2015_Fig2c_C3a.txt",',')
MEASURED_ARRAY_C5a_O3 = readdlm("./data/Shaw2015_Fig3aiii_C5a_original.txt",',')

function calculate_error_obj_1(T,X)

  obj_array = BIG*ones(2)

  # Need to interpolate the simulation onto the experimental time scale -
  number_of_measurements = 20
  IC3a = np.interp(MEASURED_ARRAY_C3a_O1[:,1],T,X[:,19])
  IC5a = np.interp(MEASURED_ARRAY_C5a_O1[:,1],T,X[:,15])

  # Compute the error values -
  MC3a = 0.11002*MEASURED_ARRAY_C3a_O1[:,2]
  MC5a = 0.00009615*MEASURED_ARRAY_C5a_O1[:,2]

  # Compute the scaled IC3a and IC5a -
  min_IC3a = minimum(IC3a)
  max_IC3a = maximum(IC3a)
  min_IC5a = minimum(IC5a)
  max_IC5a = maximum(IC5a)

  min_MC3a = minimum(MC3a)
  max_MC3a = maximum(MC3a)
  min_MC5a = minimum(MC5a)
  max_MC5a = maximum(MC5a)

  shape_C3a = (IC3a - min_IC3a)/(max_IC3a - min_IC3a)
  shape_C5a = (IC5a - min_IC5a)/(max_IC5a - min_IC5a)
  shape_MC3a = (MC3a - min_MC3a)/(max_MC3a - min_MC3a)
  shape_MC5a = (MC5a - min_MC5a)/(max_MC5a - min_MC5a)

  # Compute the scale error -
  scale_error_C3a = (max_IC3a - max_MC3a)/(max_MC3a)
  scale_error_C5a = (max_IC5a - max_MC5a)/(max_MC5a)

  error_C3a = sum((shape_C3a - shape_MC3a).^2) + 10*(scale_error_C3a)^2
  error_C5a = sum((shape_C5a - shape_MC5a).^2) + 10*(scale_error_C5a)^2

  # Compute the objective array -
  obj_array[1] = error_C3a
  obj_array[2] = error_C5a

  return (sum(obj_array))
end

function calculate_error_obj_2(T,X)

  obj_array = BIG*ones(2)

  # Need to interpolate the simulation onto the experimental time scale -
  number_of_measurements = 20
  IC3a = np.interp(MEASURED_ARRAY_C3a_O2[:,1],T,X[:,19])
  IC5a = np.interp(MEASURED_ARRAY_C5a_O2[:,1],T,X[:,15])

  # Compute the error values -
  MC3a = 0.11002*MEASURED_ARRAY_C3a_O2[:,2]
  MC5a = 0.00009615*MEASURED_ARRAY_C5a_O2[:,2]

  # error_C3a = sum((MC3a - IC3a).^2)
  # error_C5a = sum((MC5a - IC5a).^2)

  # Compute the scaled IC3a and IC5a -
  min_IC3a = minimum(IC3a)
  max_IC3a = maximum(IC3a)
  min_IC5a = minimum(IC5a)
  max_IC5a = maximum(IC5a)

  min_MC3a = minimum(MC3a)
  max_MC3a = maximum(MC3a)
  min_MC5a = minimum(MC5a)
  max_MC5a = maximum(MC5a)

  shape_C3a = (IC3a - min_IC3a)/(max_IC3a - min_IC3a)
  shape_C5a = (IC5a - min_IC5a)/(max_IC5a - min_IC5a)
  shape_MC3a = (MC3a - min_MC3a)/(max_MC3a - min_MC3a)
  shape_MC5a = (MC5a - min_MC5a)/(max_MC5a - min_MC5a)

  # Compute the scale error -
  scale_error_C3a = (max_IC3a - max_MC3a)/(max_MC3a)
  scale_error_C5a = (max_IC5a - max_MC5a)/(max_MC5a)

  error_C3a = sum((shape_C3a - shape_MC3a).^2) + 10*(scale_error_C3a)^2
  error_C5a = sum((shape_C5a - shape_MC5a).^2) + 10*(scale_error_C5a)^2

  # Compute the objective array -
  obj_array[1] = error_C3a
  obj_array[2] = error_C5a

  return (sum(obj_array))
end

function calculate_error_obj_3(T,X)

  obj_array = BIG*ones(2)

  # Need to interpolate the simulation onto the experimental time scale -
  number_of_measurements = 20
  IC3a = np.interp(MEASURED_ARRAY_C3a_O3[:,1],T,X[:,19])
  IC5a = np.interp(MEASURED_ARRAY_C5a_O3[:,1],T,X[:,15])

  # Compute the error values -
  MC3a = 0.11002*MEASURED_ARRAY_C3a_O3[:,2]
  MC5a = 0.00009615*MEASURED_ARRAY_C5a_O3[:,2]
  # error_C3a = sum((MC3a - IC3a).^2)
  # error_C5a = sum((MC5a - IC5a).^2)

  # Compute the scaled IC3a and IC5a -
  min_IC3a = minimum(IC3a)
  max_IC3a = maximum(IC3a)
  min_IC5a = minimum(IC5a)
  max_IC5a = maximum(IC5a)

  min_MC3a = minimum(MC3a)
  max_MC3a = maximum(MC3a)
  min_MC5a = minimum(MC5a)
  max_MC5a = maximum(MC5a)

  shape_C3a = (IC3a - min_IC3a)/(max_IC3a - min_IC3a)
  shape_C5a = (IC5a - min_IC5a)/(max_IC5a - min_IC5a)
  shape_MC3a = (MC3a - min_MC3a)/(max_MC3a - min_MC3a)
  shape_MC5a = (MC5a - min_MC5a)/(max_MC5a - min_MC5a)

  # Compute the scale error -
  scale_error_C3a = (max_IC3a - max_MC3a)/(max_MC3a)
  scale_error_C5a = (max_IC5a - max_MC5a)/(max_MC5a)

  error_C3a = sum((shape_C3a - shape_MC3a).^2) + (scale_error_C3a)^2
  error_C5a = sum((shape_C5a - shape_MC5a).^2) + (scale_error_C5a)^2

  # Compute the objective array -
  obj_array[1] = error_C3a
  obj_array[2] = error_C5a

  return (sum(obj_array))
end

function objective_function(parameter_array)

  # Setup the timescale -
  tStart = 0.0;
  tStop = 25.0;
  tStep = 0.1;

  # Calculate the objective function array -
  obj_array = BIG*ones(2,1)

  # Setup data dictionary -
  data_dictionary = Dict()
  data_dictionary["PARAMETER_ARRAY"] = parameter_array

  # Experiment 1 = w/o zymosan (run with default ICs)
  initial_condition_array = zeros(19)
  initial_condition_array[1] = 0.0      # 1 zymosan
  initial_condition_array[2] = 1.9      # 2 C4
  initial_condition_array[3] = 0.322    # 3 C2
  initial_condition_array[8] = 7.57     # 4 C3
  initial_condition_array[14] = 0.195   # 5 C5
  initial_condition_array[17] = 2.23    # 6 Factor H
  initial_condition_array[18] = 0.417   # 7 C4BP
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  (T1,X1) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # Experiment 2 = w 1g zymosan -
  initial_condition_array[1] = 1.0      # 1 zymosan
  data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  (T2,X2) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # # Experiment 3 = w 0.01g zymosan -
  # initial_condition_array[1] = 0.01      # 1 zymosan
  # data_dictionary["INITIAL_CONDITION_ARRAY"] = initial_condition_array;
  # (T3,X3) = SolveBalances(tStart,tStop,tStep,data_dictionary)

  # call the experiment functions -
  obj_array[1] = calculate_error_obj_1(T1,X1)
  obj_array[2] = calculate_error_obj_2(T2,X2)
  #obj_array[3] = calculate_error_obj_3(T3,X3)

  # return -
  return obj_array
end

# --------------------------------------- DO NOT EDIT BELOW THIS LINE  --------------------------------------- #
# Evaluates the objective function values -
function local_refienment_step(parameter_array)

  SIGMA = 0.05

  # initialize -
  number_of_parameters = length(parameter_array)

  # Setup the bound constraints -
  LOWER_BOUND = SMALL
  UPPER_BOUND = 10*ones(number_of_parameters)

  # calculate the starting error -
  parameter_array_best = parameter_array
  error_array = BIG*ones(4)
  error_array[1] = sum(objective_function(parameter_array_best))

  # main refinement loop -
  iteration_counter = 1
  iteration_max = 1000
  while (iteration_counter<iteration_max)

    # take a step up -
    parameter_up = parameter_array_best.*(1+SIGMA*rand(number_of_parameters))
    parameter_up = parameter_bounds_function(parameter_up,LOWER_BOUND*ones(number_of_parameters),UPPER_BOUND)

    # take a step down -
    parameter_down = parameter_array_best.*(1-SIGMA*rand(number_of_parameters))
    parameter_down = parameter_bounds_function(parameter_down,LOWER_BOUND*ones(number_of_parameters),UPPER_BOUND)

    # Evaluate the obj function -
    error_array[2] = sum(objective_function(parameter_up))
    error_array[3] = sum(objective_function(parameter_down))

    # Calculate a correction factor -
    a = error_array[2]+error_array[3] - 2.0*error_array[1]
    parameter_corrected = parameter_array_best
    if (a>0.0)
      amda = -0.5*(error_array[3] - error_array[2])/a
      parameter_corrected = parameter_array_best+amda*rand(number_of_parameters)
      parameter_corrected = parameter_bounds_function(parameter_corrected,LOWER_BOUND*ones(number_of_parameters),UPPER_BOUND)
      error_array[4] = sum(objective_function(parameter_corrected))
    end

    # Which step has the min error?
    min_index = indmin(error_array)
    if (min_index == 1)
			parameter_array_best = parameter_array_best
	 	elseif (min_index == 2)
		 	parameter_array_best = parameter_up
	 	elseif (min_index == 3)
		 	parameter_array_best = parameter_down
	 	elseif (min_index == 4)
		 parameter_array_best = parameter_corrected
	 	end

		# Update the local error
		error_array[1] = error_array[min_index]

    @show iteration_counter,error_array[min_index]

    # update local counter -
    iteration_counter = iteration_counter + 1
  end

  return parameter_array_best
end

function neighbor_function(parameter_array)

  SIGMA = 0.05
  number_of_parameters = length(parameter_array)

  # calculate new parameters -
  new_parameter_array = parameter_array.*(1+SIGMA*randn(number_of_parameters))

  # Check the bound constraints -
  LOWER_BOUND = SMALL
  UPPER_BOUND = [
        50.0  ;          # k_Initiator2_C4             = 0.1*parameter_array[1];
        15.0  ;          # Km_Initiator2_C4            = 1.0*parameter_array[2];
        50.0  ;          # k_Initiator2_C2             = 0.1*parameter_array[3];
        15.0  ;          # Km_Initiator2_C2            = 1.0*parameter_array[4];
        2.0   ;          # n_In2_C4                    = parameter_array[5];
        2.0   ;          # n_In2_C2                    = parameter_array[6];
        15.0  ;          # alpha_In2_C4                = parameter_array[7];
        15.0  ;          # alpha_In2_C2                = parameter_array[8];
        2.0   ;          # n_C5a                       = parameter_array[9];
        2.0   ;          # n_C3a                       = parameter_array[10];
        15.0  ;          # k_C4bC2a                    = parameter_array[11];
        15.0  ;          # k_C3b_basal                 = parameter_array[12];
        15.0  ;          # k_C3Convertase2             = parameter_array[13];
        15.0  ;          # k_C5Convertase2             = parameter_array[14];
        15.0  ;          # k_C5_C4bC2a_bind            = parameter_array[15];
        15.0  ;          # k_C5_conversion             = parameter_array[16];
        15.0  ;          # Km_C5_conversion            = parameter_array[17];
        15.0  ;          # k_C5_conversion_alternate   = parameter_array[18];
        15.0  ;          # Km_C5_conversion_alternate  = parameter_array[19];
        15.0  ;          # k_inhibit_C4bC2a            = parameter_array[20];
        15.0  ;          # k_inhibit_C3Convertase2     = parameter_array[21];
        15.0  ;          # k_cat_C4bC2a                = parameter_array[22];
        15.0  ;          # Km_C4bC2a                   = parameter_array[23];
        15.0  ;          # k_cat_C3convertase2         = parameter_array[24];
        15.0  ;          # Km_C3Convertase2            = parameter_array[25];
        1.0   ;          # k_degradationC3a            = parameter_array[26];
        15.0  ;          # kC5inhibit                  = parameter_array[27];
        1.0   ;          # k_degradationC5a            = parameter_array[28];
  ];

  # return the corrected parameter arrays -
  return parameter_bounds_function(new_parameter_array,LOWER_BOUND*ones(number_of_parameters),UPPER_BOUND)
end

function acceptance_probability_function(rank_array,temperature)
  return (exp(-rank_array[end]/temperature))
end

function cooling_function(temperature)

  # define my new temperature -
  alpha = 0.90
  return alpha*temperature
end

function parameter_bounds_function(parameter_array,lower_bound_array,upper_bound_array)

  # reflection_factor -
  epsilon = 0.01

  # iterate through and fix the parameters -
  new_parameter_array = copy(parameter_array)
  for (index,value) in enumerate(parameter_array)

    lower_bound = lower_bound_array[index]
    upper_bound = upper_bound_array[index]

    if (value<lower_bound)
      new_parameter_array[index] = lower_bound
    elseif (value>upper_bound)
      new_parameter_array[index] = upper_bound
    end
  end

  return new_parameter_array
end
# --------------------------------------- DO NOT EDIT ABOVE THIS LINE  --------------------------------------- #


text:
abstract type AbstractEmbed{F} end
abstract type AbstractBroadcastEmbed{F} <: AbstractEmbed{F} end

Base.eltype(::AbstractEmbed{F}) where F = F
_getF(::Type{<:AbstractEmbed{F}}) where F = F

struct CompositeEmbedding{F, T <: NamedTuple, T2 <: NamedTuple, P} <: AbstractEmbed{F}
    embeddings::T
    aggregator::T2
    postprocessor::P
    CompositeEmbedding(es::NamedTuple, as::NamedTuple, post) = new{_getF(eltype(es)), typeof(es), typeof(as), typeof(post)}(es, as, post)
end

"""
    CompositeEmbedding(;postprocessor=identity, es...)

composite several embedding into one embedding according the aggregate methods and apply `postprocessor` on it.
"""
function CompositeEmbedding(;postprocessor=identity, es...)
  eas = es.data
  emb = map(x-> x isa Tuple ? x[1] : x, eas)
  agg = map(x-> x isa Tuple ? x[2] : +, eas)
  CompositeEmbedding(emb, agg, postprocessor)
end

Flux.children(ce::CompositeEmbedding) = tuple(values(ce.embeddings)..., ce.postprocessor)
Flux.mapchildren(f, ce::CompositeEmbedding) = CompositeEmbedding(map(e->f(e), ce.embeddings), ce.aggregator, f(ce.postprocessor))

function Base.show(io::IO, e::CompositeEmbedding)
    names = keys(e.embeddings)
    firstname = first(names)

    print(io, "CompositeEmbedding(")
    print(io, firstname)
    print(io, " = ")
    print(io, e.embeddings[firstname])

    for name ∈ Base.tail(names)
        print(io, ", ")
        print(io, name)
        print(io, " = ")
        print(io, e.embeddings[name])
    end

    if e.postprocessor !== identity
        print(io, ", postprocessor = ")
        print(io, e.postprocessor)
    end

    print(io, ")")
end

get_value(e::AbstractEmbed, name::Symbol, xs::NamedTuple) = e(xs[name])
_repeatdims(bs, vs) = ntuple(i -> i > length(vs) ? bs[i] : 1, length(vs))

@inline aggregate(e::E, f::F, base::A, value::A2) where {F, T, A <: AbstractArray{T}, A2 <: AbstractArray{T}, E <: AbstractEmbed{T}} = f(base, value)
@inline aggregate(e::Eb, f::F, base::A, value::A2) where {F, T, A <: AbstractArray{T}, A2 <: AbstractArray{T}, Eb <: AbstractBroadcastEmbed{T}} = f.(base, value)
@inline aggregate(e::Eb, f::typeof(vcat), base::A, value::A2) where {T, A <: AbstractArray{T}, A2 <: AbstractArray{T}, Eb <: AbstractBroadcastEmbed{T}} = vcat(base, repeat(value, _repeatdims(size(base), size(value))...))
@inline aggregate(e::Eb, f::typeof(hcat), base::A, value::A2) where {T, A <: AbstractArray{T}, A2 <: AbstractArray{T}, Eb <: AbstractBroadcastEmbed{T}} = hcat(base, repeat(value, _repeatdims(size(base), size(value))...))

(ce::CompositeEmbedding)(;xs...) = ce(xs.data)
function (ce::CompositeEmbedding{F, T1, T2, P})(xs::NamedTuple) where {F, T1, T2, P}
    names = keys(ce.embeddings)
    firstname = first(names)
    init = get_value(ce.embeddings[firstname], firstname, xs)

    for name ∈ Base.tail(names)
        emb = ce.embeddings[name]
        agg = ce.aggregator[name]
        value = get_value(emb, name, xs)
        init = aggregate(emb, agg, init, value)
    end

    if P === typeof(identity)
        init
    else
        ce.postprocessor(init)
    end
end


text:
const KALLSYMS_CACHE = Ref{Vector{String}}()
const KALLSYMS_REGEX = r"([0-9abcdef]*) ([a-zA-Z]) ([0-9a-zA-Z_\-]*)"
function find_ksym(addr::UInt64)
    start, stop = UInt64(0), addr
    if !isassigned(KALLSYMS_CACHE)
        KALLSYMS_CACHE[] = readlines("/proc/kallsyms")
    end
    L = 1
    R = length(KALLSYMS_CACHE[])
    start_addr = start
    kind = "?"
    name = ""
    while L <= R
        M = floor(Int, (L + R) / 2)
        m = match(KALLSYMS_REGEX, KALLSYMS_CACHE[][M])
        @assert m !== nothing
        start_addr, kind, name = m.captures
        start_addr = parse(UInt64, "0x"*start_addr)
        if L == R
            return start_addr, kind, name
        elseif start_addr == stop
            return start_addr, kind, name
        elseif start_addr < stop
            L = M + 1
        else
            R = M - 1
        end
    end
    return start_addr, kind, name
end
function stack_to_string(nt::NTuple{N,UInt64}; reverse=false) where N
    iob = IOBuffer()
    order = reverse ? (N:-1:1) : (1:N)
    for i in order
        addr = nt[i]
        if addr == UInt64(0)
            continue
        end
        println(iob, "$(find_ksym(addr)[3])")
    end
    String(take!(iob))
end
function stack_to_frames(nts::Vector{Tuple{Symbol,UInt64,NTuple{N,UInt64}}}) where N
    data = UInt64[]
    lidict = Dict{UInt64, Vector{Base.StackTraces.StackFrame}}()
    for (proc,count,nt) in nts
        for i in 1:count
            for addr in nt
                if addr != 0
                    name = find_ksym(addr)[3]
                    frame = Base.StackTraces.StackFrame(Symbol(name), proc, -1, nothing, true, false, addr)
                    lidict[addr] = [frame]
                    push!(data, addr)
                end
            end
            push!(data, 0)
        end
    end
    data, lidict
end


text:
abstract type CovidState end

function statevars end

function infected end

function susceptible end

function recovered end

function init_SIRT end

function sim_step end

function unity_test_period end

function population(s::CovidState)
    return s.S + infected(s) + s.R
end

SIR(s::CovidState) = (s.S, infected(s), s.R)

"""
Action input to influence epidemic simulation dynamics
# Arguments
- `testing_prop::Real` - Proportion of population to be tested on one day
    - Simplification of typical "x-days between tests per person"  action strategy due to non agent-based model
"""
struct CovidAction
    testing_prop::Float64
end

Base.zero(CovidAction) = CovidAction(0.0)

Base.string(a::CovidAction) = string(round(a.testing_prop,sigdigits=2))

struct CovidActionSpace end

Base.rand(rng::AbstractRNG, ::CovidActionSpace) = CovidAction(rand(rng))
Base.rand(A::CovidActionSpace) = rand(Random.GLOBAL_RNG, A)

function CovidActionSpace(n::Int; zero::Bool=false)
    start = zero ? 0 : 1
    return Tuple(CovidAction(i/n) for i in start:n)
end

function CovidActionSpace(v::AbstractVector{<:Real})
    return Tuple(CovidAction(v_i) for v_i in v)
end

struct InfParams
    pos_test_probs::Vector{Float64}
    symptom_probs::Vector{Float64}
    asymptomatic_prob::Float64
    symptomatic_isolation_prob::Float64
    infectiousness::Vector{Gamma{Float64}}
end

const DEFAULT_PARAMS = InfParams(
    POS_TEST_PROBS,
    SYMPTOM_PROBS,
    ASYMPTOMATIC_PROB,
    SYMPTOMATIC_ISOLATION_PROB,
    INF_DIST
)

function cdf_step(d::Distribution, N::Int)
    return Float64[cdf(d, i) - cdf(d, i-1) for i in 1:N]
end

function init_params(asymptomatic_prob=0.10)
    infection_distributions = similar(INF_DIST)
    pos_test_probs = copy(POS_TEST_PROBS)
    for (i,d) in enumerate(INF_DIST)
        k,θ = Distributions.params(d)
        k′ = k*(rand() + 0.5)
        θ′ = θ*(rand() + 0.5)
        infection_distributions[i] = Gamma(k′,θ′)

        λ = (k′*θ′) / (k*θ)
        pos_test_probs[i] *= λ
    end

    return InfParams(
        min.(pos_test_probs, 1.0),
        cdf_step(SYMPTOM_DIST, INFECTION_HORIZON),
        asymptomatic_prob,
        SYMPTOMATIC_ISOLATION_PROB,
        infection_distributions
        )
end

const AbstractCovidPOMDP = POMDP{<:CovidState, CovidAction, Int}

function simplex_sample(N::Int, m::Float64, rng::AbstractRNG=Random.GLOBAL_RNG)
    v = rand(rng, N-1)*m
    push!(v, 0, m)
    sort!(v)
    return (v - circshift(v,1))[2:end]
end

function CovidState(pomdp::POMDP{S}, params::InfParams=DEFAULT_PARAMS) where S<:CovidState
    return S(init_SIRT(pomdp)..., params, CovidAction(0.0))
end

function rand_initialstate(pomdp::POMDP{S}) where S<:CovidState
    params = init_params()
    return S(init_SIRT(pomdp)..., params, CovidAction(0.0))
end

function rand_initialstate(pomdp::POMDP{S}, Idist::Distribution) where S<:CovidState
    inf = clamp(floor(Int, rand(Idist)), 0, pomdp.N)
    params = init_params()
    return S(init_SIRT(pomdp, inf)..., params, CovidAction(0.0))
end

function POMDPs.initialstate(pomdp::AbstractCovidPOMDP)
    return ImplicitDistribution() do rng
        rand_initialstate(pomdp)
    end
end

function POMDPs.initialstate(pomdp::AbstractCovidPOMDP, Idist::Distribution)
    return ImplicitDistribution() do rng
        rand_initialstate(pomdp, Idist)
    end
end

struct SimHist{S<:CovidState}
    S::Vector{S}
    N::Int # Total Population
    T::Int # Simulation Time
    pos_test::Vector{Int}
    actions::Vector{CovidAction}
    rewards::Vector{Float64}
    beliefs::Matrix{S}
end

infected(h::SimHist) = infected.(h.S)
susceptible(h::SimHist) = susceptible.(h.S)
recovered(h::SimHist) = recovered.(h.S)


text:
# discovering binary CUDA dependencies

using Pkg, Pkg.Artifacts
import Libdl

## global state

const __toolkit_origin = Ref{Symbol}()

"""
    toolkit_origin()

Returns the origin of the CUDA toolkit in use (either :artifact, or :local).
"""
toolkit_origin() = @after_init(__toolkit_origin[])

const __toolkit_version = Ref{VersionNumber}()

"""
    toolkit_version()

Returns the version of the CUDA toolkit in use.
"""
toolkit_version() = @after_init(__toolkit_version[])

"""
    toolkit_release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
toolkit_release() = @after_init(VersionNumber(__toolkit_version[].major, __toolkit_version[].minor))

const __nvdisasm = Ref{String}()
const __libdevice = Ref{String}()
const __libcudadevrt = Ref{String}()
const __libcupti = Ref{Union{Nothing,String}}()
const __libnvtx = Ref{Union{Nothing,String}}()
const __libcublas = Ref{String}()
const __libcusparse = Ref{String}()
const __libcusolver = Ref{String}()
const __libcufft = Ref{String}()
const __libcurand = Ref{String}()
const __libcudnn = Ref{Union{Nothing,String}}(nothing)
const __libcutensor = Ref{Union{Nothing,String}}(nothing)

nvdisasm() = @after_init(__nvdisasm[])
libdevice() = @after_init(__libdevice[])
libcudadevrt() = @after_init(__libcudadevrt[])
function libcupti()
    @after_init begin
        @assert has_cupti() "This functionality is unavailable as CUPTI is missing."
        __libcupti[]
    end
end
function libnvtx()
    @after_init begin
        @assert has_nvtx() "This functionality is unavailable as NVTX is missing."
        __libnvtx[]
    end
end

export has_cupti, has_nvtx
has_cupti() = @after_init(__libcupti[]) !== nothing
has_nvtx() = @after_init(__libnvtx[]) !== nothing

libcublas() = @after_init(__libcublas[])
libcusparse() = @after_init(__libcusparse[])
libcusolver() = @after_init(__libcusolver[])
libcufft() = @after_init(__libcufft[])
libcurand() = @after_init(__libcurand[])
function libcudnn()
    @after_init begin
        @assert has_cudnn() "This functionality is unavailabe as CUDNN is missing."
        __libcudnn[]
    end
end
function libcutensor()
    @after_init begin
        @assert has_cutensor() "This functionality is unavailabe as CUTENSOR is missing."
        __libcutensor[]
    end
end

export has_cudnn, has_cutensor
has_cudnn() = @after_init(__libcudnn[]) !== nothing
has_cutensor() = @after_init(__libcutensor[]) !== nothing

## discovery

# utilities to look up stuff in the artifact (at known locations)
artifact_binary(artifact_dir, name) = joinpath(artifact_dir, "bin", Sys.iswindows() ? "$name.exe" : name)
artifact_static_library(artifact_dir, name) = joinpath(artifact_dir, "lib", Sys.iswindows() ? "$name.lib" : "lib$name.a")
artifact_file(artifact_dir, path) = joinpath(artifact_dir, path)
function artifact_library(artifact, name, version)
    dir = joinpath(artifact, Sys.iswindows() ? "bin" : "lib")
    all_names = library_versioned_names(name, version)
    for name in all_names
        path = joinpath(dir, name)
        ispath(path) && return path
    end
    error("Could not find $name ($(join(all_names, ", ", " or "))) in $dir")
end

function artifact_cuda_library(artifact, library, toolkit_release)
    version = cuda_library_version(library, toolkit_release)
    name = get(cuda_library_names, library, library)
    artifact_library(artifact, name, version)
end

# CUDA

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_artifacts = Dict(
    v"11.0" => ()->artifact"CUDA110",
    v"10.2" => ()->artifact"CUDA102",
    v"10.1" => ()->artifact"CUDA101",
    v"10.0" => ()->artifact"CUDA100",
    v"9.2"  => ()->artifact"CUDA92",
    v"9.0"  => ()->artifact"CUDA90",
)

function use_artifact_cuda()
    @debug "Trying to use artifacts..."

    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted_version = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        @debug "Selecting artifacts based on requested version $wanted_version"
        filter!(((version,artifact),) -> version == wanted_version, cuda_artifacts)
    else
        driver_version = release()
        @debug "Selecting artifacts based on driver version $driver_version"
        filter!(((version,artifact),) -> version <= driver_version, cuda_artifacts)
    end

    # download and install
    artifact = nothing
    for release in sort(collect(keys(cuda_artifacts)); rev=true)
        try
            artifact = (release=release, dir=cuda_artifacts[release]())
            break
        catch
        end
    end
    if artifact == nothing
        @debug "Could not find a compatible artifact."
        return false
    end

    __nvdisasm[] = artifact_binary(artifact.dir, "nvdisasm")
    @assert isfile(__nvdisasm[])
    __toolkit_version[] = parse_toolkit_version(__nvdisasm[])

    __libcupti[] = artifact_cuda_library(artifact.dir, "cupti",  artifact.release)
    @assert isfile(__libcupti[])
    __libnvtx[] = artifact_cuda_library(artifact.dir, "nvtx", artifact.release)
    @assert isfile(__libnvtx[])

    __libcudadevrt[] = artifact_static_library(artifact.dir, "cudadevrt")
    @assert isfile(__libcudadevrt[])
    __libdevice[] = artifact_file(artifact.dir, joinpath("share", "libdevice", "libdevice.10.bc"))
    @assert isfile(__libdevice[])

    for library in ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CUDA, Symbol("__lib$library"))

        handle[] = artifact_cuda_library(artifact.dir, library, artifact.release)
        Libdl.dlopen(handle[])
    end

    @debug "Using CUDA $(__toolkit_version[]) from an artifact at $(artifact.dir)"
    __toolkit_origin[] = :artifact
    use_artifact_cudnn(artifact.release)
    use_artifact_cutensor(artifact.release)
    return true
end

function use_local_cuda()
    @debug "Trying to use local installation..."

    cuda_dirs = find_toolkit()

    let path = find_cuda_binary("nvdisasm", cuda_dirs)
        if path === nothing
            @debug "Could not find nvdisasm"
            return false
        end
        __nvdisasm[] = path
    end

    cuda_version = parse_toolkit_version(__nvdisasm[])
    __toolkit_version[] = cuda_version

    __libcupti[] = find_cuda_library("cupti", cuda_dirs, cuda_version)
    __libnvtx[] = find_cuda_library("nvtx", cuda_dirs, cuda_version)

    let path = find_libcudadevrt(cuda_dirs)
        if path === nothing
            @debug "Could not find libcudadevrt"
            return false
        end
        __libcudadevrt[] = path
    end
    let path = find_libdevice(cuda_dirs)
        if path === nothing
            @debug "Could not find libdevice"
            return false
        end
        __libdevice[] = path
    end

    for library in ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CUDA, Symbol("__lib$library"))

        path = find_cuda_library(library, cuda_dirs, cuda_version)
        if path === nothing
            @debug "Could not find $library"
            return false
        end
        handle[] = path
    end

    @debug "Found local CUDA $(cuda_version) at $(join(cuda_dirs, ", "))"
    __toolkit_origin[] = :local
    use_local_cudnn(cuda_dirs)
    use_local_cutensor(cuda_dirs)
    return true
end

# CUDNN

# NOTE: we currently support both CUDNN v8 and v7

const cudnn_artifacts = Dict(
    v"11.0" => ()->(artifact"CUDNN_CUDA110", v"8"),
    v"10.2" => ()->(artifact"CUDNN_CUDA102", v"8"),
    v"10.1" => ()->(artifact"CUDNN_CUDA101", v"8"),
    v"10.0" => ()->(artifact"CUDNN_CUDA100", v"7"),
    v"9.2"  => ()->(artifact"CUDNN_CUDA92",  v"7"),
    v"9.0"  => ()->(artifact"CUDNN_CUDA90",  v"7"),
)

function use_artifact_cudnn(release)
    artifact_dir, version = try
        cudnn_artifacts[release]()
    catch ex
        @debug "Could not use CUDNN from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    __libcudnn[] = artifact_library(artifact_dir, "cudnn", version)
    if version >= v"8"
        # HACK: eagerly open CUDNN v8 sublibraries to avoid dlopen discoverability issues
        for sublibrary in ("ops_infer", "ops_train",
                           "cnn_infer", "cnn_train",
                           "adv_infer", "adv_train")
            sublibrary_path = artifact_library(artifact_dir, "cudnn_$(sublibrary)", version)
            Libdl.dlopen(sublibrary_path)
        end
    end
    Libdl.dlopen(__libcudnn[])
    @debug "Using CUDNN from an artifact at $(artifact_dir)"
    return true
end

function use_local_cudnn(cuda_dirs)
    path = find_library("cudnn", v"8"; locations=cuda_dirs)
    if path !== nothing
        # HACK: eagerly open CUDNN v8 sublibraries to avoid dlopen discoverability issues
        for sublibrary in ("ops_infer", "ops_train",
                           "cnn_infer", "cnn_train",
                           "adv_infer", "adv_train")
            sublibrary_path = find_library("cudnn_$(sublibrary)", v"8"; locations=cuda_dirs)
            @assert sublibrary_path !== nothing "Could not find CUDNN sublibrary $sublibrary"
            Libdl.dlopen(sublibrary_path)
        end
    end
    if path === nothing
        path = find_library("cudnn", v"7"; locations=cuda_dirs)
    end
    if path === nothing
        path = find_library("cudnn"; locations=cuda_dirs)
    end
    path === nothing && return false

    Libdl.dlopen(path)
    __libcudnn[] = path
    @debug "Using local CUDNN at $(path)"
    return true
end

# CUTENSOR

const cutensor_artifacts = Dict(
    v"11.0" => ()->artifact"CUTENSOR_CUDA110",
    v"10.2" => ()->artifact"CUTENSOR_CUDA102",
    v"10.1" => ()->artifact"CUTENSOR_CUDA101",
)

function use_artifact_cutensor(release)
    artifact_dir = try
        cutensor_artifacts[release]()
    catch ex
        @debug "Could not use CUTENSOR from artifacts" exception=(ex, catch_backtrace())
        return false
    end

    # cutensor.dll is unversioned on Windows
    version = Sys.iswindows() ? nothing : v"1"

    __libcutensor[] = artifact_library(artifact_dir, "cutensor", version)
    Libdl.dlopen(__libcutensor[])
    @debug "Using CUTENSOR from an artifact at $(artifact_dir)"
    return true
end

function use_local_cutensor(cuda_dirs)
    path = find_library("cutensor", v"1"; locations=cuda_dirs)
    if path === nothing
        path = find_library("cutensor"; locations=cuda_dirs)
    end
    path === nothing && return false

    Libdl.dlopen(path)
    __libcutensor[] = path
    @debug "Using local CUTENSOR at $(path)"
    return true
end

function __init_dependencies__()
    found = false

    # CI runs in a well-defined environment, so prefer a local CUDA installation there
    if parse(Bool, get(ENV, "CI", "false")) && !haskey(ENV, "JULIA_CUDA_USE_BINARYBUILDER")
        found = use_local_cuda()
    end

    if !found && parse(Bool, get(ENV, "JULIA_CUDA_USE_BINARYBUILDER", "true"))
        found = use_artifact_cuda()
    end

    if !found
        found = use_local_cuda()
    end

    return found
end


text:
# # Pricing under RI w/o Endogenous Feedback

# This document goes through a couple of examples for solving pricing under rational inattention without endogenous feedback using the [DRIPs](https://github.com/afrouzi/DRIPs.jl) package.

# [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/afrouzi/DRIPs.jl/binder?filepath=examples) to run and modify the following code (no software is needed on the local machine).

# See [Afrouzi and Yang (2020)](http://www.afrouzi.com/dynamic_inattention.pdf) for background on the theory.

# ## Contents
# * [Ex. 1A: One Shock Case](@ref ex1a_1shock)
#     * [Initialization](@ref ex1a_initialize)
#     * [Solution](@ref ex1a_solution)
#     * [Measure Performance](@ref ex1a_performance)
#     * [IRFs](@ref ex1a_irfs)
# * [Ex. 1B: Two Shocks Case](@ref ex1b_twoshocks)
#     * [Initialization](@ref ex1b_initialize)
#     * [Solution](@ref ex1b_solution)
#     * [Measure Performance](@ref ex1b_performance)
#     * [IRFs](@ref ex1b_irfs)

# ## [Ex. 1A: One Shock Case](@id ex1a_1shock)
# There is a measure of firms indexed by $i\in[0,1]$. Firm $i$ chooses its price $p_{i,t}$ at time $t$ to track its ideal price $p_{i,t}^*$. Formally, her flow profit is
#     $$-(p_{i,t}-p_{i,t}^*)^2$$

# We start by assuming that $p_{i,t}^*=q_t$ where
# ```math
# \begin{aligned}
#     \Delta q_t&=\rho \Delta q_{t-1}+u_t,\quad u_t\sim \mathcal{N}(0,\sigma_u^2)
# \end{aligned}
# ```
# Here $q_t$ can be interpreted as money growth or the nominal aggregate demand. Therefore, the state-space representation of the problem is

# ```math
# \begin{aligned}
#    \vec{x}_{t}&=\left[\begin{array}{c}
#     q_{t}\\
#     \Delta q_{t}
#     \end{array}\right]
#     = \underset{\mathbf{A}}{\underbrace{\left[\begin{array}{cc}
#     1 & \rho \\
#     0 & \rho \\
#     \end{array}\right]}}\, \vec{x}_{t-1}
#     + \underset{\mathbf{Q}}{\underbrace{\left[\begin{array}{c}
#     \sigma_u  \\
#     \sigma_u \\
#     \end{array}\right]}}\, u_t, \\
#     p_{i,t}^*&=\underset{\mathbf{H}}{\underbrace{\left[\begin{array}{c}
#     1 \\
#     0 \\
#     \end{array}\right]}}'\vec{x}_t
# \end{aligned}
# ```

# ### [Initialization](@id ex1a_initialize)
# Include the package:
using DRIPs;

# Assign value to deep parameters and define the structure of the problem

ρ   = 0.6;        #persistence of money growth
σ_u = 1;          #std. deviation of shocks to money growth

# Primitives of the DRIP:

ω   = 100;
β   = 0.96^0.25;
A   = [1 ρ; 0 ρ];
Q   = σ_u*[1; 1];
H   = [1; 0];

# ### [Solution](@id ex1a_solution)
ex1 = Drip(ω,β,A,Q,H);

# ### [Measure Performance](@id ex1a_performance)

using BenchmarkTools;
@benchmark Drip(ω,β,A,Q,H) setup = (ω = 100*rand()) # solves and times the function for a random set of ω's

# ### [IRFs](@id ex1a_irfs)

ex1irfs = irfs(ex1, T = 20);

# Let's plot how the average price $p=\int_0^1 p_{i,t}di$ responds to a shock to money growth:

using Plots, LaTeXStrings; pyplot();
plot(1:ex1irfs.T,[ex1irfs.x[1,1,:],ex1irfs.a[1,1,:]],
    xlabel     = "Time",
    label      = [L"Nominal Agg. Demand ($q$)" L"Price ($p$)"],
    title      = "IRFs to 1 Std. Dev. Expansionary Shock",
    xlim       = (1,ex1irfs.T),
    lw         = 3,
    legend     = :bottomright,
    legendfont = font(12),
    tickfont   = font(12),
    framestyle = :box)

# We can also plot the IRFs of inflation $\pi_t\equiv p_t-p_{t-1}$ and output $y_t\equiv q_t-p_t$ to 1 percent expansionary shock to $q$:

p1 = plot(1:ex1irfs.T,ex1irfs.x[1,1,:]-ex1irfs.a[1,1,:],
    title  = L"Output ($y_t$)")

p2 = plot(1:ex1irfs.T,[ex1irfs.a[1,1,1];ex1irfs.a[1,1,2:end]-ex1irfs.a[1,1,1:end-1]],
    title  = L"Inflation ($\pi_t$)")

plot(p1,p2,
    layout     = (1,2),
    xlim       = (1,ex1irfs.T),
    lw         = 3,
    legend     = false,
    tickfont   = font(12),
    framestyle = :box)

# ## [Ex. 1B: Two Shocks Case](@id ex1b_twoshocks)
# Suppose now that $p_{i,t}^*=q_t-z_{t}$ where
# ```math
# \begin{aligned}
#     \Delta q_t&=\rho \Delta q_{t-1}+u_t,\quad u_t\sim \mathcal{N}(0,\sigma_u^2) \\
#     z_{t}&\sim \mathcal{N}(0,\sigma_z^2)
# \end{aligned}
# ```
# Here $q_t$ can be interpreted as money growth and $z_{i,t}$ as an idiosyncratic TFP shock. Therefore,
# ```math
# \begin{aligned}
#    \vec{x}_{t}&=\left[\begin{array}{c}
#     q_{t}\\
#     \Delta q_{t} \\
#     z_{t}
#     \end{array}\right]
#     = \underset{\mathbf{A}}{\underbrace{\left[\begin{array}{ccc}
#     1 & \rho & 0\\
#     0 & \rho & 0\\
#     0 & 0 & 0 \\
#     \end{array}\right]}}\, \vec{x}_{t-1}
#     + \underset{\mathbf{Q}}{\underbrace{\left[\begin{array}{cc}
#     \sigma_u & 0 \\
#     \sigma_u & 0 \\
#     0 & \sigma_z   \\
#     \end{array}\right]}}\, \left[\begin{array}{c}
#     u_t  \\
#     z_{t}   \\
#     \end{array}\right], \\
#     p_{i,t}^*&=\underset{\mathbf{H}}{\underbrace{\left[\begin{array}{c}
#     1     \\
#     0 \\
#     -1 \\
#     \end{array}\right]}}'\vec{x}_{t}
# \end{aligned}
# ```

# ### [Initialization](@id ex1b_initialize)
# Assign values:
ρ   = 0.6;        #persistence of money growth
σ_u = 1;          #std. deviation of shocks to money growth
σ_z = √10;      #std. deviation of idiosyncratic shock

# Specifying the primitives of the drip
ω   = 100;
β   = 0.96^0.25;
A   = [1 ρ 0; 0 ρ 0; 0 0 0];
Q   = [σ_u 0; σ_u 0; 0 σ_z];
H   = [1; 0; -1];

# ### [Solution](@id ex1b_solution)

ex2  = Drip(ω,β,A,Q,H);

# ### [Measure Performance](@id ex1b_performance)

@benchmark Drip(ω,β,A,Q,H) setup = (ω = 100*rand()) # solves and times the function for a random set of ω's

# ### [IRFs](@id ex1b_irfs)

ex2irfs = irfs(ex2, T = 20);

# To get the IRFs simply use the law of motion for actions:

p1 = plot(1:ex2irfs.T,[ex2irfs.x[1,1,:],ex2irfs.a[1,1,:]],
    title  = L"IRFs to $q$ shock");
p2 = plot(1:ex1irfs.T,[ex2irfs.x[1,2,:],ex2irfs.a[1,2,:]],
    title  = L"IRFs to $z$ shock");

plot(p1,p2, layout = (1,2),
    xlabel     = "Time",
    label      = [L"Agg. Demand ($q$)" L"Price ($p$)"],
    xlim       = (1,ex2irfs.T),
    lw         = 3,
    legend     = :bottomright,
    legendfont = font(12),
    tickfont   = font(12),
    framestyle = :box)

# More IRFs:
p1 = plot(1:ex2irfs.T,ex2irfs.x[1,1,:]-ex2irfs.a[1,1,:],
    title  = L"Output ($q_0\to y_t$)");

p2 = plot(1:ex2irfs.T,[ex2irfs.a[1,1,1];ex2irfs.a[1,1,2:end]-ex2irfs.a[1,1,1:end-1]],
    title  = L"Inflation ($q_0\to \pi_t$)")

p3 = plot(1:ex2irfs.T,ex2irfs.x[1,2,:]-ex2irfs.a[1,2,:],
    title  = L"Output ($z_0\to y_t$)");

p4 = plot(1:ex2irfs.T,[ex2irfs.a[1,2,1];ex2irfs.a[1,2,2:end]-ex2irfs.a[1,2,1:end-1]],
    title  = L"Inflation ($z_0\to \pi_t$)")

plot(p1,p2,p3,p4, layout = (2,2),
    xlim       = (1,ex2irfs.T),
    lw         = 3,
    legend     = false,
    tickfont   = font(12),
    framestyle = :box)