text:
##################################################################################
# This file is part of ModelBaseEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

export SteadyStateEquation

"""
    struct SteadyStateEquation <: AbstractEquation

Data structure representing an individual steady state equation.
"""
struct SteadyStateEquation <: AbstractEquation
    type::Symbol
    vinds::Vector{Int64}
    vsyms::Vector{Symbol}
    expr::ExtExpr
    eval_resid::Function
    eval_RJ::Function
end

########################################################

"""
    struct SteadyStateVariable

Holds the steady state solution for one variable.
"""
struct SteadyStateVariable{DATA<:AbstractVector{Float64},MASK<:AbstractVector{Bool}}
    # The corresponding entry in m.allvars. Needed for its type
    name::ModelVariable
    # Its index in the m.allvars array
    index::Int
    # A view in the SteadyStateData.values location for this variable
    data::DATA
    # A view into the mask array
    mask::MASK
end

for sym in (:lin, :log, :neglog, :steady, :exog, :shock)
    issym = Symbol("is", sym)
    eval(quote
        @inline $(issym)(s::SteadyStateVariable) = $(issym)(s.name)
    end)
end

@inline transformation(v::SteadyStateVariable) = transformation(v.name)
@inline inverse_transformation(v::SteadyStateVariable) = inverse_transformation(v.name)
@inline transform(x, v::SteadyStateVariable) = transform(x, v.name)
@inline inverse_transform(x, v::SteadyStateVariable) = inverse_transform(x, v.name)

#############################################################################
# Access to .level and .slope

Base.getproperty(v::SteadyStateVariable, name::Symbol) = begin
    if name ∈ fieldnames(SteadyStateVariable)
        return getfield(v, name)
    end
    data = getfield(v, :data)
    # we store transformed data, must invert to give back to user
    return name == :level ? inverse_transform(data[1], v) :
           name == :slope ?
           (islog(v) || isneglog(v) ? exp(data[2]) : data[2]) :
           getfield(v, name)
end

Base.setproperty!(v::SteadyStateVariable, name::Symbol, val) = begin
    data = getfield(v, :data)
    # we store transformed data, must transform user input
    if name == :level
        data[1] = transform(val, v)
    elseif name == :slope
        data[2] = (islog(v) || isneglog(v) ? log(val) : val)
    else
        setfield!(v, name, val)  # this will error (immutable)
    end
end

# use [] to get a time series of values.
# the ref= value is the t at which it equals its level
function Base.getindex(v::SteadyStateVariable, t; ref = 1)
    if eltype(t) != eltype(ref)
        throw(ArgumentError("Must provide reference time of the same type as the time index"))
    end
    int_t = convert(Int, first(t) - ref):convert(Int, last(t) - ref)
    # we must inverse transform internal data before returning to user
    return inverse_transform(v.data[1] .+ v.data[2] .* int_t, v)
end

# pretty printing

"""
    alignment5(io::IO, v::SteadyStateVariable)

Return a 5-tuple with the number of characters for the name, and the alignment 2-tuples for level and slope.

"""
function alignment5(io::IO, v::SteadyStateVariable)
    name = sprint(print, string(v.name.name), context = io, sizehint = 0)
    lvl_a = Base.alignment(io, v.level)
    if isshock(v) || issteady(v)
        slp_a = (0, 0)
    else
        slp_a = Base.alignment(io, v.slope)
    end
    (length(name), lvl_a..., slp_a...)
end

function alignment5(io::IO, vars::AbstractVector{SteadyStateVariable})
    a = (0, 0, 0, 0, 0)
    for v in vars
        a = max.(a, alignment5(io, v))
    end
    return a
end

function show_aligned5(io::IO, v::SteadyStateVariable, a = alignment5(io, v);
    mask = trues(2), sep1 = " = ",
    sep2 = islog(v) || isneglog(v) ? " * " : " + ",
    sep3 = islog(v) || isneglog(v) ? "^t" : "*t")
    name = sprint(print, string(v.name.name), context = io, sizehint = 0)
    if mask[1]
        lvl_a = Base.alignment(io, v.level)
        lvl = sprint(show, v.level, context = io, sizehint = 0)
    else
        lvl_a = (0, 1)
        lvl = "?"
    end
    if mask[2]
        slp_a = Base.alignment(io, v.slope)
        slp = sprint(show, v.slope, context = io, sizehint = 0)
    else
        slp_a = (0, 1)
        slp = "?"
    end
    print(io, "  ", repeat(' ', a[1] - length(name)), name,
        sep1, repeat(' ', a[2] - lvl_a[1]), lvl, repeat(' ', a[3] - lvl_a[2]))
    if (!issteady(v) && !isshock(v)) && !(v.data[2] + 1.0 ≈ 1.0)
        print(io, sep2, repeat(' ', a[4] - slp_a[1]), slp, repeat(' ', a[5] - slp_a[2]), sep3)
    end
end

@inline Base.show(io::IO, v::SteadyStateVariable) = show_aligned5(io, v)

########################################################

export SteadyStateData

"""
    SteadyStateData

Data structure that holds information about the steady state solution of the Model.
"""
struct SteadyStateData
    "List of steady state variables."
    vars::Vector{SteadyStateVariable}
    "Steady state solution vector."
    values::Vector{Float64}
    "`mask[i] == true if and only if `values[i]` holds the steady state value."
    mask::BitArray{1}
    "Steady state equations derived from the dynamic system."
    equations::Vector{SteadyStateEquation}
    "Steady state equations explicitly added with @steadystate."
    constraints::Vector{SteadyStateEquation}
    # default constructor
    SteadyStateData() = new([], [], [], [], [])
end

Base.push!(ssd::SteadyStateData, vars...) =
    for v in vars
        push!(ssd, v)
    end
Base.push!(ssd::SteadyStateData, var::Symbol) = push!(ssd, convert(ModelSymbol, var))
function Base.push!(ssd::SteadyStateData, var::ModelSymbol)
    for v in getfield(ssd, :vars)
        if v.name == var
            return v
        end
    end
    push!(getfield(ssd, :values), isexog(var) || isshock(var) ? 0.0 : 0.1, 0.0)
    push!(getfield(ssd, :mask), isexog(var) || isshock(var), isexog(var) || isshock(var) || issteady(var))
    ind = length(getfield(ssd, :vars)) + 1
    v = SteadyStateVariable(var, ind, @view(getfield(ssd, :values)[2ind.+(-1:0)]), @view(getfield(ssd, :mask)[2ind.+(-1:0)]))
    push!(getfield(ssd, :vars), v)
    return v
end

export alleqns
"""
    alleqns(ssd::SteadyStateData)

Return a list of all steady state equations.

The list contains all equations derived from the dynamic system and all explicitly added steady state constraints.
"""
@inline alleqns(ssd::SteadyStateData) = vcat(ssd.equations, ssd.constraints)

export neqns
"""
    neqns(ssd::SteadyStateData)

Return the total number of equations in the steady state system, including the ones derived from the dynamic system and the
ones added explicitly as steady state constraints.
"""
@inline neqns(ssd::SteadyStateData) = length(ssd.equations) + length(ssd.constraints)

export geteqn
"""
    geteqn(i, ssd::SteadyStateData)

Return the i-th steady state equation. Index i is interpreted as in the output of `alleqns`.
Calling `geteqn(i, sdd)` has the same effect as `alleqn(ssd)[i]`, but it's more efficient.

### Example
```julia
# Iterate all equations like this:
for i = 1:neqns(ssd)
    eqn = geteqn(i, ssd)
    # do something awesome with `eqn` and `i`
end
```
"""
function geteqn(i::Integer, ssd::SteadyStateData)
    ci = i - length(ssd.equations)
    return ci > 0 ? ssd.constraints[ci] : ssd.equations[i]
end

Base.show(io::IO, ::MIME"text/plain", ssd::SteadyStateData) = show(io, ssd)
Base.show(io::IO, ssd::SteadyStateData) = begin
    if issssolved(ssd)
        println(io, "Steady state solved.")
    else
        println(io, "Steady state not solved.")
    end
    len = length(ssd.constraints)
    if len == 0
        println(io, "No additional constraints.")
    else
        println(io, len, " additional constraint", ifelse(len > 1, "s.", "."))
        for c in ssd.constraints
            println(io, "    ", c)
        end
    end
end

#########
# Implement access to steady state values using dot notation and index notation

function Base.propertynames(ssd::SteadyStateData, private::Bool = false)
    if private
        return ((v.name.name for v in ssd.vars)..., fieldnames(SteadyStateData)...,)
    else
        return ((v.name.name for v in ssd.vars)...,)
    end
end

function Base.getproperty(ssd::SteadyStateData, sym::Symbol)
    if sym ∈ fieldnames(SteadyStateData)
        return getfield(ssd, sym)
    else
        for v in ssd.vars
            if v.name == sym
                return v
            end
        end
        throw(ArgumentError("Unknown variable $sym."))
    end
end

Base.getindex(ssd::SteadyStateData, ind::Int) = ssd.vars[ind]
Base.getindex(ssd::SteadyStateData, sym::ModelSymbol) = getproperty(ssd, sym.name)
Base.getindex(ssd::SteadyStateData, sym::Symbol) = getproperty(ssd, sym)
Base.getindex(ssd::SteadyStateData, sym::AbstractString) = getproperty(ssd, Symbol(sym))

@inline ss_symbol(ssd::SteadyStateData, vi::Int) = Symbol("#", ssd.vars[(1+vi)÷2].name.name, "#", (vi % 2 == 1) ? :lvl : :slp, "#")

#########################
#

export printsstate

"""
    printsstate([io::IO,] ssd::SteadyStateData)

Display steady state solution.

Steady state solution is presented in a table, where the first column is
the name of the variable, the second and third columns are the corresponding
values of the level and the slope. If the value is not determined
(as per its `mask` value) then it is displayed as "*".
"""
function printsstate(io::IO, model::AbstractModel)
    io = IOContext(io, :compact => get(io, :compact, true))
    ssd = model.sstate
    println(io, "Steady State Solution:")
    a = max.(alignment5(io, ssd.vars), (0, 0, 3, 0, 3))
    for v in ssd.vars
        show_aligned5(io, v, a, mask = v.mask)
        println(io)
    end
end
printsstate(model::AbstractModel) = printsstate(Base.stdout, model)

###########################
# Make steady state equation from dynamic equation

# Idea:
#   in the steady state equation, we assume that the variable y_ss
#   follows a linear motion expressed as y_ss[t] = y_ss#lvl + t * y_ss#slp
#   where y_ss#lvl and y_ss#slp are two unknowns we solve for.
#
#   The dynamic equation has mentions of lags and leads. We replace those
#   with the above expression.
#
#   Since we have two parameters to determine, we need two steady state equations
#   from each dynamic equation. We get this by writing the dynamic equation at
#   two different values of `t` - 0 and another one we call `shift`.
#
#   Shift is a an option in the model object, which the user can set to any integer
#   other than 0. The default is 10.

"""
    SSEqnData

Internal structure used for evaluation of the residual of the steady state equation
derived from a dynamic equation.

"""
struct SSEqnData{M<:AbstractModel}
    "Whether or not to add model.shift"
    shift::Bool
    "Reference to the model object (needed for the current value of shift"
    model::Ref{M}
    "Information needed to compute the Jacobian matrix of the transformation between steady state and dynamic unknowns"
    JT::Vector
    "The dynamic equation instance"
    eqn::Equation
    SSEqnData(s, m::M, jt, e) where {M<:AbstractModel} = new{M}(s, Ref(m), jt, e)
    SSEqnData(s, m::Ref{M}, jt, e) where {M<:AbstractModel} = new{M}(s, m, jt, e)
end

