text:
function check_simulation_chronology(
    horizons::OrderedDict{Symbol, Int},
    intervals::OrderedDict{Symbol, Dates.Millisecond},
    resolutions::OrderedDict{Symbol, Dates.Millisecond},
)
    models = collect(keys(resolutions))

    for (model, horizon) in horizons
        horizon_time = resolutions[model] * horizon
        if horizon_time < intervals[model]
            throw(IS.ConflictingInputsError("horizon ($horizon_time) is
                                shorter than interval ($interval) for $(model)"))
        end
    end

    for i in 2:length(models)
        upper_level_model = models[i - 1]
        lower_level_model = models[i]
        if horizons[lower_level_model] * resolutions[lower_level_model] >
           horizons[upper_level_model] * resolutions[upper_level_model]
            throw(
                IS.ConflictingInputsError(
                    "The lookahead length $(horizons[upper_level_model]) in model $(upper_level_model) is insufficient to syncronize with $(lower_level_model)",
                ),
            )
        end
        if (intervals[upper_level_model] % intervals[lower_level_model]) !=
           Dates.Millisecond(0)
            throw(
                IS.ConflictingInputsError(
                    "The system's intervals are not compatible for simulation. The interval in model $(upper_level_model) needs to be a mutiple of the interval $(lower_level_model) for a consistent time coordination.",
                ),
            )
        end
    end
    return
end

"""
_calculate_interval_inner_counts(intervals::OrderedDict{String,<:Dates.TimePeriod})

Calculates how many times a problem is executed for every interval of the previous problem
"""
function _calculate_interval_inner_counts(intervals::OrderedDict{Symbol, Dates.Millisecond})
    order = collect(keys(intervals))
    reverse_order = length(intervals):-1:1
    interval_run_counts = Vector{Int}(undef, length(intervals))
    interval_run_counts[1] = 1
    for k in reverse_order[1:(end - 1)]
        model_name = order[k]
        previous_model_name = order[k - 1]
        problem_interval = intervals[model_name]
        previous_problem_interval = intervals[previous_model_name]
        if Dates.Millisecond(previous_problem_interval % problem_interval) !=
           Dates.Millisecond(0)
            throw(
                IS.ConflictingInputsError(
                    "The interval configuration provided results in a fractional number of executions of problem $model_name",
                ),
            )
        end
        interval_run_counts[k] = previous_problem_interval / problem_interval
        @debug "problem $k is executed $(interval_run_counts[k]) time within each interval of problem $(k-1)"
    end
    return interval_run_counts
end

"""
Function calculates the total number of problem executions in the simulation and allocates the appropiate vector
"""
function _allocate_execution_order(interval_run_counts::Vector{Int})
    total_size_of_vector = 0
    for k in eachindex(interval_run_counts)
        mult = 1
        for i in 1:k
            mult *= interval_run_counts[i]
        end
        total_size_of_vector += mult
    end
    return -1 * ones(Int, total_size_of_vector)
end

function _fill_execution_order(
    execution_order::Vector{Int},
    interval_run_counts::Vector{Int},
)
    function _fill_problem(index::Int, problem::Int)
        last_problem = problems[end]
        if problem < last_problem
            next_problem = problem + 1
            for i in 1:interval_run_counts[next_problem]
                index = _fill_problem(index, next_problem)
            end
        end
        execution_order[index] = problem
        index -= 1
    end

    index = length(execution_order)
    problems = sort!(collect(keys(interval_run_counts)))
    _fill_problem(index, problems[1])
    return
end

function _get_execution_order_vector(intervals::OrderedDict{Symbol, Dates.Millisecond})
    length(intervals) == 1 && return [1]
    interval_run_counts = _calculate_interval_inner_counts(intervals)
    execution_order_vector = _allocate_execution_order(interval_run_counts)
    _fill_execution_order(execution_order_vector, interval_run_counts)
    @assert isempty(findall(x -> x == -1, execution_order_vector))
    return execution_order_vector
end

function _get_num_executions_by_model(
    models::SimulationModels,
    execution_order::Vector{Int},
)
    model_names = get_model_names(models)
    executions_by_model = OrderedDict(x => 0 for x in model_names)
    for number in execution_order
        executions_by_model[model_names[number]] += 1
    end
    return executions_by_model
end

function _attach_feedforwards(models::SimulationModels, feedforwards)
    names = Set(string.(get_model_names(models)))
    ff_dict = Dict{Symbol, Vector}()
    for (model_name, model_feedforwards) in feedforwards
        if model_name ∈ names
            model_name_symbol = Symbol(model_name)
            ff_dict[model_name_symbol] = model_feedforwards
            for ff in model_feedforwards
                sim_model = get_simulation_model(models, model_name_symbol)
                device_model = get_model(get_template(sim_model), get_component_type(ff))
                if device_model === nothing
                    throw(
                        IS.ConflictingInputsError(
                            "Device model $(get_component_type(ff)) not found in model $model_name",
                        ),
                    )
                end
                attach_feedforward!(device_model, ff)
            end
        else
            error("Model $model_name not present in the SimulationModels")
        end
    end
    return ff_dict
end

@doc raw"""
    SimulationSequence(
                        models::SimulationModels,
                        feedforward::Dict{Symbol, <:AbstractAffectFeedforward}
                        ini_cond_chronology::Dict{Symbol, <:FeedforwardChronology}
                        )
"""
mutable struct SimulationSequence
    horizons::OrderedDict{Symbol, Int}
    intervals::OrderedDict{Symbol, Dates.Millisecond}
    feedforwards::Dict{Symbol, Vector{<:AbstractAffectFeedforward}}
    ini_cond_chronology::InitialConditionChronology
    execution_order::Vector{Int}
    executions_by_model::OrderedDict{Symbol, Int}
    current_execution_index::Int64
    uuid::Base.UUID

    function SimulationSequence(;
        models::SimulationModels,
        feedforwards=Dict{String, Vector{<:AbstractAffectFeedforward}}(),
        ini_cond_chronology=InterProblemChronology(),
    )
        # Allow strings or symbols as keys; convert to symbols.
        intervals = determine_intervals(models)
        horizons = determine_horizons!(models)
        resolutions = determine_resolutions(models)

        if length(models.decision_models) > 1
            check_simulation_chronology(horizons, intervals, resolutions)
        end

        if length(models.decision_models) == 1
            # TODO: Not implemented yet
            # ini_cond_chronology = IntraProblemChronology()
        end

        execution_order = _get_execution_order_vector(intervals)
        executions_by_model = _get_num_executions_by_model(models, execution_order)
        sequence_uuid = IS.make_uuid()
        initialize_simulation_internals!(models, sequence_uuid)
        new(
            horizons,
            intervals,
            _attach_feedforwards(models, feedforwards),
            ini_cond_chronology,
            execution_order,
            executions_by_model,
            0,
            sequence_uuid,
        )
    end
end

get_step_resolution(sequence::SimulationSequence) = first(values(sequence.intervals))

function get_interval(sequence::SimulationSequence, problem::Symbol)
    return sequence.intervals[problem]
end

function get_interval(sequence::SimulationSequence, model::DecisionModel)
    return sequence.intervals[get_name(model)]
end

get_execution_order(sequence::SimulationSequence) = sequence.execution_order


text:
"""
`pandoc2html(formatted::AbstractString, doc::WeaveDoc)`

Convert output from pandoc markdown to html using Weave.jl template
"""
function pandoc2html(formatted::AbstractString, doc::WeaveDoc, outname::AbstractString, pandoc_options)
  weavedir = dirname(@__FILE__)
  html_template = joinpath(weavedir, "../templates/pandoc_skeleton.html")
  css_template = joinpath(weavedir, "../templates/pandoc_skeleton.css")
  css = stylesheet(MIME("text/html"), doc.highlight_theme)

  path, wsource = splitdir(abspath(doc.source))
  #wversion = string(Pkg.installed("Weave"))
  wversion = ""
  wtime =  string(Date(now()))

  #Header is inserted from displayed plots
  header_script = doc.header_script
  self_contained = (header_script ≠ "") ? [] : "--self-contained"

  if haskey(doc.header, "bibliography")
    filt = "--filter"
    citeproc = "pandoc-citeproc"
  else
    filt = []
    citeproc = []
  end

  #Change path for pandoc
  old_wd = pwd()
  cd(doc.cwd)
  html =""
  outname = basename(outname)

  open("temp.md", "w") do io
    println(io, formatted)
  end

  try
    cmd = `pandoc -f markdown+raw_html -s --mathjax=""
    $filt $citeproc $pandoc_options
    --template $html_template -H $css_template $self_contained
     -V wversion=$wversion -V wtime=$wtime -V wsource=$wsource
     -V highlightcss=$css
     -V headerscript=$header_script
     -o $outname`
     proc = open(cmd, "r+")
     println(proc.in, formatted)
     close(proc.in)
     proc_output = read(proc.out, String)
     cd(old_wd)
  catch e
    cd(old_wd)
    @warn("Error converting document to HTML")
    throw(e)
  end
end

"""
`pandoc2pdf(formatted::AbstractString, doc::WeaveDoc)`

Convert output from pandoc markdown to pdf using Weave.jl template
"""
function pandoc2pdf(formatted::AbstractString, doc::WeaveDoc, outname::AbstractString, pandoc_options)
  weavedir = dirname(@__FILE__)
  header_template = joinpath(weavedir, "../templates/pandoc_header.txt")

  path, wsource = splitdir(abspath(doc.source))
  #wversion = string(Pkg.installed("Weave"))
  wversion = ""
  wtime =  Date(now())
  outname = basename(outname)

  #Change path for pandoc
  old_wd = pwd()
  cd(doc.cwd)
  html =""

  if haskey(doc.header, "bibliography")
    filt = "--filter"
    citeproc = "pandoc-citeproc"
  else
    filt = []
    citeproc = []
  end

  @info("Done executing code. Running xelatex")
  try
    cmd = `pandoc -f markdown+raw_tex -s  --pdf-engine=xelatex --highlight-style=tango
     $filt $citeproc $pandoc_options
     --include-in-header=$header_template
     -V fontsize=12pt -o $outname`
    proc = open(cmd, "r+")
    println(proc.in, formatted)
    close(proc.in)
    proc_output = read(proc.out, String)
    cd(old_wd)
  catch e
    cd(old_wd)
    @warn("Error converting document to pdf")
    throw(e)
  end
end

function run_latex(doc::WeaveDoc, outname, latex_cmd = "xelatex")
  old_wd = pwd()
  cd(doc.cwd)
  xname = basename(outname)
  @info("Weaved code to $outname. Running $latex_cmd")
  textmp = mktempdir(".")
  try
    out = read(`$latex_cmd -shell-escape $xname -aux-directory $textmp -include-directory $(doc.cwd)`, String)
    out = read(`$latex_cmd -shell-escape $xname -aux-directory $textmp -include-directory $(doc.cwd)`, String)
    rm(xname)
    rm(textmp, recursive=true)
    cd(old_wd)
    return true
  catch e
    @warn("Error converting document to pdf. Try running latex manually")
    cd(old_wd)
    rm(textmp)
    return false
  end
end

ext: jl


text:
# This file contains the static systems of Causal.

import UUIDs: uuid4

"""
    @def_static_system ex

where `ex` is the expression to define to define a new AbstractStaticSystem component type. The usage is as follows:
```julia
@def_source struct MyStaticSystem{T1,T2,T3,...,TN,OP, RO} <: AbstractStaticSystem
    param1::T1 = param1_default     # optional field
    param2::T2 = param2_default     # optional field
    param3::T3 = param3_default     # optional field
        ⋮
    paramN::TN = paramN_default     # optional field
    input::IP = input_default       # mandatory field
    output::OP = output_default     # mandatory field
    readout::RO = readout_function  # mandatory field
end
```
Here, `MyStaticSystem` has `N` parameters, an `output` port, an `input` port and a `readout` function.

!!! warning
    `input`, `output` and `readout` are mandatory fields to define a new static system. The rest of the fields are the parameters of the system.

!!! warning
    `readout` must be a two-argument function, i.e. a function of time `t` and input value `u`.

!!! warning
    New static system must be a subtype of `AbstractStaticSystem` to function properly.

# Example
```jldoctest
julia> @def_static_system struct MyStaticSystem{IP, OP, RO} <: AbstractStaticSystem
       α::Float64 = 1.
       β::Float64 = 2.
       input::IP = Inport()
       output::OP = Outport()
       readout::RO = (t,u) -> α * u[1] + β * u[2]
       end

julia> sys = MyStaticSystem();

julia> sys.α
1.0

julia> sys.input
1-element Inport{Inpin{Float64}}:
 Inpin(eltype:Float64, isbound:false)
```
"""
macro def_static_system(ex)
    ex.args[2].head == :(<:) && ex.args[2].args[2] in [:AbstractStaticSystem, :AbstractMemory] ||
        error("Invalid usage. The type should be a subtype of AbstractStaticSystem or AbstractMemory.\n$ex")

    fields = quote
        trigger::$(TRIGGER_TYPE_SYMBOL) = Inpin()
        handshake::$(HANDSHAKE_TYPE_SYMBOL) = Outpin{Bool}()
        callbacks::$(CALLBACKS_TYPE_SYMBOL) = nothing
        name::Symbol = Symbol()
        id::$(ID_TYPE_SYMBOL) = Causal.uuid4()
    end, [TRIGGER_TYPE_SYMBOL, HANDSHAKE_TYPE_SYMBOL, CALLBACKS_TYPE_SYMBOL, ID_TYPE_SYMBOL]
    _append_common_fields!(ex, fields...)
    deftype(ex)
end

##### Define prototipical static systems.

"""
    StaticSystem(; readout, input, output)

Consructs a generic static system with `readout` function, `input` port and `output` port.

# Example
```jldoctest
julia> ss = StaticSystem(readout = (t,u) -> u[1] + u[2], input=Inport(2), output=Outport(1));

julia> ss.readout(0., ones(2))
2.0
```
"""
@def_static_system struct StaticSystem{RO, IP, OP} <: AbstractStaticSystem
    readout::RO
    input::IP
    output::OP
end

@doc raw"""
    Adder(signs=(+,+))

Construts an `Adder` with input bus `input` and signs `signs`. `signs` is a tuplle of `+` and/or `-`. The output function `g` of `Adder` is of the form,
```math
    y = g(u, t) =  \sum_{j = 1}^n s_k u_k
```
where `n` is the length of the `input`, ``s_k`` is the `k`th element of `signs`, ``u_k`` is the `k`th value of `input` and ``y`` is the value of `output`. The default value of `signs` is all `+`.

# Example
```jldoctest
julia> adder = Adder(signs=(+, +, -));

julia> adder.readout([3, 4, 5], 0.) == 3 + 4 - 5
true
```
"""
@def_static_system struct Adder{S, IP, OP, RO} <: AbstractStaticSystem
    signs::S = (+, +)
    input::IP = Inport(length(signs))
    output::OP = Outport()
    readout::RO = (u, t, signs=signs) -> sum([sign(val) for (sign, val) in zip(signs, u)])
end

@doc raw"""
    Multiplier(ops=(*,*))

Construts an `Multiplier` with input bus `input` and signs `signs`. `signs` is a tuplle of `*` and/or `/`. The output function `g` of `Multiplier` is of the form,
```math
    y = g(u, t) =  \prod_{j = 1}^n s_k u_k
```
where `n` is the length of the `input`, ``s_k`` is the `k`th element of `signs`, ``u_k`` is the `k`th value of `input` and ``y`` is the value of the `output`. The default value of `signs` is all `*`.

# Example
```jldoctest
julia> mlt = Multiplier(ops=(*, *, /));

julia> mlt.readout([3, 4, 5], 0.) == 3 * 4 / 5
true
```
"""
@def_static_system struct Multiplier{S, IP, OP, RO} <: AbstractStaticSystem
    ops::S = (*,*)
    input::IP = Inport(length(ops))
    output::OP = Outport()
    readout::RO = (u, t, ops=ops) -> begin
        ops = ops
        val = 1
        for i = 1 : length(ops)
            val = ops[i](val, u[i])
        end
        val
    end
end

@doc raw"""
    Gain(input; gain=1.)

Constructs a `Gain` whose output function `g` is of the form
```math
    y = g(u, t) =  K u
```
where ``K`` is `gain`, ``u`` is the value of `input` and `y` is the value of `output`.

# Example
```jldoctest
julia> K = [1. 2.; 3. 4.];

julia> sfunc = Gain(input=Inport(2), gain=K);

julia> sfunc.readout([1., 2.], 0.) == K * [1., 2.]
true
```
"""
@def_static_system struct Gain{G, IP, OP, RO} <: AbstractStaticSystem
    gain::G = 1.
    input::IP = Inport()
    output::OP = Outport(length(gain * zeros(length(input))))
    readout::RO = (u, t, gain=gain) -> gain * u
end

@doc raw"""
    Terminator(input::Inport)

Constructs a `Terminator` with input bus `input`. The output function `g` is eqaul to `nothing`. A `Terminator` is used just to sink the incomming data flowing from its `input`.
"""
@def_static_system struct Terminator{IP, OP, RO} <: AbstractStaticSystem
    input::IP = Inport()
    output::OP = nothing
    readout::RO = nothing
end

"""
    Memory(delay=1.; initial::AbstractVector{T}=zeros(1), numtaps::Int=5, t0=0., dt=0.01, callbacks=nothing,
        name=Symbol()) where T

Constructs a 'Memory` with input bus `input`. A 'Memory` delays the values of `input` by an amount of `numdelay`.
`initial` determines the transient output from the `Memory`, that is, until the internal buffer of `Memory` is full,
the values from `initial` is returned.

# Example
```jldoctest
julia> Memory(delay=0.1)
Memory(delay:0.1, numtaps:5, input:Inport(numpins:1, eltype:Inpin{Float64}), output:Outport(numpins:1, eltype:Outpin{Float64}))

julia> Memory(delay=0.1, numtaps=5)
Memory(delay:0.1, numtaps:5, input:Inport(numpins:1, eltype:Inpin{Float64}), output:Outport(numpins:1, eltype:Outpin{Float64}))
```
"""
@def_static_system struct Memory{D, IN, TB, DB, IP, OP, RO} <: AbstractMemory
    delay::D = 1.
    initial::IN = zeros(1)
    numtaps::Int = 5
    timebuf::TB = Buffer(numtaps)
    databuf::DB = length(initial) == 1 ? Buffer(numtaps) : Buffer(length(initial), numtaps)
    input::IP = Inport(length(initial))
    output::OP = Outport(length(initial))
    readout::RO = (u, t, delay=delay, initial=initial, numtaps=numtaps, timebuf=timebuf, databuf=databuf) -> begin
        if t <= delay
            return initial
        else
            tt = content(timebuf, flip=false)
            uu = content(databuf, flip=false)
            if length(tt) == 1
                return uu[1]
            end
            if ndims(databuf) == 1
                itp = CubicSplineInterpolation(range(tt[end], tt[1], length=length(tt)), reverse(uu), extrapolation_bc=Line())
                return itp(t - delay)
            else
                itp = map(row -> CubicSplineInterpolation(range(tt[end], tt[1], length=length(tt)), reverse(row), extrapolation_bc=Line()), eachrow(uu))
                return map(f -> f(t - delay), itp)
            end
        end
    end
end

@doc raw"""
    Coupler(conmat::AbstractMatrix, cplmat::AbstractMatrix)

Constructs a coupler from connection matrix `conmat` of size ``n \times n`` and coupling matrix `cplmat` of size ``d \times d``. The output function `g` of `Coupler` is of the form
```math
    y = g(u, t) = (E \otimes P) u
```
where ``\otimes`` is the Kronecker product, ``E`` is `conmat` and ``P`` is `cplmat`, ``u`` is the value of `input` and `y` is the value of `output`.
"""
@def_static_system struct Coupler{C1, C2, IP, OP, RO} <: AbstractStaticSystem
    conmat::C1 = [-1. 1; 1. 1.]
    cplmat::C2 = [1 0 0; 0 0 0; 0 0 0]
    input::IP = Inport(size(conmat, 1) * size(cplmat, 1))
    output::OP = Outport(size(conmat, 1) * size(cplmat, 1))
    readout::RO = typeof(conmat) <: AbstractMatrix{<:Real} ?
        ( (u, t, conmat=conmat, cplmat=cplmat) ->  kron(conmat, cplmat) * u ) :
        ( (u, t, conmat=conmat, cplmat=cplmat) ->  kron(map(f -> f(t), conmat), cplmat) * u )
end

@doc raw"""
    Differentiator(kd=1; callbacks=nothing, name=Symbol())

Consructs a `Differentiator` whose input output relation is of the form
```math
    y(t) = k_d \dot{u}(t)
```
where ``u(t)`` is the input and ``y(t)`` is the output and ``kd`` is the differentiation constant.
"""
@def_static_system struct Differentiator{IP, OP, RO} <: AbstractStaticSystem
    kd::Float64 = 1.
    t::Float64 = zeros(0.)
    u::Float64 = zeros(0.)
    input::IP = Inport()
    output::OP = Outport()
    readout::RO = (uu, tt, t=t, u=u, kd=kd) -> begin
        val = only(uu)
        sst = t[1]
        ssu = u[1]
        out = tt ≤ sst ? ssu : (val - ssu) / (tt - sst)
        t .= t
        u .= val
        kd * out
    end
end

##### Pretty-printing
show(io::IO, ss::StaticSystem) = print(io,"StaticSystem(readout:$(ss.readout), input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Adder) = print(io, "Adder(signs:$(ss.signs), input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Multiplier) = print(io, "Multiplier(ops:$(ss.ops), input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Gain) = print(io, "Gain(gain:$(ss.gain), input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Terminator) = print(io, "Terminator(input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Memory) =
    print(io, "Memory(delay:$(ss.delay), numtaps:$(length(ss.timebuf)), input:$(ss.input), output:$(ss.output))")
show(io::IO, ss::Coupler) = print(io, "Coupler(conmat:$(ss.conmat), cplmat:$(ss.cplmat))")
show(io::IO, ss::Differentiator) = print(io, "Differentiator(kd:$(ss.kd))")


text:
## make options for the stationary multichannel lengthscale GP
nminlog10λ, nmaxlog10λ = 2, 200
pnorm = 2.
Klog10λ = transD_GP.GP.Mat32()
λx,λy = 1, 1
x = 0:(0.005λx):λx
y = 0:(0.005λy):0.5λy
λlog10λ = [0.05maximum(x), 0.05maximum(y)]
demean = false
sdev_poslog10λ = [0.01maximum(y), 0.01maximum(x)]
log10bounds = [-1 -0.69; -1 -0.69]
δlog10λ = 0.1
sdev_proplog10λ = [0.1, 0.1]
xall = zeros(2,length(x)*length(y))
for i in 1:size(xall,2)
    xid, yid = Tuple(CartesianIndices((length(x),length(y)))[i])
    xall[:,i] = [x[xid]; y[yid]]
end
xbounds = zeros(Float64,size(xall, 1), 2)
for dim in 1:size(xall, 1)
    xbounds[dim,:] = [minimum(xall[dim,:]), maximum(xall[dim,:])]
end
## Initialize a model using these options
Random.seed!(12)
optlog10λ = transD_GP.OptionsStat(nmin = nminlog10λ,
                        nmax = nmaxlog10λ,
                        xbounds = xbounds,
                        fbounds = log10bounds,
                        xall = xall,
                        λ = λlog10λ,
                        δ = δlog10λ,
                        demean = demean,
                        sdev_prop = sdev_proplog10λ,
                        sdev_pos = sdev_poslog10λ,
                        pnorm = pnorm,
                        quasimultid = false,
                        K = Klog10λ,
                        timesλ = 3,
                        needλ²fromlog = true,
                        updatenonstat = true,
                        )
@time  log10λ = transD_GP.init(optlog10λ)
## make options for the nonstationary GP
nmin, nmax = 2, 200
fbounds = [-2. 2]
δ = 0.1
sdev_prop = [0.1]
sdev_pos = [0.05, 0.05]
K = transD_GP.GP.Mat32()
demean_ns = true
## Initialize model for the nonstationary GP
Random.seed!(13)
opt = transD_GP.OptionsNonstat(optlog10λ,
                        nmin = nmin,
                        nmax = nmax,
                        fbounds = fbounds,
                        δ = δ,
                        demean = demean_ns,
                        sdev_prop = sdev_prop,
                        sdev_pos = sdev_pos,
                        pnorm = pnorm,
                        K = K
                        )
@time m = transD_GP.init(opt, log10λ)
## run tests for the different McMC moves -- stationary GP changes have to update and propagate
# to the nonstationary GP
fracthresh = 0.05
@testset "GP and MCMC move do and undo state tests" begin
    @testset "init test" begin
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        for i = 1:150
            transD_GP.birth!(m, opt, log10λ)
        end
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftestm)) < 1e-12
    end
    @testset "birth tests" begin
        for i = 1:50
            transD_GP.birth!(log10λ, optlog10λ, m, opt)
        end
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
    end
    @testset "death tests" begin
        for i = 1:20
            transD_GP.death!(log10λ, optlog10λ, m, opt)
        end
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
    end
    # birth and death hold correct states if tests above passed
    @testset "undo birth" begin
        log10λold = deepcopy(log10λ)
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_birth!(log10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo_birth holds state if gotten till here
    @testset "undo multiple birth" begin
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        log10λold = deepcopy(log10λ)
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_birth!(log10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo_birth holds state as well for multiple births and deaths till here
    @testset "undo death" begin
        log10λold = deepcopy(log10λ)
        transD_GP.death!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_death!(log10λ, optlog10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo death holds state if here
    @testset "undo multiple death" begin
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        transD_GP.birth!(log10λ, optlog10λ, m, opt)
        log10λold = deepcopy(log10λ)
        transD_GP.death!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_death!(log10λ, optlog10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo_death holds state as well for multiple births and deaths till here
    @testset "property change" begin
        transD_GP.property_change!(log10λ, optlog10λ, m, opt)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
    end
    # property change works if here
    @testset "undo property change" begin
        log10λold = deepcopy(log10λ)
        transD_GP.property_change!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_property_change!(log10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo property change works if here
    @testset "position change" begin
        transD_GP.position_change!(log10λ, optlog10λ, m, opt)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
    end
    # position change works if here
    @testset "undo position change" begin
        log10λold = deepcopy(log10λ)
        transD_GP.position_change!(log10λ, optlog10λ, m, opt)
        transD_GP.undo_position_change!(log10λ, optlog10λ, m)
        transD_GP.sync_model!(log10λ, optlog10λ)
        ftest = transD_GP.testupdate(optlog10λ, log10λ)
        transD_GP.sync_model!(m, opt)
        ftestm = transD_GP.testupdate(opt, log10λ, m)
        @test mean(abs.((ftestm - m.fstar)./ftest)) < fracthresh
        @test norm(mean(ftest' - 0.5log10.(log10λ.fstar))) < 1e-12
        @test norm(mean(log10λold.fstar - log10λ.fstar)) < 1e-12
    end
    # undo position change works if here
end
# ## timing for birth
# @time for i = 1:150
#           transD_GP.birth!(log10λ, optlog10λ, m, opt)
# end
# ## timing for s birth in s model for ns mode
# log10λ = transD_GP.init(optlog10λ)
# m = transD_GP.init(opt, log10λ)
# for i = 1:148
#     transD_GP.birth!(log10λ, optlog10λ, m, opt)
# end
# @time for i = 1:148
#     transD_GP.birth!(m, opt, log10λ)
# end
# NTIMES = 20
# ntimes = 150
# T = zeros(NTIMES)
# for I = 1:NTIMES
#     T[I] = time()
#     for i = 1:ntimes
#         transD_GP.death!(log10λ, optlog10λ, m, opt)
#         transD_GP.birth!(log10λ, optlog10λ, m, opt)
#     end
#     T[I] = (time() - T[I])/2ntimes
# end
# @info "time for $ntimes birth/death is $(mean(T)) +- $(std(T)/sqrt(NTIMES))"
# ## plot
# l = log10.(log10λ.fstar.^0.5)
# fig,ax = plt.subplots(1,3, sharex=true, sharey=true)
# vmin, vmax = minimum(l[1,:]), maximum(l[1,:])
# im1 = ax[1].imshow(reshape(l[1,:],length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im1, ax=ax[1])
# ax[1].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], marker="+",c="r")
# ax[1].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], c=log10λ.ftrain[1,1:log10λ.n], alpha=0.8)
# vmin, vmax = minimum(l[2,:]), maximum(l[2,:])
# im2 = ax[2].imshow(reshape(l[2,:],length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im2, ax=ax[2])
# ax[2].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], marker="+",c="r")
# ax[2].scatter(log10λ.xtrain[2,1:log10λ.n], log10λ.xtrain[1,1:log10λ.n], c=log10λ.ftrain[2,1:log10λ.n], alpha=0.8,
#     vmin=vmin, vmax=vmax)
# vmin, vmax = minimum(m.fstar), maximum(m.fstar)
# im3 = ax[3].imshow(reshape(m.fstar,length(x), length(y)), extent=[y[1],y[end],x[end],x[1]],
#     vmin=vmin, vmax=vmax)
# colorbar(im3, ax=ax[3])
# ax[3].scatter(m.xtrain[2,1:m.n], m.xtrain[1,1:m.n], marker="+",c="r")
# ax[3].scatter(m.xtrain[2,1:m.n], m.xtrain[1,1:m.n], c=m.ftrain[1:m.n], alpha=0.8,
#         vmin=vmin, vmax=vmax)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
"""
    module Shapes

Definitions of various shapes.

All shapes have the signature: `shape(parameters::Tuple, x0, y0, θ=0)`, and some
have keyword argument `reference` which determines which point `(x0,y0)` is referring to.

The length of `parameters` depends on the shape. For example, `Circle` has just one
parameter, the radius, while `Ellipse` has two, the semi-major and semi-minor axes.
"""
module Shapes

export AbstractShape
export nsides

abstract type AbstractShape{NDIMS,NSIDES} end
Base.ndims(::AbstractShape{NDIMS,NSIDES}) where {NDIMS,NSIDES} = NDIMS
nsides(::AbstractShape{NDIMS,NSIDES}) where {NDIMS,NSIDES} = NSIDES

include("1D/Shapes.jl")
# include("2D/Shapes.jl")
# include("3D/Shapes.jl")

end #module


text:
using AsmMacro, Test, InteractiveUtils

@testset "@asm" begin
    @asm function add_loop_vec4(x::Ptr{Float64},n::Int,z::Ptr{Float64})
        movq(n, rcx)
        movapd(x[0*16], xmm0)
        movapd(x[1*16], xmm1)
        xorpd(xmm2,xmm2)
        xorpd(xmm3,xmm3)
        @loop
        addpd(xmm0,xmm2)
        addpd(xmm1,xmm3)
        dec(rcx)
        jnz(@loop)
        movapd(xmm2, z[0*16])
        movapd(xmm3, z[1*16])
    end

    x = [1.0,2.0,3.0,4.0]
    n = 10
    z = similar(x)
    add_loop_vec4(pointer(x),n,pointer(z))
    @test z == x*n
    @test occursin("\${:uid}", sprint(code_llvm, add_loop_vec4, (Ptr{Float64}, Int, Ptr{Float64})))

    @asm function add(z::Ptr{Int64}, x::Int64, y::Int64)
        addq(x, y)
        movq(y, z[0])
    end
    z = Int64[100]
    add(pointer(z), Int64(1), Int64(2))
    @test z[1] === Int64(3)

    @asm function add(z::Ptr{Int32}, x::Int32, y::Int32)
        addl(x, y)
        movl(y, z[0])
    end
    z = Int32[100]
    add(pointer(z), Int32(1), Int32(2))
    @test z[1] === Int32(3)

    @asm function add(z::Ptr{Float64}, x::Float64, y::Float64)
        vaddsd(xmm0, xmm1, xmm1)
        movq(xmm1, z[0])
    end
    z = Float64[100]
    add(pointer(z), Float64(1), Float64(2))
    @test z[1] === Float64(3.0)

    @asm function add(z::Ptr{Float32}, x::Float32, y::Float32)
        vaddss(xmm0, xmm1, xmm1)
        movq(xmm1, z[0])
    end
    z = Float32[100]
    add(pointer(z), Float32(1), Float32(2))
    @test z[1] === Float32(3.0)

    @test_throws Any (@eval @asm function add(x::Float64)
        ___
    end)

    @test_throws Any (@eval @asm function add(x::Float64)
        x :: 1
    end)

    @test_throws Any (@eval @asm function add(x::Float64)
                          foo(:: 1)
    end)

    @test_throws Any (@eval @asm function (x::Float64)
                          foo(:: 1)
    end)

    @test_throws Any (@eval @asm function f(x)
                          foo(:: 1)
    end)
end


text:
using MathOptInterface, GalacticOptim, Optim, Test, Random, Flux, ForwardDiff, Zygote
using Nonconvex

rosenbrock(x, p) =  (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
x0 = zeros(2)
_p  = [1.0, 100.0]
l1 = rosenbrock(x0, _p)

@testset "Optim, CMAEvolutionStrategy" begin
    f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
    prob = OptimizationProblem(f, x0, _p)
    Random.seed!(1234)
    sol = solve(prob, SimulatedAnnealing())
    @test 10*sol.minimum < l1

    Random.seed!(1234)
    prob = OptimizationProblem(f, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, SAMIN())
    @test 10*sol.minimum < l1

    using CMAEvolutionStrategy
    sol = solve(prob, CMAEvolutionStrategyOpt())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(rosenbrock, x0, _p)
    sol = solve(prob, Optim.NelderMead())
    @test 10*sol.minimum < l1
end

@testset "Constrained Optimisation" begin
    cons= (x,p) -> [x[1]^2 + x[2]^2]
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons)

    prob = OptimizationProblem(optprob, x0, _p)

    sol = solve(prob, Flux.ADAM(0.1), maxiters = 1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, Optim.BFGS())
    @test 10*sol.minimum < l1

    sol = solve(prob, Optim.Newton())
    @test 10*sol.minimum < l1

    sol = solve(prob, Optim.KrylovTrustRegion())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf], ucons = [Inf])
    sol = solve(prob, IPNewton())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons = [-5.0], ucons = [10.0])
    sol = solve(prob, IPNewton())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf], ucons = [Inf], lb = [-500.0,-500.0], ub=[50.0,50.0])
    sol = solve(prob, IPNewton())
    @test sol.minimum < l1

    function con2_c(x,p)
        [x[1]^2 + x[2]^2, x[2]*sin(x[1])-x[1]]
    end

    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= con2_c)
    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf,-Inf], ucons = [Inf,Inf])
    sol = solve(prob, IPNewton())
    @test 10*sol.minimum < l1

    cons_circ = (x,p) -> [x[1]^2 + x[2]^2]
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff();cons= cons_circ)
    prob = OptimizationProblem(optprob, x0, _p, lcons = [-Inf], ucons = [0.25^2])
    sol = solve(prob, IPNewton())
    @test sqrt(cons(sol.u,nothing)[1]) ≈ 0.25 rtol = 1e-6

    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())

    prob = OptimizationProblem(optprob, x0, _p)
    sol = solve(prob, Flux.ADAM(), maxiters = 1000, progress = false)
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, Optim.Fminbox())
    @test 10*sol.minimum < l1

    Random.seed!(1234)
    prob = OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, Optim.SAMIN())
    @test 10*sol.minimum < l1
