text:
struct TimeSeriesConstraintSpecInternal
    constraint_infos::Vector{DeviceTimeSeriesConstraintInfo}
    constraint_name::Symbol
    variable_name::Symbol
    bin_variable_name::Union{Nothing, Symbol}
    param_reference::Union{Nothing, UpdateRef}
end

function lazy_lb!(psi_container::PSIContainer, inputs::TimeSeriesConstraintSpecInternal)
    time_steps = model_time_steps(psi_container)
    names = [get_component_name(x) for x in inputs.constraint_infos]
    variable = get_variable(psi_container, inputs.variable_name)
    lb_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "lb")
    con_lb = add_cons_container!(psi_container, lb_name, names, time_steps)

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        if constraint_info.range.limits.min > -Inf &&
           !isempty(constraint_info.range.additional_terms_lb)
            con_lb[ci_name, :].data .= JuMP.AffExpr(0.0)
            continue
        end
        for t in time_steps
            expression_lb = JuMP.AffExpr(0.0, variable[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_lb
                JuMP.add_to_expression!(
                    expression_lb,
                    get_variable(psi_container, val)[ci_name, t],
                    -1.0,
                )
            end
            lb_val = max(0.0, constraint_info.range.limits.min)
            con_lb[ci_name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, expression_lb >= lb_val)
        end
    end
    return
end

@doc raw"""
Constructs upper bound for given variable and time series data and a multiplier.

# Constraint

```variable[name, t] <= constraint_infos[name].multiplier * ts_data[name].timeseries[t] ```

# LaTeX

`` x_t \leq r^{val} r_t, \forall t ``
"""
function device_timeseries_ub!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    names = [get_component_name(x) for x in inputs.constraint_infos]
    variable = get_variable(psi_container, inputs.variable_name)
    ub_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "ub")
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    lazy_add_lb = false

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(psi_container, val)[ci_name, t],
                )
            end
            con_ub[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub <= constraint_info.multiplier * constraint_info.timeseries[t]
            )
        end
        if constraint_info.range.limits.min > -Inf ||
           !isempty(constraint_info.range.additional_terms_lb)
            lazy_add_lb = true
        end
    end

    @debug lazy_add_lb
    lazy_add_lb && lazy_lb!(psi_container, inputs)

    return
end

@doc raw"""
Constructs lower bound for given variable subject to time series data and a multiplier.

# Constraint

``` constraint_infos[name].multiplier * ts_data[name].timeseries[t] <= variable[name, t] ```

# LaTeX

`` r^{val} r_t \leq x_t, \forall t ``

where (name, data) in range_data.
"""
function device_timeseries_lb!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    variable = get_variable(psi_container, inputs.variable_name)
    lb_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "lb")
    names = [get_component_name(x) for x in inputs.constraint_infos]
    constraint = add_cons_container!(psi_container, lb_name, names, time_steps)

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            expression_lb = JuMP.AffExpr(0.0, variable[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_lb
                JuMP.add_to_expression!(
                    expression_lb,
                    get_variable(psi_container, val)[ci_name, t],
                    -1.0,
                )
            end
            constraint[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_lb >= constraint_info.multiplier * constraint_info.timeseries[t]
            )
        end
    end
    return
end

@doc raw"""
Constructs upper bound for given variable using a parameter. The constraint is
    built with a time series data vector and a multiplier

# Constraint

``` variable[name, t] <= constraint_infos[name].multiplier * param[name, t] ```

# LaTeX

`` x^{var}_t \leq r^{val} x^{param}_t, \forall t ``
"""
function device_timeseries_param_ub!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    names = [get_component_name(x) for x in inputs.constraint_infos]
    variable = get_variable(psi_container, inputs.variable_name)
    ub_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "ub")
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    container =
        add_param_container!(psi_container, inputs.param_reference, names, time_steps)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    lazy_add_lb = false

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(psi_container, val)[ci_name, t],
                )
            end
            param[ci_name, t] =
                PJ.add_parameter(psi_container.JuMPmodel, constraint_info.timeseries[t])
            con_ub[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub <= constraint_info.multiplier * param[ci_name, t]
            )
            multiplier[ci_name, t] = constraint_info.multiplier
        end
        if constraint_info.range.limits.min > -Inf ||
           !isempty(constraint_info.range.additional_terms_lb)
            lazy_add_lb = true
        end
    end

    @debug lazy_add_lb
    lazy_add_lb && lazy_lb!(psi_container, inputs)
    return
end

@doc raw"""
Constructs lower bound for given variable using a parameter. The constraint is
    built with a time series data vector and a multiplier

# Constraint

``` constraint_infos[name].multiplier * param[name, t] <= variable[name, t] ```

# LaTeX

`` r^{val} x^{param}_t \leq x^{var}_t, \forall t ``
"""
function device_timeseries_param_lb!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    variable = get_variable(psi_container, inputs.variable_name)
    lb_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "lb")
    names = [get_component_name(x) for x in inputs.constraint_infos]
    constraint = add_cons_container!(psi_container, lb_name, names, time_steps)
    container =
        add_param_container!(psi_container, inputs.param_reference, names, time_steps)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            expression_lb = JuMP.AffExpr(0.0, variable[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_lb
                JuMP.add_to_expression!(
                    expression_lb,
                    get_variable(psi_container, val)[ci_name, t],
                    -1.0,
                )
            end
            param[ci_name, t] =
                PJ.add_parameter(psi_container.JuMPmodel, constraint_info.timeseries[t])
            constraint_info[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_lb >= constraint_info.multiplier * param[ci_name, t]
            )
            multiplier[ci_name, t] = constraint_info.multiplier
        end
    end

    return
end

@doc raw"""
Constructs upper bound for variable and time series or confines to 0 depending on binary variable.
    The upper bound is defined by a time series and a multiplier.

# constraint_infos

``` varcts[name, t] <= varbin[name, t]* constraint_infos[name].multiplier * ts_data[name].timeseries[t] ```

where (name, data) in range_data.

# LaTeX

`` x^{cts}_t \leq r^{val} r_t x^{bin}_t, \forall t ``
"""
function device_timeseries_ub_bin!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "ub")
    varcts = get_variable(psi_container, inputs.variable_name)
    varbin = get_variable(psi_container, inputs.bin_variable_name)
    names = [get_component_name(x) for x in inputs.constraint_infos]
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            forecast = constraint_info.timeseries[t]
            multiplier = constraint_info.multiplier
            expression_ub = JuMP.AffExpr(0.0, varcts[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(psi_container, val)[ci_name, t],
                )
            end
            con_ub[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub <= varbin[ci_name, t] * multiplier * forecast
            )
        end
    end
    return
end

@doc raw"""
Constructs upper bound for variable and time series and a multiplier or confines to 0 depending on binary variable.
    Uses BigM constraint type to allow for parameter since ParameterJuMP doesn't support var*parameter

# constraint_infos

``` varcts[name, t] - constraint_infos[name].multipliers * param[name, t] <= (1 - varbin[name, t]) * M_value ```

``` varcts[name, t] <= varbin[name, t]*M_value ```

# LaTeX

`` x^{cts}_t - r^{val} x^{param}_t \leq M(1 - x^{bin}_t ), forall t ``

`` x^{cts}_t \leq M x^{bin}_t, \forall t ``
"""
function device_timeseries_ub_bigM!(
    psi_container::PSIContainer,
    inputs::TimeSeriesConstraintSpecInternal,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "ub")
    key_status = middle_rename(inputs.constraint_name, PSI_NAME_DELIMITER, "status")

    varcts = get_variable(psi_container, inputs.variable_name)
    varbin = get_variable(psi_container, inputs.bin_variable_name)
    names = [get_component_name(x) for x in inputs.constraint_infos]
    con_ub = add_cons_container!(psi_container, ub_name, names, time_steps)
    con_status = add_cons_container!(psi_container, key_status, names, time_steps)
    container =
        add_param_container!(psi_container, inputs.param_reference, names, time_steps)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)

    for constraint_info in inputs.constraint_infos
        ci_name = get_component_name(constraint_info)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, varcts[ci_name, t] => 1.0)
            for val in constraint_info.range.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(psi_container, val)[ci_name, t],
                )
            end
            param[ci_name, t] =
                PJ.add_parameter(psi_container.JuMPmodel, constraint_info.timeseries[t])
            con_ub[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub - param[ci_name, t] * constraint_info.multiplier <=
                (1 - varbin[ci_name, t]) * M_VALUE
            )
            con_status[ci_name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub <= varbin[ci_name, t] * M_VALUE
            )
            multiplier[ci_name, t] = constraint_info.multiplier
        end
    end
    return
end


text:
using puMPS

D = 8
N = 32

M = rand_puMPState(ComplexF64, 2, D, N)

H = ising_local_MPO(ComplexF64)

vumps_opt!(M, H, 1e-6, maxitr=5) #Pre-optimization using DMRG-like method
minimize_energy_local!(M, H, 1000, step=0.1)

println("Computing excitations!")

ks_tocompute = [-2,-1,0,1,2]
num_states = [5,4,7,4,5]

ens, ks, exs = excitations!(M, ising_PBC_MPO_split(ComplexF64), ks_tocompute, num_states)

H1 = Hn_in_basis(M, ising_Hn_MPO_split(ComplexF64, 1, N), exs, ks_tocompute)
H2 = Hn_in_basis(M, ising_Hn_MPO_split(ComplexF64, 2, N), exs, ks_tocompute)

ind1 = argmin(real.(ens))
indT = argmax(abs.(H2[:,ind1]))

en1 = ens[ind1]
enT = ens[indT]

fac = 2.0 / (enT-en1)

c = 2*abs2(H2[indT,ind1] * fac)
@show real(c)

using PyPlot

plot(ks, real(ens), "o")
show()

ext: jl

max_stars_repo_name: FuTen/puMPS

max_issues_repo_name: FuTen/puMPS

max_forks_repo_name: FuTen/puMPS


text:
using HybridFunnelGraphs
using HPD

domain = load_domain("test/hpd/domain.hpd")
problem = load_problem("test/hpd/problem.hpd")
1


text:
"""
Low-level version of `pycall!(ret, o, ...; kwargs...)`
Sets `ret.o` to the result of the call, and returns `ret::PyObject`
"""
function _pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, args, kwargs)
    if isempty(kwargs)
        kw = C_NULL
    else
        kw = PyObject(Dict{AbstractString, Any}([Pair(string(k), v) for (k, v) in kwargs]))
    end
    _pycall!(ret, o, args, length(args), kw)
end

"""
Low-level version of `pycall!(ret, o, ...)` for when `kw` is already in python
friendly format but you don't have the python tuple to hold the arguments
(`pyargsptr`). Sets `ret.o` to the result of the call, and returns `ret::PyObject`.
"""
function _pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, args, nargs::Int=length(args),
                  kw::Union{Ptr{Cvoid}, PyObject}=C_NULL)
    pyargsptr = @pycheckn ccall((@pysym :PyTuple_New), PyPtr, (Int,), nargs)
    try
        for i = 1:nargs
            pyarg = PyObject(args[i])
            pyincref(pyarg) # PyTuple_SetItem steals the reference
            @pycheckz ccall((@pysym :PyTuple_SetItem), Cint,
                                (PyPtr,Int,PyPtr), pyargsptr, i-1, pyarg)
        end
        return __pycall!(ret, pyargsptr, o, kw) #::PyObject
    finally
        pydecref_(pyargsptr)
    end
end

"""
Lowest level version of  `pycall!(ret, o, ...)`, assumes `pyargsptr` and `kw`
have all their args set to Python values, so we can just call the function `o`.
Sets `ret.o` to the result of the call, and returns `ret::PyObject`.
"""
function __pycall!(ret::PyObject, pyargsptr::PyPtr, o::Union{PyObject,PyPtr},
  kw::Union{Ptr{Cvoid}, PyObject})
    disable_sigint() do
        retptr = @pycheckn ccall((@pysym :PyObject_Call), PyPtr, (PyPtr,PyPtr,PyPtr), o,
                        pyargsptr, kw)
        pydecref_(ret)
        setfield!(ret, :o, retptr)
    end
    return ret #::PyObject
end