function sseqn_resid_RJ(s::SSEqnData)
    buffer = Vector{Float64}(undef, length(s.JT))
    @inline lag(jt) = ifelse(s.shift, jt.tlag + s.model[].shift, jt.tlag)
    function to_dyn_pt(pt::AbstractVector{Float64})
        # This function applies the transformation from steady
        # state equation unknowns to dynamic equation unknowns
        fill!(buffer, 0.0)
        for (i, jt) in enumerate(s.JT)
            if length(jt.ssinds) == 1
                pti = pt[jt.ssinds[1]]
            else
                pti = pt[jt.ssinds[1]] + lag(jt) * pt[jt.ssinds[2]]
            end
            buffer[i] += pti
        end
        return buffer
    end
    function to_ssgrad(pt::AbstractVector{Float64}, jj::AbstractVector{Float64})
        # This function inverts the transformation. jj is the gradient of the
        # dynamic equation residual with respect to the dynamic equation unknowns.
        # Here we compute the Jacobian of the transformation and use it to compute
        ss = zeros(size(pt))
        for (i, jt) in enumerate(s.JT)
            if length(jt.ssinds) == 1
                pti = pt[jt.ssinds[1]]
                ss[jt.ssinds[1]] += jj[i]
            else
                local lag_jt = lag(jt)
                pti = pt[jt.ssinds[1]] + lag_jt * pt[jt.ssinds[2]]
                ss[jt.ssinds[1]] += jj[i]
                ss[jt.ssinds[2]] += jj[i] * lag_jt
            end
        end
        # NOTE regarding the above: The dynamic equation is F(x_t) = 0
        # Here we're solving F(u(l+t*s)) = 0
        # The derivative is dF/dl = F' * u' and dF/ds = F' * u' * t
        # F' is in jj[i]
        # u(x) = x, so u'(x) = 1
        return ss
    end
    function _resid(pt::AbstractVector{Float64})
        return s.eqn.eval_resid(to_dyn_pt(pt))
    end
    function _RJ(pt::AbstractVector{Float64})
        R, jj = s.eqn.eval_RJ(to_dyn_pt(pt))
        return R, to_ssgrad(pt, jj)
    end
    return _resid, _RJ
end

"""
    make_sseqn(model::AbstractModel, eqn::Equation; shift::Int64=0)

Create a steady state equation from the given dynamic equation for the given model.

Internal function, do not call directly.

"""
function make_sseqn(model::AbstractModel, eqn::Equation, shift::Bool)
    local allvars = model.allvars
    @inline tvalue(t) = ifelse(shift, t + model.shift, t)
    # ssind converts the dynamic index (t, v) into
    # the corresponding indexes of steady state unknowns.
    # Returned value is a list of length 0, 1, or 2.
    function ssind((ti, vi),)::Array{Int64,1}
        no_slope = isshock(allvars[vi]) || issteady(allvars[vi])
        # The level unknown has index 2*vi-1.
        # The slope unknown has index 2*vi. However:
        #  * :steady and :shock variables don't have slopes
        #  * :lin and :log variables the slope is in the equation
        #    only if the effective t-index is not 0.
        if no_slope || (!shift && ti == 0)
            return [2vi - 1]
        else
            return [2vi - 1, 2vi]
        end
    end
    local ss = model.sstate
    # The steady state indexes.
    vinds = unique(vcat(map(ssind, eqn.vinds)...))
    # The corresponding steady state symbols
    vsyms = Symbol[ss_symbol(ss, vi) for vi in vinds]
    # In the next loop we build the matrix JT which transforms
    # from the steady state values to the dynamic point values.
    JT = []
    for (i, (ti, vi)) in enumerate(eqn.vinds)
        val = (ssinds = indexin(ssind((ti, vi)), vinds), tlag = ti)
        push!(JT, val)
    end
    type = ifelse(shift == 0, :tzero, :tshift)
    let sseqndata = SSEqnData(shift, Ref(model), JT, eqn)
        return SteadyStateEquation(type, vinds, vsyms, eqn.expr, sseqn_resid_RJ(sseqndata)...)
    end
end

###########################
# Make steady state equation from user input

"""
    setss!(model::AbstractModel, expr::Expr; type::Symbol, modelmodule::Module)

Add a steady state equation to the model. Equations added by `setss!` are in
addition to the equations generated automatically from the dynamic system.

Internal function, do not call directly. Use [`@steadystate`](@ref) instead.

"""
function setss!(model::AbstractModel, expr::Expr; type::Symbol, modelmodule::Module = moduleof(model))

    if expr.head != :(=)
        error("Expected an equation, not $(expr.head)")
    end

    @assert type ∈ (:level, :slope) "Unknown steady state equation type $type. Expected either `level` or `slope`."

    local ss = sstate(model)

    local allvars = model.allvars

    ###############################################
    #     ssprocess(val)
    #
    # Process the given value to extract information about mentioned parameters and variables.
    # This function has the side effect of populating the vectors
    # `vinds`, `vsyms`, `val_params` and `source`
    #
    # Algorithm is recursive over the given expression. The bottom of the recursion is the
    # processing of a `Number`, a `Symbol`, (or a `LineNumberNode`).
    #
    # we will store indices
    local vinds = Int64[]
    local vsyms = Symbol[]
    # we will store parameters mentioned in `expr` here
    local val_params = Symbol[]
    local source = LineNumberNode[]
    # nothing to do with a number
    ssprocess(val::Number) = val
    # a symbol could be a variable (shock, auxvar), a parameter, or unknown.
    function ssprocess(val::Symbol)
        if val ∈ keys(model.parameters)
            # parameter - keep track that it's mentioned
            push!(val_params, val)
            return val
        end
        vind = indexin([val], allvars)[1]
        if vind !== nothing
            # it's a vriable of some sort: make a symbol and an index for the
            # corresponding steady state unknown
            vsym = Symbol("#", type, "#", val, "#")
            push!(vsyms, vsym)
            push!(vinds, type == :level ? 2vind - 1 : 2vind)
            if need_transform(allvars[vind])
                func = inverse_transformation(allvars[vind])
                return :($func($vsym))
            else
                return vsym
            end
        end
        # what to do with unknown symbols?
        error("unknown parameter $val")
    end
    # a sorce line information: store it and remove it from the expression
    ssprocess(val::LineNumberNode) = (push!(source, val); nothing)
    # process an expression recursively
    function ssprocess(val::Expr)
        if val.head == :(=)
            # we process the lhs and rhs seperately: there shouldn't be any equal signs
            error("unexpected equation.")
        end
        if val.head == :block
            # in a begin-end block, process each line and gather the results
            args = filter(x -> x !== nothing, map(ssprocess, val.args))
            if length(args) == 1
                # Only one thing left - no need for the begin-end anymore
                return args[1]
            else
                # reassemble the processed expressions back into a begin-end block
                return Expr(:block, args...)
            end
        elseif val.head == :call
            # in a function call, process each argument, but not the function name (args[1]) and reassemble the call
            args = filter(x -> x !== nothing, map(ssprocess, val.args[2:end]))
            return Expr(:call, val.args[1], args...)
        else
            # whatever this it, process each subexpression and reassemble it
            args = filter(x -> x !== nothing, map(ssprocess, val.args))
            return Expr(val.head, args...)
        end
    end
    # end of ssprocess() definition
    ###############################################
    #
    lhs, rhs = expr.args
    lhs = ssprocess(lhs)
    rhs = ssprocess(rhs)
    expr.args .= MacroTools.unblock.(expr.args)
    #
    nargs = length(vinds)
    # In case there's no source information, add a dummy one
    push!(source, LineNumberNode(0))
    # create the resid and RJ functions for the new equation
    # To do this, we use `makefuncs` from evaluation.jl
    resid, RJ = let mparams = parameters(model)
        # create a list of expressions that assign the values of model parameters to
        # variables of the same name
        param_assigments = Expr(:block)
        for p in unique(val_params)
            push!(param_assigments.args, :(local $(p) = $(mparams).$(p)))
        end
        residual = Expr(:block, source[1], :($(lhs) - $(rhs)))
        funcs_expr = makefuncs(residual, vsyms, param_assigments; mod = modelmodule)
        modelmodule.eval(funcs_expr)
    end
    # We have all the ingredients to create the instance of SteadyStateEquation
    for i = 1:2
        # remove blocks with line numbers from expr.args[i]
        a = expr.args[i]
        if Meta.isexpr(a, :block)
            args = filter(x -> !isa(x, LineNumberNode), a.args)
            if length(a.args) == 1
                expr.args[i] = args[1]
            end
        end
    end
    sscon = SteadyStateEquation(type, vinds, vsyms, expr, resid, RJ)
    if nargs == 1
        # The equation involves only one variable. See if there's already an equation
        # with just that variable and, if so, remove it.
        for (i, ssc) in enumerate(ss.constraints)
            if ssc.type == type && length(ssc.vinds) == 1 && ssc.vinds[1] == sscon.vinds[1]
                ss.constraints[i] = sscon
                return sscon
            end
        end
    end
    push!(ss.constraints, sscon)
    return sscon
end

export @steadystate

"""
    @steadystate model [type] equation

Add a steady state equation to the model.

The steady state system of the model is automatically derived from the dynamic
system. Use this macro to define additional equations for the steady state.
This is particularly useful in the case of a non-linear model that might have
multiple steady state, or the steady state might be difficult to solve for,
to help the steady state solver find the one you want to use.

  * `model` is the model instance you want to update
  * `type` (optional) is the type of constraint you want to add. This can be `level`
  or `slope`. If missing, the default is `level`
  * `equation` is the expression defining the steady state constraint. In the
  equation, use variables and shocks from the model, but without any t-references.

"""
macro steadystate(model, type::Symbol, equation::Expr)
    thismodule = @__MODULE__
    modelmodule = __module__
    return esc(:($(thismodule).setss!($(model), $(Meta.quot(equation)); type = $(QuoteNode(type)))))  # , modelmodule=$(modelmodule))))
end

macro steadystate(model, equation::Expr)
    thismodule = @__MODULE__
    return esc(:($(thismodule).setss!($(model), $(Meta.quot(equation)); type = :level))) # , modelmodule=$(modelmodule))))
end

"""
    initssdata!(m::AbstractModel)

Initialize the steady state data structure of the given model.

Do not call directly. This is an internal function, called during
[`@initialize`](@ref)

"""
function initssdata!(model::AbstractModel)
    ss = sstate(model)
    empty!(ss.vars)
    empty!(ss.values)
    empty!(ss.mask)
    for var in model.allvars
        push!(ss, var)
    end
    empty!(ss.equations)
    for eqn in alleqns(model)
        push!(ss.equations, make_sseqn(model, eqn, false))
    end
    if !model.flags.ssZeroSlope
        for eqn in alleqns(model)
            push!(ss.equations, make_sseqn(model, eqn, true))
        end
    end
    empty!(ss.constraints)
    return nothing
end

export issssolved
"""
    issssolved(sstate::SteadyStateData)

Return `true` if the steady state has been solved, or `false` otherwise.

"""
@inline issssolved(ss::SteadyStateData) = all(ss.mask)


text:
using BenchmarkTools
using Ripserer

include(joinpath(@__DIR__, "utils.jl"))
include(joinpath(@__DIR__, "benchmarks.jl"))

function run_benchmarks(benchmark_list; alg)
    @info "setting up $alg"
    group = BenchmarkGroup()
    for arg in benchmark_list
        data = load_data(joinpath(@__DIR__, "../datasets", arg.filename))
        group[arg.filename] = @benchmarkable ripserer(
            $data; $(arg.args)..., $(arg.best)..., alg=Symbol($alg)
        ) gcsample=true samples=5 seconds=1000
    end
    @info "running $alg"
    @time run(group)
end