end

@testset "MaxSense test" begin
    optprob = OptimizationFunction((x,p) -> -rosenbrock(x,p), GalacticOptim.AutoZygote())
    prob = OptimizationProblem(optprob, x0, _p; sense = GalacticOptim.MaxSense)

    import Ipopt
    sol = solve(prob, Ipopt.Optimizer())
    @test 10*sol.minimum < l1

    import NLopt
    sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
    @test 10*sol.minimum < l1

    sol = solve(prob, NelderMead())
    @test 10*sol.minimum < l1

    sol = solve(prob, BFGS())
    @test 10*sol.minimum < l1

    function g!(G, x)
        G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
        G[2] = 200.0 * (x[2] - x[1]^2)
    end
    optprob = OptimizationFunction((x,p) -> -rosenbrock(x,p), GalacticOptim.AutoZygote(), grad = g!)
    prob = OptimizationProblem(optprob, x0, _p; sense = GalacticOptim.MaxSense)
    sol = solve(prob, BFGS())
    @test 10*sol.minimum < l1
end

@testset "MinSense test" begin
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
    prob = OptimizationProblem(optprob, x0, _p; sense = GalacticOptim.MinSense)

    sol = solve(prob, Ipopt.Optimizer())
    @test 10*sol.minimum < l1

    sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(Ipopt.Optimizer, "max_cpu_time" => 60.0))
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0)
    sol = solve(prob, NLopt.Opt(:LN_BOBYQA, 2))
    @test 10*sol.minimum < l1

    sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LN_BOBYQA))
    @test 10*sol.minimum < l1

    sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
    @test 10*sol.minimum < l1

    sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LD_LBFGS))
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, NLopt.Opt(:LD_LBFGS, 2))
    @test 10*sol.minimum < l1

    sol = solve(prob, GalacticOptim.MOI.OptimizerWithAttributes(NLopt.Optimizer, "algorithm" => :LD_LBFGS))
    @test 10*sol.minimum < l1

    sol = solve(prob, NLopt.Opt(:G_MLSL_LDS, 2), local_method = NLopt.Opt(:LD_LBFGS, 2), maxiters=10000)
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0)
    sol = solve(prob, NLopt.LN_BOBYQA())
    @test 10*sol.minimum < l1

    sol = solve(prob, NLopt.LD_LBFGS())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(optprob, x0, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, NLopt.LD_LBFGS())
    @test 10*sol.minimum < l1

    sol = solve(prob, NLopt.G_MLSL_LDS(), local_method = NLopt.LD_LBFGS(), local_maxiters=10000, maxiters=10000, population=10)
    @test 10*sol.minimum < l1