"""
```
pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, returntype::Type, args...; kwargs...)
```
Set `ret` to the result of `pycall(o, returntype, args...; kwargs)` and return
`convert(returntype, ret)`.
Avoids allocating an extra PyObject for `ret`. See `pycall` for other details.
"""
pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...) =
    return convert(returntype, _pycall!(ret, o, args, kwargs))

function pycall!(ret::PyObject, o::T, returntype::Type{PyObject}, args...; kwargs...) where
  T<:Union{PyObject,PyPtr}
    return _pycall!(ret, o, args, kwargs)
end

pycall!(ret::PyObject, o::Union{PyObject,PyPtr}, ::Type{PyAny}, args...; kwargs...) =
    return convert(PyAny, _pycall!(ret, o, args, kwargs))

"""
```
pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...)
```
Call the given Python function (typically looked up from a module) with the
given args... (of standard Julia types which are converted automatically to the
corresponding Python types if possible), converting the return value to
returntype (use a returntype of PyObject to return the unconverted Python object
reference, or of PyAny to request an automated conversion)
"""
pycall(o::Union{PyObject,PyPtr}, returntype::TypeTuple, args...; kwargs...) =
    return convert(returntype, _pycall!(PyNULL(), o, args, kwargs))::returntype

pycall(o::Union{PyObject,PyPtr}, ::Type{PyAny}, args...; kwargs...) =
    return convert(PyAny, _pycall!(PyNULL(), o, args, kwargs))

(o::PyObject)(args...; kwargs...) =
    return convert(PyAny, _pycall!(PyNULL(), o, args, kwargs))

PyAny(o::PyObject) = convert(PyAny, o)

"""
    @pycall func(args...)::T

Convenience macro which turns `func(args...)::T` into pycall(func, T, args...)
"""
macro pycall(ex)
    if !(isexpr(ex,:(::)) && isexpr(ex.args[1],:call))
        throw(ArgumentError("Usage: @pycall func(args...)::T"))
    end
    func = ex.args[1].args[1]
    args, kwargs = ex.args[1].args[2:end], []
    if isexpr(args[1],:parameters)
        kwargs, args = args[1], args[2:end]
    end
    T = ex.args[2]
    :(pycall($(map(esc,[kwargs; func; T; args])...)))
end


text:
#=

https://discourse.julialang.org/t/how-to-replace-consume-and-produce-with-channels/5125/4

Of course, you can specify larger capacity for channels to avoid to frequent switching between tasks.

To summarize:

Create a channel available for both - producer and consumer tasks.

Use put!(chnl, x) instead of produce() and take!(chnl) instead of consume(), or just iterate over the channel.
Call close(chnl) to finish.

=#

fib(n::Integer) = n < 3 ? 1 : fib(n-1) + fib(n-2)

n = 15;
@async for i=1:n
  put!(chnl, fib(i))
end

# chnl.data

OR

fib(n::Integer) = n < 3 ? 1 : fib(n-1) + fib(n-2)

n = 15;
@async begin
  for i=1:n
    put!(chnl, fib(i))
  end
  put!(chnl,-1)
end

for x in chnl
    x < 0 ? break : println(x)
end

close(chnl)


text:
macro latexrun(expr)
    return quote
        $(esc(expr))
        latexify($(string(expr)))
    end
end

macro latexify(expr)
    return quote
        latexify($(string(expr)))
    end
end


text:
const RD = RobotDynamics
using TrajectoryOptimization: QuatLQRCost, ErrorQuadratic, QuatVecEq

function YakProblems(;
        integration=RD.RK4,
        N = 101,
        vecstate=false,
        scenario=:barrellroll,
        costfun=:Quadratic,
        termcon=:goal,
        quatnorm=:none,
        kwargs...
    )
    model = RobotZoo.YakPlane(UnitQuaternion)

    opts = SolverOptions(
        cost_tolerance_intermediate = 1e-1,
        penalty_scaling = 10.,
        penalty_initial = 10.;
        kwargs...
    )

    s = RD.LieState(model)
    n,m = RD.dims(model)
    rsize = RD.state_dim(model) - 9
    vinds = SA[1,2,3,8,9,10,11,12,13]

    # Discretization
    tf = 1.25
    dt = tf/(N-1)

    if scenario == :barrellroll
        ey = @SVector [0,1,0.]

        # Initial and final condition
        p0 = MRP(0.997156, 0., 0.075366) # initial orientation
        pf = MRP(0., -0.0366076, 0.) # final orientation (upside down)

        x0 = RD.build_state(model, [-3,0,1.5], p0, [5,0,0], [0,0,0])
        utrim  = @SVector  [41.6666, 106, 74.6519, 106]
        xf = RD.build_state(model, [3,0,1.5], pf, [5,0,0], [0,0,0])

        # Xref trajectory
        x̄0 = RBState(model, x0)
        x̄f = RBState(model, xf)
        Xref = map(1:N) do k
            x̄0 + (x̄f - x̄0)*((k-1)/(N-1))
        end

        # Objective
        Qf_diag = RD.fill_state(model, 100, 500, 100, 100.)
        Q_diag = RD.fill_state(model, 0.1, 0.1, 0.1, 0.1)
        Qf = Diagonal(Qf_diag)
        Q = Diagonal(Q_diag)
        R = Diagonal(@SVector fill(1e-3,4))
        if quatnorm == :slack
            m += 1
            R = Diagonal(push(R.diag, 1e-6))
            utrim = push(utrim, 0)
        end
        if costfun == :Quadratic
            costfuns = map(Xref) do xref
                LQRCost(Q*dt, R*dt, xf, utrim)
            end
            costfun = LQRCost(Q*dt, R*dt, xf, utrim)
            costterm = LQRCost(Qf, R, xf, utrim)
            costfuns[end] = costterm
        elseif costfun == :QuatLQR
            costfuns = map(Xref) do xref
                QuatLQRCost(Q*dt, R*dt, xf, utrim, w=0.1)
            end
            costfun = QuatLQRCost(Q*dt, R*dt, xf, utrim; w=0.1)
            costterm = QuatLQRCost(Qf, R, xf, utrim; w=200.0)
            costfuns[end] = costterm
        elseif costfun == :LieLQR
            costfun = LieLQR(s, Q*dt, R*dt, xf, utrim)
            costterm = LieLQR(s, Qf, R*dt, xf, utrim)
        elseif costfun == :ErrorQuadratic
            costfuns = map(Xref) do xref
                ErrorQuadratic(model, Q*dt, R*dt, xref, utrim)
            end
            costfun = ErrorQuadratic(model, Q*dt, R*dt, xf, utrim)
            costterm = ErrorQuadratic(model, Qf*10, R, xf, utrim)
            costfuns[end] = costterm
        end
        obj = Objective(costfuns)

        # Constraints
        conSet = ConstraintList(n,m,N)
        vecgoal = GoalConstraint(xf, vinds)
        rot_diffmethod = RD.UserDefined()
        if termcon == :goal
            rotgoal = GoalConstraint(xf, SA[4,5,6,7])
        elseif termcon == :quatvec
            rotgoal = QuatVecEq(n, m, pf)
            rot_diffmethod = RD.ForwardAD()
        elseif termcon == :quaterr
            rotgoal = QuatErr(n, UnitQuaternion(pf), SA[4,5,6,7])
        else
            throw(ArgumentError("$termcon is not a known option for termcon. Options are :goal, :quatvec, :quaterr"))
        end
        add_constraint!(conSet, vecgoal, N)
        add_constraint!(conSet, rotgoal, N, diffmethod=rot_diffmethod)

    else
        throw(ArgumentError("$scenario isn't a known scenario"))
    end

    # Initialization
    U0 = [copy(utrim) for k = 1:N-1]

    # Use a standard model (no special handling of rotation states)
    if quatnorm == :renorm
        model = QuatRenorm(model)
    elseif quatnorm == :slack
        model = QuatSlackModel(model)
        slackcon = UnitQuatConstraint(model)
        add_constraint!(conSet, slackcon, 1:N-1)
    end
    if vecstate
        model = VecModel(model)
    end

    # Build problem
    prob = Problem(model, obj, x0, tf, xf=xf, constraints=conSet, integration=integration(model))
    initial_controls!(prob, U0)
    prob, opts
end


text:
using MusicXML

@MX begin
## Part List
### Piano
midiinstrument_piano = MidiInstrument(channel= 1, program =1, volume = 100, pan =0, id = "P1-I1")
scorepart_piano = ScorePart(name = "Piano", midiinstrument = midiinstrument_piano, id = "P1")

### Acoustic Guitar
midiinstrument_guitar = MidiInstrument(channel= 1, program =26, volume = 100, pan =0, id = "P2-I1")
scorepart_guitar = ScorePart(name = "Acoustic Guitar", midiinstrument = midiinstrument_guitar, id = "P2")

###
partlist = PartList(scoreparts = [scorepart_piano, scorepart_guitar])

## Part
### Piano

attributes1_piano = Attributes(
   time = Time(beats = 4, beattype = 4), # 4/4
   divisions = 4, # we want to use 16th notes at minimum
   clef = [Clef(number = 1, sign = "G", line = 2), Clef(number = 2, sign = "F", line = 4)], # Piano clefs
   staves = 2, # Piano staves
   key = Key(fifths = 0, mode = "major"), # no accidentals, major key
)

measure1_notes_piano = [
Note(pitch = Pitch(step = "C", alter = 0, octave = 4), duration =  4),
Note(pitch = Pitch(step = "D", alter = 0, octave = 4), duration =  4),
Note(pitch = Pitch(step = "E", alter = 0, octave = 4), duration =  4),
Note(pitch = Pitch(step = "F", alter = +1, octave = 4), duration =  4),
]

measure2_notes_piano = [
Note(pitch = Pitch(step = "G", alter = 0, octave = 5), duration =  1),
Note(pitch = Pitch(step = "G", alter = +1, octave = 5), duration =  1),
Note(pitch = Pitch(step = "B", alter = 0, octave = 5), duration =  1),
Note(pitch = Pitch(step = "A", alter = +1, octave = 5), duration =  1),
Note(rest = Rest(), duration =  4), # Rest
Note(pitch = Pitch(step = "A", alter = 0, octave = 5), duration =  4),
Note(pitch = Pitch(step = "B", alter = 0, octave = 5), duration =  4),
]

measures_piano = [
Measure(attributes = attributes1_piano, notes = measure1_notes_piano) # measure 1 has attributes
# Measure(notes = measure2_notes_piano)
]

part_piano = Part(measures = measures_piano, id = "P1")

### Guitar

attributes1_guitar = Attributes(
   time = Time(beats = 4, beattype = 4), # 4/4
   divisions = 4, # we want to use 16th notes at minimum
   clef = [Clef(number = 1, sign = "G", line = 2), Clef(number = 2, sign = "TAB", line = 6)], # Guitar Clefs
   staves = 2, # Guitar staves
   key = Key(fifths = 0, mode = "major"), # no accidentals, major key
)

measure1_notes_guitar = [
# G Major chord for a bar
Note(pitch = Pitch(step = "G", alter = 0, octave = 2), duration =  16),
Note(pitch = Pitch(step = "B", alter = 0, octave = 2), duration =  16, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 3), duration =  16, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 3), duration =  16, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 4), duration =  16, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 4), duration =  16, chord = Chord()),
]

measure2_notes_guitar = [
# G Major chord for half a bar
Note(pitch = Pitch(step = "G", alter = 0, octave = 2), duration =  8),
Note(pitch = Pitch(step = "B", alter = 0, octave = 2), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 3), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 3), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 4), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 4), duration =  8, chord = Chord()),

# G Major chord for half a bar
Note(pitch = Pitch(step = "G", alter = 0, octave = 2), duration =  8),
Note(pitch = Pitch(step = "B", alter = 0, octave = 2), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 3), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 3), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "D", alter = 0, octave = 4), duration =  8, chord = Chord()),
Note(pitch = Pitch(step = "G", alter = 0, octave = 4), duration =  8, chord = Chord()),
]

measures_guitar = [
Measure(attributes = attributes1_guitar, notes = measure1_notes_guitar) # measure 1 has attributes
Measure(notes = measure2_notes_guitar)
]

part_guitar = Part(measures = measures_guitar, id = "P2")

##
score = ScorePartwise(
partlist = partlist,
parts =  [part_piano, part_guitar],
)

writemusicxml("myscore.musicxml", score)

end


text:
include("poincare_invariant_2nd_module.jl")