"""
    run_and_save_all()

Run all involuted homology benchmarks.
"""
function run_and_save_all()
    timings(f) = joinpath(@__DIR__, "..", "results", "involuted", f)
    @info "cubical"
    BenchmarkTools.save(
        timings("cubical-cohomology.json"), run_benchmarks(cubical_benchmarks; alg="involuted")
    )
    BenchmarkTools.save(
        timings("cubical-involuted.json"), run_benchmarks(cubical_benchmarks; alg="cohomology")
    )
    BenchmarkTools.save(
        timings("cubical-homology.json"), run_benchmarks(cubical_benchmarks; alg="homology")
    )
    @info "rips"
    BenchmarkTools.save(
        timings("rips-cohomology.json"), run_benchmarks(rips_benchmarks; alg="involuted")
    )
    BenchmarkTools.save(
        timings("rips-involuted.json"), run_benchmarks(rips_benchmarks; alg="cohomology")
    )
    BenchmarkTools.save(
        timings("rips-homology.json"), run_benchmarks(rips_benchmarks; alg="homology")
    )
end

ext: jl


text:
abstract type AbstractDefault end

struct DataBasedDefault <: AbstractDefault end

function _set_defaults end

requires_defaults(::AbstractDefault) = true
requires_defaults(::Type) = false

function requires_defaults(obj)
    if nfields(obj) == 0
        return false
    else
        _fields = getfield.(Ref(obj), fieldnames(typeof(obj)))
        any(requires_defaults.(_fields))
    end
end

function set_defaults(obj, data...; hints = Dict())
    hints = Dict{Any,Any}(hints)
    get!(hints, :field_parents, [])

    if !requires_defaults(obj)
        return obj
    end

    _fieldnames = fieldnames(typeof(obj))
    _fields = getfield.(Ref(obj), _fieldnames)

    if any(isa.(_fields, AbstractDefault))
        updated_obj = _set_defaults(obj, data...; hints = hints)
        all(.!isa.(getfield.(Ref(updated_obj), _fieldnames), AbstractDefault)) ||
            throw("_set_defaults did not set all DataBasedDefaults. Please specify manually.")
        return set_defaults(updated_obj, data...; hints = hints)
    end

    idx = findfirst(requires_defaults.(_fields))
    subfield_hints = deepcopy(hints)
    push!(subfield_hints[:field_parents], obj)
    subfield = set_defaults(_fields[idx], data...; hints = subfield_hints)

    _lens = Setfield.PropertyLens{_fieldnames[idx]}()
    set_defaults(set(obj, _lens, subfield), data...; hints = hints)
end


text:
using COBREXA, Test

using Clp
using Distributed
using Downloads
using GLPK
using JSON
using JuMP
using LinearAlgebra
using MAT
using OrderedCollections
using OSQP
using SHA
using SparseArrays
using Statistics
using Tulip

# tolerance for comparing analysis results (should be a bit bigger than the
# error tolerance in computations)
TEST_TOLERANCE = 10 * COBREXA._constants.tolerance

function run_test_file(path...)
    fn = joinpath(path...)
    t = @elapsed include(fn)
    @info "$(fn) done in $(round(t; digits = 2))s"
end

function run_test_dir(dir, comment = "Directory $dir/")
    @testset "$comment" begin
        run_test_file.(joinpath.(dir, filter(fn -> endswith(fn, ".jl"), readdir(dir))))
    end
end

function check_data_file_hash(path, expected_checksum)
    actual_checksum = bytes2hex(sha256(open(path)))
    if actual_checksum != expected_checksum
        @error "The downloaded data file `$path' seems to be different from the expected one. Tests will likely fail." actual_checksum expected_checksum
    end
end

function download_data_file(url, path, hash)
    if isfile(path)
        check_data_file_hash(path, hash)
        @info "using cached `$path'"
        return path
    end

    Downloads.download(url, path)
    check_data_file_hash(path, hash)
    return path
end

# load the test models
run_test_file("data", "test_models.jl")

# import base files
@testset "COBREXA test suite" begin
    run_test_dir(joinpath("base", "types", "abstract"), "Abstract types")
    run_test_dir(joinpath("base", "types"), "Base model types")
    run_test_dir(joinpath("base", "logging"), "Logging")
    run_test_dir("base", "Base functionality")
    run_test_dir("io", "I/O functions")
    run_test_dir("reconstruction")
    run_test_dir("analysis")
    run_test_dir(joinpath("analysis", "sampling"), "Sampling")
end


text:
include("factors_sandbox.jl")

using Optim

# TODO 2 variables x1 and x2 of type SE2 one factor and one prior on x1

# X = hat(M, ϵSE2, [1.1, 0.1, pi+0.1])
X = hat(M, ϵSE2, [1.1, 0.1, pi + 0.05])
priorx1 = compose(M, ϵSE2, exp(M, ϵSE2, X))

# X = hat(M, ϵSE2, [1.1, 0.1, pi-0.1])
X = hat(M, ϵSE2, [0.9, -0.1, pi + 0.15])
priorx2 = compose(M, ϵSE2, exp(M, ϵSE2, X))

measX = hat(M, ϵSE2, [1, 0, pi/4])
measx1x2 = compose(M, ϵSE2, exp(M, ϵSE2, measX))

X2 = hat(M, ϵSE2, [0, 0, pi + pi/4])
x2 = compose(M, ϵSE2, exp(M, ϵSE2, X2))

## ======================================================================================
## with Optim "point" as LieAlgebra coordinates and only one variable at a time as non-parametric
## ======================================================================================
using Manifolds
using Optim

M = SpecialEuclidean(2)

representation_size.(M.manifold.manifolds)

# algorithm=Optim.BFGS
algorithm=Optim.NelderMead
alg = algorithm()
manifold = ManifoldsMani(SpecialEuclidean(2))
alg = algorithm(;manifold, algorithmkwargs...)

test_retract = false
function Optim.retract!(MM::ManifoldsMani, X)
    test_retract && (X[3] = rem2pi(X[3], RoundNearest))
    return X
end

function Optim.project_tangent!(MM::ManifoldsMani, G, x)
    return G
end

options = Optim.Options(allow_f_increases=true,
                        iterations = 200,
                        time_limit = 100,
                        show_trace = true,
                        show_every = 10,
                        )

# NOTE
# Only for Lie Groups
# Input: Lie Algebra coordinates
# f: exp to Group and calcuate residual
# Output: sum of squares residual

function cost(X)
    x = exp(M, ϵSE2, hat(M, ϵSE2, X))
    return PriorPose2(priorx1, x)^2 + PriorPose2(priorx2, x)^2 + Pose2Pose2(measx1x2, x, x2)^2
end

initValues = @MVector [0.,0.,0.0]

@time result = Optim.optimize(cost, initValues, alg, options)

rv = Optim.minimizer(result)

## result on group:

rv_G = exp(M, ϵSE2, hat(M, ϵSE2, rv))

##
autodiff = :forward
# autodiff = :finite

initValues = [0.,0.,0.1]

tdtotalCost = Optim.TwiceDifferentiable((x)->cost(x), initValues; autodiff)

@time result = Optim.optimize(tdtotalCost, initValues, alg, options)

rv = Optim.minimizer(result)

##
H = Optim.hessian!(tdtotalCost, rv)
Σ = pinv(H)

## ======================================================================================
## with Optim flattened ProductRepr and only one variable at a time as non-parametric
## ======================================================================================
using Manifolds
using Optim

M = SpecialEuclidean(2)

representation_size.(M.manifold.manifolds)

struct ManifoldsMani <: Optim.Manifold
  mani::AbstractManifold
end

# flatten should be replaced by a view or @cast
function flatten(M::SpecialEuclidean{2}, p)
    return mapreduce(vec, vcat, p.parts)
end

fp = flatten(M, p)

function unflatten(M::SpecialEuclidean{2}, fp::MVector)
    ProductRepr(MVector{2}(fp[1:2]), MMatrix{2,2}(fp[3:6]))
end

_p = unflatten(M, fp)

function Optim.retract!(MM::ManifoldsMani, fx)
    M = MM.mani
    x = unflatten(M, fx)
    project!(M, x, x)
    fx .= flatten(M, x)
    return fx
end

function Optim.project_tangent!(MM::ManifoldsMani, fG, fx)
    M = MM.mani
    x = unflatten(M, fx)
    G = unflatten(M, fG)
    project!(M, G, x, G)
    fG .= flatten(M, G)
    return fG
end

# autodiff = :forward
algorithm=Optim.BFGS
autodiff = :finite
# algorithm=Optim.NelderMead # does not work with manifolds
algorithmkwargs=() # add manifold to overwrite computed one
options = Optim.Options(allow_f_increases=true,
                        iterations = 100,
                        time_limit = 100,
                        show_trace = true,
                        show_every = 1,
                        )
##

function cost(fx)
    x = unflatten(M, fx)
    return PriorPose2(priorx1, x)^2 + PriorPose2(priorx2, x)^2
end

manifold = ManifoldsMani(SpecialEuclidean(2))
alg = algorithm(;manifold, algorithmkwargs...)
# alg = algorithm(; algorithmkwargs...)

initValues = MVector(flatten(M, ϵSE2))
# tdtotalCost = Optim.TwiceDifferentiable((x)->cost(x), initValues; autodiff)
# result = Optim.optimize(tdtotalCost, initValues, alg, options)

result = Optim.optimize(cost, initValues, alg, options)

rv = Optim.minimizer(result)
unflatten(M, rv)


text:
using HTTP

"""
    HTTPStore

A basic HTTP store without any credentials. The underlying data is supposed to be
consolidated and only read operations are supported. This store is compatible to
datasets being served through the [xpublish](https://xpublish.readthedocs.io/en/latest/)
python package.
"""
struct HTTPStore <: AbstractStore
    url::String
end

function Base.getindex(s::HTTPStore, k::String)
r = HTTP.request("GET",joinpath(s.url,k),status_exception = false)
if r.status >= 300
    if r.status == 404
        nothing
    else
        error("Error connecting to $(s.url) :", String(r.body))
    end
else
    r.body
end
end
getsub(s::HTTPStore,n) = HTTPStore(joinpath(s.url,n))
zname(s::HTTPStore) = split(s.url,"/")[end]

## This is a server implementation for Zarr datasets

function zarr_req_handler(s::AbstractStore)
  if s[".zmetadata"] === nothing
    consolidate_metadata(s)
  end
  request -> begin
    k = request.target
    while startswith(k,"/")
      k = k[2:end]
    end
    k_split = filter(!isequal(".."),split(k,"/"))
    storenew = if length(k_split)>1
      foldl((ss,kk)->getsub(ss,kk),k_split[1:end-1],init = s)
    else
      s
    end
    # @show zname(storenew), k_split
    # @show keys(storenew), k_split[end]
    r = storenew[k_split[end]]

    try
      if r ===  nothing
        return HTTP.Response(404, "Error: Key $k not found")
      else
        return HTTP.Response(200, r)
      end
    catch e
      return HTTP.Response(404, "Error: $e")
    end
  end
end

HTTP.serve(s::AbstractStore, args...; kwargs...) = HTTP.serve(zarr_req_handler(s),args...;kwargs...)


text:
struct FiniteDiscreteMeasure{X<:AbstractVector,P<:AbstractVector}
    support::X
    p::P

    function FiniteDiscreteMeasure{X,P}(support::X, p::P) where {X,P}
        length(support) == length(p) || error("length of `support` and `p` must be equal")
        isprobvec(p) || error("`p` must be a probability vector")
        return new{X,P}(support, p)
    end
end

"""
    discretemeasure(
        support::AbstractVector,
        probs::AbstractVector{<:Real}=FillArrays.Fill(inv(length(support)), length(support)),
    )

Construct a finite discrete probability measure with `support` and corresponding
`probabilities`. If the probability vector argument is not passed, then
equal probability is assigned to each entry in the support.

# Examples
```julia
using KernelFunctions
# rows correspond to samples
μ = discretemeasure(RowVecs(rand(7,3)), normalize!(rand(10),1))

# columns correspond to samples, each with equal probability
ν = discretemeasure(ColVecs(rand(3,12)))
```

!!! note
    If `support` is a 1D vector, the constructed measure will be sorted,
    e.g. for `mu = discretemeasure([3, 1, 2],[0.5, 0.2, 0.3])`, then
    `mu.support` will be `[1, 2, 3]` and `mu.p` will be `[0.2, 0.3, 0.5]`.
    Also, avoid passing 1D distributions as `RowVecs(rand(3))` or `[[1],[3],[4]]`,
    since this will be dispatched to the multivariate case instead
    of the univariate case for which the algorithm is more efficient.

!!! warning
    This function and in particular its return values are not stable and might be changed in future releases.
"""
function discretemeasure(
    support::AbstractVector{<:Real},
    probs::AbstractVector{<:Real}=Fill(inv(length(support)), length(support)),
)
    return DiscreteNonParametric(support, probs)