end

# using QuadDIRECT
# sol = solve(prob, QuadDirect(); splits = ([-0.5, 0.0, 0.5],[-0.5, 0.0, 0.5]))
# @test 10*sol.minimum < l1

@testset "Evolutionary, BlackBoxOptim, Metaheuristics, Nonconvex, GCMAES, SpeedMapping, MultistartOptimization, NOMAD" begin
    optprob = OptimizationFunction(rosenbrock, GalacticOptim.AutoZygote())
    using Evolutionary
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p)
    sol = solve(prob, CMAES(μ =40 , λ = 100),abstol=1e-15)
    @test 10*sol.minimum < l1

    using BlackBoxOptim
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[0.8, 0.8])
    sol = solve(prob, BBO_adaptive_de_rand_1_bin_radiuslimited())
    @test 10*sol.minimum < l1

    using Metaheuristics
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p, lb=[-1.0, -1.0], ub=[1.5, 1.5])
    sol = solve(prob, ECA())
    @test 10*sol.minimum < l1

    sol = solve(prob, Metaheuristics.DE())
    @test 10*sol.minimum < l1

    sol = solve(prob, PSO())
    @test 10*sol.minimum < l1

    sol = solve(prob, ABC())
    @test 10*sol.minimum < l1

    sol = solve(prob, CGSA(N=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, SA())
    @test 10*sol.minimum < l1

    sol = solve(prob, WOA())
    @test 10*sol.minimum < l1

    sol = solve(prob, ECA())
    @test 10*sol.minimum < l1

    sol = solve(prob, Metaheuristics.DE(), use_initial=true)
    @test 10*sol.minimum < l1

    sol = solve(prob, PSO(), use_initial=true)
    @test 10*sol.minimum < l1

    sol = solve(prob, ABC(), use_initial=true)
    @test 10*sol.minimum < l1

    sol = solve(prob, CGSA(N=100), use_initial=true)
    @test 10*sol.minimum < l1

    sol = solve(prob, SA(), use_initial=true)
    @test 10*sol.minimum < l1

    sol = solve(prob, WOA(), use_initial=true)
    @test 10*sol.minimum < l1

    using ModelingToolkit
    f = OptimizationFunction(rosenbrock,ModelingToolkit.AutoModelingToolkit())
    prob = OptimizationProblem(f, x0, _p)
    @test_broken sol = solve(prob,Optim.Newton())

    #=
    ### Nonconvex test
    Nonconvex.@load MMA
    Nonconvex.@load Ipopt
    Nonconvex.@load NLopt
    Nonconvex.@load BayesOpt
    Nonconvex.@load Juniper
    Nonconvex.@load Pavito
    Nonconvex.@load Hyperopt
    Nonconvex.@load MTS
    prob = GalacticOptim.OptimizationProblem(optprob, x0, _p, lb = [-1.0,-1.0], ub = [1.5,1.5])

    sol = solve(prob, MMA02())
    @test 10*sol.minimum < l1

    sol = solve(prob, GCMMA())
    @test 10*sol.minimum < l1

    sol = solve(prob, IpoptAlg())
    @test 10*sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LN_NELDERMEAD))
    @test 10*sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LD_LBFGS))
    @test 10*sol.minimum < l1

    sol = solve(prob, MTSAlg())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler())
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    sol = solve(prob, MMA02(), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, GCMMA(), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, IpoptAlg(), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LN_NELDERMEAD), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, NLoptAlg(:LD_LBFGS), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, MTSAlg(), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=1000)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), maxiters=100)
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)))
    @test 10*sol.minimum < l1

    ### suboptions
    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=LHSampler(), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=CLHSampler(), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(NLoptAlg(:LN_NELDERMEAD)), sampler=GPSampler(), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=LHSampler(), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=CLHSampler(), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, HyperoptAlg(IpoptAlg()), sampler=GPSampler(), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(IpoptAlg()), sub_options=(;max_iter=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1

    sol = solve(prob, BayesOptAlg(NLoptAlg(:LN_NELDERMEAD)), sub_options=(;maxeval=100))
    @test 10*sol.minimum < l1
    =#

    using NOMAD
    f = OptimizationFunction(rosenbrock)

    prob = OptimizationProblem(f, x0, _p)
    sol = GalacticOptim.solve(prob,NOMADOpt())
    @test 10*sol.minimum < l1

    prob = OptimizationProblem(f, x0, _p; lb = [-1.0,-1.0], ub = [1.5,1.5])
    sol = GalacticOptim.solve(prob,NOMADOpt())
    @test 10*sol.minimum < l1

    using GCMAES
    f_ad = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
    f_noad = OptimizationFunction(rosenbrock)

    prob = GalacticOptim.OptimizationProblem(f_ad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
    sol = solve(prob, GCMAESOpt(), maxiters=1000)
    @test 10*sol.minimum < l1

    prob = GalacticOptim.OptimizationProblem(f_noad, x0, _p, lb=[-1.0, -1.0], ub=[1.0, 1.0])
    sol = solve(prob, GCMAESOpt(), maxiters=1000)
    @test 10*sol.minimum < l1

    using SpeedMapping
    f = OptimizationFunction(rosenbrock, GalacticOptim.AutoForwardDiff())
    prob = OptimizationProblem(f, x0, _p)
    sol = solve(prob,SpeedMappingOpt())

    prob = OptimizationProblem(f, x0, _p;lb = [-1.0,-1.0], ub = [1.5,1.5])
    sol = solve(prob,SpeedMappingOpt())

    f = OptimizationFunction(rosenbrock)
    prob = OptimizationProblem(f, x0, _p)
    sol = solve(prob,SpeedMappingOpt())

    prob = OptimizationProblem(f, x0, _p;lb = [-1.0,-1.0], ub = [1.5,1.5])
    sol = solve(prob,SpeedMappingOpt())

    using MultistartOptimization
    f = OptimizationFunction(rosenbrock,GalacticOptim.AutoForwardDiff())
    prob = GalacticOptim.OptimizationProblem(f, x0, _p, lb = [-1.0,-1.0], ub = [1.5,1.5])
    sol = solve(prob, MultistartOptimization.TikTak(100), local_method = NLopt.LD_LBFGS())
    @test 10*sol.minimum < l1
end


text:
# This file is part of IntegerTriangles.jl.
# Copyright Peter Luschny. License is MIT.

# Builds the module IntegerTriangles, runtests and perftests from the src modules.
#
# The source files are scanned and the source modules are joined together,
# not merely included. This architecture has three advantages:
# * Source and test modules can be placed in the same file.
# * A uniform layout of the source modules is enforced.
# * The build can be easily customized by filtering source files.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module BuildTriangles
using Dates, UUIDs
export build_all

cloudpath = "https://github.com/OpenLibMathSeq/IntegerTriangles.jl/blob/master/src/"

#  Version information
version() = v"0.3.0"
print("IntegerTriangles version $(version())\n")

pkgdir = dirname((@__DIR__))
@info("Package directory is: " * pkgdir)

srcdir = dirname((@__FILE__))
cd(srcdir)
@info("Working directory is: " * pwd())

tstdir = joinpath(pkgdir, "test")
@info("Test files are in:    " * tstdir)

docdir = joinpath(pkgdir, "docs")
@info("Docs directory is:    " * docdir)

docsrcdir = joinpath(docdir, "src")
@info("Docs sources are in:  " * docsrcdir)

const exclude = [
    "BuildTriangles.jl",
    "IntegerTriangles.jl",
    "TrianglesTypes.jl",
    "SequenceID.jl",
    "_EXPORT.jl",
    "_TEMP.jl",
    "_INDEX.jl",
    "_SINDEX.jl",
    "_IS.jl",
    "tempCodeRunnerFile.jl"
]

#     Plots  fails
const references =
    "Nemo,
    HTTP,
    DocStringExtensions,
    OrderedCollections,
    PrettyTables"

function header(f)
    println(f, "# This file is part of IntegerTriangles.")
    println(f, "# Copyright Peter Luschny. License is MIT.")
    println(f)
    println(f, "# Version of: UTC ", Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))
    println(f, "# ", UUIDs.uuid1())
    println(f)
    println(f, "# Do not edit this file, it is generated from the modules and will be overwritten!")
    println(f, "# Edit the modules in the src directory and build this file with BuildTriangles.jl!")
    println(f)
end

function sortnames()

    index = open("_INDEX.jl", "r")
    sindex = open("_SINDEX.jl", "w")

    dict = Dict{Int64,Array{String}}()

    for l ∈ eachline(index, keep=true)

        i = 0
        c = 2
        if occursin(r"^[ACFGILMPRSTV][0-9]{6},$", l)
            c = 0
        end
        if occursin(r"^is[0-9]{6},$", l)
            # || occursin(r"^TL[0-9]{6},$", l)
            # || occursin(r"^TF[0-9]{6},$", l)
            i = 2
            c = 0
        end

        c == 2 && (print(sindex, l); continue)

        Anum = l[1:end - 2]
        N = Meta.parse(Anum[2 + i:end - c])

        if !haskey(dict, N)
            dict[N] = [Anum]
        else
            push!(dict[N], Anum)
        end
    end

    # We have to avoid the final comma in the export list
    K = sort(collect(keys(dict)))
    for key ∈ K[1:end - 1]
        for s ∈ dict[key]
            print(sindex, s * ",")
        end
        println(sindex)
    end

    E = dict[K[end]]
    for s ∈ E[1:end - 1]
        print(sindex, s * ",")
    end
    println(sindex, E[end])

    close(index)
    close(sindex)
end

# Nota bene: we make use of the convention that an Triangles module
# is closed by "end # module" (and not merely by "end")!
function build_seq()

    tmp = open("_TEMP.jl", "w")
    exp = open("_EXPORT.jl", "w")

    seq_modules = filter!(s -> occursin(r"\.jl$", s), readdir(srcdir))

    for filename ∈ seq_modules
        filename ∈ exclude && continue
        path = joinpath(srcdir, filename)
        mod = open(path, "r")
        println(" + ", filename)
        println(tmp, "# *** ", splitdir(path)[2], " ****************")
        doc = false
        mlcomment = false

        for l ∈ eachline(mod, keep=true)

            n = lstrip(l)
            if mlcomment && startswith(n, "=#")
                mlcomment = false
                continue
            end
            if mlcomment || startswith(n, "#=")
                mlcomment = true
                continue
            end
            startswith(n, "# START-TEST") && break
            startswith(n, "end # module") && break

            startswith(n, '#') && continue
            startswith(n, "using") && continue
            startswith(n, "module") && continue
            startswith(n, "(@__DIR__)") && continue
            startswith(n, "__precompile__") && continue

            startswith(n, "\"\"\"") && (doc = !doc)
            if doc
                print(tmp, l)
                continue
            end

            if startswith(n, "export")
                print(exp, n)
                continue
            end

            print(tmp, n)
        end
        close(mod)
    end

    flush(exp)
    close(exp)
    flush(tmp)
    close(tmp)

    sleep(0.1)
    exp = open("_EXPORT.jl", "r")
    sor = open("_INDEX.jl", "w")

    s = ""
    for l ∈ eachline(exp, keep=true)
        n = lstrip(l)
        if startswith(n, "export")
            n = n[7:end]
        end
        s *= n
    end
    R = replace(s, ',' => ' ')
    T = sort(split(R))

    println(sor, "export ")
    for t ∈ T[1:end]
        println(sor, t, ",")
    end

    close(sor)
    close(exp)
    sleep(0.01)

    sortnames()

    tmp = open("_TEMP.jl", "r")
    sor = open("_SINDEX.jl", "r")
    target = joinpath(srcdir, "_IS.jl")
    olm = open(target, "w")
    header(olm)
    println(olm, "__precompile__()")

    println(olm, "module IntegerTriangles")
    println(olm, "using " * references)

    for l ∈ eachline(sor, keep=true)
        print(olm, l)
    end

    for l ∈ eachline(tmp, keep=true)
        print(olm, l)
    end
    print(olm, "end")

    close(sor)
    close(tmp)
    close(olm)
end

# Builds the file test/runtests.jl from the Triangles modules.
# Nota bene: we make use of the convention that a Triangles module
# has three functions: test(), demo() and perf(), in this order.
function build_test()
    path = joinpath(tstdir, "runtests.jl")
    o = open(path, "w")
    header(o)
    println(o, "module runtests")
    println(o, "tstdir = realpath(joinpath(dirname(@__FILE__)))")
    println(o, "srcdir = joinpath(dirname(tstdir), \"src\")")
    println(o, "tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)")
    println(o, "srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)")
    println(o, "using Nemo, Test, IntegerTriangles, PrettyTables")

    path = joinpath(tstdir, "runtests.jl")
    i = open(path, "r")
    buff = Array{String,1}()
    for l ∈ eachline(i, keep=true)
        n = lstrip(l)
        startswith(n, '#') && continue
        startswith(n, "module") && continue
        startswith(n, "using") && continue
        startswith(n, "export") && continue
        push!(buff, n)
    end
    j = length(buff)
    while j > 0 && buff[j] == ""
        j -= 1
    end

    for k ∈ 1:j - 1
        print(o, buff[k])
    end
    close(i)

    seq_modules = filter!(s -> occursin(r"\.jl$", s), readdir(srcdir))
    for filename ∈ seq_modules
        filename ∈ exclude && continue
        path = joinpath(srcdir, filename)
        i = open(path, "r")
        inside = false
        println(o, "# *** ", splitdir(path)[2], " *********")

        buff = Array{String,1}()
        for l ∈ eachline(i, keep=true)
            n = lstrip(l)
            startswith(n, '#') && continue
            b = startswith(n, "function test")
            if b
                inside = true
                continue
            end

            if inside
                c = startswith(n, "function") && n[10] ≠ 't'
                if c
                    break
                end
                push!(buff, n)
            end
        end

        j = length(buff)
        while j > 0 && buff[j] == ""
            j -= 1
        end

        for k ∈ 1:j - 1
            print(o, buff[k])
        end
        close(i)
    end
    print(o, "end # module")
    close(o)
end

# Builds the file test/perftests.jl.
function build_perf()
    path = joinpath(tstdir, "perftests.jl")
    o = open(path, "w")

    header(o)
    println(o, "tstdir = realpath(joinpath(dirname(@__FILE__)))")
    println(o, "srcdir = joinpath(dirname(tstdir), \"src\")")
    println(o, "tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)")
    println(o, "srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)")

    println(o, "module perftests")
    println(o, "using IntegerTriangles, Dates, InteractiveUtils")

    println(o, "InteractiveUtils.versioninfo()")
    println(o, "start = Dates.now()")

    seq_modules = filter!(s -> occursin(r"\.jl$", s), readdir(srcdir))
    for filename ∈ seq_modules
        filename ∈ exclude && continue
        path = joinpath(srcdir, filename)
        i = open(path, "r")
        inside = false
        println(o, "# +++ ", filename, " +++")

        buff = Array{String,1}()
        s = ""
        for l ∈ eachline(i, keep=true)
            n = lstrip(l)
            b = startswith(n, "function perf()")
            if b
                inside = true
                continue
            end
            if inside
                c = startswith(n, "function main")
                if c
                    break
                end
                startswith(n, '#') && continue
                if startswith(n, "@time")
                    s = chomp(n[7:end]) * "\")\n"
                    push!(buff, "println(\"\\nTEST: ", s)
                end
                push!(buff, n)
            end
        end

        j = length(buff)
        while j > 0 && buff[j] == ""
            j -= 1
        end

        for k ∈ 1:j - 1
            print(o, buff[k])
        end
        close(i)
    end

    println(o, "stop = Dates.now()")
    println(o, "tdiff = stop - start")
    println(o, "println(\"\\nJulia version: \" * string(VERSION) )")
    println(o, "println(start)")
    println(o, "println(\"Total test time: \", tdiff)")
    println(o, "end # module")
    close(o)
end

# Builds the file index.md in docs/src.
function make_index()
    path = joinpath(docsrcdir, "index.md")
    ind = open(path, "w")
    tind = open("_SINDEX.jl", "r")

    first = true
    for l ∈ eachline(tind, keep=false)
        if first
            first = false
            println(ind, "# Library")
            println(ind, "")
            continue
        end
        endswith(l, ',') && (l = chop(l))
        for f ∈ split(l, ",")
            println(ind, "```@docs")
            println(ind, f)
            println(ind, "```")
        end
    end

    close(ind)
    close(tind)
end

function make_modules()

    mdpath = joinpath(docsrcdir, "modules.md")
    mod = open(mdpath, "w")
    seq_modules = filter!(s -> occursin(r"\.jl$", s), readdir(srcdir))

    for filename ∈ seq_modules
        filename ∈ exclude && continue
        path = joinpath(srcdir, filename)
        src = open(path, "r")
        indoc = false

        for l ∈ eachline(src, keep=false)
            n = lstrip(l)
            if startswith(n, "\"\"\"")
                indoc && break
                indoc = true
                path = joinpath(cloudpath, filename)
                name = splitext(filename)
                println(mod, "\n##   🔶  ", '[', name[1], "](", path, ")\n")
            else
                n = replace(n, "\\\\" => "\\")
                indoc && println(mod, n)
            end
        end
        close(src)
    end

    close(mod)
end

function addsig(srcfile, docfile)

    nextline(srcfile) =
        !eof(srcfile) ? (return readline(srcfile)) : (return nothing)
    n = nextline(srcfile)

    while true

        n === nothing && return
        while !startswith(n, "\"\"\"")
            println(docfile, n)
            n = nextline(srcfile)
            n === nothing && return
        end

        if startswith(n, "\"\"\"")
            println(docfile, n)
            n = nextline(srcfile)
            n === nothing && return
            while !startswith(n, "\"\"\"")
                println(docfile, n)
                n = nextline(srcfile)
                n === nothing && return
            end

            nn = nextline(srcfile)
            # Do you want signatures?
            #if !startswith(nn, "const Module")
            #    println(docfile, "\$(SIGNATURES)")
            #end
            println(docfile, n)
            n = nn
        end
    end
end

function build_all()

    @info("Uses the following external modules:")
    println(references)

    @info("The local modules included in IntegerTriangles.jl:")
    build_seq()

    srcdir = dirname(@__FILE__)
    srcfile = open(joinpath(srcdir, "_IS.jl"), "r")
    profile = open(joinpath(srcdir, "IntegerTriangles.jl"), "w")
    addsig(srcfile, profile)
    close(srcfile)
    close(profile)

    build_test()
    build_perf()
    make_index()
    make_modules()

    rm("_EXPORT.jl")
    rm("_INDEX.jl")
    rm("_IS.jl")
    rm("_SINDEX.jl")
    rm("_TEMP.jl")
end

build_all()

end # module


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

using .Meshes
using .GeoStatsBase
import .GeoStatsBase: preprocess, solvesingle

export TPS, PARAMS1

const PARAMS1 = [
  Params(2,  4,   0.01),
  Params(5,  10,  0.02),
  Params(10, 20,  0.03),
  Params(20, 40,  0.04),
  Params(50, 100, 0.05)
]

"""
    TPS(var₁=>param₁, var₂=>param₂, ...)

Turing pattern simulation.

## Parameters

* `params` - basic parameters (default to `PARAMS1`)
* `blur`   - blur algorithm (default to `BoxBlur`)
* `edge`   - edge condition (default to `Clamp`)
* `iter`   - number of iterations (default to `100`)

### References

Turing 1952. *The chemical basis of morphogenesis.*
"""
@simsolver TPS begin
  @param params = PARAMS1
  @param blur = BoxBlur
  @param edge = Clamp
  @param iter = 100
end

function preprocess(problem::SimulationProblem, solver::TPS)
  # retrieve domain size
  sz = size(domain(problem))

  # result of preprocessing
  preproc = Dict{Symbol,NamedTuple}()

  for covars in covariables(problem, solver)
    for var in covars.names
      # get user parameters
      varparams = covars.params[(var,)]

      # determine simulation parameters
      params = varparams.params
      blur = varparams.blur(sz)
      edge = varparams.edge()
      iter = varparams.iter

      # construct patterns from parameters
      patterns = [SimplePattern(param, sz) for param in params]

      preproc[var] = (patterns=patterns,blur=blur,edge=edge,iter=iter)
    end
  end

  preproc
end

function solvesingle(problem::SimulationProblem, covars::NamedTuple, ::TPS, preproc)
  # retrieve domain size
  sz = size(domain(problem))

  mactypeof = Dict(name(v) => mactype(v) for v in variables(problem))

  varreal = map(covars.names) do var
    # unpack preprocessed parameters
    patterns, blur, edge, iter = preproc[var]

    # determine value type
    V = mactypeof[var]

    # perform simulation
    sim = Sim(rand(V, sz), patterns, edge, blur)
    for i in 1:iter
      step!(sim)
    end
    real = scale01(sim.fluid)

    # flatten result
    var => vec(real)
  end

  Dict(varreal)
end


text:
module TestHistoricalWords

import Wordle: play
import Test: @test

# https://screenrant.com/wordle-answers-updated-word-puzzle-guide/

# 01/21 - #216 - PRICK
# 01/20 - #215 - ROBOT
# 01/19 - #214 - POINT
# 01/18 - #213 - PROXY
# 01/17 - #212 - SHIRE
# 01/16 - #211 - SOLAR
# 01/15 - #210 - PANIC
# 01/14 - #209 - TANGY
# 01/13 - #208 - ABBEY
# 01/12 - #207 - FAVOR
# 01/11 - #206 - DRINK
# 01/10 - #205 - QUERY
# 01/09 - #204 - GORGE
# 01/08 - #203 - CRANK
# 01/07 - #202 - SLUMP
# 01/06 - #201 - BANAL
# 01/05 - #200 - TIGER
# 01/04 - #199 - SIEGE
# 01/03 - #198 - TRUSS
# 01/02 - #197 - BOOST
# 01/01 - #196 - REBUS

words = (
    "prick",
    "robot",
    "point",
    "proxy",
    "shire",
    "solar",
    "panic",
    "tangy",
    "abbey",
    "favor",
    "drink",
    "query",
    "gorge",
    "crank",
    "slump",
    "banal",
    "tiger",
    "siege",
    "truss",
    "boost",
    "rebus",
)

for word in words
    @test play(word) > 0
end

end # module


text:
function gradient(O::Operation, vars::AbstractVector{Operation}; simplify = true)
    out = [expand_derivatives(Differential(v)(O)) for v in vars]
    simplify ? simplify_constants.(out) : out
end

function jacobian(ops::AbstractVector{Operation}, vars::AbstractVector{Operation}; simplify = true)
    out = [expand_derivatives(Differential(v)(O)) for O in ops, v in vars]
    simplify ? simplify_constants.(out) : out
end

function hessian(O::Operation, vars::AbstractVector{Operation}; simplify = true)
    out = [expand_derivatives(Differential(v2)(Differential(v1)(O))) for v1 in vars, v2 in vars]
    simplify ? simplify_constants.(out) : out
end

function simplified_expr(O::Operation)
  if O isa Constant
    return O.value
  elseif isa(O.op, Differential)
    return :(derivative($(simplified_expr(O.args[1])),$(simplified_expr(O.op.x))))
  elseif isa(O.op, Variable)
    isempty(O.args) && return O.op.name
    return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
  end
  return Expr(:call, Symbol(O.op), simplified_expr.(O.args)...)
end

function simplified_expr(c::Constant)
    c.value
end

function simplified_expr(eq::Equation)
    Expr(:(=), simplified_expr(eq.lhs), simplified_expr(eq.rhs))
end

macro I(ex)
    name = :ICompile
    ret = return quote
        macro $(esc(name))()
            esc($ex)
        end
    end
end


text:
using LogNormals
using SimpleTraits
using MissingStrategies

macro handlemissings_stub2(
    # omit the forwarding method
    fun,pos_missing=1, pos_strategy=pos_missing +1, type_missing=Any,
    defaultstrategy=nothing, argname_strategy=:ms, suffix = "_hm",
    )
    #gens = (mgen.forwarder, mgen.missingstrategy_notsuperofeltype)
    gens = (mgen.missingstrategy_notsuperofeltype,)
    #gens = (mgen.forwarder, )
    #gens = ()
    esc(handlemissings(fun, pos_missing, pos_strategy, type_missing,
        defaultstrategy, gens, argname_strategy, suffix))
 end
# #tmp2 = @macroexpand
# macro handlemissings_stub2(
#     fun,pos_missing=1, pos_strategy=pos_missing +1, type_missing=Any,
#     defaultstrategy=nothing,suffix = "_hm",
#     )
#     #gens = (mgen.forwarder, mgen.missingstrategy_notsuperofeltype)
#     gens = (mgen.missingstrategy_notsuperofeltype,)
#     esc(handlemissings(fun, pos_missing, pos_strategy, type_missing,
#         defaultstrategy, gens, suffix))
#   end
# @handlemissings_stub2(
#     autocor(x::AbstractVector{<:Real}; demean::Bool=true) = 0,
#     1,2,AbstractVector{<:Union{Missing,Real}}, nothing,
# )
@handlemissings_stub2(
    autocor(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) = 0,
    1,2,AbstractVector{<:Union{Missing,Real}}, nothing,
)
#autocor(x::AbstractVector{<:Union{Missing, Real}}, ms::MissingStrategies.MissingStrategy; kwargs...) = 1

# @handlemissings_stub(
#     autocor(x::AbstractVector{<:Real}; demean::Bool=true) = 0,
#     1,3,AbstractVector{<:Union{Missing,Real}}, PassMissing(),
# )

# tmp3 = @macroexpand @traitfn function autocor_hmnolag(ms::MissingStrategies.MissingStrategy, x::::IsEltypeSuperOfMissing;
#     demean::Bool=true
#     )
#     "Super of missing defined in LogNormals"
#     # lags not given -> forward to method with standard lag size
#     #autocor_hmlag(ms, x, StatsBase.default_autolags(size(x,1)); demean)
# end
# @traitfn function autocor_hm(ms::MissingStrategy, x::::IsEltypeSuperOfMissing;
#     demean::Bool=true
#     )
#     "Super of missing defined in LogNormals"
#     # lags not given -> forward to method with standard lag size
#     #autocor_hmlag(ms, x, StatsBase.default_autolags(size(x,1)); demean)
# end
@traitfn function autocor_hm(ms::MissingStrategy, x::::IsEltypeSuperOfMissing,
    lags::AbstractVector{<:Integer};
    demean::Bool=true
    )
    "Super of missing defined in LogNormals"
    # lags not given -> forward to method with standard lag size
    #autocor_hmlag(ms, x, StatsBase.default_autolags(size(x,1)); demean)
end

# module Test1
#     using SimpleTraits, MissingStrategies
#     @traitfn function autocor_hmnolag(ms::MissingStrategies.MissingStrategy, x::::IsEltypeSuperOfMissing;
#         demean::Bool=true
#         )
#         "Super of missing defined in Test1"
#         # lags not given -> forward to method with standard lag size
#         #autocor_hmlag(ms, x, StatsBase.default_autolags(size(x,1)); demean)
#     end
# end

# MissingStrategies.@testdefinetrait_esc_opt()
# @traitfn function f1_esc_opt(x::::!(SimpleTraits.BaseTraits.IsBits); demean::Bool=true)
#     "f1_esc_opt !IsBist from LogNormals"
# end

using Missings
m = allowmissing(rand(Float64,(3,3)))
convert.(Float64,m) # convert preserves dimensions


text:
"""
   Traits for the NEID spectrograph.
   https://neid.psu.edu/
Author: Eric Ford and collaborators
Created: August 2020
Updated: March 2021 for DRP v0.7
         May 2021 for DRP v1.0
         November 2021 for DRP v1.1
"""

""" Delegates loading of code specifying types essential to the package.  """

import RvSpectMLBase: min_order, max_order, min_pixel_in_order, max_pixel_in_order, min_pixel, max_pixel

#= Pre-ship
min_order(::NEID2D) = 1
max_order(::NEID2D) = 90
=#
# At KPNO
#min_order(::NEID2D) = 1
min_order(::NEID2D) = 4  # DRP 0.7
#max_order(::NEID2D) = 126 # DRS 0.5
#max_order(::NEID2D) = 122 # DRS 0.6
#max_order(::NEID2D) = 118 # DRS 0.7
max_order(::NEID2D) = 122 # DRS 1.0

min_pixel_in_order(inst::NEID2D) = 1
max_pixel_in_order(inst::NEID2D) = 9216

min_pixel(::NEID1D) = 1
max_pixel(::NEID1D) = (max_order(NEID1D())-min_order(NEID1D())+1)*9216 # TODO: Update once know size of NEID's 1d extracted spectra

import RvSpectMLBase: bad_col_ranges
#bad_col_ranges(inst::NEID2D, ord::Int) = [439:449, 1934:1943, 6714:6714] # DRP v0.6
function bad_col_ranges(inst::NEID2D, ord::Int)   # DRP v0.7
    if ord == 1      return [1:1712, 1938:1938]
    elseif ord == 2  return [1:866, 1938:1938]
    elseif ord == 3  return [1:152, 438:450, 1938:1938]
    elseif 4 <= ord <= 44   return [438:450, 1938:1938]
    elseif 45 <= ord <= 46   return [438:450, 1938:1938, 6714:6714]
    elseif 47 <= ord <= 122   return [438:450, 1932:1946, 6714:6714]
    end
end

import RvSpectMLBase: orders_to_use_default, min_col_default, max_col_default

#= Pre-ship
#orders_to_use_default(inst::NEID2D) = 1:52   # Avoiding redder orders due to tellurics
orders_to_use_default(inst::NEID2D) = 1:71   # Avoiding 72 because of NaNs in solar data
min_col_default(::NEID2D, ord::Integer) = 451              # Avoiding smaller columns due to NaNs
#min_col_default(::NEID2D) = 2000              # Avoiding smaller columns due to lower flux and distortions
max_col_default(::NEID2D, ord::Integer) = 9216 - (min_col_default(NEID2D(),ord)-1)   # Avoiding larger columns for symmetry
=#
# At KPNO
#orders_to_use_default(::NEID2D) = 55:113 # everything plaussibly usable as of DRS 0.5
#orders_to_use_default(::NEID2D) = 56:111 # everything plaussibly usable as of DRS 0.6
#orders_to_use_default(::NEID2D) = 60:90  # relatively safe
orders_to_use_default(::NEID2D) = 56:98   # avoid worst of tellurics
#=
function min_col_default(::NEID2D, ord::Integer)
#    return 1445 +500 # DRS 0.6, avoiding NaN in cols 1435-1444
    if ord == 55
        return 787
    else
        return 500
    end
end
=#
#max_col_default(::NEID2D, ord::Integer) = 8429  # DRS 0.5
#max_col_default(::NEID2D, ord::Integer) = 6214  # DRS 0.6, avoiding NaN in col 6215
#max_col_default(::NEID2D, ord::Integer) = 9215  # DRS 0.6, 0.7
min_col_default(::NEID2D, ord::Integer) = 1500  # Avoid where continuum normalization effected by edges/scattered light
max_col_default(::NEID2D, ord::Integer) = 8000  # Avoid where continuum normalization effected by edges/scattered light

import RvSpectMLBase: get_pixel_range
function get_pixel_range(inst::NEID2D, ord::Integer)
    minc = max(min_col_default(inst, ord)) #, min_col_excalibur(inst,order), min_col_nonnan(inst,order))
    maxc = min(max_col_default(inst, ord)) #, max_col_excalibur(inst,order), max_col_nonnan(inst,order))
    return minc:maxc
end

import RvSpectMLBase: metadata_symbols_default, metadata_strings_default
#= Pre-ship headers
#metadata_symbols_default(::AnyNEID) = Symbol[:bjd, :target, :ssbz]
#metadata_strings_default(::AnyNEID) = String["OBSJD", "SKY-OBJ", "SSBZ000"]
=#
# Headers at KPNO
metadata_symbols_default(::AnyNEID) = Symbol[:bjd, :target, :exptime, :airmass, :ssbz,
                            :wavecal, :driftfun, :drpversion, :drpextsnr,
                            :envotmp, :envohum, :envwinds, :envwindd ]
metadata_strings_default(::AnyNEID) = String["OBSJD", "OBJECT", "EXPTIME", "AIRMASS",  "SSBZ100",
                            "WAVECAL", "DRIFTFUN", "E_VER", "EXTSNR",
                            "ENVOTMP", "ENVOHUM", "ENVWINDS", "ENVWINDD" ]

import RvSpectMLBase: default_ccf_mask_v_width
default_ccf_mask_v_width(::AnyNEID) = 620.953

import RvSpectMLBase: get_inst_module
get_inst_module(::AnyNEID) = NEID

import RvSpectMLBase: get_λ_range
function get_λ_range(data::CLT) where { T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2},
                                       IT<:AnyNEID, CLT<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT} }
   (λmin, λmax) = NaNMath.extrema(data.λ)
   return (min=λmin, max=λmax)