using BenchmarkTools
using GeometricIntegrators
using Requires
using Test

using .PoincareInvariant2ndTest

set_config(:verbosity, 1)

const nx = 141
const ny = 142

println("\nSecond Canonical Euler-Poincaré Integral Invariant (ApproxFun):")
@btime compute_canonical_invariant_approxfun(nx, ny)

println("\nSecond Euler-Poincaré Integral Invariant (ApproxFun):")
@btime compute_invariant_approxfun(nx, ny)

# @require ClassicalOrthogonalPolynomials = "b30e2e7b-c4ee-47da-9d5f-2c5c27239acd" begin
#     @require FastTransforms = "057dd010-8810-581a-b7be-e3fc3b93f78c" begin
        println("\nSecond Euler-Poincaré Integral Invariant (OrthogonalPolynomials):")
        @btime compute_invariant_cop(nx, ny)
#     end
# end

println("\nSecond Euler-Poincaré Integral Invariant (Trapezoidal):")
@btime compute_invariant_trapezoidal(nx, ny)


text:
using Oceananigans: cell_advection_timescale
using Oceananigans.TurbulenceClosures: cell_diffusion_timescale

"""
    CFL{D, S}

An object for computing the Courant-Freidrichs-Lewy (CFL) number.
"""
struct CFL{D, S}
           Δt :: D
    timescale :: S
end

"""
    CFL(Δt [, timescale=Oceananigans.cell_advection_timescale])

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and `timescale`.

See also `AdvectiveCFL` and `DiffusiveCFL`.
"""
CFL(Δt) = CFL(Δt, cell_advection_timescale)

(c::CFL{<:Number})(model) = c.Δt / c.timescale(model)
(c::CFL{<:TimeStepWizard})(model) = c.Δt.Δt / c.timescale(model)

"""
    AdvectiveCFL(Δt)

Returns an object for computing the Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale
for advection across a cell.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(8, 8, 8)));

julia> cfl = AdvectiveCFL(1.0);

julia> data(model.velocities.u) .= π;

julia> cfl(model)
6.283185307179586
```
"""
AdvectiveCFL(Δt) = CFL(Δt, cell_advection_timescale)

"""
    DiffusiveCFL(Δt)

Returns an object for computing the diffusive Courant-Freidrichs-Lewy (CFL) number
associated with time step or `TimeStepWizard` `Δt` and the time scale for diffusion
across a cell associated with `model.closure`.

The maximum diffusive CFL number among viscosity and all tracer diffusivities is
returned.

Example
=======
```julia
julia> model = Model(grid=RegularCartesianGrid(size=(16, 16, 16), length=(1, 1, 1)));

julia> dcfl = DiffusiveCFL(0.1);

julia> dcfl(model)
2.688e-5
```
"""
DiffusiveCFL(Δt) = CFL(Δt, cell_diffusion_timescale)


text:
# does MCMC using the NN estimate as the statistic
using LinearAlgebra, Statistics #, Optim

# the main MCMC routine, does several short chains to tune proposal
# then a longer final chain
# covreps: replications used to compute weight matrix (the R in the paper, eqn. 5)
function MCMC(θnn, length, model::SNMmodel, nnmodel, nninfo; covreps = 1000, tuningloops = 10, verbosity = false, do_cue = false, burnin = 0, nthreads=1, tuning = 1.0) #, rt=0.25)
    nParams = size(model.lb,1)
    # make sure these are defined at this scope
    Σ = EstimateΣ(θnn, covreps, model, nnmodel, nninfo)
    # set up the proposal
    P = ((cholesky(Σ)).U)' # transpose it here
    reps = 1.0
    Σinv = 1.0
    lnL = θ -> 1.0
    # define things for MCMC
    if !do_cue
        reps = 10 # fewer for the MC chain
        Σinv = inv((1.0+1/reps).*Σ)
        lnL = θ -> H(θ, θnn, reps, model, nnmodel, nninfo, Σinv)
    else
        reps = 100
        lnL = θ -> H(θ, θnn, reps, model, nnmodel, nninfo, true)
    end
    ChainLength = 200
    # loops to tune proposal
    chain = 0.0
    @inbounds for j = 1:tuningloops+1
        Proposal = θ -> θ + tuning*P*randn(size(θ))
        if j == tuningloops + 1
            ChainLength = length
        end
        if j > 1
            θinit = mean(chain[:,1:nParams], dims=1)[:]
        else
            θinit = θnn
        end
        chain = mcmc(θinit, ChainLength, burnin, model.prior, lnL, Proposal, verbosity, nthreads)
        # adjust tuning to try to keep acceptance rate between 0.25 - 0.35
        if j <= tuningloops
            accept = mean(chain[:,end])
            if accept > 0.3
                tuning *= 1.1
            elseif accept < 0.2
                tuning /= 1.1
            end
        end
    end
    return chain[:,1:nParams], P, tuning
end


text:
init_project() = init_project(pwd())
function init_project(path)
    config = joinpath(path, "config.json")
    if !isfile(config)
        throw(ArgumentError("cannot find $(normpath(config))"))
    end

    GAMEENV["PROJECT"] = path
    GAMEENV["CONFIG"] = config

    CACHE["config"] = ConfigData(joinpath(GAMEENV["PROJECT"], "config.json"))

    # Game specific JuliaModule implementation
    if haskey(GAMEENV, "JULIAMODULE")
        if !in(GAMEENV["JULIAMODULE"], LOAD_PATH)
            push!(LOAD_PATH, GAMEENV["JULIAMODULE"])
        end
    end

    # @info "\"$(CACHE["config"]["name"])\" Project has loaded successfully!"
end