end
function discretemeasure(
    support::AbstractVector,
    probs::AbstractVector{<:Real}=Fill(inv(length(support)), length(support)),
)
    return FiniteDiscreteMeasure{typeof(support),typeof(probs)}(support, probs)
end

Distributions.support(d::FiniteDiscreteMeasure) = d.support
Distributions.probs(d::FiniteDiscreteMeasure) = d.p


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitorIcmpConfiguration <: SwaggerModel
    disableTraceRoute::Any # spec type: Union{ Nothing, Bool } # spec name: disableTraceRoute

    function ConnectionMonitorIcmpConfiguration(;disableTraceRoute=nothing)
        o = new()
        validate_property(ConnectionMonitorIcmpConfiguration, Symbol("disableTraceRoute"), disableTraceRoute)
        setfield!(o, Symbol("disableTraceRoute"), disableTraceRoute)
        o
    end
end # type ConnectionMonitorIcmpConfiguration

const _property_map_ConnectionMonitorIcmpConfiguration = Dict{Symbol,Symbol}(Symbol("disableTraceRoute")=>Symbol("disableTraceRoute"))
const _property_types_ConnectionMonitorIcmpConfiguration = Dict{Symbol,String}(Symbol("disableTraceRoute")=>"Bool")
Base.propertynames(::Type{ ConnectionMonitorIcmpConfiguration }) = collect(keys(_property_map_ConnectionMonitorIcmpConfiguration))
Swagger.property_type(::Type{ ConnectionMonitorIcmpConfiguration }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorIcmpConfiguration[name]))}
Swagger.field_name(::Type{ ConnectionMonitorIcmpConfiguration }, property_name::Symbol) =  _property_map_ConnectionMonitorIcmpConfiguration[property_name]

function check_required(o::ConnectionMonitorIcmpConfiguration)
    true
end

function validate_property(::Type{ ConnectionMonitorIcmpConfiguration }, name::Symbol, val)
end


text:
using Documenter, FlexibilityAnalysis, JuMP

makedocs(modules = [FlexibilityAnalysis],
         format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
         sitename = "FlexibilityAnalysis.jl - A Framework for Flexibility Analysis",
         authors = "Joshua Pulsipher",
         pages = ["Home" => "index.md",
                  "Background" => "background.md",
                  "User Guide" => "guide.md",
                  "Examples" => "examples.md",
                  "Library" => "api.md"])

 deploydocs(
     repo   = "github.com/pulsipher/FlexibilityAnalysis.jl.git",
     target = "build",
     deps   = nothing,
     make   = nothing
 )


text:
using ADFPEPS
using ADFPEPS: parity_conserving,double_ipeps_energy
using CUDA
using LinearAlgebra
using Random
using Test
using VUMPS:num_grad
using Zygote
CUDA.allowscalar(false)

@testset "parity_conserving" for atype in [Array], dtype in [ComplexF64], Ni = [2], Nj = [2]
    Random.seed!(100)
    Nv = 1
    D = 2^Nv
    T = atype(rand(dtype,D,D,4,D,D,4))
    function foo(T)
        ipeps = reshape([parity_conserving(T[:,:,:,:,:,i]) for i = 1:4], (2, 2))
        norm(ipeps)
    end
    @test Zygote.gradient(foo, T)[1] ≈ num_grad(foo, T) atol = 1e-8
end

@testset "energy" for atype in [Array], dtype in [ComplexF64], Ni = [1], Nj = [1]
    Random.seed!(100)
    folder = "E:/1 - research/4.9 - AutoDiff/data/ADFPEPS/"
    model = Hubbard(1.0,4.0,2.0)
    ipeps, key = init_ipeps(model; Ni = Ni, Nj = Nj, atype = atype, folder = folder, D=2, χ=4, tol=1e-10, maxiter=10)
    function foo(ipeps)
        double_ipeps_energy(ipeps, model; Ni=Ni,Nj=Nj,χ=χ,maxiter=10,infolder=folder,outfolder=folder)
    end
    folder, model, Ni, Nj, atype, D, χ, tol, maxiter = key
    @test Zygote.gradient(foo, ipeps)[1] ≈ num_grad(foo, ipeps) atol = 1e-3
end


text:
module TableFrames

export DataFrame, record, columns, rows

using NamedTuples
using DataStructures

const Store{T<:AbstractVector} = Associative{Symbol, T}

struct DataFrame{S<:Store{<:AbstractVector}, cols}
    data::S
end

function DataFrame(data::Store{<:AbstractVector})
    samelen(values(data)) || throw(ArgumentError("All columns of the DataFrame must be of the same length"))
    DataFrame{typeof(data), tuple(keys(data)...)}(data)
end

function DataFrame(data::Pair{Symbol, <:AbstractVector}...)
    T = typejoin((typeof(d.second) for d in data)...)
    return DataFrame(OrderedDict{Symbol, T}(data...))
end

Base.eltype(df::DataFrame) = eltype(df.data)
Base.isempty(df::DataFrame) = isempty(df.data)
Base.size(df::DataFrame) = isempty(df) ? (0, 0) : (length(df.data), length(first(values(df.data))))
Base.size(df::DataFrame, dim::Int) = size(df)[dim]

#####################
# Indexing Protocol
#####################
# Single-column based indexing
Base.getindex(df::DataFrame, col::Symbol) = df.data[col]
Base.getindex(df::DataFrame, col::Symbol, row::Int) = df.data[col][row]
Base.getindex(df::DataFrame, col::Symbol, rows::AbstractVector{Int}) = view(df.data[col], rows)

# Multi-column based indexing
function Base.getindex(df::DataFrame, cols::Vector{Symbol})
    selected = map(cols) do col
        col => df.data[col]
    end |> typeof(df.data)

    return DataFrame(selected)
end

function Base.getindex(df::DataFrame, cols::Vector{Symbol}, row::Int)
    return map(cols) do col
        col => df.data[col][row]
    end |> record
end

function Base.getindex(df::DataFrame, cols::Vector{Symbol}, rows::AbstractVector{Int})
    selected = map(cols) do col
        col => view(df.data[col], rows)
    end |> typeof(df.data)

    return DataFrame(selected)
end

# Row-only based indexing
Base.getindex(df::DataFrame, row::Int) = getindex(df, collect(keys(df.data)), row)
Base.getindex(df::DataFrame, rows::AbstractVector{Int}) = getindex(df, collect(keys(df.data)), rows)

# Set single value
Base.setindex!(df::DataFrame, val::Any, col::Symbol, row::Int)= df[col][row] = val

# Set entire column
function Base.setindex!(df::DataFrame, val::AbstractVector, col::Symbol)
    samelen(values(df), length(val)) || throw(ArgumentError("All columns of the DataFrame must be of the same length"))
    df[col] = val
end

# Set a row
function Base.setindex!(df::DataFrame, val::Union{NamedTuple, Associative}, row::Int)
    for col in keys(val)
        df[col, row] = val[col]
    end
end

#######################################
# Iteration Protocols
#######################################
Base.keys(df::DataFrame) = keys(df.data)
Base.values(df::DataFrame) = values(df.data)
names(df::DataFrame) = keys(df)
columns(df::DataFrame) = values(df)

struct RowIterator
    df::DataFrame
end

rows(df::DataFrame) = RowIterator(df)
Base.start(itr::RowIterator) = 1
Base.done(itr::RowIterator, i::Int) = i > size(itr.df, 2)
Base.next(itr::RowIterator, i::Int) = (itr.df[i], i + 1)
Base.length(itr::RowIterator) = size(itr.df, 2)
Base.eltype(itr::RowIterator) = NamedTuple

################################################
# Methods for creating NamedTuple row records.
###############################################
record(fields::Pair...) = record([fields...])

function record(fields::Vector{<:Pair})
    unzipped = collect(zip(fields...))
    return record([unzipped[1]...], [unzipped[2]...])
end

function record(names::Vector{<:Symbol}, vals::Vector)
    nt = NamedTuples.create_namedtuple_type(names)
    return nt(vals...)
end

#=
Simply checks that all columns are the same length
or that all columns match the default value.
=#
function samelen(cols, default=-1)
    result = true
    len = default

    for col in cols
        if len < 0
            len = length(col)
        elseif length(col) != len
            result = false
            break
        end
    end

    return result
end

include("show.jl")

end  # module


text:
##### Generic methods #####

## sampling

function rand!(d::MultivariateDistribution, x::AbstractVector)
    length(x) == length(d) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(d, x)
end

function rand!(d::MultivariateDistribution, A::AbstractMatrix)
    size(A,1) == length(d) ||
        throw(DimensionMismatch("Output size inconsistent with sample length."))
    _rand!(sampler(d), A)
end

rand(d::MultivariateDistribution) = _rand!(d, Array(eltype(d), length(d)))
rand(d::MultivariateDistribution, n::Int) = _rand!(sampler(d), Array(eltype(d), length(d), n))

## domain

function insupport!{D<:MultivariateDistribution}(r::AbstractArray, d::Union{D,Type{D}}, X::AbstractMatrix)
    n = length(r)
    size(X) == (length(d),n) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    for i in 1:n
        @inbounds r[i] = insupport(d, view(X, :, i))
    end
    return r
end

insupport{D<:MultivariateDistribution}(d::Union{D,Type{D}}, X::AbstractMatrix) =
    insupport!(BitArray(size(X,2)), d, X)

## statistics

entropy(d::MultivariateDistribution, b::Real) = entropy(d) / log(b)

function cor(d::MultivariateDistribution)
    C = cov(d)
    n = size(C, 1)
    @assert size(C, 2) == n
    R = Array(eltype(C), n, n)

    for j = 1:n
        for i = 1:j-1
            @inbounds R[i, j] = R[j, i]
        end
        R[j, j] = 1.0
        for i = j+1:n
            @inbounds R[i, j] = C[i, j] / sqrt(C[i, i] * C[j, j])
        end
    end

    return R
end

# pdf and logpdf

_pdf(d::MultivariateDistribution, X::AbstractVector) = exp(_logpdf(d, X))

function logpdf(d::MultivariateDistribution, X::AbstractVector)
    length(X) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf(d, X)
end

function pdf(d::MultivariateDistribution, X::AbstractVector)
    length(X) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf(d, X)
end

function _logpdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    for i in 1 : size(X,2)
        @inbounds r[i] = logpdf(d, view(X,:,i))
    end
    return r
end

function _pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    for i in 1 : size(X,2)
        @inbounds r[i] = pdf(d, view(X,:,i))
    end
    return r
end

function logpdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    size(X) == (length(d), length(r)) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _logpdf!(r, d, X)
end

function pdf!(r::AbstractArray, d::MultivariateDistribution, X::AbstractMatrix)
    size(X) == (length(d), length(r)) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _pdf!(r, d, X)
end