end

#=
# Pre-ship
default_λmin = 3950.0  # Based on HD solar data from PSU, should generalize
default_λmax = 9500.0  #
=#
# At KPNO
#default_λmin = 5040.04045429369  # Based on one FITS file, should generalize
#default_λmax = 9810.539033567005 # Based on one FITS file, should generalize
# default_λmin = 5145.0 # DRP v0.6
# default_λmax = 9822.0 # DRP v0.6
default_λmin = 3571.0 # DRP v0.7 rounded down
default_λmax = 11048.0 # DRP v0.7

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using MathOptInterface
const MOI = MathOptInterface

export DSOSCone, SDSOSCone, SOSCone

export CoDSOSCone, CoSDSOSCone, CoSOSCone
export SOSMatrixCone
export getslack, certificate_monomials, lagrangian_multipliers

struct DSOSCone <: PolyJuMP.PolynomialSet end
struct CoDSOSCone <: PolyJuMP.PolynomialSet end
_varconetype(::DSOSCone) = DSOSPoly
_nococone(::CoDSOSCone) = DSOSCone()

struct SDSOSCone <: PolyJuMP.PolynomialSet end
struct CoSDSOSCone <: PolyJuMP.PolynomialSet end
_varconetype(::SDSOSCone) = SDSOSPoly
_nococone(::CoSDSOSCone) = SDSOSCone()