text:
function diffeq(du,u,p,t)

    # extracellular volume to cytoplasmic volume ratio
    VeVc = 33.3

    # fraction definitions
    if u[V.sigmaS] + u[V.sigmaSP] + u[V.sigmaSP_G] > 0.0
        fsigmaS = u[V.sigmaS]/(u[V.sigmaS] + u[V.sigmaSP] + u[V.sigmaSP_G])
    else
        fsigmaS = 0.0
    end

    if u[V.sigmaG] + u[V.sigmaG_A] + u[V.sigmaG_O] + u[V.A_sigmaG_O] > 0.0
        fsigmaG = u[V.sigmaG]/(u[V.sigmaG] + u[V.sigmaG_A] + u[V.sigmaG_O] + u[V.A_sigmaG_O])
    else
        fsigmaG = 0.0
    end

    if u[V.sigmaA] + u[V.sigmaAP] + u[V.sigmaAP_S] + u[V.sigmaAP_R] + u[V.sigmaAP_I] + u[V.sigmaAP_T] > 0.0
        fsigmaA = u[V.sigmaA]/(u[V.sigmaA] + u[V.sigmaAP] + u[V.sigmaAP_S] + u[V.sigmaAP_R] + u[V.sigmaAP_I] + u[V.sigmaAP_T])
    else
        fsigmaA = 0.0
    end

    if u[V.sigmaR] + u[V.sigmaRP] > 0.0
        fsigmaR = u[V.sigmaR]/(u[V.sigmaR] + u[V.sigmaRP])
    else
        fsigmaR = 0.0
    end

    sigmaEP = u[V.E11P] + u[V.E12P] + u[V.E23P] + u[V.E24P] + u[V.E34P] + u[V.E44P] + u[V.E13P] + u[V.E14P]
    if sigmaEP > 0.0
        f11 = u[V.E11P]/sigmaEP
    else
        f11 = 0.0
    end

    # E11 degradation
    v10d = p[C.kdeg]*u[V.E11P]
    v16d = p[C.kdeg]*u[V.E11G]
    v17d = p[C.kdeg]*u[V.E11S]
    v18d = p[C.kdeg]*u[V.E11R]
    v73d = p[C.kdeg]*u[V.E11T]

    v = Dict{Int64,Float64}()

    v[1] = (p[C.kon1]*u[V.E]*u[V.E1] - p[C.EGF_off]*u[V.E_E1])
    v[2] = (p[C.kon2]*u[V.H]*u[V.E3] - p[C.HRGoff_3]*u[V.H_E3])
    v[3] = (p[C.kon3]*u[V.H]*u[V.E4] - p[C.HRGoff_4]*u[V.H_E4])
    v[4] = (p[C.kon4]*u[V.E_E1]*u[V.E_E1] - p[C.koff4]*u[V.E11])
    v[5] = (p[C.kon5]*u[V.E_E1]*u[V.E2] - p[C.koff5]*u[V.E12])
    v[6] = (p[C.kon6]*u[V.H_E3]*u[V.E2] - p[C.koff6]*u[V.E23])
    v[7] = (p[C.kon7]*u[V.H_E3]*u[V.H_E4] - p[C.koff7]*u[V.E34])
    v[8] = (p[C.kon8]*u[V.H_E4]*u[V.E2] - p[C.koff8]*u[V.E24])
    v[9] = (p[C.kon9]*u[V.H_E4]*u[V.H_E4] - p[C.koff9]*u[V.E44])
    v[10] = (p[C.kf10]*u[V.E11] - p[C.VmaxPY]*u[V.E11P]/(p[C.KmPY] + u[V.E11P]) - p[C.kPTP10]*u[V.sigmaT]*u[V.E11P])
    v[11] = (p[C.kf11]*u[V.E12] - p[C.VmaxPY]*u[V.E12P]/(p[C.KmPY] + u[V.E12P]) - p[C.kPTP11]*u[V.sigmaT]*u[V.E12P])
    v[12] = (p[C.kf12]*u[V.E23] - p[C.VmaxPY]*u[V.E23P]/(p[C.KmPY] + u[V.E23P]) - p[C.kPTP12]*u[V.sigmaT]*u[V.E23P])
    v[13] = (p[C.kf13]*u[V.E34] - p[C.VmaxPY]*u[V.E34P]/(p[C.KmPY] + u[V.E34P]) - p[C.kPTP13]*u[V.sigmaT]*u[V.E34P])
    v[14] = (p[C.kf14]*u[V.E24] - p[C.VmaxPY]*u[V.E24P]/(p[C.KmPY] + u[V.E24P]) - p[C.kPTP14]*u[V.sigmaT]*u[V.E24P])
    v[15] = (p[C.kf15]*u[V.E44] - p[C.VmaxPY]*u[V.E44P]/(p[C.KmPY] + u[V.E44P]) - p[C.kPTP15]*u[V.sigmaT]*u[V.E44P])
    v[16] = (4*p[C.kon16]*u[V.E11P]*u[V.G] - p[C.koff16]*fsigmaG*u[V.E11G])
    v[17] = (8*p[C.kon17]*u[V.E11P]*u[V.S] - p[C.koff17]*fsigmaS*u[V.E11S])
    v[18] = (2*p[C.kon18]*u[V.E11P]*u[V.R] - p[C.koff18]*fsigmaR*u[V.E11R])
    v[19] = (3*p[C.kon19]*u[V.E12P]*u[V.G] - p[C.koff19]*fsigmaG*u[V.E12G])
    v[20] = (6*p[C.kon20]*u[V.E12P]*u[V.S] - p[C.koff20]*fsigmaS*u[V.E12S])
    v[21] = (2*p[C.kon21]*u[V.E12P]*u[V.R] - p[C.koff21]*fsigmaR*u[V.E12R])
    v[22] = (3*p[C.kon22]*u[V.E23P]*u[V.G] - p[C.koff22]*fsigmaG*u[V.E23G])
    v[23] = (3*p[C.kon23]*u[V.E23P]*u[V.S] - p[C.koff23]*fsigmaS*u[V.E23S])
    v[24] = (3*p[C.kon24]*u[V.E23P]*u[V.I] - p[C.koff24]*u[V.E23I])
    v[25] = (2*p[C.kon25]*u[V.E23P]*u[V.R] - p[C.koff25]*fsigmaR*u[V.E23R])
    v[26] = (4*p[C.kon26]*u[V.E34P]*u[V.G] - p[C.koff26]*fsigmaG*u[V.E34G])
    v[27] = (3*p[C.kon27]*u[V.E34P]*u[V.S] - p[C.koff27]*fsigmaS*u[V.E34S])
    v[28] = (4*p[C.kon28]*u[V.E34P]*u[V.I] - p[C.koff28]*u[V.E34I])
    v[29] = (2*p[C.kon29]*u[V.E34P]*u[V.R] - p[C.koff29]*fsigmaR*u[V.E34R])
    v[30] = (3*p[C.kon30]*u[V.E24P]*u[V.G] - p[C.koff30]*fsigmaG*u[V.E24G])
    v[31] = (4*p[C.kon31]*u[V.E24P]*u[V.S] - p[C.koff31]*fsigmaS*u[V.E24S])
    v[32] = (1*p[C.kon32]*u[V.E24P]*u[V.I] - p[C.koff32]*u[V.E24I])
    v[33] = (2*p[C.kon33]*u[V.E24P]*u[V.R] - p[C.koff33]*fsigmaR*u[V.E24R])
    v[34] = (4*p[C.kon34]*u[V.E44P]*u[V.G] - p[C.koff34]*fsigmaG*u[V.E44G])
    v[35] = (4*p[C.kon35]*u[V.E44P]*u[V.S] - p[C.koff35]*fsigmaS*u[V.E44S])
    v[36] = (2*p[C.kon36]*u[V.E44P]*u[V.I] - p[C.koff36]*u[V.E44I])
    v[37] = (2*p[C.kon37]*u[V.E44P]*u[V.R] - p[C.koff37]*fsigmaR*u[V.E44R])
    v[38] = (p[C.kf38]*u[V.sigmaS]*sigmaEP - p[C.VmaxPY]*u[V.sigmaSP]/(p[C.KmPY] + u[V.sigmaSP]) - p[C.kPTP38]*u[V.sigmaT]*u[V.sigmaSP])
    v[39] = (p[C.kf39]*u[V.sigmaA]*sigmaEP - p[C.VmaxPY]*u[V.sigmaAP]/(p[C.KmPY] + u[V.sigmaAP]) - p[C.kPTP39]*u[V.sigmaT]*u[V.sigmaAP])
    v[40] = (p[C.kon40]*u[V.sigmaG]*u[V.O] - p[C.koff40]*u[V.sigmaG_O])
    v[41] = (p[C.kon41]*u[V.sigmaG]*u[V.A] - p[C.koff41]*u[V.sigmaG_A]*fsigmaA)
    v[42] = (p[C.kon42]*u[V.sigmaSP]*u[V.G] - p[C.koff42]*u[V.sigmaSP_G]*fsigmaG)
    v[43] = (3*p[C.kon43]*u[V.sigmaAP]*u[V.S] - p[C.koff43]*u[V.sigmaAP_S]*fsigmaS)
    v[44] = (3*p[C.kon44]*u[V.sigmaAP]*u[V.I] - p[C.koff44]*u[V.sigmaAP_I])
    v[45] = (2*p[C.kon45]*u[V.sigmaAP]*u[V.R] - p[C.koff45]*u[V.sigmaAP_R]*fsigmaR)
    v[46] = (p[C.kon46]*u[V.P3]*u[V.A] - p[C.koff46]*u[V.P3_A]*fsigmaA)
    v[47] = (p[C.kf47]*u[V.P3]*u[V.Akt]/(p[C.Kmf47] + u[V.Akt]) - p[C.Vmaxr47]*u[V.Aktstar]/(p[C.Kmr47] + u[V.Aktstar]))
    v[48] = (p[C.kf48]*(1 - u[V.fint]*f11)*u[V.sigmaI]*u[V.P2]/(p[C.Kmf48] + u[V.P2]) - 3*p[C.PTEN]*u[V.P3]/(p[C.Kmr48] + u[V.P3]))
    v[49] = (p[C.kf49]*u[V.sigmaO]*u[V.RsD]/(p[C.Kmf49] + u[V.RsD]) - p[C.kr49]*u[V.sigmaR]*u[V.RsT]/(p[C.Kmr49] + u[V.RsT]) - p[C.kr49b]*u[V.sigmaRP]*u[V.RsT]/(p[C.Kmr49b] + u[V.RsT]) - p[C.kcon49]*u[V.RsT])
    v[50] = (p[C.kf50]*u[V.sigmaR]*sigmaEP - p[C.VmaxPY]*u[V.sigmaRP]/(p[C.KmPY] + u[V.sigmaRP]) - p[C.kPTP50]*u[V.sigmaT]*u[V.sigmaRP])
    v[51] = (p[C.kf51]*u[V.RsT]*u[V.Raf]/(p[C.Kmf51] + u[V.Raf]) - p[C.Vmaxr51]*u[V.Rafstar]/(p[C.Kmrb51] + u[V.Rafstar]))
    v[52] = (p[C.kf52]*u[V.Rafstar]*u[V.MEK]/(p[C.Kmf52] + u[V.MEK]) - p[C.Vmaxr52]*u[V.MEKstar]/(p[C.Kmr52] + u[V.MEKstar]))
    # v[53]: No reaction
    v[54] = (p[C.kf54]*u[V.O]*u[V.ERKstar]/(p[C.Kmf54] + u[V.O]) - p[C.Vmaxr54]*u[V.OP]/(p[C.Kmr54] + u[V.OP]))
    v[55] = (p[C.kf55]*u[V.A]*u[V.ERKstar]/(p[C.Kmf55] + u[V.A]) - p[C.Vmaxr55]*u[V.AP]/(p[C.Kmr55] + u[V.AP]))
    # v[56]: No reaction
    v[57] = (p[C.kon57]*u[V.P3_A]*u[V.G] - p[C.koff57]*u[V.sigmaA_G])
    v[58] = (p[C.kon58]*u[V.sigmaA_G]*u[V.O] - p[C.koff58]*u[V.sigmaA_G_O])
    v[59] = (p[C.kon59]*u[V.sigmaG_O]*u[V.A] - p[C.koff59]*u[V.A_sigmaG_O]*fsigmaA)
    v[60] = (p[C.kon60]*u[V.sigmaG_A]*u[V.O] - p[C.koff60]*u[V.A_sigmaG_O])
    v[61] = (p[C.kon61]*u[V.H_E3]*u[V.E_E1] - p[C.koff61]*u[V.E13])
    v[62] = (p[C.kon62]*u[V.H_E4]*u[V.E_E1] - p[C.koff62]*u[V.E14])
    v[63] = (p[C.kf63]*u[V.E13] - p[C.VmaxPY]*u[V.E13P]/(p[C.KmPY] + u[V.E13P]) - p[C.kPTP63]*u[V.sigmaT]*u[V.E13P])
    v[64] = (p[C.kf64]*u[V.E14] - p[C.VmaxPY]*u[V.E14P]/(p[C.KmPY] + u[V.E14P]) - p[C.kPTP64]*u[V.sigmaT]*u[V.E14P])
    v[65] = (4*p[C.kon65]*u[V.E13P]*u[V.G] - p[C.koff65]*fsigmaG*u[V.E13G])
    v[66] = (5*p[C.kon66]*u[V.E13P]*u[V.S] - p[C.koff66]*fsigmaS*u[V.E13S])
    v[67] = (3*p[C.kon67]*u[V.E13P]*u[V.I] - p[C.koff67]*u[V.E13I])
    v[68] = (2*p[C.kon68]*u[V.E13P]*u[V.R] - p[C.koff68]*fsigmaR*u[V.E13R])
    v[69] = (4*p[C.kon69]*u[V.E14P]*u[V.G] - p[C.koff69]*fsigmaG*u[V.E14G])
    v[70] = (6*p[C.kon70]*u[V.E14P]*u[V.S] - p[C.koff70]*fsigmaS*u[V.E14S])
    v[71] = (1*p[C.kon71]*u[V.E14P]*u[V.I] - p[C.koff71]*u[V.E14I])
    v[72] = (2*p[C.kon72]*u[V.E14P]*u[V.R] - p[C.koff72]*fsigmaR*u[V.E14R])
    v[73] = (4*p[C.kon73]*u[V.E11P]*u[V.T] - p[C.koff73]*u[V.E11T])
    v[74] = (3*p[C.kon74]*u[V.E12P]*u[V.T] - p[C.koff74]*u[V.E12T])
    v[75] = (2*p[C.kon75]*u[V.E23P]*u[V.T] - p[C.koff75]*u[V.E23T])
    v[76] = (2*p[C.kon76]*u[V.E34P]*u[V.T] - p[C.koff76]*u[V.E34T])
    v[77] = (2*p[C.kon77]*u[V.E24P]*u[V.T] - p[C.koff77]*u[V.E24T])
    v[78] = (2*p[C.kon78]*u[V.E44P]*u[V.T] - p[C.koff78]*u[V.E44T])
    v[79] = (3*p[C.kon79]*u[V.E13P]*u[V.T] - p[C.koff79]*u[V.E13T])
    v[80] = (3*p[C.kon80]*u[V.E14P]*u[V.T] - p[C.koff80]*u[V.E14T])
    v[81] = (p[C.kf81]*u[V.E1]*u[V.ERKstar]/(p[C.Kmf81] + u[V.E1]) - p[C.Vmaxr81]*u[V.E1_PT]/(p[C.Kmr81] + u[V.E1_PT]))
    v[82] = (p[C.kf82]*u[V.E2]*u[V.ERKstar]/(p[C.Kmf82] + u[V.E2]) - p[C.Vmaxr82]*u[V.E2_PT]/(p[C.Kmr82] + u[V.E2_PT]))
    v[83] = (p[C.kf83]*u[V.E4]*u[V.ERKstar]/(p[C.Kmf83] + u[V.E4]) - p[C.Vmaxr83]*u[V.E4_PT]/(p[C.Kmr83] + u[V.E4_PT]))
    v[84] = (p[C.kf84]*u[V.E_E1]*u[V.ERKstar]/(p[C.Kmf84] + u[V.E_E1]) - p[C.Vmaxr84]*u[V.E_E1_PT]/(p[C.Kmr84] + u[V.E_E1_PT]))
    v[85] = (p[C.kf85]*u[V.H_E4]*u[V.ERKstar]/(p[C.Kmf85] + u[V.H_E4]) - p[C.Vmaxr85]*u[V.H_E4_PT]/(p[C.Kmr85] + u[V.H_E4_PT]))
    v[86] = (p[C.kon86]*u[V.E]*u[V.E1_PT] - p[C.EGF_off]*u[V.E_E1_PT])
    v[87] = (p[C.kon87]*u[V.H]*u[V.E4_PT] - p[C.HRGoff_4]*u[V.H_E4_PT])
    v[88] = (2*p[C.kon88]*u[V.sigmaAP]*u[V.T] - p[C.koff88]*u[V.sigmaAP_T])
    v[89] = (p[C.kon89]*u[V.ERK]*u[V.MEKstar] - p[C.koff89]*u[V.ERK_MEKstar])
    v[90] = p[C.kcat90]*u[V.ERK_MEKstar]
    v[91] = (p[C.kon91]*u[V.pERK]*u[V.MEKstar] - p[C.koff91]*u[V.pERK_MEKstar])
    v[92] = p[C.kcat92]*u[V.pERK_MEKstar]
    v[93] = (p[C.kon93]*u[V.ERKstar]*u[V.ERKPpase] - p[C.koff93]*u[V.ERKstar_ERKPpase])
    v[94] = p[C.kcat94]*u[V.ERKstar_ERKPpase]
    v[95] = (p[C.kon95]*u[V.pERK]*u[V.ERKPpase] - p[C.koff95]*u[V.pERK_ERKPpase])
    v[96] = p[C.kcat96]*u[V.pERK_ERKPpase]

    du[V.E1] = -v[1] - v[81]
    du[V.E2] = -v[5] - v[6] - v[8] - v[82]
    du[V.E3] = -v[2]
    du[V.E4] = -v[3] - v[83]
    du[V.E_E1] = v[1] - v[4] - v[4] - v[5] - v[61] - v[62] - v[84]
    du[V.H_E3] = v[2] - v[6] - v[7] - v[61]
    du[V.H_E4] = v[3] - v[7] - v[8] - v[9] - v[9] - v[62] - v[85]
    du[V.E11] = v[4] - v[10]
    du[V.E12] = v[5] - v[11]
    du[V.E23] = v[6] - v[12]
    du[V.E34] = v[7] - v[13]
    du[V.E24] = v[8] - v[14]
    du[V.E44] = v[9] - v[15]
    du[V.E11P] = v[10] - v[16] - v[17] - v[18] - v[73] - v10d
    du[V.E12P] = v[11] - v[19] - v[20] - v[21] - v[74]
    du[V.E23P] = v[12] - v[22] - v[23] - v[24] - v[25] - v[75]
    du[V.E34P] = v[13] - v[26] - v[27] - v[28] - v[29] - v[76]
    du[V.E24P] = v[14] - v[30] - v[31] - v[32] - v[33] - v[77]
    du[V.E44P] = v[15] - v[34] - v[35] - v[36] - v[37] - v[78]
    du[V.G] = -v[16] - v[19] - v[22] - v[26] - v[30] - v[34] - v[42] - v[57] - v[65] - v[69] + v16d
    du[V.S] = -v[17] - v[20] - v[23] - v[27] - v[31] - v[35] - v[43] - v[66] - v[70] + v17d
    du[V.I] = -v[24] - v[28] - v[32] - v[36] - v[44] - v[67] - v[71]
    du[V.R] = -v[18] - v[21] - v[25] - v[29] - v[33] - v[37] - v[45] - v[68] - v[72] + v18d
    du[V.O] = -v[40] - v[54] - v[58] - v[60]
    du[V.A] = -v[41] - v[46] - v[55] - v[59]
    du[V.E11G] = v[16] - v16d
    du[V.E11S] = v[17] - v17d
    du[V.E11R] = v[18] - v18d
    du[V.E12G] = v[19]
    du[V.E12S] = v[20]
    du[V.E12R] = v[21]
    du[V.E23G] = v[22]
    du[V.E23S] = v[23]
    du[V.E23I] = v[24]
    du[V.E23R] = v[25]
    du[V.E34G] = v[26]
    du[V.E34S] = v[27]
    du[V.E34I] = v[28]
    du[V.E34R] = v[29]
    du[V.E24G] = v[30]
    du[V.E24S] = v[31]
    du[V.E24I] = v[32]
    du[V.E24R] = v[33]
    du[V.E44G] = v[34]
    du[V.E44S] = v[35]
    du[V.E44I] = v[36]
    du[V.E44R] = v[37]
    du[V.sigmaG] = v[16] + v[19] + v[22] + v[26] + v[30] + v[34] - v[40] - v[41] + v[42] + v[65] + v[69] - v16d
    du[V.sigmaS] = v[17] + v[20] + v[23] + v[27] + v[31] + v[35] - v[38] + v[43] + v[66] + v[70] - v17d
    du[V.sigmaI] = v[24] + v[28] + v[32] + v[36] + v[44] + v[67] + v[71]
    du[V.sigmaR] = v[18] + v[21] + v[25] + v[29] + v[33] + v[37] + v[45] - v[50] + v[68] + v[72] - v18d
    du[V.sigmaA] = -v[39] + v[41] + v[46] + v[59]
    du[V.sigmaSP] = v[38] - v[42]
    du[V.sigmaAP] = v[39] - v[43] - v[44] - v[45] - v[88]
    du[V.sigmaG_O] = v[40] - v[59]
    du[V.sigmaG_A] = v[41] - v[60]
    du[V.sigmaSP_G] = v[42]
    du[V.sigmaAP_S] = v[43]
    du[V.sigmaAP_I] = v[44]
    du[V.sigmaAP_R] = v[45]
    du[V.P3_A] = v[46] - v[57]
    du[V.P2] = -v[48]
    du[V.P3] = -v[46] + v[48]
    du[V.Akt] = -v[47]
    du[V.RsD] = -v[49]
    du[V.RsT] = v[49]
    du[V.sigmaRP] = v[50]
    du[V.Raf] = -v[51]
    du[V.Rafstar] = v[51]
    du[V.MEK] = -v[52]
    du[V.MEKstar] = v[52] - v[89] + v[90] - v[91] + v[92]
    du[V.ERK] = -v[89] + v[96]
    du[V.ERKstar] = v[92] - v[93]
    du[V.OP] = v[54]
    du[V.AP] = v[55]
    du[V.A_sigmaG_O] = v[59] + v[60]
    du[V.sigmaA_G] = v[57] - v[58]
    du[V.sigmaA_G_O] = v[58]
    du[V.sigmaO] = v[40] + v[58] + v[60]
    du[V.E13] = v[61] - v[63]
    du[V.E14] = v[62] - v[64]
    du[V.E13P] = v[63] - v[65] - v[66] - v[67] - v[68] - v[79]
    du[V.E14P] = v[64] - v[69] - v[70] - v[71] - v[72] - v[80]
    du[V.E13G] = v[65]
    du[V.E13S] = v[66]
    du[V.E13I] = v[67]
    du[V.E13R] = v[68]
    du[V.E14G] = v[69]
    du[V.E14S] = v[70]
    du[V.E14I] = v[71]
    du[V.E14R] = v[72]
    du[V.T] = -v[73] - v[74] - v[75] - v[76] - v[77] - v[78] - v[79] - v[80] - v[88] + v73d
    du[V.E11T] = v[73] - v73d
    du[V.E12T] = v[74]
    du[V.E23T] = v[75]
    du[V.E34T] = v[76]
    du[V.E24T] = v[77]
    du[V.E44T] = v[78]
    du[V.E13T] = v[79]
    du[V.E14T] = v[80]
    du[V.sigmaT] = v[73] + v[74] + v[75] + v[76] + v[77] + v[78] + v[79] + v[80] + v[88] - v73d
    du[V.E1_PT] = v[81] - v[86]
    du[V.E2_PT] = v[82]
    du[V.E4_PT] = v[83] - v[87]
    du[V.E_E1_PT] = v[84] + v[86]
    du[V.H_E4_PT] = v[85] + v[87]
    du[V.Aktstar] = v[47]
    du[V.sigmaAP_T] = v[88]
    du[V.E] = (-v[1] - v[86])/VeVc
    du[V.H] = (-v[2] - v[3] - v[87])/VeVc
    du[V.fint] = p[C.a98]*(-u[V.fint] + p[C.b98])
    du[V.pERK] = v[90] - v[91] + v[94] - v[95]
    du[V.ERK_MEKstar] = v[89] - v[90]
    du[V.pERK_MEKstar] = v[91] - v[92]
    du[V.pERK_ERKPpase] = v[95] - v[96]
    du[V.ERKPpase] = -v[93] + v[94] - v[95] + v[96]
    du[V.ERKstar_ERKPpase] = v[93] - v[94]