function logpdf(d::MultivariateDistribution, X::AbstractMatrix)
    size(X, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    T = promote_type(partype(d), eltype(X))
    _logpdf!(Array(T, size(X,2)), d, X)
end

function pdf(d::MultivariateDistribution, X::AbstractMatrix)
    size(X, 1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    T = promote_type(partype(d), eltype(X))
    _pdf!(Array(T, size(X,2)), d, X)
end

## log likelihood

function _loglikelihood(d::MultivariateDistribution, X::AbstractMatrix)
    ll = 0.0
    for i in 1:size(X, 2)
        ll += _logpdf(d, view(X,:,i))
    end
    return ll
end

function loglikelihood(d::MultivariateDistribution, X::AbstractMatrix)
    size(X,1) == length(d) ||
        throw(DimensionMismatch("Inconsistent array dimensions."))
    _loglikelihood(d, X)
end

##### Specific distributions #####

for fname in ["dirichlet.jl",
              "multinomial.jl",
              "mvnormal.jl",
              "mvnormalcanon.jl",
              "mvlognormal.jl",
              "mvtdist.jl",
              "vonmisesfisher.jl"]
    include(joinpath("multivariate", fname))
end


text:
function solve_coefficients!(
            bands::AbstractMatrix{Float64},
            bb::AbstractVector{T},
            coefs::AbstractVector{T}
    ) where T

    M = size(coefs,1) - 2

    # Solve interpolating equations
    # First and last rows are different

    bands[1, 2] /= bands[1, 1]
    bands[1, 3] /= bands[1, 1]
    bb[1] /= bands[1, 1]
    bands[1, 1] = 1.0
    bands[2, 2] -= bands[2, 1] * bands[1, 2]
    bands[2, 3] -= bands[2, 1] * bands[1, 3]
    bb[2] -= bands[2, 1] * bb[1]
    bands[1, 1] = 0.0
    bands[2, 3] /= bands[2, 2]
    bb[2] /= bands[2, 2]
    bands[2, 2] = 1.0

    # Now do rows 2 through M+1
    for row = 3:(M+1)
        bands[row, 2] -= bands[row, 1] * bands[row - 1, 3]
        bb[row] -= bands[row, 1] * bb[row - 1]
        bands[row, 3] /= bands[row, 2]
        bb[row] /= bands[row, 2]
        bands[row, 1] = 0.0
        bands[row, 2] = 1.0
    end
    # Do last row
    bands[M + 2, 2] -= bands[M + 2, 1] * bands[M, 3]
    bb[M + 2] -= bands[M + 2, 1] * bb[M]
    bands[M + 2, 3] -= bands[M + 2, 2] * bands[M + 1, 3]
    bb[M + 2] -= bands[M + 2, 2] * bb[M + 1]
    bb[M + 2] /= bands[M + 2, 3]
    bands[M + 2, 3] = 1.0

    coefs[M + 2] = bb[M + 2]
    # Now back substitute up
    # for row in (M+1):2
    for ii in 1:M
        row = M+2 - ii
        coefs[row] = bb[row] - bands[row, 3] * coefs[row + 1]
    end
    # Finish with first row
    coefs[1] = bb[1] - bands[1, 2] * coefs[2] - bands[1, 3] * coefs[3]

end

function fill_bands!(M::Int, bands::AbstractMatrix, bb::AbstractVector{T}, data::AbstractVector{T}) where T

    # boundary conditions
    bands[1, 1] = 1.0
    bands[1, 2] = -2.0
    bands[1, 3] = 1.0
    bands[M+2, 1] = 1.0
    bands[M+2, 2] = -2.0
    bands[M+2, 3] = 1.0

    bands[2:end-1,1] .= 1.0/6.0
    bands[2:end-1,2] .= 2.0/3.0
    bands[2:end-1,3] .= 1.0/6.0

    bb[:] = data[:]
    bb[1] = zero(T)
    bb[M+2] = zero(T)

end

function prefilter!(data::AbstractVector{T}) where T

    # data is 1d array of length M
    # data is 1d array of length M+2
    M = length(data) - 2
    bands = zeros(M+2, 3)
    bb = zeros(T, M+2)
    fill_bands!(M, bands, bb, data)
    prefilter!(data, bands, bb)

end

function prefilter!(data, bands, bb)

    M = length(data) - 2
    fill_bands!(M, bands, bb, data)
    solve_coefficients!(bands, bb, data)

end

function prefilter!(data::Array{T,1}) where T
    M = size(data,1)-2
    bands = zeros(M+2,3)
    bb = zeros(T,M+2)
    fill_bands!(M, bands, bb, data)
    prefilter!(data, bands, bb)
end

function prefilter!(data::Array{T,2}) where T
    I,J = size(data)

    M = J-2
    bands = zeros(M+2,3)
    bb = zeros(T,M+2)
    for i=1:I
            dat = view(data, i, :)
            # prefilter!(dat)
            fill_bands!(M, bands, bb, dat)
            prefilter!(dat, bands, bb)
    end
    M = I-2
    bands = zeros(M+2, 3)
    bb = zeros(T,M+2)
    for j=1:J
        dat = view(data, :, j)
        # prefilter!(dat)
        #
        fill_bands!(M, bands, bb, dat)
        prefilter!(dat, bands, bb)
    end
end

function prefilter!(data::Array{T,3}) where T
    I,J,K = size(data)

    M = K-2
    bands = zeros(M+2,3)
    bb = zeros(T,M+2)
    for i=1:I
        for j=1:J
            dat = view(data, i,j, :)
            fill_bands!(M, bands, bb, dat)
            prefilter!(dat, bands, bb)
        end
    end
    M = J-2
    bands = zeros(M+2, 3)
    bb = zeros(T,M+2)
    for i=1:I
        for k=1:K
            dat = view(data, i,:, k)
            fill_bands!(M, bands, bb, dat)
            prefilter!(dat, bands, bb)
        end
    end
    M = I-2
    bands = zeros(M+2, 3)
    bb = zeros(T,M+2)
    for j=1:J
        for k=1:K
            dat = view(data, :, j,k)
            fill_bands!(M, bands, bb, dat)
            prefilter!(dat, bands, bb)
        end
    end
end

function prefilter!(data::Array{T,4}) where T
    I,J,K,L = size(data)

    M = L-2
    bands = zeros(M+2,3)
    bb = zeros(T,M+2)
    for i=1:I
        for j=1:J
            for k=1:K
                dat = view(data,i,j,k,:)
                fill_bands!(M, bands, bb, dat)
                prefilter!(dat, bands, bb)
            end
        end
    end
    M = K-2
    bands = zeros(M+2,3)
    bb = zeros(T,M+2)
    for i=1:I
        for j=1:J
            for l=1:L
                dat = view(data,i,j,:,l)
                fill_bands!(M, bands, bb, dat)
                prefilter!(dat, bands, bb)
            end
        end
    end
    M = J-2
    bands = zeros(M+2, 3)
    bb = zeros(T,M+2)
    for i=1:I
        for k=1:K
            for l=1:L
                dat = view(data, i,:, k,l)
                fill_bands!(M, bands, bb, dat)
                prefilter!(dat, bands, bb)
            end
        end
    end
    M = I-2
    bands = zeros(M+2, 3)
    bb = zeros(T,M+2)
    for j=1:J
        for k=1:K
            for l=1:L
                dat = view(data, :, j,k,l)
                fill_bands!(M, bands, bb, dat)
                prefilter!(dat, bands, bb)
            end
        end
    end
end


text:
# This is the run function used for commandline access
# for the DECaPS2 survey error bar correction. It can be
# used as a model for other survey specific run structures.

using CloudCovErr
using ArgParse

"""
    parse_commandline()

A simple commandline parser that determines which files are run and makes most
code parameters user-facing.
"""
function parse_commandline()
    s=ArgParseSettings()
    @add_arg_table s begin
        # required file/dir name parameters
        "base"
            required=true
            arg_type=String
            help="base directory containing raw exposure image files"
        "date"
            required=true
            arg_type=String
            help="file name date (exposure datetime)"
        "filt"
            required=true
            arg_type=String
            help="file name filter"
        "vers"
            required=true
            arg_type=String
            help="file name version"
        "basecat"
            required=true
            arg_type=String
            help="base directory containing crowdsource cat files"
        # method parameter values
        "--thr", "-t"
            help = "threshold for psf-based masking of the residuals (int for ease)"
            arg_type = Int
            default = 20
        "--outthr"
            help = "threshold for residual-based masking (int for ease)"
            arg_type = Int
            default = 20000
        "--Np", "-p"
            help = "sidelength size of spatial covariance matrix (must be int, odd)"
            arg_type = Int
            default = 33
        "--wx", "-w"
            help = "sidelength size of sample averaging (must be int, odd)"
            arg_type = Int
            default = 129
        "--wy"
            help = "sidelength size of sample averaging (must be int, odd)"
            arg_type = Int
            default = -1
        "--tilex"
            help = "number of tiles to divide image into along x"
            arg_type = Int
            default = 1
        "--tiley"
            help = "number of tiles to divide image into along y"
            arg_type = Int
            default = -1
        "--cS7"
            help = "when set attempts to correct the amplifier offset in S7 fields"
            action = :store_true
        "--ftype", "-f"
            help = "use float n (32 or 64)"
            arg_type = Int
            default = 32
        # run settings
        "--resume", "-r"
            help = "resume partially completed file where leftoff"
            action = :store_true
        # limit run scope/specific subsets for testing
        "--ccdlist", "-c"
            help = "list of ccd_CAT names to be run"
            nargs = '+'
            default = String[]
    end
    return parse_args(s)
end

"""
    run_wrapper()

Reads commandline arguments and passes them to `proc_all` or start
processing all or a subset of a given exposure.
"""
function run_wrapper()
    parg = parse_commandline()

    # handle default x==y conditions
    if parg["wy"] == -1
        wy = parg["wx"]
    else
        wy = parg["wy"]
    end
    if parg["tiley"] == -1
        tiley = parg["tilex"]
    else
        tiley = parg["tiley"]
    end

    CloudCovErr.proc_all(parg["base"],parg["date"],parg["filt"],parg["vers"],
        parg["basecat"],ccdlist=parg["ccdlist"],resume=parg["resume"],corrects7=parg["cS7"],
        thr=parg["thr"],outthr=parg["outthr"],Np=parg["Np"],widx=parg["wx"],widy=wy,tilex=parg["tilex"],
        tiley=tiley,ftype=parg["ftype"])
end

run_wrapper()


text:
# Some simple tests on the entropy/information code

using SubmodularMaximization

grid_size = 10

grid = Grid(grid_size, grid_size)

sensor = RangingSensor()
histogram_filter = Filter(grid)

# Test entropy on uniform prior
@show grid_entropy = entropy(histogram_filter)
real_entropy = -log(2, 1/num_states(grid))
if !isapprox(grid_entropy, real_entropy; rtol=1e-3)
  println("Incorrect entropy for uniform prior")
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
using Readables, Test

# using the default settings

@test readablestring(0) == "0"
@test readablestring(12) == "12"
@test readablestring(123) == "123"
@test readablestring(1234) == "1,234"
@test readablestring(12345) == "12,345"
@test readablestring(123456) == "123,456"
@test readablestring(1234567) == "1,234,567"

@test readablestring(0.0) == "0.0"
@test readablestring(12.0) == "12.0"
@test readablestring(123.0) == "123.0"
@test readablestring(1234.0) == "1,234.0"
@test readablestring(12345.0) == "12,345.0"
@test readablestring(123456.0) == "123,456.0"
@test readablestring(1234567.0) == "1.23456_7e+6"

@test readablestring(0.12345) == "0.12345"
@test readablestring(12.12345) == "12.12345"
@test readablestring(123.12345) == "123.12345"
@test readablestring(1234.12345) == "1,234.12345"
@test readablestring(12345.12345) == "12,345.12345"
@test readablestring(123456.12345) == "123,456.12345"
@test readablestring(1234567.12345) == "1.23456_71234_5e+6"

@test readablestring(0.12345678) == "0.12345_678"
@test readablestring(12.12345678) == "12.12345_678"
@test readablestring(123.12345678) == "123.12345_678"
@test readablestring(1234.12345678) == "1,234.12345_678"
@test readablestring(12345.12345678) == "12,345.12345_678"
@test readablestring(123456.12345678) == "123,456.12345_678"

@test readablestring(1234567.12325675) == "1.23456_71232_5675e+6"
@test readablestring(BigFloat("1234567.12325675")) == "1.23456_71232_5675e+06"


text:
## ShortRateTree methods ##
get_state_prices!(tree::ShortRateTree, i::Int) = get_state_prices!(tree.treeLattice, i)

struct PrivateConstraint{P <: Parameter} <: Constraint
  arguments::Vector{P}
end

# type ShortRateModelCommon
#   observers::Vector
# end
#
# ShortRateModelCommon() = ShortRateModelCommon([])

function QuantLib.Math.test(c::PrivateConstraint, x::Vector{Float64})
  k = 1
  for i = 1:length(c.arguments)
    sz = length(c.arguments[i].data)
    testParams = zeros(sz)
    for j = 1:sz
      testParams[j] = x[k]
      k += 1
    end
    if !test_params(c.arguments[i], testParams)
      return false
    end
  end

  return true
end

mutable struct CalibrationFunction{M <: ShortRateModel, C <: CalibrationHelper} <: CostFunction
  model::M
  helpers::Vector{C}
  weights::Vector{Float64}
  projection::Projection
end

struct SolvingFunction <: Function
  lambda::Vector{Float64}
  Bb::Vector{Float64}
end

function (solvFunc::SolvingFunction)(y::Float64)
  value_ = 1.0
  for i = 1:length(solvFunc.lambda)
    value_ -= solvFunc.lambda[i] * exp(-solvFunc.Bb[i] * y)
  end
  return value_
end

# function operator(solvFunc::SolvingFunction)
#   function _inner(y::Float64)
#     value_ = 1.0
#     for i = 1:length(solvFunc.lambda)
#       value_ -= solvFunc.lambda[i] * exp(-solvFunc.Bb[i] * y)
#     end
#     return value_
#   end
#
#   return _inner
# end

function func_values(calibF::CalibrationFunction, params::Vector{Float64})
  set_params!(calibF.model, include_params(calibF.projection, params))
  values = zeros(length(calibF.helpers))
  for i = 1:length(values)
    values[i] = calibration_error(calibF.helpers[i].calibCommon.calibrationErrorType, calibF.helpers[i]) * sqrt(calibF.weights[i])
  end

  return values
end

function value(calibF::CalibrationFunction, params::Vector{Float64})
  set_params!(calibF.model, include_params(calibF.projection, params))
  _value = 0.0
  for i = 1:length(calibF.helpers)
    diff = calibration_error(calibF.helpers[i].calibCommon.calibrationErrorType, calibF.helpers[i])
    _value += diff * diff * calibF.weights[i]
  end

  return sqrt(_value)
end

# accessor methods ##
get_a(m::ShortRateModel) = m.a.data[1]
get_sigma(m::ShortRateModel) = m.sigma.data[1]
get_b(m::ShortRateModel) = m.b.data[1]

check_params_equal(m::ShortRateModel, params::Vector{Float64}) = get_params(m) == params

function set_params!(model::ShortRateModel, params::Vector{Float64})
  if check_params_equal(model, params)
    return model
  end
  paramCount = 1
  args = model.privateConstraint.arguments
  for i = 1:length(args)
    arg = model.privateConstraint.arguments[i]
    for j = 1:length(arg.data)
       arg.data[j] = params[paramCount]
       paramCount += 1
     end
   end
   generate_arguments!(model)

   notify_observers!(model)
   # model.phi = G2FittingParameter(get_a(model), get_sigma(model), get_b(model), get_eta(model), get_rho(model), model.ts)
   return model
 end

## Short Rate Model calibration function #
function calibrate!(model::ShortRateModel, instruments::Vector{C}, method::OptimizationMethod, endCriteria::EndCriteria,
                    constraint::Constraint = model.privateConstraint, weights::Vector{Float64} = ones(length(instruments)),
                    fixParams::BitArray{1} = BitArray(undef, 0)) where {C <: CalibrationHelper}

  w = length(weights) == 0 ? ones(length(instruments)) : weights
  prms = get_params(model)
  # println("model params: ", prms)
  all = falses(length(prms))
  proj = Projection(prms, length(fixParams) > 0 ? fixParams : all)
  calibFunc = CalibrationFunction(model, instruments, w, proj)
  pc = ProjectedConstraint(constraint, proj)
  prob = Problem(calibFunc, pc, project(proj, prms))

  # minimization
  minimize!(method, prob, endCriteria)
  res = prob.currentValue
  set_params!(model, include_params(proj, res))

  return model
end


text:
# ##############################################################################
# ZonedDateTime
# ##############################################################################
# the structure
# utc_datetime : isbits
# timezone : VariableTimeZone
#   name : string
#   transitions : Vector{TimeZones.Transitions}
#     utc_datetime : isbits
#     zone: FixedTimeZone
#       name : string
#       offset : isbits
#   cutoff : isbits (Union{Nothing, DateTime})
# zone
#   name : string
#   offset : isbits
compress_then_write(b::Vector{ZonedDateTime}, io) = begin
    utc_datetime_meta = compress_then_write([b.utc_datetime for b in b], io)
    timezone_meta = compress_then_write([b.timezone for b in b], io)
    zone_meta = compress_then_write([b.zone for b in b], io)

    (
        type = ZonedDateTime,
        len = utc_datetime_meta.len + timezone_meta.len + zone_meta.len,
        utc_datetime_meta = utc_datetime_meta,
        timezone_meta = timezone_meta,
        zone_meta = zone_meta,
    )
end

# load a ZonedDateTime column
column_loader(::Type{ZonedDateTime}, io, metadata) = begin
    v1 = column_loader(metadata.utc_datetime_meta.type, io, metadata.utc_datetime_meta)
    v2 = column_loader(metadata.timezone_meta.type, io, metadata.timezone_meta)
    v3 = column_loader(metadata.zone_meta.type, io, metadata.zone_meta)
    [ZonedDateTime(v1, v2, v3) for (v1, v2, v3) in zip(v1, v2, v3)]
end

column_loader!(_, ::Type{ZonedDateTime}, io, metadata) = begin
    column_loader(ZonedDateTime, io, metadata)
end

some_elm(::Type{ZonedDateTime}) = ZonedDateTime(1984, 1, 5, tz"Asia/Shanghai")

# ##############################################################################
# TimeZones.Transitions
# ##############################################################################
compress_then_write(b::Vector{TimeZones.Transition}, io) = begin
    # the structure of Transition is like
    # (utc_datetime::isbits, zone::string)
    utc_datetime_meta = compress_then_write([b.utc_datetime for b in b], io)
    zone_meta = compress_then_write([b.zone for b in b], io)

    (
        type = TimeZones.Transition,
        len = utc_datetime_meta.len + zone_meta.len,
        utc_datetime_meta = utc_datetime_meta,
        zone_meta = zone_meta,
    )
end

column_loader(::Type{TimeZones.Transition}, io, metadata) = begin
    v1 = column_loader(metadata.utc_datetime_meta.type, io, metadata.utc_datetime_meta)
    v2 = column_loader(metadata.zone_meta.type, io, metadata.zone_meta)
    [TimeZones.Transition(v1, v2) for (v1, v2) in zip(v1, v2)]
end

column_loader!(::Type{TimeZones.Transition}, io, metadata) = begin
    column_loader(TimeZones.Transition, io, metadata)
end

# ##############################################################################
# VariableTimeZone
# ##############################################################################
compress_then_write(b::Vector{VariableTimeZone}, io) = begin
    # the structure of Transition is like
    # (name::string, transitions, cutoff::string)
    name_meta = compress_then_write([b.name for b in b], io)

    len_transitions = [length(b.transitions) for b in b]
    len_transitions_meta = compress_then_write(len_transitions, io)

    transitions = Vector{TimeZones.Transition}(undef, sum(len_transitions))
    i = 1
    for b in b
        for t in b.transitions
            transitions[i] = t
            i += 1
        end
    end
    transitions_meta = compress_then_write(transitions, io)
    cutoff_meta = compress_then_write([b.cutoff for b in b], io)

    (
        type = VariableTimeZone,
        len = name_meta.len +
              len_transitions_meta.len +
              transitions_meta.len +
              cutoff_meta.len,
        name_meta = name_meta,
        len_transitions_meta = len_transitions_meta,
        transitions_meta = transitions_meta,
        cutoff_meta = cutoff_meta,
    )
end

column_loader(::Type{VariableTimeZone}, io, metadata) = begin
    name = column_loader(metadata.name_meta.type, io, metadata.name_meta)

    len_transitions = column_loader(
        metadata.len_transitions_meta.type,
        io,
        metadata.len_transitions_meta,
    )
    hi = cumsum(len_transitions)
    lo = vcat(0, hi[1:end-1]) .+ 1

    transitions =
        column_loader(metadata.transitions_meta.type, io, metadata.transitions_meta)

    metadata.cutoff_meta.type
    cutoff = column_loader(metadata.cutoff_meta.type, io, metadata.cutoff_meta)

    [
        VariableTimeZone(name, transitions[lo:hi], cutoff)
        for (name, lo, hi, cutoff) in zip(name, lo, hi, cutoff)
    ]
end

column_loader!(::Type{VariableTimeZone}, io, metadata) = begin
    column_loader(VariableTimeZone, io, metadata)
end

# ##############################################################################
# FixedTimeZone
# ##############################################################################
compress_then_write(b::Vector{FixedTimeZone}, io) = begin
    # the structure of Transition is like
    # (name, offset)
    m1 = compress_then_write([b.name for b in b], io)
    m2 = compress_then_write([b.offset for b in b], io)

    (type = FixedTimeZone, len = m1.len + m2.len, name_meta = m1, offset_meta = m2)
end

column_loader(::Type{FixedTimeZone}, io, metadata) = begin
    v1 = column_loader(metadata.name_meta.type, io, metadata.name_meta)
    v2 = column_loader(metadata.offset_meta.type, io, metadata.offset_meta)
    [FixedTimeZone(v1, v2) for (v1, v2) in zip(v1, v2)]
end

column_loader!(::Type{FixedTimeZone}, io, metadata) = begin
    column_loader(FixedTimeZone, io, metadata)
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF


text:
filename = Pkg.dir("DataFrames", "test", "data", "sample_data.csv")

ds = DataStream(filename)
next_df = start(ds)
(df, next_df) = next(ds, next_df)
@assert done(ds, next_df) == false
(df, next_df) = next(ds, next_df)
@assert done(ds, next_df) == false
(df, next_df) = next(ds, next_df)
@assert done(ds, next_df) == true

filename = Pkg.dir("DataFrames", "test", "data", "big_data.csv")

ds = DataStream(filename, 100)
next_df = start(ds)
(df, next_df) = next(ds, next_df)
@assert done(ds, next_df) == false

ds = DataStream(filename, 5)

for minibatch in ds
  @assert isa(minibatch, DataFrame)
end

means = colmeans(ds)
@assert abs(means[1, 4] - (-0.005686449)) < 10e-4
@assert abs(means[1, 5] - 19.01197) < 10e-4

vars = colvars(ds)
@assert abs(vars[1, 4] - 0.98048) < 10e-4
@assert abs(vars[1, 5] - 0.989416) < 10e-4

(mins, maxs) = colranges(ds)
@assert abs(mins[1, 4] - (-4.33635)) < 10e-4
@assert abs(mins[1, 5] - 15.6219) < 10e-4
@assert abs(maxs[1, 4] - 3.86857) < 10e-4
@assert abs(maxs[1, 5] - 22.574) < 10e-4

covariances = cov(ds)
for i in 1:ncol(covariances)
  if isna(covariances[i, i])
    @assert isna(vars[1, i])
  else
    @assert abs(covariances[i, i] - vars[1, i]) < 10e-4
  end
end
@assert abs(covariances[4, 4] - 0.980479916) < 10e-4
@assert abs(covariances[4, 5] - 0.009823644) < 10e-4
@assert abs(covariances[5, 4] - 0.009823644) < 10e-4
@assert abs(covariances[5, 5] - 0.989415811) < 10e-4

correlations = cor(ds)
for i in ncol(correlations)
  for j in ncol(correlations)
    true_value = covariances[i, j] / sqrt(covariances[i, i] * covariances[j, j])
    if isna(true_value)
      @assert isna(correlations[i, j])
    else
      @assert abs(correlations[i, j] - true_value) < 10e-4
    end
  end
end

# Deal with different delimiters
filename = Pkg.dir("DataFrames", "test", "data", "sample_data.csv")
ds = DataStream(filename)
for row in ds
  @assert nrow(row) <= 1
end

filename = Pkg.dir("DataFrames", "test", "data", "sample_data.tsv")
ds = DataStream(filename)
for row in ds
  @assert nrow(row) <= 1
end

filename = Pkg.dir("DataFrames", "test", "data", "sample_data.wsv")
ds = DataStream(filename)
for row in ds
  @assert nrow(row) <= 1
end

# DataFrame to DataStream conversion
df = DataFrame(quote A = 1:25 end)

ds = DataStream(df, 5)

for mini in ds
	@assert nrow(mini) <= 5
end

colmeans(ds)


text:
# module for calculating standard deviation parameters
# robot moves for 4m at 0.1m/s on straight trajectory
# including random noise
# this is iterated for 100 times

module MotionTestForward
  using DataFrames, Statistics

  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function main(time_interval=0.1)
    w = World(-5.0, 5.0, -5.0, 5.0)

    rs = []
    init_pose = [0.0, 0.0, 0.0]
    straight = Agent(0.1, 0.0)
    for i in 1:100
      r = DifferentialWheeledRobot(init_pose, 0.2, "black",
                                   straight, time_interval,
                                   noise_per_meter=5,
                                   noise_std=pi/60)
      push!(rs, r)
    end

    for t in 0:time_interval:40
      draw(w)
      annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
      for r in rs
        draw!(r)
      end
    end

    poses = DataFrame(range=[sqrt(r.pose[1]^2 + r.pose[2]^2) for r in rs],
                      theta=[r.pose[3] for r in rs])
    theta_var = var(poses.theta, corrected=false)
    range_mean = mean(poses.range)
    println("Variance of range and theta: ", poses)
    println("Variance of theta: ", theta_var)
    println("Mean of range: ", range_mean)
    println("Standard deviation per m of theta: ", sqrt(theta_var/range_mean))
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
import MultivariatePolynomials
import MultivariatePolynomials: AbstractPolynomialLike

for (Z, AZ) in ((:ContinuousIdentitySystem, :AbstractContinuousSystem),
                (:DiscreteIdentitySystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z) <: $(AZ)
            statedim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    ContinuousIdentitySystem <: AbstractContinuousSystem

Trivial identity continuous-time system of the form:

```math
    x' = 0.
```

### Fields

- `statedim` -- number of state variables
"""
ContinuousIdentitySystem

@doc """
    DiscreteIdentitySystem <: AbstractDiscreteSystem

Trivial identity discrete-time system of the form:
```math
    x_{k+1} = x_k.
```

### Fields

- `statedim` -- number of state variables
"""
DiscreteIdentitySystem

for (Z, AZ) in ((:ConstrainedContinuousIdentitySystem, :AbstractContinuousSystem),
                (:ConstrainedDiscreteIdentitySystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){ST} <: $(AZ)
            statedim::Int
            X::ST
        end
        statedim(s::$Z) = s.statedim
        stateset(s::$Z) = s.X
        inputdim(::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedContinuousIdentitySystem <: AbstractContinuousSystem

Trivial identity continuous-time system with state constraints of the form:

```math
    x' = 0, x(t) ∈ \\mathcal{X}.
```

### Fields

- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedContinuousIdentitySystem

@doc """
    ConstrainedDiscreteIdentitySystem <: AbstractDiscreteSystem

Trivial identity discrete-time system with state constraints of the form:

```math
    x_{k+1} = x_k, x_k ∈ \\mathcal{X}.
```

### Fields

- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedDiscreteIdentitySystem

for (Z, AZ) in ((:LinearContinuousSystem, :AbstractContinuousSystem),
                (:LinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}} <: $(AZ)
            A::MT
        end
        statedim(s::$Z) = checksquare(s.A)
        inputdim(::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    LinearContinuousSystem

Continuous-time linear system of the form:

```math
    x' = A x.
```

### Fields

- `A` -- square matrix
"""
LinearContinuousSystem

@doc """
    LinearDiscreteSystem

Discrete-time linear system of the form:

```math
    x_{k+1} = A x_k.
```

### Fields

- `A` -- square matrix
"""
LinearDiscreteSystem

for (Z, AZ) in ((:AffineContinuousSystem, :AbstractContinuousSystem),
                (:AffineDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}} <: $(AZ)
            A::MT
            b::VT
            function $(Z)(A::MT, b::VT) where {T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}}
                @assert checksquare(A) == length(b)
                return new{T, MT, VT}(A, b)
            end
        end
        statedim(s::$Z) = length(s.b)
        inputdim(::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = true
    end
end

@doc """
    AffineContinuousSystem

Continuous-time affine system of the form:

```math
    x' = A x + b.
```

### Fields

- `A` -- square matrix
- `b` -- vector
"""
AffineContinuousSystem

@doc """
    AffineDiscreteSystem

Discrete-time affine system of the form:

```math
    x_{k+1} = A x_k + b.
```

### Fields

- `A` -- square matrix
- `b` -- vector
"""
AffineDiscreteSystem

for (Z, AZ) in ((:LinearControlContinuousSystem, :AbstractContinuousSystem),
                (:LinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            B::MTB
            function $(Z)(A::MTA, B::MTB) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}}
                @assert checksquare(A) == size(B, 1)
                return new{T, MTA, MTB}(A, B)
            end
        end
        statedim(s::$Z) = checksquare(s.A)
        inputdim(s::$Z) = size(s.B, 2)
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    LinearControlContinuousSystem

Continuous-time linear control system of the form:

```math
    x' = A x + B u.
```

### Fields

- `A` -- square matrix
- `B` -- matrix
"""
LinearControlContinuousSystem

@doc """
    LinearControlDiscreteSystem

Discrete-time linear control system of the form:

```math
    x_{k+1} = A x_k + B u_k.
```

### Fields

- `A` -- square matrix
- `B` -- matrix
"""
LinearControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, ST} <: $(AZ)
            A::MT
            X::ST
        end
        statedim(s::$Z) = checksquare(s.A)
        stateset(s::$Z) = s.X
        inputdim(::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedLinearContinuousSystem

Continuous-time linear system with state constraints of the form:

```math
    x' = A x, x(t) ∈ \\mathcal{X}.
```

### Fields

- `A` -- square matrix
- `X` -- state constraints
"""
ConstrainedLinearContinuousSystem

@doc """
    ConstrainedLinearDiscreteSystem

Discrete-time linear system with state constraints of the form:

```math
    x_{k+1} = A x_k, x_k ∈ \\mathcal{X}.
```

### Fields

- `A` -- square matrix
- `X` -- state constraints
"""
ConstrainedLinearDiscreteSystem

for (Z, AZ) in ((:ConstrainedAffineContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedAffineDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST} <: $(AZ)
            A::MT
            b::VT
            X::ST
            function $(Z)(A::MT, b::VT, X::ST) where {T, MT <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST}
                @assert checksquare(A) == length(b)
                return new{T, MT, VT, ST}(A, b, X)
            end
        end
        statedim(s::$Z) = length(s.b)
        stateset(s::$Z) = s.X
        inputdim(::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedAffineContinuousSystem

Continuous-time affine system with state constraints of the form:

```math
    x' = A x + b, x(t) ∈ \\mathcal{X}.
```

### Fields

- `A` -- square matrix
- `b` -- vector
- `X` -- state constraints
"""
ConstrainedAffineContinuousSystem

@doc """
    ConstrainedAffineDiscreteSystem

Discrete-time affine system with state constraints of the form:

```math
    x_{k+1} = A x_k + b, x_k ∈ \\mathcal{X} \\text{ for all } k.
```

### Fields

- `A` -- square matrix
- `b` -- vector
- `X` -- state constraints
"""
ConstrainedAffineDiscreteSystem

for (Z, AZ) in ((:ConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST, UT} <: $(AZ)
            A::MTA
            B::MTB
            c::VT
            X::ST
            U::UT
            function $(Z)(A::MTA, B::MTB, c::VT, X::ST, U::UT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, VT <: AbstractVector{T}, ST, UT}
                @assert checksquare(A) == length(c) == size(B, 1)
                return new{T, MTA, MTB, VT, ST, UT}(A, B, c, X, U)
            end
        end
        statedim(s::$Z) = length(s.c)
        stateset(s::$Z) = s.X
        inputdim(s::$Z) = size(s.B, 2)
        inputset(s::$Z) = s.U
        islinear(::$Z) = false
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedAffineControlContinuousSystem

Continuous-time affine control system with state constraints of the form:

```math
    x' = A x + B u + c, x(t) ∈ \\mathcal{X}, u(t) ∈ \\mathcal{U} \\text{ for all } t,
```
and ``c`` a vector.

### Fields

- `A` -- square matrix
- `B` -- matrix
- `c` -- vector
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedAffineControlContinuousSystem

@doc """
    ConstrainedAffineControlDiscreteSystem

Continuous-time affine control system with state constraints of the form:

```math
    x_{k+1} = A x_k + B u_k + c, x_k ∈ \\mathcal{X}, u_k ∈ \\mathcal{U} \\text{ for all } k,
```
and ``c`` a vector.

### Fields

- `A` -- square matrix
- `B` -- matrix
- `c` -- vector
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedAffineControlDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST, UT} <: $(AZ)
            A::MTA
            B::MTB
            X::ST
            U::UT
            function $(Z)(A::MTA, B::MTB, X::ST, U::UT) where {T, MTA <: AbstractMatrix{T}, MTB <: AbstractMatrix{T}, ST, UT}
                @assert checksquare(A) == size(B, 1)
                return new{T, MTA, MTB, ST, UT}(A, B, X, U)
            end
        end
        statedim(s::$Z) = checksquare(s.A)
        stateset(s::$Z) = s.X
        inputdim(s::$Z) = size(s.B, 2)
        inputset(s::$Z) = s.U
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedLinearControlContinuousSystem

Continuous-time linear control system with state constraints of the form:
```math
    x' = A x + B u, x(t) ∈ \\mathcal{X}, u(t) ∈ \\mathcal{U} \\text{ for all } t.
```

### Fields

- `A` -- square matrix
- `B` -- matrix
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedLinearControlContinuousSystem

@doc """
    ConstrainedLinearControlDiscreteSystem

Discrete-time linear control system with state constraints of the form:

```math
    x_{k+1} = A x_k + B u_k, x_k ∈ \\mathcal{X}, u_k ∈ \\mathcal{U} \\text{ for all } k.
```

### Fields

- `A` -- square matrix
- `B` -- matrix
- `X` -- state constraints
- `U` -- input constraints
"""
ConstrainedLinearControlDiscreteSystem

for (Z, AZ) in ((:LinearAlgebraicContinuousSystem, :AbstractContinuousSystem),
                (:LinearAlgebraicDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}} <: $(AZ)
            A::MTA
            E::MTE
            function $(Z)(A::MTA, E::MTE) where {T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}}
                @assert size(A) == size(E)
                return new{T, MTA, MTE}(A, E)
            end
        end
        statedim(s::$Z) = size(s.A, 1)
        inputdim(::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    LinearAlgebraicContinuousSystem

Continuous-time linear algebraic system of the form:

```math
    E x' = A x.
```

### Fields

- `A` -- matrix
- `E` -- matrix, same size as `A`
"""
LinearAlgebraicContinuousSystem

@doc """
    LinearAlgebraicDiscreteSystem

Discrete-time linear algebraic system of the form:

```math
    E x_{k+1} = A x_k.
```

### Fields

- `A` -- matrix
- `E` -- matrix, same size as `A`
"""
LinearAlgebraicDiscreteSystem

for (Z, AZ) in ((:ConstrainedLinearAlgebraicContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedLinearAlgebraicDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, ST} <: $(AZ)
            A::MTA
            E::MTE
            X::ST
            function $(Z)(A::MTA, E::MTE, X::ST) where {T, MTA <: AbstractMatrix{T}, MTE <: AbstractMatrix{T}, ST}
                @assert size(A) == size(E)
                return new{T, MTA, MTE, ST}(A, E, X)
            end
        end
        statedim(s::$Z) = size(s.A, 1)
        stateset(s::$Z) = s.X
        inputdim(::$Z) = 0
        islinear(::$Z) = true
        isaffine(::$Z) = true
    end
end

@doc """
    ConstrainedLinearAlgebraicContinuousSystem

Continuous-time linear system with state constraints of the form:

```math
    E x' = A x, x(t) ∈ \\mathcal{X}.
```

### Fields

- `A` -- matrix
- `E` -- matrix, same size as `A`
- `X` -- state constraints
"""
ConstrainedLinearAlgebraicContinuousSystem

@doc """
    ConstrainedLinearAlgebraicDiscreteSystem

Discrete-time linear system with state constraints of the form:

```math
    E x_{k+1} = A x_k, x_k ∈ \\mathcal{X}.
```

### Fields

- `A` -- matrix
- `E` -- matrix, same size as `A`
- `X` -- state constraints
"""
ConstrainedLinearAlgebraicDiscreteSystem

for (Z, AZ) in ((:PolynomialContinuousSystem, :AbstractContinuousSystem),
                (:PolynomialDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}} <: $(AZ)
            p::VPT
            statedim::Int
            function $(Z)(p::VPT, statedim::Int) where {T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}}
                @assert statedim == MultivariatePolynomials.nvariables(p) "the state dimension $(statedim) does not match the number of state variables"
                return new{T, PT, VPT}(p, statedim)
            end
        end
        statedim(s::$Z) = s.statedim
        inputdim(::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = false

        MultivariatePolynomials.variables(s::$Z) = MultivariatePolynomials.variables(s.p)
        MultivariatePolynomials.nvariables(s::$Z) = s.statedim

        $(Z)(p::AbstractVector{<:AbstractPolynomialLike}) = $(Z)(p, MultivariatePolynomials.nvariables(p))
        $(Z)(p::AbstractPolynomialLike) = $(Z)([p])
    end
end

@doc """
    PolynomialContinuousSystem

Continuous-time polynomial system of the form:

```math
    x' = p(x).
```

### Fields

- `p`        -- polynomial vector field
- `statedim` -- number of state variables
"""
PolynomialContinuousSystem

@doc """
    PolynomialDiscreteSystem

Discrete-time polynomial system of the form:

```math
    x_{k+1} = p(x_k), x_k ∈ \\mathcal{X}.
```

### Fields

- `p`        -- polynomial vector field
- `statedim` -- number of state variables
"""
PolynomialDiscreteSystem

for (Z, AZ) in ((:ConstrainedPolynomialContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedPolynomialDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}, ST} <: $(AZ)
            p::VPT
            statedim::Int
            X::ST
            function $(Z)(p::VPT, statedim::Int, X::ST) where {T, PT <: AbstractPolynomialLike{T}, VPT <: AbstractVector{PT}, ST}
                @assert statedim == MultivariatePolynomials.nvariables(p) "the state dimension $(statedim) does not match the number of state variables"
                return new{T, PT, VPT, ST}(p, statedim, X)
            end
        end
        statedim(s::$Z) = s.statedim
        stateset(s::$Z) = s.X
        inputdim(::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = false

        MultivariatePolynomials.variables(s::$Z) = MultivariatePolynomials.variables(s.p)
        MultivariatePolynomials.nvariables(s::$Z) = s.statedim

        $Z(p::AbstractVector{<:AbstractPolynomialLike}, X::ST) where {ST} = $(Z)(p, MultivariatePolynomials.nvariables(p), X)
        $Z(p::AbstractPolynomialLike, X::ST) where {ST} = $(Z)([p], X)
    end
end

@doc """
    ConstrainedPolynomialContinuousSystem

Continuous-time polynomial system with state constraints:

```math
    x' = p(x), x(t) ∈ \\mathcal{X}
```

### Fields

- `p`        -- polynomial vector field
- `X`        -- constraint set
- `statedim` -- number of state variables
"""
ConstrainedPolynomialContinuousSystem

@doc """
    ConstrainedPolynomialDiscreteSystem

Discrete-time polynomial system with state constraints:

```math
    x_{k+1} = p(x_k), x_k ∈ \\mathcal{X}.
```

### Fields

- `p`        -- polynomial
- `X`        -- constraint set
- `statedim` -- number of state variables
"""
ConstrainedPolynomialDiscreteSystem

for (Z, AZ) in ((:BlackBoxContinuousSystem, :AbstractContinuousSystem),
                (:BlackBoxDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT} <: $(AZ)
            f::FT
            statedim::Int
        end
        statedim(s::$Z) = s.statedim
        inputdim(s::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = false
    end
end

@doc """
    BlackBoxContinuousSystem <: AbstractContinuousSystem

Continuous-time system defined by a right-hand side of the form:

```math
    x' = f(x(t))
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
"""
BlackBoxContinuousSystem

@doc """
    BlackBoxDiscreteSystem <: AbstractDiscreteSystem

Discrete-time system defined by a right-hand side of the form:

```math
    x_{k+1} = f(x_k)
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
"""
BlackBoxDiscreteSystem

for (Z, AZ) in ((:ConstrainedBlackBoxContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedBlackBoxDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT, ST} <: $(AZ)
            f::FT
            statedim::Int
            X::ST
        end
        statedim(s::$Z) = s.statedim
        stateset(s::$Z) = s.X
        inputdim(s::$Z) = 0
        islinear(::$Z) = false
        isaffine(::$Z) = false
    end
end

@doc """
    ConstrainedBlackBoxContinuousSystem <: AbstractContinuousSystem

Continuous-time system defined by a right-hand side with state constraints of the
form:

```math
    x' = f(x(t)), x(t) ∈ \\mathcal{X}.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedBlackBoxContinuousSystem

@doc """
    ConstrainedBlackBoxDiscreteSystem <: AbstractDiscreteSystem

Discrete-time system defined by a right-hand side with state constraints
of the form:

```math
    x_{k+1} = f(x_k), x_k ∈ \\mathcal{X}.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `X`        -- state constraints
"""
ConstrainedBlackBoxDiscreteSystem

for (Z, AZ) in ((:ConstrainedBlackBoxControlContinuousSystem, :AbstractContinuousSystem),
                (:ConstrainedBlackBoxControlDiscreteSystem, :AbstractDiscreteSystem))
    @eval begin
        struct $(Z){FT, ST, UT} <: $(AZ)
            f::FT
            statedim::Int
            inputdim::Int
            X::ST
            U::UT
        end
        statedim(s::$Z) = s.statedim
        stateset(s::$Z) = s.X
        inputdim(s::$Z) = s.inputdim
        inputset(s::$Z) = s.U
        islinear(::$Z) = false
        isaffine(::$Z) = false
    end
end

@doc """
    ConstrainedBlackBoxControlContinuousSystem <: AbstractContinuousSystem

Continuous-time control system defined by a right-hand side with state constraints
of the form:

```math
    x' = f(x(t), u(t)), x(t) ∈ \\mathcal{X}, u(t) ∈ \\mathcal{U}.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `X`        -- state constraints
- `U`        -- input constraints
"""
ConstrainedBlackBoxControlContinuousSystem

@doc """
    ConstrainedBlackBoxControlDiscreteSystem <: AbstractDiscreteSystem

Discrete-time control system defined by a right-hand side with state constraints
of the form:

```math
    x_{k+1} = f(x_k), x_k ∈ \\mathcal{X}, u_k ∈ \\mathcal{U}.
```

### Fields

- `f`        -- function that holds the right-hand side
- `statedim` -- number of state variables
- `inputdim` -- number of input variables
- `X`        -- state constraints
- `U`        -- input constraints
"""
ConstrainedBlackBoxControlDiscreteSystem


text:
using NeutralLandscapes, Test

algorithms = (
   DiamondSquare(),
   DiamondSquare(; H=0.2),
   DiamondSquare(0.4),
   DiscreteVoronoi(),
   DiscreteVoronoi(; n=2),
   DiscreteVoronoi(3),
   DistanceGradient(),
   DistanceGradient(; sources=[2]),
   DistanceGradient([1, 2]),
   EdgeGradient(),
   EdgeGradient(; direction=120.0),
   EdgeGradient(90.0),
   MidpointDisplacement(),
   MidpointDisplacement(;),
   NearestNeighborElement(),
   NearestNeighborElement(; k=1),
   NearestNeighborElement(2),
   NearestNeighborCluster(),
   # NearestNeighborCluster(; n=:queen),
   NearestNeighborCluster(; n=:diagonal),
   NearestNeighborCluster(0.2, :diagonal),
   NoGradient(),
   PerlinNoise(),
   PerlinNoise(; octaves=3, lacunarity=3, persistance=0.3),
   PerlinNoise((1, 2), 3, 2, 0.2),
   PlanarGradient(),
   PlanarGradient(; direction=120.0),
   PlanarGradient(90),
   RectangularCluster(),
   RectangularCluster(; maximum=5),
   RectangularCluster(2),
   WaveSurface(),
   WaveSurface(90.0),
   WaveSurface(; periods=2),
)

sizes = (50, 200), (100, 100), (301, 87)

for alg in algorithms, sze in sizes
    A = rand(alg, sze)
    @test size(A) == sze
end


text:
module ActiveAppearanceModels

export AAModel,
       train,
       fit,
       Shape,  # convenient alias
       viewshape

include("core.jl")

end


text:
# Reproject Bearing Range factors to see what is happening.

struct TagkTl
	tagID::Symbol
	kTl::Vector{Float64}
end

"""
    $(SIGNATURES)

Draw tags on pose.
""" #TODO: confirm xTc or cTx
function visTagsOnPose!(vis::Visualizer,
                        tagsOnPoses::Dict{Symbol,Vector{TagkTl}},
                        sessionId::Symbol )
    #
	lmpoint = HyperSphere(Point(0.,0,0), 0.05)
	blueMat = MeshPhongMaterial(color=RGBA(0, 0, 1, 0.5))

	for x = keys(tagsOnPoses)
	    for tag in tagsOnPoses[x]
			setobject!(botvis.vis[sessionId][:poses][x][tag.tagID], lmpoint, blueMat)
			trans = CTs.Translation(tag.kTl)
			settransform!(botvis.vis[sessionId][:poses][x][tag.tagID], trans)
	    end
	end
end

"""
    $(SIGNATURES)

Main plugin callback function
"""
function reprojectBearingRange(vis::MeshCat.Visualizer,
                               params::Dict{Symbol, Any},
                               rose_fgl)
    #

    # create the necessary cache Dict
    tagsOnPoses = Dict{Symbol,Vector{TagkTl}}()
    if !haskey(params, :tagsOnPoses)
        # point to the required memory
        params[:tagsOnPoses] = tagsOnPoses
    else
        # convenience variable avoiding Dict look-up in later loops
        tagsOnPoses = params[:tagsOnPoses]
    end

    # populate / update the tag cache
    for (xx, val) in params[:cachevars]
        # add projections only once and ensure pose is already drawn in viewer
        if string(xx)[1] == 'x' && val[2][1] && !haskey(tagsOnPose, xx)
            # @show xx, val[1] == :Pose2  # assume this is a pose

            # get tags and measurements
            poseTags = getBearingRangesOnPose(rose_fgl, xx)

            # cache the tags to be projected from this pose
            # if haskey(tagsOnPoses, xx)
            #     push!(tagsOnPoses[xx], TagkTl(???) )
            # else
            #     tagsOnPoses[xx] = TagkTl[TagkTl(???);]
            #     @info "reprojectBearingRange plugin needs to update an already existing tagOnPose ???"
            # end
        end
    end

    # visualize the contents of the cache
    visTagsOnPose!(vis, params[:tagsOnPoses], params[:sessionId])

    nothing
end