struct SOSCone <: PolyJuMP.PolynomialSet end
struct CoSOSCone <: PolyJuMP.PolynomialSet end
_varconetype(::SOSCone) = SOSPoly
_nococone(::CoSOSCone) = SOSCone()

struct SOSMatrixCone <: PolyJuMP.PolynomialSet end

const SOSLikeCones = Union{DSOSCone, SDSOSCone, SOSCone}
const CoSOSLikeCones = Union{CoDSOSCone, CoSDSOSCone, CoSOSCone}
const SOSSubCones = Union{CoSOSLikeCones, SOSLikeCones}

mutable struct SOSConstraint{MT <: AbstractMonomial, MVT <: AbstractVector{MT}, JS<:JuMP.AbstractJuMPScalar, F<:MOI.AbstractVectorFunction} <: PolyJuMP.ConstraintDelegate
    # JS is AffExpr for CoSOS and is Variable for SOS
    slack::MatPolynomial{JS, MT, MVT}
    zero_constraint::PolyJuMP.ZeroConstraint{MT, MVT, F}
    lagrangian_multipliers::Vector{MatPolynomial{JS, MT, MVT}}
end
function SOSConstraint(slack::MatPolynomial{JS, MT, MVT},
                       zero_constraint::PolyJuMP.ZeroConstraint{MT, MVT, F}) where {MT <: AbstractMonomial, MVT <: AbstractVector{MT}, JS<:JuMP.AbstractJuMPScalar, F<:MOI.AbstractVectorFunction}
    return SOSConstraint(slack, zero_constraint, MatPolynomial{JS, MT, MVT}[])