end

function param_values()

    p::Vector{Float64} = zeros(C.NUM)

    p[C.VmaxPY] = 223.8776
    p[C.KmPY] = 486.1398
    p[C.kdeg] = 0.0259
    p[C.kf47] = 24.6048
    p[C.Vmaxr47] = 590.5058
    p[C.Kmf47] = 698.6022
    p[C.Kmr47] = 483.8622
    p[C.kf48] = 16.833
    p[C.Kmf48] = 715.5688
    p[C.Kmr48] = 324.9294
    p[C.PTEN] = 693.5786
    p[C.kf49] = 44.3501
    p[C.kr49] = 552.6746
    p[C.Kmf49] = 343.2483
    p[C.Kmr49] = 753.1667
    p[C.Kmr49b] = 381.2208
    p[C.kr49b] = 640.8212
    p[C.kf51] = 3.6515
    p[C.Vmaxr51] = 16.737
    p[C.Kmf51] = 599.7076
    p[C.Kmr51] = 346.4779
    p[C.Kmrb51] = 988.4496
    p[C.kf52] = 0.7745
    p[C.Vmaxr52] = 199.2773
    p[C.Kmf52] = 545.4408
    p[C.Kmr52] = 675.2994
    p[C.kf54] = 0.0538
    p[C.Vmaxr54] = 588.2671
    p[C.Kmf54] = 457.9645
    p[C.Kmr54] = 336.183
    p[C.kf55] = 0.2256
    p[C.Vmaxr55] = 646.9003
    p[C.Kmf55] = 460.9446
    p[C.Kmr55] = 643.9247
    p[C.kf38] = 279.9929
    p[C.kf39] = 385.7428
    p[C.kf50] = 389.1061
    p[C.a98] = 0.0849
    p[C.b98] = 0.1833
    p[C.koff46] = 0.5194
    p[C.EGF_off] = 0.0175
    p[C.HRGoff_3] = 9.0E-4
    p[C.HRGoff_4] = 0.0973
    p[C.koff4] = 0.1717
    p[C.koff5] = 4.3985
    p[C.koff6] = 2.6619
    p[C.koff7] = 8.0557
    p[C.koff8] = 9.1034
    p[C.koff9] = 5.5425
    p[C.koff57] = 0.4526
    p[C.koff16] = 0.5737
    p[C.koff17] = 4.6874
    p[C.koff18] = 2.2768
    p[C.koff19] = 2.3361
    p[C.koff20] = 0.6761
    p[C.koff21] = 4.7291
    p[C.koff22] = 3.6962
    p[C.koff23] = 2.3619
    p[C.koff24] = 4.4226
    p[C.koff25] = 2.225
    p[C.koff26] = 0.0103
    p[C.koff27] = 1.8922
    p[C.koff28] = 4.6432
    p[C.koff29] = 2.0148
    p[C.koff30] = 4.9936
    p[C.koff31] = 1.2204
    p[C.koff32] = 3.8752
    p[C.koff33] = 1.2817
    p[C.koff34] = 3.2036
    p[C.koff35] = 1.8696
    p[C.koff36] = 1.2567
    p[C.koff37] = 0.4059
    p[C.koff65] = 0.1185
    p[C.koff66] = 2.2988
    p[C.koff67] = 1.6142
    p[C.koff40] = 3.1051
    p[C.koff41] = 7.0487
    p[C.koff42] = 3.5195
    p[C.koff43] = 0.5441
    p[C.koff44] = 0.4265
    p[C.koff45] = 3.9967
    p[C.koff58] = 6.3059
    p[C.koff59] = 9.172
    p[C.koff68] = 2.8871
    p[C.kPTP10] = 29.8531
    p[C.kPTP11] = 78.204
    p[C.kPTP12] = 11.4211
    p[C.kPTP13] = 55.2104
    p[C.kPTP14] = 57.7506
    p[C.kPTP15] = 60.2628
    p[C.kPTP63] = 7.4766
    p[C.kPTP64] = 48.6335
    p[C.koff73] = 3.0048
    p[C.koff74] = 1.2496
    p[C.koff75] = 1.4323
    p[C.koff76] = 2.1542
    p[C.koff77] = 1.2237
    p[C.koff78] = 0.2007
    p[C.koff79] = 1.1852
    p[C.koff80] = 2.9373
    p[C.kPTP38] = 83.4465
    p[C.kPTP39] = 79.6132
    p[C.koff88] = 3.9255
    p[C.kPTP50] = 96.5716
    p[C.kf81] = 1.361
    p[C.Vmaxr81] = 242.6034
    p[C.Kmf81] = 485.2626
    p[C.Kmr81] = 323.4012
    p[C.kf82] = 6.9987
    p[C.Vmaxr82] = 398.1931
    p[C.Kmf82] = 781.4374
    p[C.Kmr82] = 595.8395
    p[C.kf83] = 1.7628
    p[C.Vmaxr83] = 534.0531
    p[C.Kmf83] = 609.4766
    p[C.Kmr83] = 653.5184
    p[C.kf84] = 4.6894
    p[C.Vmaxr84] = 634.1626
    p[C.Kmf84] = 622.3847
    p[C.Kmr84] = 258.4637
    p[C.kf85] = 6.7591
    p[C.Vmaxr85] = 369.2261
    p[C.Kmf85] = 179.6486
    p[C.Kmr85] = 290.7667
    p[C.kcon49] = 9.9783
    p[C.kon1] = 1.0086E-4
    p[C.kon86] = 0.0038
    p[C.kon2] = 0.0059
    p[C.kon3] = 0.0334
    p[C.kon87] = 8.0E-4
    p[C.kon4] = 0.5005
    p[C.kon5] = 2.5427
    p[C.kon6] = 0.2283
    p[C.kon7] = 1.0606
    p[C.kon8] = 1.0259
    p[C.kon9] = 2.2868
    p[C.kon57] = 0.0039
    p[C.kf10] = 0.6496
    p[C.kf11] = 0.3721
    p[C.kf12] = 1.8012
    p[C.kf13] = 0.8875
    p[C.kf14] = 6.1726
    p[C.kf15] = 1.3565
    p[C.kf63] = 0.9297
    p[C.kf64] = 1.2083
    p[C.kon16] = 0.0097
    p[C.kon17] = 0.0049
    p[C.kon18] = 0.0117
    p[C.kon73] = 0.0116
    p[C.kon19] = 0.0896
    p[C.kon20] = 0.0478
    p[C.kon21] = 0.0114
    p[C.kon74] = 0.0133
    p[C.kon22] = 7.0E-4
    p[C.kon23] = 0.0138
    p[C.kon24] = 0.005
    p[C.kon25] = 0.0995
    p[C.kon75] = 0.0137
    p[C.kon26] = 0.0355
    p[C.kon27] = 0.0201
    p[C.kon28] = 0.0074
    p[C.kon29] = 0.0346
    p[C.kon76] = 0.0053
    p[C.kon30] = 0.002
    p[C.kon31] = 0.0032
    p[C.kon32] = 9.0E-4
    p[C.kon33] = 0.0335
    p[C.kon77] = 0.0101
    p[C.kon34] = 1.0E-4
    p[C.kon35] = 0.0602
    p[C.kon36] = 0.0043
    p[C.kon37] = 0.0791
    p[C.kon78] = 0.0076
    p[C.kon79] = 0.0078
    p[C.kon65] = 0.0123
    p[C.kon66] = 1.9264E-4
    p[C.kon67] = 6.6667E-5
    p[C.kon80] = 2.0E-4
    p[C.kon40] = 0.0191
    p[C.kon41] = 0.0051
    p[C.kon42] = 0.0023
    p[C.kon43] = 0.0127
    p[C.kon44] = 0.0122
    p[C.kon45] = 0.0028
    p[C.kon88] = 0.0108
    p[C.kon46] = 0.0148
    p[C.kon58] = 0.0215
    p[C.kon59] = 0.0077
    p[C.kon60] = 1.1994E-4
    p[C.koff60] = 4.9981
    p[C.koff61] = 5.229
    p[C.kon61] = 0.8048
    p[C.kon62] = 1.782
    p[C.koff62] = 5.5142
    p[C.kon68] = 0.0045
    p[C.kon69] = 0.0084
    p[C.koff69] = 3.97
    p[C.kon70] = 0.0116
    p[C.koff70] = 2.6069
    p[C.kon71] = 0.0078
    p[C.koff71] = 2.2988
    p[C.kon72] = 0.0355
    p[C.koff72] = 0.907
    p[C.kon89] = 0.1997
    p[C.koff89] = 99.9637
    p[C.kcat90] = 20.0037
    p[C.kon91] = 0.1966
    p[C.koff91] = 99.9983
    p[C.kcat92] = 0.2004
    p[C.kon93] = 0.2003
    p[C.koff93] = 100.0037
    p[C.kcat94] = 0.9966
    p[C.kon95] = 0.1993
    p[C.koff95] = 100.0023
    p[C.kcat96] = 19.9851

    return p