end

certificate_monomials(c::PolyJuMP.PolyConstraintRef) = certificate_monomials(PolyJuMP.getdelegate(c))
certificate_monomials(c::SOSConstraint) = c.slack.x
lagrangian_multipliers(c::PolyJuMP.PolyConstraintRef) = lagrangian_multipliers(PolyJuMP.getdelegate(c))
lagrangian_multipliers(c::SOSConstraint) = c.lagrangian_multipliers

JuMP.dual(c::SOSConstraint) = JuMP.dual(c.zero_constraint)

PolyJuMP.getslack(c::SOSConstraint) = JuMP.value(c.slack)

function PolyJuMP.addpolyconstraint!(m::JuMP.Model, P::Matrix{PT}, ::SOSMatrixCone, domain::AbstractBasicSemialgebraicSet, basis) where PT <: APL
    n = LinearAlgebra.checksquare(P)
    if !issymmetric(P)
        throw(ArgumentError("The polynomial matrix constrained to be SOS must be symmetric"))
    end
    y = [similarvariable(PT, gensym()) for i in 1:n]
    p = dot(y, P * y)
    PolyJuMP.addpolyconstraint!(m, p, SOSCone(), domain, basis)
end

function _createslack(m, x, set::SOSLikeCones)
    createpoly(m, _varconetype(set)(x), false, false)