end

function initial_values()

    u0::Vector{Float64} = zeros(V.NUM)

    u0[V.E1] = 274.0
    u0[V.E2] = 158.0
    u0[V.E3] = 294.0
    u0[V.E4] = 399.0
    u0[V.G] = 82.4
    u0[V.S] = 11.5
    u0[V.I] = 46.4
    u0[V.R] = 93.6
    u0[V.T] = 500.0
    u0[V.O] = 82.3
    u0[V.A] = 43.1
    u0[V.P2] = 197
    u0[V.Akt] = 444.0
    u0[V.RsD] = 95.7
    u0[V.Raf] = 743.0
    u0[V.MEK] = 772.0
    u0[V.ERK] = 750.0
    u0[V.ERKPpase] = 35.0

    return u0
end


text:
using Test

include("adaptive_arithmetic.jl")
include("geometry_predicates.jl")
include("triangulation.jl")
include("delaunay.jl")
include("biased_random_order.jl")

max_stars_repo_name: PaulMatlashewski/DelaunayTriangulations

max_issues_repo_name: PaulMatlashewski/DelaunayTriangulations

max_forks_repo_name: PaulMatlashewski/DelaunayTriangulations


text:
## BEFORE RUNNING THIS SCRIPT, MAKE SURE ros is in the environment

## Get user commands

# Populates `parsed_args`
include(joinpath(dirname(@__DIR__),"parsecommands.jl"))

# assume in Atom editor (not scripted use)
if length(ARGS) == 0
  parsed_args["folder_name"] = "labrun8"
  parsed_args["remoteprocs"] = 10
  parsed_args["vis2d"] = true
end

## load main process processes only

if parsed_args["vis2d"]
  using Cairo
  using Gadfly
  using RoMEPlotting
  Gadfly.set_default_plot_size(35cm, 20cm)
end

# create context for vis3d
vis = if parsed_args["vis3d"]
  using MeshCat, GeometryTypes, ColorTypes, CoordinateTransformations
  vis = Visualizer()
  open(vis)
  vis
end

## bring in all the required source code

# bring up the multiprocess cluster
include(joinpath(@__DIR__, "CarFeedCommonSetup.jl"))

# a few required utils
include(joinpath(@__DIR__, "CarSlamUtilsMono.jl"))

##

gui, canvases = if parsed_args["imshow"]
  gui = imshow_gui((1000, 200), (1, 2))  # 2 columns, 1 row of images (each initially 300×300)
  canvases = gui["canvas"]
  # allow gui time to load
  sleep(1.0)
  gui, canvases
else
  nothing, nothing
end

tools = RacecarTools(detector)

## SLAM portion

slam = SLAMWrapperLocal()
getSolverParams(slam.dfg).drawtree = true
getSolverParams(slam.dfg).showtree = false

dfg_datafolder = getLogPath(slam.dfg)
datastore = FileDataStore("$dfg_datafolder/bigdata")

# also store parsed_args used in this case
fid = open(joinLogPath(slam.dfg,"args.json"),"w")
println(fid, JSON2.write(parsed_args))
close(fid)

# TODO add to results.log
fid = open(joinpath(dirname(getLogPath(slam.dfg)),"results.log"),"a")
resdirname = splitpath(getLogPath(slam.dfg))[end]
thisfile = splitpath(@__FILE__)[end]
println(fid, "$resdirname -- $thisfile, $(parsed_args["folder_name"]), $ARGS")
close(fid)

bagSubscriber = RosbagSubscriber(bagfile)

syncz = SynchronizeCarMono(30,syncList=[:leftFwdCam;:camOdo])
fec = FrontEndContainer(slam,bagSubscriber,syncz,tools,datastore)

# callbacks via Python
bagSubscriber(leftimgtopic, leftImgHdlr, fec)
bagSubscriber(zedodomtopic, odomHdlr, fec, WEIRDOFFSET)
bagSubscriber(joysticktopic, joystickHdlr, fec)

defaultFixedLagOnTree!(slam.dfg, 50, limitfixeddown=true)
# getSolverParams(slam.dfg).dbg = true
getSolverParams(slam.dfg).drawtree = true
getSolverParams(slam.dfg).showtree = false

##

# (fec.synchronizer.leftFwdCam |> last)[2] |> collect |> imshow

# consume a few messages to find initialization timestamp
for i in 1:100
  @show i, bagSubscriber.nextMsgChl
  0 == length(fec.synchronizer.leftFwdCam) || 0 == length(fec.synchronizer.camOdo) ? nothing : break
  loop!(bagSubscriber)
end

## initialize the factor graph

T0 = nanosecond2datetime(fec.synchronizer.leftFwdCam[1][2])

addVariable!(slam.dfg, :x0, Pose2, timestamp=T0)
addFactor!(slam.dfg, [:x0], PriorPose2(MvNormal(zeros(3),diagm([0.1,0.1,0.01].^2))), timestamp=T0)

## start solver

ST = manageSolveTree!(slam.dfg, slam.solveSettings, dbg=false)

##

sleep(0.01)  # allow gui sime time to setup
while loop!(bagSubscriber)
# for i in 1:2000
#   loop!(bagSubscriber)
  blockProgress(slam) # required to prevent duplicate solves occuring at the same time
end

## close all

stopManageSolveTree!(slam)

parsed_args["vis3d"] ? delete!(vis) : nothing

# wiat for previous and then ensure all data is solved at least once
@info "CarFeedMono.jl is waiting on penultimate trigger solve to finalize."
blockSolvingInProgress(fec.slam)
checkSolveStrideTrigger!(fec.slam, force=true)
sleep(10)
@info "CarFeedMono.jl is waiting on last trigger solve to finalize."
blockSolvingInProgress(fec.slam)

## interpose results

allD = jsonResultsSLAM2D(fec)

allStr = JSON2.write(allD)

fid = open(joinLogPath(fec.slam.dfg, "$(runfile)_results.json"),"w")
println(fid, allStr)
close(fid)

## save the factor graph

saveDFG(fec.slam.dfg, joinLogPath(fec.slam.dfg, "fg_$(slam.poseCount)"))

## draw results

if parsed_args["vis2d"]

# drawPoses(slam.dfg, spscale=0.3, drawhist=false)
pl = drawPosesLandms(fec.slam.dfg, spscale=0.3, drawhist=false)
pl |> PDF(joinLogPath(fec.slam.dfg,"fg_$(slam.poseCount).pdf"))
# drawPosesLandms(fg4, spscale=0.3, drawhist=false)

if parsed_args["report_factors"]
  reportFactors(fec.slam.dfg, Pose2Pose2, show=false)
end

end

## batch resolve

if parsed_args["batch_resolve"]

# fg2 = deepcopy(fec.slam.dfg)

dontMarginalizeVariablesAll!(fec.slam.dfg)
foreach(x->setSolvable!(fec.slam.dfg, x, 1), ls(fec.slam.dfg))
foreach(x->setSolvable!(fec.slam.dfg, x, 1), lsf(fec.slam.dfg))

tree, smt, hist = solveTree!(fec.slam.dfg)

saveDFG(fec.slam.dfg, joinLogPath(fec.slam.dfg, "fg_$(slam.poseCount)_resolve"))

end

## after resolve interpose results

allD = jsonResultsSLAM2D(fec)

allStr = JSON2.write(allD)

fid = open(joinLogPath(fec.slam.dfg, "$(runfile)_results_batchresolve.json"),"w")
println(fid, allStr)
close(fid)

## draw a second time

if parsed_args["batch_resolve"] && parsed_args["vis2d"]

pl = drawPosesLandms(fec.slam.dfg, spscale=0.3, drawhist=false)
pl |> PDF(joinLogPath(slam.dfg,"fg_$(slam.poseCount)_resolve.pdf"))

end

#


text:
import ..BifurcationsBase

BifurcationsBase.aliasof(::Codim2Sweep) = "Codim2Sweep"
BifurcationsBase.aliasof(::Codim2Solution) = "Codim2Solution"
BifurcationsBase.aliasof(::Codim2Solver) = "Codim2Solver"


text:
module Scanf

using Base.Ryu

export @scanf, @sscanf, scanf

const EOF = -1
const ARGNUM_TYPE = UInt16
const WIDTH_TYPE = UInt32
# whitespace characters in format strings
const WHITESPACE = b" \n\t\r\f\v"
# format escape character
const ESC = UInt8('%')
# format specifier to skip whitespace
const SKIP = UInt8('_')
# format indicator no-assign
const NOASSIGN = UInt8('*')
# format specifiers for character sets
const CSOPEN = UInt8('[')
const CSCLOSE = UInt8(']')
const CSNEG = UInt8('^')
const CSMINUS = UInt8('-')

# format specifier categories
const Ints = Union{Val{'d'},Val{'i'},Val{'u'},Val{'x'},Val{'X'},Val{'o'}}
const Floats = Val{'f'}
const Chars = Val{'c'}
const Strings = Val{'s'}
const Pointer = Val{'p'}
const HexBases = Union{Val{'x'},Val{'X'},Val{'a'},Val{'A'}}
const DecBases = Union{Val{'d'},Val{'u'}}
const OctBases = Val{'o'}
const IntBases = Val{'i'}
const PositionCounter = Val{'n'}
const Whitespace = Val{Char(SKIP)}
const CharSets = Union{Val{Char(CSOPEN)},Val{Char(CSNEG)}}
const EscapeChar = Val{Char(ESC)}

const DECIMAL = b"0123456789"
const OCTAL = b"01234567"
const HEXADECIMAL = b"0123456789ABCDEFabcdef"
abstract type AbstractSpec{T} end

"""
    Spec{T}

Typed representation of a format specifier.

`T` is a `Val{'X'}`, where `X` is a valid format specifier character.

Special case `X == '_'` represents whitespace (not assigning)

See also `CharsetSpec` and `LiteralSpec`.
"""
struct Spec{T<:Val} <: AbstractSpec{T} # T => %type => Val{'type'}
    assign::ARGNUM_TYPE
    width::WIDTH_TYPE
end

"""
    CharsetSpec{T,S}

Format specifier representing a character set (inclusive or exclusive)
"""
struct CharsetSpec{T<:Val,S} <: AbstractSpec{T}
    assign::ARGNUM_TYPE
    width::WIDTH_TYPE
    set::S
end

"""
    LiteralSpec{S<:AbstractString}

Non-assigning Format specifier representing a literal matching string.
"""
struct LiteralSpec{S} <: AbstractSpec{0}
    string::S
end

# replace double %% by single %
function LiteralSpec(str::T) where {T<:AbstractString}
    if contains(str, "%%")
        str = replace(str, "%%" => "%")
        S = typeof(str)
    else
        S = T
    end
    LiteralSpec{S}(str)
end

"""
    Scanf.Format(format_str)

Create a C scanf-compatible format object that can be used for parse formatted texts.

The input `format_str` can include any valid format specifier character and modifiers.

A `Format` object can be passed to `scanf(str::String, f::Format, args...)` to parse a
formatted string, or `scanf(io::IO, f::Format, args...)` to parse the
formatted string directly from `io`.

For convenience, the `Scanf.format"..."` string macro form can be used for building
a `Scanf.Format` object at macro-expansion-time.

The `@scanf` macro converts a literal string into a `Scanf.Format` implicitly.
"""
struct Format{S,T}
    str::S # original full format string as CodeUnits
    formats::T # Tuple of Specs (including whitespace specs and literal specs)
end

"""
    Scanf.scanf(b::String, f::Scanf.Format, args...)
    Scanf.scanf([io::IO,] f::Scanf.Format, args...)

Apply a Scanf format object `f` to provided String
(1st method), or scan directly from an `io` object (2nd method). See [`@scanf`](@ref)
for more details on C `scanf` support.
`io` defaults to `stdin`.
The args determine the type and default value for the output data.

Return the number of assigned arguments, followed by the assigned output data.
"""
function scanf end

scanf(io::IO, f::Format, args...) = scanner(io, f, args...)
scanf(f::Format, args...) = scanf(stdin, f, args...)
scanf(str::AbstractString, f::Format, args...) = scanf(IOBuffer(str), f, args...)

"""
    @scanf([io:Union{IO,String}, ], "%Fmt", args::...)

Scan input stream or string of using C `scanf` style format specification string and assign values
to arguments.

The format string (not a variable) is analyzed at macro expansion time.
Equivalent to `scanf(io, Scanf.format"%Fmt", args...)`.

# Examples
```jldoctest
r, a... = @scanf("%f%f%f%f", zeros(4)...)

julia> r, a, b = @scanf "23.4 text" "%f %2s" Float64 String
(2, 23.4, "te")
```
"""
macro scanf(arg, args...)
    n = length(args)
    if arg isa String && !(n > 0 && args[1] isa String)
        fmt = arg
        io = :stdin
    elseif n >= 1 && args[1] isa String
        fmt = args[1]
        io = arg
        args = Base.tail(args)
    else
        throwa("invalid macro arguments: @scanf $arg $(join(args, ' '))")
    end
    f = Format(unescape_string(fmt))
    return esc(:($Scanf.scanf($io, $f, $(args...))))
end

"""
    Scanf.format"..."

Convenience string macro to call `Scanf.Format("...")
"""
macro format_str(str)
    str = unescape_string(str)
    esc(:(Scanf.Format($str)))
end

# Implementation details

# construct Format object by parsing the format string
function Format(f::AbstractString)
    isempty(f) && throwa("empty format string")
    bytes = codeunits(f)
    len = length(bytes)
    fmts = AbstractSpec[]
    pos, b = pushliteral!(fmts, f, bytes, 1)
    ac = 0 # count number of assigning specifiers

    while pos <= len || b == SKIP
        assign = true
        width = 0
        charset = nothing
        if b == ESC
            b = bytes[pos]
            pos += 1
            # positioned at start of first format str %
            # parse flags
            if b == NOASSIGN
                assign = false
                pos > len && throwa("incomplete format string: '$f'")
                b = bytes[pos]
                pos += 1
            end
            # parse width
            while 0x00 <= b - UInt8('0') < 0x0a
                width = 10 * width + (b - UInt8('0'))
                b = bytes[pos]
                pos += 1
                pos > len && break
            end
            if !(0 <= width <= typemax(WIDTH_TYPE))
                throwi("\"$f\", width not in 1..$(typemax(Int32))")
            end
            # type modifier characters (ignored)
            if b == UInt8('h') || b == UInt8('l')
                prev = b
                pos > len && throwa("incomplete format string: \"$f\"")
                b = bytes[pos]
                pos += 1
                if b == prev # cases "ll" and "hh"
                    pos > len && throwa("incomplete format string: \"$f\"")
                    b = bytes[pos]
                    pos += 1
                end
            elseif b in b"Ljqtz" # more type modifiers (ignored)
                b = bytes[pos]
                pos += 1
            end
            # parse conversion specifier
            if b in b"diouxXcspn" # uper case variants are invalid
                type = Val{lowercase(Char(b))}
            elseif b in b"eEfFgGaA"
                type = Val{Char('f')}
            elseif b == CSOPEN
                txt = "["
                start = pos
                xpos = findnext(isequal(CSCLOSE), bytes, start)
                if start <= len && bytes[start] == CSNEG
                    b = CSNEG
                    txt = "[^"
                    start += 1
                end
                if xpos === nothing
                    throwi("\"$f\", after '$txt' a Char(CSCLOSE) is missing")
                end
                if bytes[start] == CSCLOSE # first character after "[" or "[^"
                    xpos = findnext(isequal(CSCLOSE), bytes, start + 1)
                end
                if (xpos === nothing || start >= xpos)
                    throwi("\"$f\", after '$txt]' a Char(CSCLOSE) is missing")
                end
                charset = view(bytes, start:xpos-1)
                pos = xpos + 1
                type = Val{Char(b)}
            else
                throwi("\"$f\", invalid conversion specifier: '$(Char(b))'")
            end
        else # format contains a WS character
            type = Whitespace
            assign = false
        end
        ac += assign
        ass = ifelse(assign, ac, 0)
        if ass > typemax(ARGNUM_TYPE)
            throwi("\"$f\", too many assignable conversion specifiers")
        end
        push!(fmts, make_spec(type, ass, width, charset))
        pos, b = pushliteral!(fmts, f, bytes, pos)
    end
    return Format(bytes, Tuple(fmts))
end

# consume characters in format string up to next % or whitespace and insert LiteralSpec
@inline function pushliteral!(fmts, f, bytes, pos)
    len = length(bytes)
    start = pos
    b = 0x00
    j = -1
    while pos <= len
        b = bytes[pos]
        pos += 1
        if b == ESC
            pos > len && throwa("incomplete format string: \"$f\"")
            if bytes[pos] == ESC # "%%" will be removed in LiteralSpec
                pos += 1
            else
                j = pos - 2
                break
            end
        elseif b in WHITESPACE
            j = pos - 2
            while pos <= len && bytes[pos] in WHITESPACE
                pos += 1
            end
            b = SKIP
            break
        end
    end
    j = j < 0 ? pos - 1 : j
    start <= j && push!(fmts, LiteralSpec(SubString(f, start:j)))
    return pos, b
end

Format(f::Format) = f

# extract text items from input stream

# match literal strings
@inline function fmt(io, pos, spec::LiteralSpec)
    str = spec.string
    len = sizeof(str)
    ix = 1
    while ix <= len && !eof(io)
        c = str[ix]
        peek(io, Char) != c && break
        skip(io, ncodeunits(c))
        ix = nextind(str, ix)
    end
    l = ix - 1
    IOBuffer(), l == len, pos + l
end

# skip whitespace spec
@inline function fmt(io, pos, ::Spec{Whitespace})
    pos = skip_ws(io, pos)
    IOBuffer(), true, pos
end

# string spec
@inline function fmt(io, pos, spec::Spec{T}) where {T<:Strings}
    pos = skip_ws(io, pos)
    width = spec.width
    l = 0
    out = IOBuffer()
    while (width == 0 || l < width) && !eof(io)
        c = peek(io, Char)
        (isspace(c) || !isvalid(Char, c)) && break
        write(out, c)
        n = ncodeunits(c)
        skip(io, n)
        l += 1
    end
    out, l > 0, pos + l
end

# charset and chars specs
@inline function fmt(io, pos, spec::S) where {S<:Union{CharsetSpec,Spec{<:Chars}}}
    width = spec.width
    width = ifelse(width == 0, S <: Chars ? 1 : typemax(Int), width)
    l = 0
    out = IOBuffer()
    while l < width && !eof(io)
        c = peek(io, Char)
        n = ncodeunits(c)
        check_set(c, spec) || break
        skip(io, n)
        write(out, c)
        l += 1
        pos += 1
    end
    out, l > 0, pos
end

# integer specs
@inline function fmt(io, pos, spec::Spec{T}) where {T<:Ints}
    pos = skip_ws(io, pos)
    out = IOBuffer()
    eof(io) && return out, false
    width = spec.width
    width = ifelse(width == 0, typemax(Int), width)
    l = 0
    sig = false
    b = peek(io)
    negate = false
    if b in b"+-"
        l = writeskip(io, out, b, l)
        negate = b == UInt('-')
        sig = true
    end
    digits = DECIMAL
    base, digits = basespec(T)
    ndig = 0
    if base === nothing || base == 16
        b = peek(io)
        if b == UInt('0')
            ndig += 1
            l = writeskip(io, out, b, l)
            if !eof(io) && (b = peek(io)) in b"xX"
                digits = HEXADECIMAL
                l = writeskip(io, out, b, l)
                ndig = 0
            elseif base === nothing
                digits = OCTAL
            end
        end
    end
    while l < width && !eof(io)
        b = peek(io)
        if b in digits
            ndig += 1
            l = writeskip(io, out, b, l)
        else
            break
        end
    end
    succ = ndig > 0
    out, succ, pos + l
end

# pointer spec
@inline function fmt(io, pos, spec::Spec{T}) where {T<:Pointer}
    pos = skip_ws(io, pos)
    out = IOBuffer()
    eof(io) && return out, false
    width = spec.width
    width = ifelse(width == 0, typemax(Int), width)
    digits = HEXADECIMAL
    l = 0
    while l < width && !eof(io)
        b = peek(io)
        if l == 0 && b == UInt8('0')
            l = writeskip(io, out, b, l)
            if !eof(io) && (b = peek(io)) in b"xX"
                l = writeskip(io, out, b, l)
            end
        elseif b in digits
            l = writeskip(io, out, b, l)
        else
            break
        end
    end
    out, l > 2, pos + l
end