end
function _matposynomial(m, x)
    p = _matpolynomial(m, x, false, false)
    for q in p.Q
        JuMP.set_lower_bound(q, 0)
    end
    p
end
function _createslack(m, x, set::CoSOSLikeCones)
    _matplus(_createslack(m, x, _nococone(set)), _matposynomial(m, x))
end

function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::AbstractAlgebraicSet, basis)
    r = rem(p, ideal(domain))
    X = getmonomialsforcertificate(monomials(r))
    slack = _createslack(m, X, set)
    q = r - slack
    zero_constraint = PolyJuMP.addpolyconstraint!(m, q, ZeroPoly(), domain, basis)
    SOSConstraint(slack, zero_constraint)
end

function lagrangian_multiplier(model::JuMP.Model, p, set::SOSSubCones, q, mindegree::Integer, maxdegree::Integer)
    mindegree_q, maxdegree_q = extdegree(q)
    # extdegree's that s^2 should have so that s^2 * p has degrees between mindegree and maxdegree
    mindegree_s2 = mindegree - mindegree_q
    maxdegree_s2 = maxdegree - maxdegree_q
    # extdegree's for s
    mindegree_s = max(0, div(mindegree_s2, 2))
    # If maxdegree_s2 is odd, div(maxdegree_s2,2) would make s^2 have degree up to maxdegree_s2-1
    # for this reason, we take div(maxdegree_s2+1,2) so that s^2 have degree up to maxdegree_s2+1
    maxdegree_s = div(maxdegree_s2 + 1, 2)
    # FIXME handle the case where `p`, `q_i`, ...  do not have the same variables
    # so instead of `variable(p)` we would have the union of them all
    @assert variables(q) ⊆ variables(p)
    return createpoly(model, _varconetype(set)(monomials(variables(p), mindegree_s:maxdegree_s)), false, false)
end

function PolyJuMP.addpolyconstraint!(m::JuMP.Model, p, set::SOSSubCones, domain::BasicSemialgebraicSet, basis;
                                     mindegree=MultivariatePolynomials.mindegree(p),
                                     maxdegree=MultivariatePolynomials.maxdegree(p))
    λ = lagrangian_multiplier.(Ref(m), p, Ref(set), domain.p, mindegree, maxdegree)
    p -= dot(λ, domain.p)
    constraint = PolyJuMP.addpolyconstraint!(m, p, set, domain.V, basis)
    constraint.lagrangian_multipliers = λ
    return constraint
end


text:
module CollabUtils2FakeMiniHeart
# The code for fake mini hearts is in miniHeart.jl.
# This file is only here to show up in the custom entity catalog.
end


text:
export toOrthoNormal, flatten, activation, segmentplot, rms, modelfit