# floating point spec
@inline function fmt(io, pos, spec::Spec{T}) where {T<:Floats}
    pos = skip_ws(io, pos)
    width = spec.width
    width = ifelse(width == 0, typemax(Int), width)
    digits = DECIMAL
    expch = b"eEfF"
    x_sign = 0x01
    x_sep = 0x02
    x_exp = 0x04
    x_base = 0x08
    x_inexp = 0x10
    x_nexp = 0x20
    x_mdigits = 0x40
    x_edigits = 0x80
    l = 0
    out = IOBuffer()
    status = x_sign | x_sep | x_base | x_edigits
    while l < width && !eof(io)
        b = peek(io)
        if status & x_base != 0 && b == UInt8('0')
            l = writeskip(io, out, b, l)
            if !eof(io) && (b = peek(io)) in b"xX"
                digits = HEXADECIMAL
                expch = b"pP"
                l = writeskip(io, out, b, l)
            else
                status |= x_mdigits
            end
            status = (status | x_exp | x_edigits) & ~(x_base | x_sign)
            continue
        elseif b in digits
            status = (status | x_exp) & ~(x_base | x_sign)
            status |= (status & x_inexp != 0) ? x_edigits : x_mdigits
        elseif b == UInt8('-') && (status & x_sign) != 0
            status |= ifelse(status & x_base == 0, x_nexp, 0)
            status = status & ~x_sign
        elseif b == UInt8('+') && (status & x_sign) != 0
            status = status & ~x_sign
        elseif b == UInt8('.') && (status & x_sep) != 0
            status = status & ~(x_base | x_sign | x_sep)
        elseif b in expch && (status & x_exp) != 0
            status = (status & ~(x_base | x_exp | x_edigits | x_sep)) | x_sign | x_inexp
            digits = DECIMAL
        elseif b in b"iI"
            n = expect(io, out, b"INFINITY", 3)
            l = (n == 3 || n == 8) ? l + n : 0
            status |= x_mdigits | x_edigits
            break
        elseif b in b"nN"
            n = expect(io, out, b"NAN", 3)
            l = n == 3 ? l + n : 0
            status |= x_mdigits | x_edigits
            break
        else
            break
        end
        l = writeskip(io, out, b, l)
    end
    succ = l > 0 && (status & x_mdigits != 0) && (status & x_edigits != 0)
    out, succ, pos + l
end

# helpers for matching nan and inf
@inline function expect(io::IO, out::IO, bytes::AbstractVector{UInt8}, n)
    l = 1
    while l <= length(bytes) && !eof(io)
        b = peek(io)
        if b == bytes[l] || b == bytes[l] + (UInt8('a') - UInt8('A'))
            l <= n && write(out, b)
            skip(io, 1)
            l += 1
        else
            break
        end
    end
    return l - 1
end

# position counter spec
@inline function fmt(_, pos, ::Spec{PositionCounter})
    IOBuffer(string(pos - 1)), true, pos
end

# convert to output type specializations
function toout(::Type{<:Val}, out, arg::AbstractString)
    r = String(take!(out))
    oftype(arg, r)
end

function toout(::Type{<:Val}, out, arg::AbstractChar)
    r = String(take!(out))
    length(r) >= 1 ? oftype(arg, r[1]) : arg
end

function toout(::Type{<:Val}, out, arg::AbstractVector{<:AbstractChar})
    r = String(take!(out))
    resize!(arg, length(r))
    arg .= collect(r)
    arg
end

function toout(::Type{T}, out, arg::R) where {T<:Union{Ints,Val{'p'},Val{'n'}},R<:Union{Real,Ptr}}
    S = inttype(typeof(arg))
    r = String(take!(out))
    x = xparse(S, r, base = basespec(T)[1])
    x !== nothing ? convert(R, x) : nothing
end

function toout(::Type{T}, out, arg::R) where {T<:Floats,R<:Real}
    r = String(take!(out))
    A = floattype(typeof(arg))
    s = xparse(A, r)
    s !== nothing ? convert(R, s) : nothing
end

function toout(::Type{<:Strings}, out, arg::R) where R<:Real
    r = String(take!(out))
    xparse(R, r, base = R <: Integer ? 10 : nothing)
end

# fmt helpers

# try parse float - assume input string is syntaxtically correct
@inline function xparse(::Type{A}, r::AbstractString; base = nothing) where {A<:AbstractFloat}
    if 'f' in r || 'F' in r
        r = replace(r, r"[fF]" => 'e')
    end
    s = tryparse(A, r)
    if s === nothing
        expos = findfirst(x -> x in "eEfF", r)
        negx = expos !== nothing && expos < length(r) && r[nextind(r, expos)] == '-'
        s = ifelse(negx, zero(A), typemax(A))
        r[1] == '-' && (s = -s)
    end
    s
end

# try parse integers - assume input string is syntaxtically correct
@inline function xparse(::Type{S}, r::AbstractString; base = nothing) where {S<:Integer}
    sig = r[1] in "+-"
    negate = r[1] == '-'
    n = length(r)
    if n >= 1 && r[1+sig] == '0'
        base = n >= 2+sig && r[2+sig] in "xX" ? nothing : base === nothing ? 8 : base
    end
    if sig && S <: Unsigned
        r = SubString(r, 2:n)
    end
    x = tryparse(S, r, base = base)
    if x !== nothing &&  negate && S <: Unsigned
        x = -x
    end
    x
end

# skip WS characters in io stream
@inline function skip_ws(io, pos)
    while !eof(io)
        b = peek(io)
        n = _ncodeunits(b)
        if n == 1
            if b in WHITESPACE
                skip(io, 1)
                pos += 1
            else
                break
            end
        else
            c = peek(io, Char)
            if isspace(c)
                skip(io, n)
                pos += 1
            else
                break
            end
        end
    end
    pos
end

# consume 1 byte
@inline function writeskip(io, out, b, l)
    write(out, b)
    skip(io, 1)
    l + 1
end

itemtype(::AbstractSpec{T}) where {T<:Union{Strings,CharSets}} = Val{'s'}
itemtype(::AbstractSpec{T}) where {T} = T

# call the format specifiers aligned with the IO stream
@inline function scanner(io::IO, f::Format, args...)
    n = length(args)
    m = countspecs(f)
    n == m || argmismatch(m, n)
    res = Any[valuefor(x) for x in args]
    eof(io) && return EOF, res...
    pos = 1

    # for each format, scan arg and next substring
    # unroll up to 8 formats
    formats = f.formats
    N = length(formats)
    j = 0
    UNROLL_UPTO = 8
    Base.@nexprs 8 i ->
    begin if N >= i
            fi = formats[i]
            out, succ, pos = fmt(io, pos, fi)
            succ || @goto BREAK
            assign, argno = assignnr(fi)
            if assign
                r = toout(itemtype(fi), out, res[argno])
                r === nothing && @goto BREAK
                res[argno] = r
                j += 1
            end
        end
    end
    if N > UNROLL_UPTO
        for i = UNROLL_UPTO+1:N
            fi = formats[i]
            out, succ, pos = fmt(io, pos, fi)
            succ || break
            assign, argno = assignnr(fi)
            if assign
                r = toout(itemtype(fi), out, res[argno])
                r === nothing && break
                res[argno] = r
                j += 1
            end
        end
    end
    @label BREAK
    return tuple(j, res...)
end

# utility functions

# default value for types
valuefor(v::AbstractVector{T}) where {T<:AbstractChar} = resize!(v, 0)
valuefor(::Type{T}) where {T<:Union{Real,Ptr,Char}} = T(0)
valuefor(::Type{T}) where {T<:AbstractString} = T("")
valuefor(a::T) where {T<:Union{Real,AbstractChar,AbstractString,Ptr}} = a

# accessor functions
assignnr(::LiteralSpec) = false, 0
assignnr(spec::AbstractSpec) = !iszero(spec.assign), Int(spec.assign)
assign(spec::AbstractSpec) = assignnr(spec)[1]

# showing formats

# recreate the format specifier string from a typed Spec
Base.string(f::Spec{Val{T}}) where {T} = string(string_header(f), T)
function Base.string(f::CharsetSpec{T}) where T<:CharSets
    string(string_header(f),
        Char(CSOPEN),
        T <: Val{Char(CSNEG)} ? Char(CSNEG) : "",
        showset(f.set),
        Char(CSCLOSE),
    )
end
Base.string(f::LiteralSpec) = string('"', f.string, '"')

function string_header(f::AbstractSpec)
    string(Char(ESC), !assign(f) ? "*" : "", f.width > 0 ? f.width : "")
end

Base.show(io::IO, f::AbstractSpec) = print(io, string(f))
Base.show(io::IO, f::Format) = print(io, string("Scanf.format\"", String(f.str), '"'))

# reconstruct internals of %[...]
showset(s::UnitRange) = "$(Char(first(s)))-$(Char(last(s)))"

function showset(t::Tuple{String,Vararg})
    s = t[1]
    m, s = special(Char(CSMINUS), s)
    b, s = special(Char(CSCLOSE), s)
    string(b, s, showset.(t[2:end])..., m)
end

showset(t::Tuple) = string(showset.(t)...)

showset(s::AbstractString) = isempty(s) ? "1-0" : s

function special(x::Char, s::AbstractString)
    if occursin(x, s)
        string(x), filter(!isequal(x), s)
    else
        "", s
    end
end

# internal representation of % and whitespace format specifiers
function make_spec(T, assign, width, cs)
    if cs === nothing || isempty(cs)
        Spec{T}(assign, width)
    else
        charset = make_charset(cs)
        CharsetSpec{T,typeof(charset)}(assign, width, charset)
    end
end

# details of internal representation of %[...] format specifiers
make_charset(cs::AbstractVector{UInt8}) = make_charset(String(cs))
function make_charset(cs::String)
    M = Char(CSMINUS)
    ch = String[]
    ra = Any[]
    start = 1
    stop = lastindex(cs)
    mind = findall(isequal(M), cs)
    a = start
    for i in mind
        if i > start && i < stop
            j = prevind(cs, i)
            k = nextind(cs, i)
            cj, ck = cs[j], cs[k]
            if cj < ck
                push!(ra, UInt32(cj):UInt32(ck))
                b = prevind(cs, j)
            elseif cj == ck
                b = j
            else
                b = prevind(cs, j)
            end
            push!(ch, SubString(cs, a:b))
            a = nextind(cs, k)
        end
    end
    push!(ch, SubString(cs, a:stop))
    str = String(unique(collect(string(ch...))))
    isempty(ra) ? str : isempty(str) ? (length(ra) == 1 ? ra[1] : (ra...,)) : (str, ra...)
end

# bases for integer specs
basespec(::Type{<:HexBases}) = 16, HEXADECIMAL
basespec(::Type{<:OctBases}) = 8, OCTAL
basespec(::Type{<:DecBases}) = 10, DECIMAL
basespec(::Type) = nothing, DECIMAL

# type conversion hints for integer specs
inttype(::Type{<:Float64}) = Int64
inttype(::Type{<:Float32}) = Int32
inttype(::Type{<:Float16}) = Int16
inttype(::Type{<:BigFloat}) = BigInt
inttype(::Type{<:Ptr}) = UInt
inttype(::Type{T}) where {T} = T
floattype(::Type{T}) where {T<:Integer} = float(T)
floattype(::Type{T}) where {T} = T

# check if character is in or not in charset
@inline check_set(c::Char, spec::Spec{<:Chars}) = true
@inline check_set(c::Char, spec::CharsetSpec{Val{Char(CSOPEN)}}) = check_in(c, spec.set)
@inline check_set(c::Char, spec::CharsetSpec{Val{Char(CSNEG)}}) = check_!in(c, spec.set)
@inline check_in(c, set::AbstractString) = occursin(c, set)
@inline check_!in(c, set::AbstractString) = !occursin(c, set)
@inline check_in(c, set::UnitRange) = in(UInt32(c), set)
@inline check_!in(c, set::UnitRange) = !in(UInt32(c), set)
@inline check_in(c, set::Tuple) = any(x -> check_in(c, x), set)
@inline check_!in(c, set::Tuple) = all(x -> check_!in(c, x), set)

# count number of assigning format specifiers in format
countspecs(f::Format) = count(assign, f.formats)

import Base: ==
function ==(f::T, g::T) where {T<:Format}
    f.formats == g.formats
end

# length of utf8 encoding of character starting with this byte
@inline _ncodeunits(a::UInt8) = a < 0xc0 ? 1 : a < 0xe0 ? 2 : a < 0xf0 ? 3 : 4

# error
@noinline function argmismatch(a, b)
    throwa("mismatch between # of format specifiers and provided args: $a != $b")
end

@noinline throwa(msg) = throw(ArgumentError(msg))
@noinline throwi(msg) = throwa("invalid format string: " * msg)

end # module


text:
# This is a temporary module to validate the ideas in
# https://github.com/JuliaImages/ImagesAPI.jl/pull/3
module NoiseAPI

"""
    AbstractImageAlgorithm

The root of image algorithms type system
"""
abstract type AbstractImageAlgorithm end

"""
    AbstractImageFilter <: AbstractImageAlgorithm

Filters are image algorithms whose input and output are both images
"""
abstract type AbstractImageFilter <: AbstractImageAlgorithm end

include("./utils.jl")
include("./ImageNoise.jl")

end