function toOrthoNormal(Ti)
    local T = deepcopy(Ti)
    U_,S_,V_ = svd(T[1:3,1:3])
    local R = U_*diagm(0=>[1,1,sign(det(U_*V_'))])*V_'
    T[1:3,1:3] = R
    return T
end

const ⟂ = toOrthoNormal

function getD(D,T)
    if D == 3
        return sparse(toeplitz([-1; zeros(T-4)],[-1 3 -3 1 zeros(1,T-4)]))
    elseif D == 2
        return sparse(toeplitz([1; zeros(T-3)],[1 -2 1 zeros(1,T-3)]))
    elseif D == 1
        return sparse(toeplitz([-1; zeros(T-2)],[-1 1 zeros(1,T-2)]))
    end
    error("Can not handle your choice of D: $D")
end

function matrices(model::Union{SimpleLTVModel, KalmanModel}, d::AbstractIdData)
    T = length(d)
    T -= 1
    n = nstates(d)
    m = ninputs(d)
    x = oftype(Matrix, state(d))
    u = oftype(Matrix, input(d))
    FT = eltype(x)
    A = spzeros(FT, T*n, n^2+n*m)
    Is = sparse(FT(1.0)*I,n,n)
    for i = 1:T
        ii = (i-1)*n+1
        ii2 = ii+n-1
        A[ii:ii2,1:n^2] = kron(Is,x[:,i]')
        A[ii:ii2,n^2+1:end] = kron(Is,u[:,i]')
    end
    vec(x[:,2:end]),A
end

function matrices(m::LTVAutoRegressive, d::AbstractIdData)
    y = output(d)
    getARregressor(vec(y),m.na)
end

flatten(A) = reshape(A,prod(size(A)[1:2]),size(A,3))'
flatten(model::LTVStateSpaceModel) = [flatten(model.At) flatten(model.Bt)]
flatten(model::LTVAutoRegressive) = copy(model.θ')
decayfun(iters, reduction) = reduction^(1/iters)

function ABfromk(k,n,m,T)
    At = reshape(k[:,1:n^2]',n,n,T)
    At = permutedims(At, [2,1,3])
    Bt = reshape(k[:,n^2+1:end]',m,n,T)
    Bt = permutedims(Bt, [2,1,3])
    At,Bt
end

@views model2statevec(model,t) = [model.At[:,:,t]  model.Bt[:,:,t]]' |> vec

function model2statevec(model) # dispatch happens in flatten
    k = flatten(model)
    if model.extended
        k = k[1:end-1,:]
    end
    k
end

function statevec2model(::Type{<:LTVStateSpaceModel}, k,n,m,extend)
    At,Bt = ABfromk(k,n,m,size(k,1))
    SimpleLTVModel(At,Bt,extend)
end

function statevec2model(::Type{<:LTVAutoRegressive}, k,n,m,extend)
    LTVAutoRegressive(k',nothing,extend,0)
end

"""
    At,Bt = segments2full(parameters,breakpoints,n,m,T)

Takes a parameter vecor and breakpoints and returns At and Bt with length T
"""
function segments2full(parameters,breakpoints,n,m,T)
    At,Bt = zeros(n,n,T), zeros(n,m,T)
    i = 1
    for t = 1:T
        t ∈ breakpoints && (i+=1)
        At[:,:,t] = reshape(parameters[i][1:n^2],n,n)'
        Bt[:,:,t] = reshape(parameters[i][n^2+1:end],m,n)'
    end
    At,Bt
end

activation(model::LTVStateSpaceModel; kwargs...) = activation(model.At,model.Bt; kwargs...)
function activation(At,Bt; normalize=false)
    diffparams = (diff([flatten(At) flatten(Bt)],dims=1)).^2
    if normalize
        diffparams .-= minimum(diffparams,1)
        diffparams ./= maximum(diffparams,1)
    end
    activation = sqrt.(sum(diffparams,dims=2)[:])
    activation
end

segmentplot(model::AbstractModel, state;  kwargs...) = segmentplot(activation(model),state;kwargs...)
function segmentplot(act, state; filterlength=10, doplot=false, kwargs...)
    plot(act, lab="Activation")
    ds = diff(state)
    ma = mean(act)
    ds[ds.==0] = NaN
    segments = findpeaks(act; filterlength=filterlength, doplot=doplot, kwargs...)
    doplot || scatter!(segments,3ma*ones(segments), lab="Automatic segments", m=(10,:xcross))
    scatter!(2ma*ds, lab="Manual segments", xlabel="Time index", m=(10,:cross))
    segments
end

"""
    A,B,x,u,n,m,N = testdata(T=10000, σ_state_drift=0.001, σ_param_drift=0.001, σ_control=1)

Create an LTVModel with Brownian-walk A[t] and B[t]
"""
function testdata(;T=10000, σ_state_drift=0.001, σ_param_drift=0.001, seed=1, σ_control=1)
    Random.seed!(seed)
    n           = 3
    m           = 2
    N           = n*(n+m)
    A           = zeros(n,n,T)
    B           = zeros(n,m,T)
    x           = zeros(n,T)
    u           = σ_control*randn(m,T)
    U,S,V       = ⟂(randn(n,n)), diagm(0=>0.4rand(n)), ⟂(randn(n,n))
    A[:,:,1]    = U*S*V'
    B[:,:,1]    = 0.5randn(n,m)
    x[:,1]      = 0.1randn(n)

    for t = 1:T-1
        x[:,t+1]   = A[:,:,t]*x[:,t] + B[:,:,t]*u[:,t] + σ_state_drift*randn(n)
        A[:,:,t+1] = A[:,:,t] + σ_param_drift*randn(n,n)
        B[:,:,t+1] = B[:,:,t] + σ_param_drift*randn(n,m)
    end
    d = iddata(x,u,x)
    A,B,d,n,m,N
end

"""
    x,xm,u,n,m = testdata(T)

Create an LTVModel that changes from `A = [0.95 0.1; 0 0.95]` to `A = [0.5 0.05; 0 0.5]`
at T÷2
"""
function testdata(T_)
    Random.seed!(1)

    n,m      = 2,1
    At_      = [0.95 0.1; 0 0.95]
    Bt_      = reshape([0.2; 1],2,1)
    u        = randn(1,T_)
    x        = zeros(n,T_)
    for t = 1:T_-1
        if t == T_÷2
            At_ = [0.5 0.05; 0 0.5]
        end
        x[:,t+1] = At_*x[:,t] + Bt_*u[:,t] + 0.2randn(n)
    end
    xm = x + 0.2randn(size(x));
    x,xm,u,n,m
end

# Optimizers ===================================================================
#
#
# struct RMSpropOptimizer{VecType}
#     α::Float64
#     rmspropfactor::Float64
#     momentum::Float64
#     Θ::VecType
#     dΘs::VecType
#     dΘs2::VecType
# end
#
# function RMSpropOptimizer(Θ, α, rmspropfactor=0.8, momentum=0.1)
#     RMSpropOptimizer(α, rmspropfactor, momentum, Θ, ones(size(Θ)), zeros(size(Θ)))
# end
#
# function (opt::RMSpropOptimizer)(dΘ)
#     opt.dΘs .= opt.rmspropfactor.*opt.dΘs .+ (1-opt.rmspropfactor).*dΘ.^2
#     ΔΘ = -opt.α * dΘ
#     ΔΘ ./= sqrt.(opt.dΘs.+1e-10) # RMSprop
#
#     # ΔΘ = dΘ./sqrt(dΘs+1e-10).*sqrt(dΘs2) # ADAdelta
#     # dΘs2 .= 0.9dΘs2 + 0.1ΔΘ.^2 # ADAdelta
#     opt.dΘs2 .= opt.momentum.*opt.dΘs2 .+ ΔΘ # Momentum + RMSProp
#     opt.Θ .+= opt.dΘs2
# end
#
#
# mutable struct ADAMOptimizer{T, VecType <: AbstractArray}
#     Θ::VecType
#     α::T
#     β1::T
#     β2::T
#     ɛ::T
#     m::VecType
#     v::VecType
#     expdecay::Float64
# end
#
# ADAMOptimizer(Θ::VecType; α::T = 0.005,  β1::T = 0.9, β2::T = 0.999, ɛ::T = 1e-8, m=zeros(size(Θ)), v=zeros(size(Θ)), expdecay=0) where {T,VecType <: AbstractArray} = ADAMOptimizer{T,VecType}(Θ, α,  β1, β2, ɛ, m, v, expdecay)
#
# """
#     (a::ADAMOptimizer{T,VecType})(g::VecType, t::Integer)
# Applies the gradient `g` to the parameters `a.Θ` (mutating) at iteration `t`
# ADAM GD update http://sebastianruder.com/optimizing-gradient-descent/index.html#adam
# """
# function (a::ADAMOptimizer)(g, t::Integer)
#     mul = (1-a.β1)
#     mul2 = (1-a.β2)
#     div  = 1/(1 - a.β1 ^ t)
#     div2 = 1/(1 - a.β2 ^ t)
#     α,β1,β2,ɛ,m,v,Θ,γ = a.α,a.β1,a.β2,a.ɛ,a.m,a.v,a.Θ,a.expdecay
#     γ *= t > 5
#     Base.Threads.@threads for i = eachindex(g)
#         @inbounds m[i] = β1 * m[i] + mul * (g[i] + γ*Θ[i])
#         @inbounds v[i] = β2 * v[i] + mul2 * g[i]^2
#         @inbounds Θ[i] -= α * m[i] * div / (sqrt(v[i] * div2) + ɛ)
#     end
#     Θ
# end
# function (a::ADAMOptimizer)(g, t::Integer)
#     a.m .= a.β1 .* a.m .+ (1-a.β1) .* g
#     a.v .= a.β2 .* a.v .+ (1-a.β2) .* g.^2
#     div  = 1/(1 - a.β1 ^ t)
#     div2 = 1/(1 - a.β2 ^ t)
#     @. a.Θ  -= a.α * a.m * div / (sqrt(a.v * div2) + a.ɛ)
# end


text:
import JLD, FileIO
import JSON
import PyPlot
using LightGraphs, MetaGraphs

include("calc_xq.jl")
include("make_roadnet_figs.jl")
include("utilities.jl")
include("plot_rwd_dists.jl")

# make experiment dict
experiment_name = "mturk_fast"
f_loc = "logs"
include("experiment_utilities.jl")

i_loc = "imgs"
myval = []

run = Dict(:figs=>true,:xq=>true,:json=>true)

if run[:figs]
    info("Making Figures")
    make_figs(exp_dict,f_loc,i_loc)
end
if run[:xq]
    info("Calculating xQ")
    # run this second, because we create a new file based on solver to which reference is compared
    calculate_sq(exp_dict)
end

function get_exit_node(g::MetaGraphs.MetaGraph)
    e_nodes = g.gprops[:exit_nodes]
    if length(e_nodes) > 1
        error("something is wrong, should only be length 1 in this experiment")
    else
        return e_nodes[1]
    end
end

if run[:json]
    info("## Making .json file for experiment data")
    srand(5555) #for consistent results
    for expr in keys(exp_dict[:conditions])
        net_type = exp_dict[:name]
        inpts = exp_dict[:xQ][:inpts]

        fname =  net_type*"_$(expr)_"*make_label_from_keys(inpts)*"_solver"
        f = joinpath("logs",fname)

        info("processing $f")
        d = JLD.load(f*".jld")["problem_dict"]

        exp_json = Dict()

        for net in d
            # need to use a HACK because I accidentally placed 13 in the :exit_loc key, when the networks
            # actually have 8 as the exit in the :graph properties. the graph properties are used for
            # the mdp solution, so unless I want to re-create the data (which is expensive in time and money
            # I'm just going to write some hack code to get the exit node from the graph
            exit_loc = get_exit_node(net[2][:graph])
            dist_to_exit = length(a_star(net[2][:graph],net[2][:evader_start],exit_loc))
            dist_to_pursuer = length(a_star(net[2][:graph],net[2][:evader_start],net[2][:pursuer_start]))
            edge_to_node_ratio = ne(net[2][:graph])/nv(net[2][:graph])
            etnr = 2.5
            dte = 2.0
            dtp = 2.0
            if dist_to_exit < dte || dist_to_pursuer < dtp || edge_to_node_ratio > etnr
                # we don't want networks where exit or pursuer is too close at the beginning
                # or where the edge/node ratio is too high (i.e. too many edges)
                if dist_to_exit < dte
                    println("Exit too close (d=$dist_to_exit): skipping $(net[1])")
                end
                if dist_to_pursuer < dtp
                    println("Pursuer too close (d=$dist_to_pursuer): skipping $(net[1])")
                end
                if edge_to_node_ratio > etnr
                    println("Too many edges ($edge_to_node_ratio): skipping $(net[1])")
                end
                continue
            end

            n = exp_json[net[1]] = Dict()
            net_vals = net[2]

            n[:xQ] = net_vals[:solver_quality][:X_Q]
            n[:xP] = net_vals[:training_data][:X4]
            n[:dist_to_exit] = dist_to_exit
            n[:dist_to_pursuer] = dist_to_pursuer
            n[:edge_to_node_ratio] = edge_to_node_ratio

            # calculate whether experiment succeded or not, take reward distribution and draw at random
            # this will have a higher probability of drawing most repeated numbers, so it is like sampling
            # a parameterized distribution, with the limitaion that it will only draw an outcome that has
            # already been seen
            # If that value is greater than the mean expected reward of the trusted distribution (:X3_1) then the experiment is a success
            # otherwise it is a failure
            r_dist = net_vals[:training_data][:r_dist]
            myval = r_dist
            o = rand(r_dist)

            #  r_star_mean = net_vals[:training_data][:X3_1] # Old, WRONG way, leave hear for illustrative purposes
            r_star_mean = net_vals[:solver_quality][:R_star][:X3_1]
            #  if all(r_star_mean .> o) #old WRONG way
            #  if all(r_star_mean .> o) && all(o .< 0.0)
            if all(o .< 0.0)
                if n[:xQ] < 0.5 && n[:xP] > 0.5
                    # this is an unexpected outcome
                    plot_rwd_dists(net[1],n[:xQ],n[:xP],net_vals[:training_data][:r_dist],net_vals[:solver_quality][:R_star],o,fldr="figs/exp_data/"*net_type*"_$(expr)",fname="$(fname)_$(net[1])_bigfail",outcome=:fail)
                end
                if n[:xQ] > 1.15 && n[:xP] > 0.25
                    # this is an unexpected outcome
                    plot_rwd_dists(net[1],n[:xQ],n[:xP],net_vals[:training_data][:r_dist],net_vals[:solver_quality][:R_star],o,fldr="figs/exp_data/"*net_type*"_$(expr)",fname="$(fname)_$(net[1])_surprisefail",outcome=:fail)
                end
                outcome = "fail"
            else
                if n[:xQ] < 0.75 && n[:xP] < -0.1
                    # this is an unexpected outcome
                    plot_rwd_dists(net[1],n[:xQ],n[:xP],net_vals[:training_data][:r_dist],net_vals[:solver_quality][:R_star],o,fldr="figs/exp_data/"*net_type*"_$(expr)",fname="$(fname)_$(net[1])_surprisesuccess",outcome=:success)
                end
                outcome = "success"
            end
            n[:outcome] = outcome
            n[:num_nodes] = net_vals[:num_nodes]
            n[:trans_prob] = net_vals[:trans_prob]
            n[:rwd_draw] = o
            n[:trusted_mean] = r_star_mean
            outcome == "success" ? oc_bool = 1 : oc_bool = 0
            #  println("rand draw: $o")
            #  println("trusted mean: $r_star_mean")
            #  println("oucome: $outcome")
            #  readline()
            #  info(n[:outcome])

            n[:image_file] = joinpath(net_vals[:img_name])
            n[:mcts_depth] = net_vals[:mcts_depth]
        end

        json_fname = "logs/experiment_data_$(fname).json"
        info("Writing to $json_fname")
        open(json_fname,"w") do f
            dat = JSON.json(exp_json,4)
            write(f,dat)
        end
    end
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using JavaCall
DBtableTypeorString = Union{DBtableType, AbstractString}

function adjbfs(A::DBtableTypeorString, v0::AbstractString, numsteps::Number, Rtable::AbstractString, RtableTranspose::AbstractString; minDegree = 0::Number, maxDegree = 2^31 - 1::Number, ADegtable = ""::AbstractString, degColumn = ""::AbstractString, degInColQ = false::Bool)

#    (String Atable, String v0, int k, String Rtable, String RtableTranspose,
#    String ADegtable, String degColumn, boolean degInColQ, int minDegree, int maxDegree)

    if isa(A, DBtablePair)
        Aname = A.name1
    elseif isa(A, DBtable)
        Aname = A.name
    else
        Aname = A
    end

    jcall(A.DB.Graphulo, "AdjBFS", JString,
        (JString, JString, jint, JString, JString, JString, JString, jboolean, jint, jint,),
        Aname, v0, numsteps, Rtable, RtableTranspose, ADegtable, degColumn, degInColQ, minDegree, maxDegree)

end

# Rows are edges, columns are nodes
function edgebfs(E::DBtableTypeorString, v0::AbstractString, numsteps::Number, Rtable::AbstractString, RtableTranspose::AbstractString; minDegree = 0::Number, maxDegree = 2^31 - 1::Number, EDegtable = ""::AbstractString, degColumn = ""::AbstractString, degInColQ= false::Bool)

    #    (String Atable, String v0, int k, String Rtable, String RtableTranspose,
    #    String ADegtable, String degColumn, boolean degInColQ, int minDegree, int maxDegree)

    if isa(E, DBtablePair)
        Ename = E.name1
    elseif isa(E, DBtable)
        Ename = E.name
    else
        Ename = E
    end

    JAuthorizations = @jimport "org.apache.accumulo.core.security.Authorizations"
    emptyAuth = JAuthorizations((), )
    JIteratorSetting = @jimport "org.apache.accumulo.core.client.IteratorSetting"
    emptyIteratorSetting = JIteratorSetting((jint, JString, JString), 500, "c", "cname")
    # We want to pass a null IteratorSetting as an argument; JavaCall hates this
    # here's a workaround:
    emptyIteratorSetting.ptr = C_NULL
    JMutableLong = @jimport "org.apache.commons.lang.mutable.MutableLong"
    emptyMutableLong = JMutableLong((), )

    startPrefixes = ","
    endPrefixes = ","
    plusOp = emptyIteratorSetting
    EScanIteratorPriority = -1
    Eauthorizations = emptyAuth
    EDegauthorizations = emptyAuth
    newVisibility = "" # TODO use null or empty string as default?
    useNewTimestamp = true
    outputUnion = true
    numEntriesWritten = emptyMutableLong

    #=
    JString, JString, jint, JString, JString, JString, JString, JString, JString, jboolean, jint, jint, IteratorSetting, jint,
    Authorizations, Authorizations, JString, jboolean, jboolean, MutableLong

    Etable, v0, k, Rtable, RTtable, startPrefixes, endPrefixes, ETDegtable, degColumn, degInColQ, minDegree,
    maxDegree, plusOp, EScanIteratorPriority, Eauthorizations, EDegauthorizations, newVisibility, useNewTimestamp, outputUnion, numEntriesWritten
    =#

    return jcall(A.DB.Graphulo, "EdgeBFS", JString,
            (JString, JString, jint, JString, JString,
                JString, JString, JString, JString,
                jboolean, jint, jint, JIteratorSetting,
                jint, JAuthorizations, JAuthorizations,
                JString, jboolean, jboolean, JMutableLong,),
            Ename, v0, numsteps, Rtable, RtableTranspose,
            startPrefixes, endPrefixes, EDegtable, degColumn,
            degInColQ, minDegree, maxDegree, plusOp,
            EScanIteratorPriority, Eauthorizations, EDegauthorizations,
            newVisibility, useNewTimestamp, outputUnion, numEntriesWritten)

end

function singlebfs(A::DBtableType, v0::AbstractString, numsteps::Number, Rtable = ""::AbstractString, RtableTranspose = ""::AbstractString; minDegree = 0::Number, maxDegree = 2^31 - 1::Number, ADegtable = ""::AbstractString, degColumn = ""::AbstractString, degInColQ = false)

    #    (String Atable, String v0, int k, String Rtable, String RtableTranspose,
    #    String ADegtable, String degColumn, boolean degInColQ, int minDegree, int maxDegree)

    if isa(A, DBtablePair)
        Atable = A.name1
    else
        Atable = A.name
    end

    return jcall(A.DB.Graphulo, "AdjBFS", JString,
            (JString, JString, jint, JString, JString, JString, JString, jboolean, jint, jint,),
            Atable, v0, numsteps, Rtable, RtableTranspose, ADegtable, degColumn, degInColQ, minDegree, maxDegree)

end