text:
# Gurobi model attributes

############################################
#
#   Low level attribute getter/setters
#
############################################

function get_intattr(model::Model, name::String)
    @assert isascii(name)
    a = Ref{Cint}()
    ret = @grb_ccall(getintattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ref{Cint}),
        model, name, a);
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    Int(a[])
end

function get_dblattr(model::Model, name::String)
    @assert isascii(name)
    a = Ref{Float64}()
    ret = @grb_ccall(getdblattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ref{Float64}),
        model, name, a);
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    a[]
end

function get_strattr(model::Model, name::String)
    @assert isascii(name)
    a = Ref{Ptr{UInt8}}()
    ret = @grb_ccall(getstrattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ref{Ptr{UInt8}}),
        model, name, a)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    unsafe_string(a[])
end

# array element

function get_intattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Ref{Cint}()
    ret = @grb_ccall(getintattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ref{Cint}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    Int(a[])
end

function get_dblattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Ref{Float64}()
    ret = @grb_ccall(getdblattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ref{Float64}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    a[]
end

function get_charattrelement(model::Gurobi.Model, name::String, element::Int)
    @assert isascii(name)
    a = Ref{Cchar}()
    ret = @grb_ccall(getcharattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ref{Cchar}),
        model, name, element - 1, a
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    Char(a[])
end

# Note: in attrarray API, the start argument is one-based (following Julia convention)

function get_intattrarray!(r::Array{Cint}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getintattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Cint}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    r
end

function get_intattrlist!(r::Array{Cint}, model::Model, name::String, inds::Vector{I}) where I<:Integer
    @assert isascii(name)
    @assert length(r) == length(inds)
    ret = @grb_ccall(getintattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Cint}),
        model, name, Cint(length(inds)), ivec(inds .- 1), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function get_intattrlist(model::Model, name::String, inds::Vector{I}) where I<:Integer
    r = Array{Cint}(undef, length(inds))
    get_intattrlist!(r::Array{Cint}, model::Model, name::String, inds)
    return r
end

function get_intattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_intattrarray!(Array{Cint}(len), model, name, start)
end

function get_dblattrarray!(r::Array{Float64}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getdblattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Float64}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    r
end

function get_dblattrlist!(r::Array{Float64}, model::Model, name::String, inds::Vector{I}) where I<:Integer
    @assert isascii(name)
    @assert length(r) == length(inds)
    ret = @grb_ccall(getdblattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Float64}),
        model, name, Cint(length(inds)), ivec(inds .- 1), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function get_dblattrlist(model::Model, name::String, inds::Vector{I}) where I<:Integer
    r = Array{Float64}(undef, length(inds))
    get_dblattrlist!(r::Array{Float64}, model::Model, name::String, inds)
    return r
end

function get_dblattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_dblattrarray!(Array{Float64}(undef, len), model, name, start)
end

function get_charattrarray!(r::Array{Cchar}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getcharattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Cchar}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    map(Char,r)
end

function get_charattrlist!(r::Array{Cchar}, model::Model, name::String, inds::Vector{I}) where I<:Integer
    @assert isascii(name)
    @assert length(r) == length(inds)
    ret = @grb_ccall(getcharattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Cchar}),
        model, name, Cint(length(inds)), ivec(inds .- 1), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function get_charattrlist(model::Model, name::String, inds::Vector{I}) where I<:Integer
    r = Array{Cchar}(undef, length(inds))
    get_charattrlist!(r::Array{Cchar}, model::Model, name::String, inds)
    return r
end

function get_charattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_charattrarray!(Array{Cchar}(undef, len), model, name, start)
end

function get_strattrarray(model::Model, name::String, start::Integer, len::Integer)
    @assert isascii(name)
    get_strattrarray!(Array{Ptr{UInt8}}(undef, len), model, name, start)
end

function get_strattrarray!(r::Array{Ptr{UInt8}}, model::Model, name::String, start::Integer)
    @assert isascii(name)
    ret = @grb_ccall(getstrattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Ptr{UInt8}}),
        model, name, start - 1, length(r), r)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    map(unsafe_string, r)
end

# setters

function set_intattr!(model::Model, name::String, v::Integer)
    @assert isascii(name)
    ret = @grb_ccall(setintattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattr!(model::Model, name::String, v::Real)
    @assert isascii(name)
    ret = @grb_ccall(setdblattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Float64), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_strattr!(model::Model, name::String, v::String)
    @assert isascii(name)
    @assert isascii(v)
    ret = @grb_ccall(setstrattr, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}), model, name, v)
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# array element

function set_intattrelement!(model::Gurobi.Model, name::String, element::Int, v::Integer)
    @assert isascii(name)
    ret = @grb_ccall(setintattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattrelement!(model::Gurobi.Model, name::String, element::Int, v::Real)
    @assert isascii(name)
    ret = @grb_ccall(setdblattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Float64),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_charattrelement!(model::Gurobi.Model, name::String, element::Int, v::Char)
    @assert isascii(name)
    ret = @grb_ccall(setcharattrelement, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cchar),
        model, name, element - 1, v
    )
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

# reminder: start is one-based

function set_intattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Cint},values)
    ret = @grb_ccall(setintattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Cint}), model, name, start-1, len, ivec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_intattrlist!(model::Model, name::String, inds::Vector, values::Vector)
    @assert isascii(name)
    @assert length(inds) == length(values)
    ret = @grb_ccall(setintattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Cint}), model, name, Cint(length(inds)), ivec(inds.-1), ivec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Float64},values)
    ret = @grb_ccall(setdblattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Float64}), model, name, start-1, len, fvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_dblattrlist!(model::Model, name::String, inds::Vector, values::Vector)
    @assert isascii(name)
    @assert length(inds) == length(values)
    ret = @grb_ccall(setdblattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Float64}), model, name, Cint(length(inds)), ivec(inds.-1), fvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_charattrarray!(model::Model, name::String, start::Integer, len::Integer, values::Vector)
    @assert isascii(name)
    values = convert(Vector{Cchar},values)
    ret = @grb_ccall(setcharattrarray, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Cint, Ptr{Cchar}), model, name, start-1, len, cvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

function set_charattrlist!(model::Model, name::String, inds::Vector, values::Vector)
    @assert isascii(name)
    @assert length(inds) == length(values)
    ret = @grb_ccall(setcharattrlist, Cint,
        (Ptr{Cvoid}, Ptr{UInt8}, Cint, Ptr{Cint}, Ptr{Cchar}), model, name, Cint(length(inds)), ivec(inds.-1), cvec(values))
    if ret != 0
        throw(GurobiError(model.env, ret))
    end
    nothing
end

############################################
#
#   Macros for array definition
#
############################################

macro grb_int_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_intattr(model, $(attrname))
end

macro grb_dbl_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_dblattr(model, $(attrname))
end

macro grb_str_attr(fun, attrname)
    @eval $(fun)(model::Model) = get_strattr(model, $(attrname))
end

############################################
#
#   Model attributes
#
############################################

# basic attributes

@grb_str_attr model_name  "ModelName"

@grb_int_attr num_vars     "NumVars"
@grb_int_attr num_constrs  "NumConstrs"
@grb_int_attr num_sos      "NumSOS"
@grb_int_attr num_qconstrs "NumQConstrs"
@grb_int_attr num_cnzs     "NumNZs"
@grb_int_attr num_qnzs     "NumQNZs"
@grb_int_attr num_qcnzs    "NumQCNZs"

@grb_int_attr num_intvars  "NumIntVars"
@grb_int_attr num_binvars  "NumBinVars"

# derived attribute functions

model_sense(model::Model) = get_intattr(model, "ModelSense") > 0 ? (:minimize) : (:maximize)

is_mip(model::Model) = get_intattr(model, "IsMIP") != 0
is_qp(model::Model)  = get_intattr(model, "IsQP") != 0
is_qcp(model::Model) = get_intattr(model, "IsQCP") != 0

function model_type(model::Model)
    is_qp(model)  ? (:QP)  :
    is_qcp(model) ? (:QCP) : (:LP)
end

function set_sense!(model::Model, sense::Symbol)
    v = sense == :maximize ? -1 :
        sense == :minimize ? 1 :
        throw(ArgumentError("Invalid model sense."))

    set_intattr!(model, "ModelSense", v)
end

# variable related attributes

lowerbounds(model::Model) = get_dblattrarray(model, "LB", 1, num_vars(model))
upperbounds(model::Model) = get_dblattrarray(model, "UB", 1, num_vars(model))
objcoeffs(model::Model) = get_dblattrarray(model, "Obj", 1, num_vars(model))

# note: this takes effect only after update_model! is called:
function set_objcoeffs!(model::Model, c::Vector)
    n = num_vars(model)
    length(c) == n || error("Inconsistent argument dimensions.")
    set_dblattrarray!(model, "Obj", 1, n, c)
end

############################################
#
#   The show method for model
#
#   - Based on attributes
#
############################################

function show(io::IO, model::Model)
    if model.ptr_model != C_NULL
        println(io, "Gurobi Model: $(model_name(model))")
        if is_mip(model)
            println(io, "    type   : $(model_type(model)) (MIP)")
        else
            println(io, "    type   : $(model_type(model))")
        end
        println(io, "    sense  : $(model_sense(model))")
        println(io, "    number of variables             = $(num_vars(model))")
        println(io, "    number of linear constraints    = $(num_constrs(model))")
        println(io, "    number of quadratic constraints = $(num_qconstrs(model))")
        println(io, "    number of sos constraints       = $(num_sos(model))")
        println(io, "    number of non-zero coeffs       = $(num_cnzs(model))")
        println(io, "    number of non-zero qp objective terms  = $(num_qnzs(model))")
        println(io, "    number of non-zero qp constraint terms = $(num_qcnzs(model))")
    else
        println(io, "Gurobi Model: NULL")
    end
end


text:
module DocumentFormat
using CSTParser
using CSTParser.Tokenize
import CSTParser.Tokenize.Tokens

mutable struct FormatOptions
    indent::Int
    indents::Bool
    ops::Bool
    tuples::Bool
    curly::Bool
    calls::Bool
    iterOps::Bool
    comments::Bool
    docs::Bool
    lineends::Bool
end
FormatOptions() = FormatOptions(4, true, true, true, true, true, true, true, true, true)

struct Edit{T}
    loc::T
    text::String
end

mutable struct State{T}
    offset::Int
    edits::T
    opts::FormatOptions
end

function format(text, formatopts::FormatOptions = FormatOptions())
    original_ast = CSTParser.remlineinfo!(Meta.parse(string("begin\n", text, "\nend")))
    state = State(0, Edit[], formatopts)
    x = CSTParser.parse(text, true)
    if formatopts.ops
        pass(x, state, operator_pass)
    end
    if formatopts.tuples
        state.offset = 0
        pass(x, state, tuple_pass)
    end
    if formatopts.curly
        state.offset = 0
        pass(x, state, curly_pass)
    end
    if formatopts.calls
        state.offset = 0
        pass(x, state, call_pass)
    end
    if formatopts.iterOps
        state.offset = 0
        pass(x, state, forloop_pass)
    end
    if formatopts.comments
        comments_pass(text, state)
    end
    if formatopts.docs
        state.offset = 0
        pass(x, state, doc_pass)
    end
    if formatopts.lineends
        lineends_pass(text, state)
    end
    sort!(state.edits, lt = (a, b)->first(a.loc) < first(b.loc), rev = true)

    for i = 1:length(state.edits)
        text = apply(text, state.edits[i])
    end
    if formatopts.indents
        text = indents(text, state.opts)
    end
    new_ast = CSTParser.remlineinfo!(Meta.parse(string("begin\n", text, "\nend")))
    original_ast != new_ast && error("Mismatch between AST of original and formatted text")
    return text
end

function pass(x, state, f = (x, state)->nothing)
    f(x, state)
    if x.args isa Vector
        for a in x.args
            pass(a, state, f)
        end
    else
        state.offset += x.fullspan
    end
    state
end

function ensure_single_space_after(x, state, offset)
    if x.fullspan == x.span
        push!(state.edits, Edit(offset + x.fullspan, " "))

    end
end

function ensure_no_space_after(x, state, offset)
    if x.fullspan != x.span
        push!(state.edits, Edit(offset .+ (x.span + 1:x.fullspan), ""))
    end
end

function apply(text, edit::Edit{Int})
    v = Vector{UInt8}(deepcopy(text))
    String(vcat(v[1:edit.loc], Vector{UInt8}(edit.text), v[edit.loc + 1:end]))
end
function apply(text, edit::Edit{UnitRange{Int}})
    v = Vector{UInt8}(deepcopy(text))
    String(vcat(v[1:first(edit.loc) - 1], Vector{UInt8}(edit.text), v[last(edit.loc) + 1:end]))
end
include("passes.jl")
include("indents.jl")
end

ext: jl


text:
module FrostHelperNoDashArea

using ..Ahorn, Maple

@mapdef Entity "FrostHelper/NoDashArea" NoDashAreaFrostHelper(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, fastMoving::Bool=false)

const placements = Ahorn.PlacementDict(
    "No Dash Area (Frost Helper)" => Ahorn.EntityPlacement(
        NoDashAreaFrostHelper,
        "rectangle"
    ),
)

Ahorn.nodeLimits(entity::NoDashAreaFrostHelper) = 0, 1
Ahorn.minimumSize(entity::NoDashAreaFrostHelper) = 8, 8
Ahorn.resizable(entity::NoDashAreaFrostHelper) = true, true

function Ahorn.selection(entity::NoDashAreaFrostHelper)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    nodes = get(entity.data, "nodes", ())
    if isempty(nodes)
        return Ahorn.Rectangle(x, y, width, height)

    else
        nx, ny = Int.(nodes[1])
        return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx, ny, width, height)]
    end
end

edgeColor = (38, 0, 0, 255) ./ 255
centerColor = (64, 0, 0, 102) ./ 255

function renderSpaceJam(ctx::Ahorn.Cairo.CairoContext, x::Number, y::Number, width::Number, height::Number)
    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)

    Ahorn.drawRectangle(ctx, x, y, width, height, edgeColor, centerColor)

    Ahorn.restore(ctx)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::NoDashAreaFrostHelper)
    x, y = Ahorn.position(entity)
    nodes = get(entity.data, "nodes", ())

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    if !isempty(nodes)
        nx, ny = Int.(nodes[1])

        cox, coy = floor(Int, width / 2), floor(Int, height / 2)

        renderSpaceJam(ctx, nx, ny, width, height)
        Ahorn.drawArrow(ctx, x + cox, y + coy, nx + cox, ny + coy, Ahorn.colors.selection_selected_fc, headLength=6)
    end
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::NoDashAreaFrostHelper, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderSpaceJam(ctx, 0, 0, width, height)
end

end

max_stars_repo_name: JaThePlayer/FrostHelper

max_issues_repo_name: JaThePlayer/FrostHelper

max_forks_repo_name: JaThePlayer/FrostHelper


text:
function bifurcate_unassigned_scheduler(env)
    # Split a problem into two sub problems to increase speed.
    # Hypthetically could be called recursively.
    child_ceil = deepcopy(env)
    child_floor = deepcopy(env)
    for day in 1:length(env["coverage"])
        for t in 1:length(env["coverage"][day])
            child_ceil["coverage"][day][t] = ceil(env["coverage"][day][t] / 2)

            child_floor["coverage"][day][t] = floor(env["coverage"][day][t] / 2)
            if child_floor["coverage"][day][t] == 0 && env["coverage"][day][t] != 0
                child_floor["coverage"][day][t] = 1
            end
        end
    end

    # Run separately
    Logging.info("Starting ceiling bifurcation child")
    ok, shifts_ceil = run_schedule_unassigned(child_ceil)
    if !ok
        return ok, Dict()
    end

    Logging.info( "Starting floor bifurcation child")
    ok, shifts_floor = run_schedule_unassigned(child_floor)
    if !ok
        return ok, Dict()
    end

    # Merge shifts
    output_shifts = deepcopy(shifts_ceil)
    for name in keys(shifts_floor)
        if !(name in keys(output_shifts))
            output_shifts[name] = Any[]
        end
        for shift in shifts_floor[name]
            push!(output_shifts[name], shift)
        end
    end
    return true, output_shifts
end

max_stars_repo_name: shdann/autoshift

max_issues_repo_name: shdann/autoshift

max_forks_repo_name: shdann/autoshift


text:
abstract type AbstractBasis{dim,T<:Number} <: AbstractMatrix{T} end

@pure Base.size(::AbstractBasis{dim}) where {dim} = (dim, dim)
Base.getindex(ℬ::AbstractBasis, i::Integer, j::Integer) = getindex(vecbasis(ℬ, :cov), i, j)
@pure Base.eltype(::Type{AbstractBasis{dim,T}}) where {dim,T} = T
@pure getdim(::AbstractBasis{dim}) where {dim} = dim

function subscriptnumber(i::Integer)
    if i < 0
        c = [Char(0x208B)]
    else
        c = []
    end
    for d in reverse(digits(abs(i)))
        push!(c, Char(0x2080+d))
    end
    return join(c)
end

function superscriptnumber(i::Integer)
    if i < 0
        c = [Char(0x207B)]
    else
        c = []
    end
    for d in reverse(digits(abs(i)))
        if d == 0 push!(c, Char(0x2070)) end
        if d == 1 push!(c, Char(0x00B9)) end
        if d == 2 push!(c, Char(0x00B2)) end
        if d == 3 push!(c, Char(0x00B3)) end
        if d > 3 push!(c, Char(0x2070+d)) end
    end
    return join(c)
end

"""
    Basis(v::AbstractMatrix{T}, ::Val{:cov})
    Basis{dim, T<:Number}()
    Basis(θ::T<:Number, ϕ::T<:Number, ψ::T<:Number)

Basis built from a square matrix `v` where columns correspond either to
- primal vectors ie `eᵢ=v[:,i]` if `var=:cov` as by default
- dual vectors ie `eⁱ=v[:,i]` if `var=:cont`.

Basis without any argument refers to the canonical basis (`CanonicalBasis`) in `Rᵈⁱᵐ` (by default `dim=3` and `T=Sym`)

Basis can also be built from Euler angles (`RotatedBasis`) `θ` in 2D and `(θ, ϕ, ψ)` in 3D

The attributes of this object can be obtained by
- `vecbasis(ℬ, :cov)`: square matrix defining the primal basis `eᵢ=e[:,i]`
- `vecbasis(ℬ, :cont)`: square matrix defining the dual basis `eⁱ=E[:,i]`
- `metric(ℬ, :cov)`: square matrix defining the covariant components of the metric tensor `gᵢⱼ=eᵢ⋅eⱼ=g[i,j]`
- `metric(ℬ, :cont)`: square matrix defining the contravariant components of the metric tensor `gⁱʲ=eⁱ⋅eʲ=gⁱʲ[i,j]`

# Examples
```julia
julia> v = Sym[1 0 0; 0 1 0; 0 1 1] ; ℬ = Basis(v)
Basis{3, Sym}
# basis: 3×3 Tensor{2, 3, Sym, 9}:
 1  0  0
 0  1  0
 0  1  1
# dual basis: 3×3 Tensor{2, 3, Sym, 9}:
 1  0   0
 0  1  -1
 0  0   1
# covariant metric tensor: 3×3 SymmetricTensor{2, 3, Sym, 6}:
 1  0  0
 0  2  1
 0  1  1
# contravariant metric tensor: 3×3 SymmetricTensor{2, 3, Sym, 6}:
 1   0   0
 0   1  -1
 0  -1   2

julia> θ, ϕ, ψ = symbols("θ, ϕ, ψ", real = true) ; ℬʳ = Basis(θ, ϕ, ψ) ; display(vecbasis(ℬʳ, :cov))
3×3 Tensor{2, 3, Sym, 9}:
 -sin(ψ)⋅sin(ϕ) + cos(θ)⋅cos(ψ)⋅cos(ϕ)  -sin(ψ)⋅cos(θ)⋅cos(ϕ) - sin(ϕ)⋅cos(ψ)  sin(θ)⋅cos(ϕ)
  sin(ψ)⋅cos(ϕ) + sin(ϕ)⋅cos(θ)⋅cos(ψ)  -sin(ψ)⋅sin(ϕ)⋅cos(θ) + cos(ψ)⋅cos(ϕ)  sin(θ)⋅sin(ϕ)
                        -sin(θ)⋅cos(ψ)                          sin(θ)⋅sin(ψ)         cos(θ)
```
"""
struct Basis{dim,T} <: AbstractBasis{dim,T}
    eᵢ::Matrix{T} # Primal basis `eᵢ=eᵢ[:,i]`
    eⁱ::Matrix{T} # Dual basis `eⁱ=eⁱ[:,i]`
    gᵢⱼ::Symmetric{T,Matrix{T}} # Metric tensor `gᵢⱼ=eᵢ⋅eⱼ=gᵢⱼ[i,j]`
    gⁱʲ::Symmetric{T,Matrix{T}} # Inverse of the metric tensor `gⁱʲ=eⁱ⋅eʲ=gⁱʲ[i,j]`
    function Basis(
        eᵢ::AbstractMatrix{T},
        eⁱ::AbstractMatrix{T},
        gᵢⱼ::AbstractMatrix{T},
        gⁱʲ::AbstractMatrix{T},
    ) where {T}
        dim = size(eᵢ, 1)
        @assert dim == size(eᵢ, 2) "v should be a square matrix"
        if isidentity(eᵢ)
            return CanonicalBasis{dim,T}()
        elseif isidentity(gᵢⱼ)
            return RotatedBasis(eᵢ)
        elseif isdiag(gᵢⱼ)
            χ = sqrt.(diag(gᵢⱼ))
            return OrthogonalBasis(RotatedBasis(eᵢ .* transpose(inv.(χ))), χ)
        else
            eᵢ = Matrix(eᵢ)
            gᵢⱼ = Symmetric(gᵢⱼ)
            gⁱʲ = Symmetric(gⁱʲ)
            eⁱ = Matrix(eⁱ)
            new{dim,T}(eᵢ, eⁱ, gᵢⱼ, gⁱʲ)
        end
    end
    function Basis(ℬ::AbstractBasis{dim,T}, χᵢ::V) where {dim,T,V}
        Χ = collect(χᵢ)
        invΧ = inv.(Χ)
        if Χ == one.(Χ)
            return ℬ
        else
            if ℬ isa OrthonormalBasis
                return OrthogonalBasis(ℬ, Χ)
            else
                eᵢ = vecbasis(ℬ, :cov) .* transpose(Χ)
                eⁱ = vecbasis(ℬ, :cont) .* transpose(invΧ)
                gᵢⱼ = Symmetric(Χ .* metric(ℬ, :cov) .* transpose(Χ))
                gⁱʲ = Symmetric(invΧ .* metric(ℬ, :cont) .* transpose(invΧ))
                new{dim,T}(eᵢ, eⁱ, gᵢⱼ, gⁱʲ)
            end
        end
    end
    function Basis(eᵢ::AbstractMatrix{T}, ::Val{:cov}) where {T}
        dim = size(eᵢ, 1)
        @assert dim == size(eᵢ, 2) "v should be a square matrix"
        if isidentity(eᵢ)
            return CanonicalBasis{dim,T}()
        else
            eᵢ = Matrix(eᵢ)
            gᵢⱼ = simplifyif(Symmetric(transpose(eᵢ) * eᵢ))
            if isidentity(gᵢⱼ)
                return RotatedBasis(eᵢ)
            elseif isdiagonal(gᵢⱼ)
                χ = sqrt.(diag(gᵢⱼ))
                return OrthogonalBasis(RotatedBasis(eᵢ .* transpose(inv.(χ))), χ)
            else
                if T == Sym
                    gⁱʲ = simplifyif(Symmetric(inv(Matrix(gᵢⱼ))))
                else
                    gⁱʲ = simplifyif(inv(gᵢⱼ))
                end
                eⁱ = simplifyif(eᵢ * transpose(gⁱʲ))
                new{dim,T}(eᵢ, eⁱ, gᵢⱼ, gⁱʲ)
            end
        end
    end
    function Basis(eⁱ::AbstractMatrix{T}, ::Val{:cont}) where {T}
        dim = size(eⁱ, 1)
        @assert dim == size(eⁱ, 2) "v should be a square matrix"
        if isidentity(eⁱ)
            return CanonicalBasis{dim,T}()
        else
            eⁱ = Matrix(eⁱ)
            gⁱʲ = simplifyif(Symmetric(transpose(eⁱ) * eⁱ))
            if isidentity(gⁱʲ)
                return RotatedBasis(eⁱ)
            elseif isdiagonal(gⁱʲ)
                uχ = inv.(sqrt.(diag(gⁱʲ)))
                return OrthogonalBasis(RotatedBasis(eⁱ .* transpose(uχ)), uχ)
            else
                if T == Sym
                    gᵢⱼ = simplifyif(Symmetric(inv(Matrix(gⁱʲ))))
                else
                    gᵢⱼ = simplifyif(inv(gⁱʲ))
                end
                eᵢ = simplifyif(eⁱ * transpose(gᵢⱼ))
                new{dim,T}(eᵢ, eⁱ, gᵢⱼ, gⁱʲ)
            end
        end
    end
    Basis(v::AbstractMatrix{T}, var::Symbol) where {T} = Basis(v, Val(var))
    Basis(v::AbstractMatrix{T}) where {T} = Basis(v, :cov)
    Basis(θ::T1, ϕ::T2, ψ::T3 = 0) where {T1,T2,T3} = RotatedBasis(θ, ϕ, ψ)
    Basis(θ::T) where {T} = RotatedBasis(θ)
    Basis{dim,T}() where {dim,T} = CanonicalBasis{dim,T}()
    Basis() = CanonicalBasis()
end

"""
    CanonicalBasis{dim, T}

Canonical basis of dimension `dim` (default: 3) and type `T` (default: Sym)

The attributes of this object can be obtained by
- `vecbasis(ℬ, :cov)`: square matrix defining the primal basis `eᵢ=e[:,i]=δᵢⱼ`
- `vecbasis(ℬ, :cont)`: square matrix defining the dual basis `eⁱ=E[:,i]=δᵢⱼ`
- `metric(ℬ, :cov)`: square matrix defining the covariant components of the metric tensor `gᵢⱼ=eᵢ⋅eⱼ=g[i,j]=δᵢⱼ`
- `metric(ℬ, :cont)`: square matrix defining the contravariant components of the metric tensor `gⁱʲ=eⁱ⋅eʲ=gⁱʲ[i,j]=δᵢⱼ`

# Examples
```julia
julia> ℬ = CanonicalBasis()
CanonicalBasis{3, Sym}
# basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# dual basis: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# covariant metric tensor: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1
# contravariant metric tensor: 3×3 TensND.LazyIdentity{3, Sym}:
 1  0  0
 0  1  0
 0  0  1

julia> ℬ₂ = CanonicalBasis{2, Float64}()
CanonicalBasis{2, Float64}
# basis: 2×2 TensND.LazyIdentity{2, Float64}:
 1.0  0.0
 0.0  1.0
# dual basis: 2×2 TensND.LazyIdentity{2, Float64}:
 1.0  0.0
 0.0  1.0
# covariant metric tensor: 2×2 TensND.LazyIdentity{2, Float64}:
 1.0  0.0
 0.0  1.0
# contravariant metric tensor: 2×2 TensND.LazyIdentity{2, Float64}:
 1.0  0.0
 0.0  1.0
```
"""
struct CanonicalBasis{dim,T} <: AbstractBasis{dim,T} end
CanonicalBasis() = CanonicalBasis{3,Sym}()

"""
    RotatedBasis(θ::T<:Number, ϕ::T<:Number, ψ::T<:Number)
    RotatedBasis(θ::T<:Number)

Orthonormal basis of dimension `dim` (default: 3) and type `T` (default: Sym) built from Euler angles `θ` in 2D and `(θ, ϕ, ψ)` in 3D

# Examples
```julia
julia> θ, ϕ, ψ = symbols("θ, ϕ, ψ", real = true) ; ℬʳ = RotatedBasis(θ, ϕ, ψ) ; display(vecbasis(ℬʳ, :cov))
3×3 Tensor{2, 3, Sym, 9}:
 -sin(ψ)⋅sin(ϕ) + cos(θ)⋅cos(ψ)⋅cos(ϕ)  -sin(ψ)⋅cos(θ)⋅cos(ϕ) - sin(ϕ)⋅cos(ψ)  sin(θ)⋅cos(ϕ)
  sin(ψ)⋅cos(ϕ) + sin(ϕ)⋅cos(θ)⋅cos(ψ)  -sin(ψ)⋅sin(ϕ)⋅cos(θ) + cos(ψ)⋅cos(ϕ)  sin(θ)⋅sin(ϕ)
                        -sin(θ)⋅cos(ψ)                          sin(θ)⋅sin(ψ)         cos(θ)
```
"""
struct RotatedBasis{dim,T} <: AbstractBasis{dim,T}
    eᵢ::Matrix{T} # Primal basis `eᵢ=e[:,i]`
    eⁱ::Matrix{T} # Dual basis `eⁱ=E[:,i]`
    angles::NamedTuple
    function RotatedBasis(R::AbstractMatrix{T}) where {T<:Number}
        dim = size(R, 1)
        if isidentity(R)
            return CanonicalBasis{dim,T}()
        else
            eᵢ = eⁱ = Matrix(R)
            return new{dim,T}(eᵢ, eⁱ, angles(R))
        end
    end
    function RotatedBasis(θ::T1, ϕ::T2, ψ::T3 = 0) where {T1<:Number,T2<:Number,T3<:Number}
        T = promote_type(T1, T2, T3)
        dim = 3
        R = RotZYZ(ϕ, θ, ψ)
        if isidentity(R)
            return CanonicalBasis{dim,T}()
        else
            eᵢ = eⁱ = Matrix(R)
            return new{dim,T}(eᵢ, eⁱ, angles(R))
        end
    end
    function RotatedBasis(θ::T) where {T<:Number}
        dim = 2
        cθ = cos(θ)
        sθ = sin(θ)
        eᵢ = eⁱ = [cθ -sθ; sθ cθ]
        if isidentity(eᵢ)
            return CanonicalBasis{dim,T}()
        else
            return new{dim,T}(eᵢ, eⁱ, angles(eᵢ))
        end
    end
    function RotatedBasis(θ::Sym)
        dim = 2
        cθ = cos(θ)
        sθ = sin(θ)
        eᵢ = eⁱ = [cθ -sθ; sθ cθ]
        if isidentity(eᵢ)
            return CanonicalBasis{dim,Sym}()
        else
            return new{dim,Sym}(eᵢ, eⁱ, (θ = θ,))
        end
    end
end

const OrthonormalBasis{dim,T} = Union{CanonicalBasis{dim,T},RotatedBasis{dim,T}}

struct OrthogonalBasis{dim,T} <: AbstractBasis{dim,T}
    parent::OrthonormalBasis{dim,T}
    λ::Vector{T}
    eᵢ::Matrix{T}
    eⁱ::Matrix{T}
    gᵢⱼ::Diagonal{T,Vector{T}}
    gⁱʲ::Diagonal{T,Vector{T}}
    function OrthogonalBasis(parent::OrthonormalBasis{dim,T}, λ::Vector) where {dim,T}
        λ = T.(λ)
        if λ == one.(λ)
            return parent
        else
            eᵢ = [λ[j] * parent[i, j] for i ∈ 1:dim, j ∈ 1:dim]
            eⁱ = [parent[i, j] / λ[j] for i ∈ 1:dim, j ∈ 1:dim]
            return new{dim,T}(parent, λ, eᵢ, eⁱ, Diagonal(λ .^ 2), Diagonal(inv.(λ) .^ 2))
        end
    end
end

relevant_OrthonormalBasis(ℬ::OrthogonalBasis) = ℬ.parent
relevant_OrthonormalBasis(ℬ::OrthonormalBasis) = ℬ
relevant_OrthonormalBasis(::Basis{dim,T}) where {dim,T} = CanonicalBasis{dim,T}()

@inline CylindricalBasis(θ) = RotatedBasis(0, θ, 0)

@inline SphericalBasis(θ, ϕ) = RotatedBasis(θ, ϕ, 0)

const AllOrthogonalBasis{dim,T} = Union{OrthonormalBasis{dim,T},OrthogonalBasis{dim,T}}

angles(M::AbstractMatrix{T}, ::Val{2}) where {T} =
    (θ = atan(M[2, 1] - M[1, 2], M[1, 1] + M[2, 2]),)
function angles(M::AbstractMatrix{T}, ::Val{3}) where {T}
    R = RotZYZ(M)
    return (θ = R.theta2, ϕ = R.theta1, ψ = R.theta3)
end

"""
    angles(M::AbstractMatrix{T})

Determines the Euler angles corresponding to the input matrix supposed to be a rotation matrix or at least a similarity

# Examples
```julia
julia> θ, ϕ, ψ = symbols("θ, ϕ, ψ", real = true) ; ℬʳ = RotatedBasis(θ, ϕ, ψ) ; display(vecbasis(ℬʳ, :cov))
3×3 Tensor{2, 3, Sym, 9}:
 -sin(ψ)⋅sin(ϕ) + cos(θ)⋅cos(ψ)⋅cos(ϕ)  -sin(ψ)⋅cos(θ)⋅cos(ϕ) - sin(ϕ)⋅cos(ψ)  sin(θ)⋅cos(ϕ)
  sin(ψ)⋅cos(ϕ) + sin(ϕ)⋅cos(θ)⋅cos(ψ)  -sin(ψ)⋅sin(ϕ)⋅cos(θ) + cos(ψ)⋅cos(ϕ)  sin(θ)⋅sin(ϕ)
                        -sin(θ)⋅cos(ψ)                          sin(θ)⋅sin(ψ)         cos(θ)

julia> angles(ℬʳ)
(θ = θ, ϕ = ϕ, ψ = ψ)
```
"""
angles(M::AbstractMatrix{T}) where {T} = angles(M, Val(size(M, 1)))
angles(ℬ::RotatedBasis) = ℬ.angles

angles(v::AbstractVector{T}, ::Val{2}) where {T} = (θ = atan(v[2], v[1]),)
angles(v::AbstractVector{T}, ::Val{3}) where {T} =
    (θ = atan(√(v[1]^2 + v[2]^2), v[3]), ϕ = atan(v[2], v[1]))
angles(v::AbstractVector{T}) where {T} = angles(v, Val(size(v, 1)))

invvar(::Val{:cov}) = :cont
invvar(::Val{:cont}) = :cov
invvar(var) = invvar(Val(var))

"""
    vecbasis(ℬ::AbstractBasis, var = :cov)

Returns the primal (if `var = :cov`) or dual (if `var = :cont`) basis
"""
vecbasis(ℬ::AbstractBasis, ::Val{:cov}) = ℬ.eᵢ
vecbasis(ℬ::AbstractBasis, ::Val{:cont}) = ℬ.eⁱ
vecbasis(::CanonicalBasis{dim,T}, ::Val{:cov}) where {dim,T} = Id2{dim,T}()
vecbasis(::CanonicalBasis{dim,T}, ::Val{:cont}) where {dim,T} = Id2{dim,T}()

vecbasis(ℬ::AbstractBasis, var) = vecbasis(ℬ, Val(var))
vecbasis(ℬ::AbstractBasis) = vecbasis(ℬ, :cov)
vecbasis(ℬ::AbstractBasis, i::Integer, j::Integer, var = :cov) = vecbasis(ℬ, Val(var))[i, j]
vecbasis(ℬ::AbstractBasis, i::Integer, var = :cov) = vecbasis(ℬ, Val(var))[:, i]

strvecbasis(::AbstractBasis, i::Integer, ::Val{:cov} ; vec = "𝐞") = vec * subscriptnumber(i)
strvecbasis(::AbstractBasis, i::Integer, ::Val{:cont} ; vec = "𝐞") = vec * superscriptnumber(i)
strvecbasis(ℬ::AbstractBasis, i::Integer, var = :cov ; vec = "𝐞") = strvecbasis(ℬ, i, Val(var) ; vec = vec)

const dsubscriptchar = Dict(
    "a" => "ₐ",
    "e" => "ₑ",
    "h" => "ₕ",
    "i" => "ᵢ",
    "j" => "ⱼ",
    "k" => "ₖ",
    "l" => "ₗ",
    "m" => "ₘ",
    "n" => "ₙ",
    "o" => "ₒ",
    "p" => "ₚ",
    "r" => "ᵣ",
    "s" => "ₛ",
    "t" => "ₜ",
    "u" => "ᵤ",
    "v" => "ᵥ",
    "x" => "ₓ",
    "β" => "ᵦ",
    "γ" => "ᵧ",
    "ρ" => "ᵨ",
    "ϕ" => "ᵩ",
    "χ" => "ᵪ",
)

const dsuperscriptchar = Dict(
    "a" => "ᵃ",
    "b" => "ᵇ",
    "c" => "ᶜ",
    "d" => "ᵈ",
    "e" => "ᵉ",
    "f" => "ᶠ",
    "g" => "ᵍ",
    "h" => "ʰ",
    "i" => "ⁱ",
    "j" => "ʲ",
    "k" => "ᵏ",
    "l" => "ˡ",
    "m" => "ᵐ",
    "n" => "ⁿ",
    "o" => "ᴼ",
    "p" => "ᵖ",
    "r" => "ʳ",
    "s" => "ˢ",
    "t" => "ᵗ",
    "u" => "ᶸ",
    "v" => "ᵛ",
    "w" => "ʷ",
    "x" => "ˣ",
    "y" => "ʸ",
    "z" => "ᶻ",
    "β" => "ᵝ",
    "γ" => "ᵞ",
    "ε" => "ᵋ",
    "θ" => "ᶿ",
    "ι" => "ᶥ ",
    "ϕ" => "ᵠ",
    "χ" => "ᵡ",
)

subscriptchar(s::String) = s ∈ keys(dsubscriptchar) ? dsubscriptchar[s] : s ∈ keys(dsuperscriptchar) ? dsuperscriptchar[s] : s
superscriptchar(s::String) = s ∈ keys(dsuperscriptchar) ? dsuperscriptchar[s] : s ∈ keys(dsubscriptchar) ? dsubscriptchar[s] : s

strvecbasis(::AbstractBasis, i::AbstractString, ::Val{:cov} ; vec = "𝐞") = vec * subscriptchar(i)
strvecbasis(::AbstractBasis, i::AbstractString, ::Val{:cont} ; vec = "𝐞") = vec * superscriptchar(i)
strvecbasis(ℬ::AbstractBasis, i::AbstractString, var = :cov ; vec = "𝐞") = strvecbasis(ℬ, i, Val(var) ; vec = vec)

"""
    metric(ℬ::AbstractBasis, var = :cov)

Returns the covariant (if `var = :cov`) or contravariant (if `var = :cont`) metric matrix
"""
metric(ℬ::AbstractBasis, ::Val{:cov}) = ℬ.gᵢⱼ
metric(ℬ::AbstractBasis, ::Val{:cont}) = ℬ.gⁱʲ
metric(::OrthonormalBasis{dim,T}, ::Val{:cov}) where {dim,T} = Id2{dim,T}()
metric(::OrthonormalBasis{dim,T}, ::Val{:cont}) where {dim,T} = Id2{dim,T}()

metric(ℬ::AbstractBasis, var) = metric(ℬ, Val(var))
metric(ℬ::AbstractBasis) = metric(ℬ, :cov)
metric(ℬ::AbstractBasis, i::Integer, j::Integer, var = :cov) = metric(ℬ, Val(var))[i, j]

"""
    normalize(ℬ::AbstractBasis, var = cov)

Builds a basis after normalization of column vectors of input matrix `v` where columns define either
- primal vectors ie `eᵢ=v[:,i]/norm(v[:,i])` if `var = :cov` as by default
- dual vector ie `eⁱ=v[:,i]/norm(v[:,i])` if `var = :cont`.
"""
function LinearAlgebra.normalize(ℬ::AbstractBasis, var = :cov)
    w = copy(vecbasis(ℬ, var))
    for i = 1:size(w, 2)
        w[:, i] /= norm(w[:, i])
    end
    return Basis(w, var)
end

"""
    isorthogonal(ℬ::AbstractBasis)

Checks whether the basis `ℬ` is orthogonal
"""
isorthogonal(ℬ::AbstractBasis) = isdiagonal(metric(ℬ))

isorthogonal(::OrthonormalBasis) = true

isorthogonal(::OrthogonalBasis) = true

"""
    isorthonormal(ℬ::AbstractBasis)

Checks whether the basis `ℬ` is orthonormal
"""
isorthonormal(ℬ::AbstractBasis) = isidentity(metric(ℬ))

isorthonormal(::OrthonormalBasis) = true

#####################
# Display Functions #
#####################
for OP in (:show, :print, :display)
    @eval begin
        function Base.$OP(ℬ::AbstractBasis)
            $OP(typeof(ℬ))
            print("→ basis: ")
            $OP(vecbasis(ℬ, :cov))
            print("→ dual basis: ")
            $OP(vecbasis(ℬ, :cont))
            print("→ covariant metric tensor: ")
            $OP(metric(ℬ, :cov))
            print("→ contravariant metric tensor: ")
            $OP(metric(ℬ, :cont))
        end
    end
end

export Basis, CanonicalBasis, RotatedBasis, CylindricalBasis, SphericalBasis, OrthogonalBasis, AllOrthogonalBasis
export getdim, vecbasis, metric, angles, isorthogonal, isorthonormal


text:
using Revise
using POMDPs
using BlinkPOMDPSimulator
using POMDPModels # requires POMDPModels with SimpleGridWorld (v0.3.3 or higher?)
using POMDPPolicies

bs = BlinkSimulator()

mdp = SimpleGridWorld()
p = FunctionPolicy(s->:up)

simulate(bs, mdp, p)


text:
using Revise

using BNPRegressionGGA2021
using Cairo
using CSV
using DataFrames
using Distributions
using Gadfly
using Random

const BNP = BNPRegressionGGA2021

#======================================#
# Replicate the 2nd simulation example #
#======================================#

function simulate_sample(N0)
    dy1 = LogNormal(0, 0.3)
    dy2 = MixtureModel(LogNormal, [(-1.0, 1.0), (0.3, 1.0)], [0.6, 0.4]);
    X0 = [ones(N0) rand([0, 1], N0, 5)]
    X1 = ones(1, 6)
    y0 = [X0[i, end] == 1 ? rand(dy2) : rand(dy1) for i in 1:N0]
    c0 = rand(Exponential(10), N0)
    y1 = zeros(1)
    z0 = deepcopy(y0)
    return dy1, dy2, z0, c0, y0, X0, y1, X1
end

function simulate_test_results(N0)
    _, _, _, c0, y0, X0, y1, X1 = simulate_sample(N0)
    m = BNP.DGPMErlang(; c0, y0, X0, y1, X1);
    _, _, chaing = BNP.sample!(m; mcmcsize = 4000, burnin = 2000)
    return mean(chaing)[2:end]
end

Random.seed!(1);
niter = 100;
results = zeros(100, 5);
for i in 1:niter
    results[i, :] = simulate_test_results(1000);
end
df = DataFrame(results, :auto)
CSV.write("data/simulation-example-02_bnp.csv", df);


text:
module CountDownLatches

export CountDownLatch, await, count_down, get_count

"""
    CountDownLatch(count::I) where {I<:Integer}

A `CountDownLatch` that starts at some initial non-negative `count`. This latch can be counted down and waited on.
"""
mutable struct CountDownLatch
    count::Threads.Atomic{Int64}
    condition::Threads.Condition

    function CountDownLatch(count::I) where {I<:Integer}
        if count < 0
            throw(ArgumentError("Count ($count) must not be negative."))
        end

        new(Threads.Atomic{Int64}(count), Threads.Condition(ReentrantLock()))
    end
end

"""
    await(latch::CountDownLatch)

Waits until the `latch` has counted down to `<= 0`.
"""
function await(latch::CountDownLatch)
    # Check the count early to possibly avoid a call to lock
    if get_count(latch) <= 0
        return
    end

    lock(latch.condition) do
        while get_count(latch) > 0
            wait(latch.condition)
        end
    end
end

"""
    count_down(latch::CountDownLatch)

Counts down the `latch` once. This may cause the count to become negative. Any tasks that were blocked in `await` will
be woken if the count becomes `<= 0`.
"""
function count_down(latch::CountDownLatch)
    # I think doing something like:
    #
    # 1: count = get_count(latch)
    # 2: count == 0 && return
    # 3: next_count = count - 1
    # 4: if Threads.atomic_cas!(latch.count, count, next_count) == count
    # 5:     lock(latch.condition) do
    # 6:         notify(latch.condition, nothing; all = true, error = false)
    # 7:     end
    # 8: end
    #
    # is wrong here because of the case where a thread decrements count after line 1 and before line 4. In that case,
    # two calls to `notify` must still happen, but only one would.

    Threads.atomic_sub!(latch.count, 1)
    lock(latch.condition) do
        notify(latch.condition, nothing; all = true, error = false)
    end
    nothing
end

"""
    get_count(latch::CountDownLatch)

Returns the current count of the `latch`. This count may be negative.
"""
function get_count(latch::CountDownLatch)
    latch.count[]
end

end


text:
struct IndexablePoint2D <: VoronoiDelaunay.AbstractPoint2D
    _x::Float64
    _y::Float64
    _index::Int64
end
IndexablePoint2D(x::Float64, y::Float64) = IndexablePoint2D(x, y, -1)

getx(p::IndexablePoint2D) = p._x
gety(p::IndexablePoint2D) = p._y
Base.getindex(p::IndexablePoint2D) = p._index

const AbstractPoints2D = Vector{VoronoiDelaunay.AbstractPoint2D}
const IndexablePoints2D = Vector{IndexablePoint2D}
const Points2D = Vector{VoronoiDelaunay.Point2D}
const Tessellation = Dict{Int64, Points2D}

# Edges of the bounding box
const LEFT = VoronoiDelaunay.min_coord
const RIGHT = VoronoiDelaunay.max_coord
const LOWER = VoronoiDelaunay.min_coord
const UPPER = VoronoiDelaunay.max_coord

const MIDDLEy = 0.5*(LOWER + UPPER)
const MIDDLEx = 0.5*(LEFT + RIGHT)

const LL = VoronoiDelaunay.Point2D(LEFT, LOWER)
const RL = VoronoiDelaunay.Point2D(RIGHT, LOWER)
const RU = VoronoiDelaunay.Point2D(RIGHT, UPPER)
const LU = VoronoiDelaunay.Point2D(LEFT, UPPER)

# Corners of the bounding box ordered by quadrant
const BoxCorners = [RU; LU; LL; RL]


text:
module Mathieu
using LinearAlgebra
using LinearAlgebra: BlasReal
using Base.Iterators: product, filter

const Order = Union{Integer, AbstractVector{T} where T<:Integer}

include("integerorder_mat.jl")
include("integerorder_evals.jl")
include("integerorder_coeffs.jl")
include("integerorder_funs.jl")
include("integerorder_derivs.jl")
include("integerorder_special.jl")
include("utils.jl")

end # module


text:
module TransformerModule

using LinearAlgebra: norm, Transpose, mul!
using FinEtools

"""
    QTEQTransformer

QTEQTransformer of element matrices.

A callable object: computes `Q^T * E * Q`, where `E` the element stiffness
matrix, and `Q` is the transformation matrix.

Buffers the intermediate result. Hence no allocation is incurred.
"""
struct QTEQTransformer
    _buff::FFltMat
    function QTEQTransformer(Q)
        _buff = fill(0.0, size(Q)...);
        return new(_buff)
    end
end

"""
    (o::QTEQTransformer)(E, Q)

Perform the transformation on the matrix `E` with the transformation matrix
`Q`: `Ebar = Q^T * E * Q`.
"""
(o::QTEQTransformer)(E, Q) = begin
    @assert size(o._buff) == size(Q)
    @assert size(E) == size(Q)
    mul!(o._buff, E, Q)
    mul!(E, Transpose(Q), o._buff)
    return E
end

"""
    QEQTTransformer

QEQTTransformer of element matrices.

A callable object: computes `Q * E * Q^T`, where `E` the element stiffness
matrix, and `Q` is the transformation matrix.

Buffers the intermediate result. Hence no allocation is incurred.
"""
struct QEQTTransformer
    _buff::FFltMat
    function QEQTTransformer(Q)
        _buff = fill(0.0, size(Q)...);
        return new(_buff)
    end
end

"""
    (o::QEQTTransformer)(E, Q)

Perform the transformation on the matrix `E` with the transformation matrix
`Q`: `Ebar = Q * E * Q^T`.
"""
(o::QEQTTransformer)(E, Q) = begin
    @assert size(o._buff) == size(Q)
    @assert size(E) == size(Q)
    mul!(o._buff, E, Transpose(Q))
    mul!(E, Q, o._buff)
    return E
end

struct Layup2ElementAngle
    M::FFltMat
    M2D::FFltMat
end

function Layup2ElementAngle()
    Layup2ElementAngle(fill(0.0, 3, 3), fill(0.0, 2, 2))
end

function (o::Layup2ElementAngle)(E_G, lcsmat)
    mul!(o.M, E_G', lcsmat)
    o.M2D[1, 1] = o.M[1, 1]
    o.M2D[1, 2] = o.M[1, 2]
    o.M2D[2, 1] = o.M[2, 1]
    o.M2D[2, 2] = o.M[2, 2]
    o.M2D[:, 1] ./= norm(@view o.M2D[:, 1])
    o.M2D[:, 2] ./= norm(@view o.M2D[:, 2])
    m = (o.M2D[1, 1] + o.M2D[2, 2]) / 2
    n = (o.M2D[1, 2] - o.M2D[2, 1]) / 2
    sn = n >= 0.0 ? +1.0 : -1.0
    n = sn * sqrt(1 - m^2)
    return m, n
end

end # module


text:
# struct ParentsBeforeChildrenIterator
#     ops::Vector{Operation}
#     visited::Vector{Bool}
# end
# function iterate(pbci::ParentsBeforeChildrenIterator)
#     for (i,op) ∈ enumerate(pbci.ops)
#         if iszero(length(parents(op)))
#             pbci.visited[i] = true
#             return op, pbci
#         end
#     end
#     nothing
# end
# function iterate()

# end

function dependent_outer_reducts(ls::LoopSet, op)
    for i ∈ ls.outer_reductions
        search_tree(parents(operations(ls)[i]), name(op)) && return true
    end
    false
end

function isnopidentity(ls::LoopSet, op::Operation, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, suffix)
    parents_op = parents(op)
    if iscompute(op) && instruction(op).instr === :identity && name(first(parents_op)) === name(op) && isone(length(parents_op))
        mvar, u₁unrolledsym, u₂unrolledsym = variable_name_and_unrolled(op, u₁loop, u₂loop, suffix)
        parents_u₁syms, parents_u₂syms = parent_unroll_status(op, u₁loop, u₂loop, suffix)
        if (u₁unrolledsym == first(parents_u₁syms)) && ((!isnothing(suffix)) == parents_u₂syms[1])
            #TODO: identifer(first(parents_op)) ∉ ls.outer_reductions is going to miss a lot of cases
            #Should probably replace that with `DVec` (demoting Vec) types, that demote to scalar.
            if (isvectorized(first(parents_op)) && !isvectorized(op)) && !dependent_outer_reducts(ls, op)
                op.instruction = reduction_to_scalar(instruction(first(parents_op)))
                op.mangledvariable = gensym(op.mangledvariable)
                false
            else
                true
            end
        else
            false
        end
    else
        false
    end
end

function set_upstream_family!(adal::Vector{T}, op::Operation, val::T, ld::Vector{Symbol}, id::Int) where {T}
    adal[identifier(op)] == val && return # must already have been set
    if ld != loopdependencies(op) || id == identifier(op)
        (adal[identifier(op)] = val)
    end
    for opp ∈ parents(op)
        identifier(opp) == identifier(op) && continue
        set_upstream_family!(adal, opp, val, ld, id)
    end
end

function addoptoorder!(
    ls::LoopSet, included_vars::Vector{Bool}, place_after_loop::Vector{Bool}, op::Operation,
    loopsym::Symbol, _n::Int, u₁loop::Symbol, u₂loop::Symbol, vectorized::Symbol, loopistiled::Bool
)
    lo = ls.loop_order
    id = identifier(op)
    included_vars[id] && return nothing
    loopsym ∈ loopdependencies(op) || return nothing
    for opp ∈ parents(op) # ensure parents are added first
        addoptoorder!(ls, included_vars, place_after_loop, opp, loopsym, _n, u₁loop, u₂loop, vectorized, loopistiled)
    end
    included_vars[id] && return nothing
    included_vars[id] = true
    isunrolled = (isu₁unrolled(op)) + 1
    istiled = isu₂unrolled(op)
    # optype = Int(op.node_type) + 1
    after_loop = place_after_loop[id] + 1
    if !isloopvalue(op)
        if istiled
            isnopidentity(ls, op, u₁loop, u₂loop, vectorized, 0) || push!(lo[isunrolled,2,after_loop,_n], op)
        else
            isnopidentity(ls, op, u₁loop, u₂loop, vectorized, nothing) || push!(lo[isunrolled,1,after_loop,_n], op)
        end
    end
    # @show op, after_loop
    # isloopvalue(op) || push!(lo[isunrolled,istiled,after_loop,_n], op)
    # all(opp -> iszero(length(reduceddependencies(opp))), parents(op)) &&
    set_upstream_family!(place_after_loop, op, false, loopdependencies(op), identifier(op)) # parents that have already been included are not moved, so no need to check included_vars to filter
    nothing
end

function fillorder!(ls::LoopSet, order::Vector{Symbol}, u₁loop::Symbol, u₂loop::Symbol, loopistiled::Bool, vectorized::Symbol)
    lo = ls.loop_order
    resize!(lo, length(ls.loopsymbols))
    ro = lo.loopnames # reverse order; will have same order as lo
    nloops = length(order)
    ops = operations(ls)
    nops = length(ops)
    included_vars = fill!(resize!(ls.included_vars, nops), false)
    place_after_loop = fill!(resize!(ls.place_after_loop, nops), true)
    # to go inside out, we just have to include all those not-yet included depending on the current sym
    empty!(lo)
    for _n ∈ 1:nloops
        n = 1 + nloops - _n
        ro[_n] = loopsym = order[n]
        #loopsym = order[n]
        for op ∈ ops
            addoptoorder!( ls, included_vars, place_after_loop, op, loopsym, _n, u₁loop, u₂loop, vectorized, loopistiled )
        end
    end
end


text:
# --- Coupling Model Creation --- #

"""
    Coupling(::QuasiSteady, ::Section)

Coupling model for coupling a quasi-steady aerodynamic model based on thin airfoil theory
(see [`QuasiSteady`](@ref)) and a two-degree of freedom typical section model
(see [`Section()`]).  This model introduces the freestream velocity ``U``, air density
``\\rho``, and air speed of sound ``c`` as additional parameters.
"""
function Coupling(::QuasiSteady, ::Section)

    # coupling function
    g = quasisteady_section_inputs

    # number of states, inputs, and parameters (use Val(N) to use inferrable dimensions)
    nx = Val(4)
    ny = Val(2)
    np = Val(14)

    # number of parameters introduced by the coupling (use Val(N) to use inferrable dimensions)
    npc = Val(3)

    # jacobian definitions
    ratejac = Nonlinear() # TODO: define rate jacobian function
    statejac = Nonlinear() # TODO: define state jacobian function
    paramjac = Nonlinear() # TODO: define parameter jacobian function
    tgrad = Zeros()

    # convenience function for setting coupling parameters
    setparam = quasisteady_section_set_parameters!

    # convenience function for separating coupling parameters
    sepparam = quasisteady_section_separate_parameters

    # return resulting coupling
    return Coupling{false}(g, nx, ny, np, npc;
        ratejac = ratejac,
        statejac = statejac,
        paramjac = paramjac,
        tgrad = tgrad,
        setparam = setparam,
        sepparam = sepparam)
end

# --- Internal Methods --- #

# coupling function
function quasisteady_section_inputs(dx, x, p, t)
    # extract state variables
    dh, dθ, dhdot, dθdot = dx
    # extract state variables
    h, θ, hdot, θdot = x
    # extract aerodynamic, structural, and aerostructural parameters
    a, b, a0, α0, cd0, cm0, kh, kθ, m, Sθ, Iθ, U, ρ, c = p
    # local freestream velocity components
    u, v, ω = section_velocities(U, θ, hdot, θdot)
    # local freestream accelerations
    udot, vdot, ωdot = section_accelerations(dhdot, dθdot)
    # calculate aerodynamic loads
    N, A, M = quasisteady_loads(a, b, ρ, c, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # lift is approximately equal to the normal force
    L = N
    # return inputs
    return SVector(L, M)
end

# convenience function for defining the coupling function parameters
function quasisteady_section_set_parameters!(p; U, rho, c)
    p[1] = U
    p[2] = rho
    p[3] = c
    return p
end

# convenience function for separating the coupling function parameters
quasisteady_section_separate_parameters(p) = (U = p[1], rho = p[2], c = p[3])


text:
using PSSFSS
import PSSFSS.Sheets: SV2, recttri, combine, read_sheet_data, write_sheet_data
using Test

bl = SV2([0.,0.])
tr = SV2([1.,1.])
nx = ny = 1
sh1 = recttri(bl, tr, nx, ny)
@testset "recttri" begin
    @test length(sh1.ρ) == 4
    @test sh1.ρ[1] == [0.,0.]
    @test sh1.ρ[2] == [1.,0.]
    @test sh1.ρ[3] == [0.,1.]
    @test sh1.ρ[4] == [1.,1.]
    @test sh1.e1 == [1,3,1,2,1]
    @test sh1.e2 == [2,4,3,4,4]
    @test sh1.fe == reshape([2,3,5, 4,5,1], 3,2)
    @test sh1.fv == reshape([1,4,3, 1,2,4], 3,2)
end

bl = SV2([1.,0.])
tr = SV2([2.,1.])
nx = ny = 1
sh2 = recttri(bl, tr, nx, ny)
dup_coor = 'x'
dup_coor_value = 1
sh3 = combine(sh1, sh2, dup_coor, dup_coor_value)
@testset "combine right" begin
    @test sh3.ρ == [[0,0], [1,0], [0,1], [1,1], [2,0], [2,1]]
    @test sh3.e1 == [1, 3, 1, 2, 1, 2, 4, 5, 2]
    @test sh3.e2 == [2, 4, 3, 4, 4, 5, 6, 6, 6]
    @test sh3.fv == reshape([1,4,3,  1,2,4,  2,6,4,  2,5,6], 3,4)
    @test sh3.fe == reshape([2,3,5,  4,5,1,  7,4,9,  8,9,6], 3,4)
end

bl = SV2([0.,1.])
tr = SV2([1.,2.])
nx = ny = 1
sh2 = recttri(bl, tr, nx, ny)
dup_coor = 'y'
dup_coor_value = 1.0
sh3 = combine(sh1, sh2, dup_coor, dup_coor_value)
@testset "combine top" begin
    @test sh3.ρ == [[0,0], [1,0], [0,1], [1,1], [0,2], [1,2]]
    @test sh3.e1 == [1, 3, 1, 2, 1, 5, 3, 4, 3]
    @test sh3.e2 == [2, 4, 3, 4, 4, 6, 5, 6, 6]
    @test sh3.fe == [2 4 6 8; 3 5 7 9; 5 1 9 2]
    @test sh3.fv == [1 1 3 3; 4 2 6 4; 3 4 5 6]
end

bl = SV2([1.,1.])
tr = SV2([2.,2.])
nx = ny = 1
sh2 = recttri(bl, tr, nx, ny)
dup_coor = 'y'
dup_coor_value = 1.0
sh3 = combine(sh1, sh2, dup_coor, dup_coor_value)
@testset "combine topright" begin
    @test sh3.ρ == [[0, 0], [1, 0], [0, 1],
                     [1, 1], [2, 1], [1, 2], [2, 2]]
    @test sh3.e1 == [1, 3, 1, 2, 1, 4, 6, 4, 5, 4]
    @test sh3.e2 == [2, 4, 3, 4, 4, 5, 7, 6, 7, 7]
    @test sh3.fe == [2 4 7 9; 3 5 8 10; 5 1 10 6]
    @test sh3.fv == [1 1 4 4; 4 2 7 5; 3 4 6 7]
end

bl = SV2([2.,0.])
tr = SV2([3.,1.])
nx = ny = 1
sh2 = recttri(bl, tr, nx, ny)
dup_coor = ' '
dup_coor_value = NaN
sh3 = combine(sh1, sh2, dup_coor, dup_coor_value)
@testset "combine disjoint" begin
    @test sh3.ρ == [[0,0], [1,0], [0,1], [1,1], [2,0], [3,0], [2,1], [3,1]]
    @test sh3.e1 == [1, 3, 1, 2, 1, 5, 7, 5, 6, 5]
    @test sh3.e2 == [2, 4, 3, 4, 4, 6, 8, 7, 8, 8]
    @test sh3.fe == [2 4 7 9; 3 5 8 10; 5 1 10 6]
    @test sh3.fv == [1 1 5 5; 4 2 8 6; 3 4 7 8]
end

@testset "read and write" begin
    fname = tempname()
    write_sheet_data(fname, sh3)
    sh4 = read_sheet_data(fname)
    @test sh3 == sh4
end


text:
# Add PATH (adapt!)
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Leaf/");
push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/Utils/");
#push!(LOAD_PATH, "/Volumes/cfranken/code/gitHub/LSM-SPAM/src/tools/");

# Use Plots:
using Plots
using BenchmarkTools
using Statistics
using Parameters

using Revise
using CanopyRTMod
using LeafPhotosynthesis

@unpack wl,wle,wlf, soil = CanopyRTMod;
# Create an array of standard leaves (needs to be in Module later on:
arrayOfLeaves = Array{leafbio{FT,length(wl), length(wle), length(wlf),length(wle)*length(wlf)}, 1}(undef, CanopyRTMod.canopy.nlayers)
for i = 1:CanopyRTMod.canopy.nlayers
    #@show i
    arrayOfLeaves[i] = leafbio{FT, length(wl), length(wle), length(wlf),length(wle)*length(wlf)}()
    CanopyRTMod.fluspect!(arrayOfLeaves[i], CanopyRTMod.optis)
end

# 4 Different steps to compute Short-Wave RT
@time CanopyRTMod.computeCanopyGeomProps!(canopy, angles,canOpt)
@time CanopyRTMod.computeCanopyMatrices!(arrayOfLeaves,canOpt);
@time CanopyRTMod.RTM_SW!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil);
@time CanopyRTMod.deriveCanopyFluxes!(canopy, canOpt, canRad,sunRad, CanopyRTMod.soil, arrayOfLeaves);
# Compute Long Wave (Last term is LW incoming in W m^-2)
@time CanopyRTMod.computeThermalFluxes!(arrayOfLeaves, canOpt, canRad, canopy, soil, [Float32(400.0)]);

#@show arrayOfLeaves[1].kChlrel
# Layer Temperatures are here:
@show canRad.T_sun;
@show canRad.T_shade;

# Net Energy fluxes
@show canRad.intNetLW_shade;
@show canRad.intNetLW_sunlit;
@show canRad.intNetSW_shade;
@show canRad.intNetSW_sunlit;
@show canRad.RnSoilLW;
@show canRad.RnSoil;
#@show 1e6 * canRad.absPAR_sunCab;
@show 1e6 * canRad.absPAR_shadeCab;

plot(wl, canRad.netSW_sunlit)

iLAI = canopy.LAI/canopy.nlayers
plot(wl,1 ./canOpt.Ps[10] / iLAI * canRad.netSW_sunlit[:,10])
plot!(wl, sunRad.E_direct)

using Leaf
l = leaf_params{Float32}();
l2 = leaf_params{Float32}();
# Create a Flux structure
f = LeafPhotosynthesis.fluxes{Float32}();
l.vcmax25=120
l.jmax25=l.vcmax25*1.8

# initialize some reasonable values
f.je = 100;
f.gbc = 100;
f.gbv = 100;
f.ceair=1500;
f.eair = 1500;
f.APAR = 1.0e6 * canRad.absPAR_shadeCab[1];
1e6*mean(canRad.absPAR_sunCab[:,:,1])

A = similar(canRad.absPAR_sunCab)
Ashade = similar(canRad.absPAR_shadeCab)
I = CartesianIndices(A)
IShade = CartesianIndices(Ashade)
for i in I
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,Float32(298.0), Float32(1.0e6) * canRad.absPAR_sunCab[i]);
    A[i]=f.an;
end
for i in IShade
    LeafPhotosynthesis.LeafPhotosynthesis(f,l,Float32(298.0), Float32(1.0e6) * canRad.absPAR_shadeCab[i]);
    Ashade[i]=f.an;
end

LeafPhotosynthesis.LeafPhotosynthesis(f,l,Float32(298.0), Float32(1.0e6) * mean(canRad.absPAR_sunCab[:,:,1]));

using Plots

plot(Float32(1.0e6) * canRad.absPAR_sunCab[:,:,1][:], A[:,:,1][:], seriestype = :scatter,legend=:bottomright, label="Individual sunlit leaves in one layer  (A(APAR))")
plot!([Float32(1.0e6) * mean(canRad.absPAR_sunCab[:,:,1])], [f.an], seriestype = :scatter,  label="A(average(APAR))")
plot!([Float32(1.0e6) * mean(canRad.absPAR_sunCab[:,:,1])], [mean(A)], seriestype = :scatter, label="average(A((APAR)))")
plot!(Float32(1.0e6) * canRad.absPAR_shadeCab[:], Ashade, seriestype = :scatter,legend=:bottomright, label="Individual shaded layers")

f

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
require("Options")

# Note that the two calls to using Stats in this file are
# strictly required: one pulls the names into DataFrames
# for easy access within the module, whereas the other call
# pushes those names into Main.
using Stats

module DataFrames

##############################################################################
##
## Dependencies
##
##############################################################################

using Base.Intrinsics
using OptionsMod
using Stats

##############################################################################
##
## Global constants
##
##############################################################################

const DEFAULT_COLUMN_TYPE = Float64
const POOLED_DATA_VEC_REF_TYPE = Uint16
const POOLED_DATA_VEC_REF_CONVERTER = uint16

##############################################################################
##
## Extend methods in Base by default
##
##############################################################################

importall Base
importall Stats
import Base.Sort.Algorithm, Base.Sort.By, Base.Sort.Forward
import Base.Sort.Ordering, Base.Sort.Perm, Base.Sort.lt
import Base.Sort.sort, Base.Sort.sort!, Base.Sort.sortby, Base.Sort.sortby!

##############################################################################
##
## Exported methods and types
##
##############################################################################

export # reconcile_groups,
       @DataFrame,
       @transform,
       AbstractDataArray,
       AbstractDataFrame,
       AbstractDataVector,
       AbstractIndex,
       anyna,
       allna,
       array,
       based_on,
       by,
       cbind,
       clean_colnames!,
       colffts,
       colmaxs,
       colmeans,
       colmedians,
       colmins,
       colnames!,
       colnames,
       colnorms,
       colprods,
       colranges,
       colstds,
       colsums,
       coltypes,
       colvars,
       colwise,
       combine,
       complete_cases,
       complete_cases!,
       cut,
       DataArray,
       databool,
       datadiagm,
       dataeye,
       datafalses,
       datafloat,
       DataFrame,
       dataint,
       DataMatrix,
       dataones,
       DataStream,
       datatrues,
       DataVector,
       datazeros,
       describe,
       dict,
       drop_duplicates!,
       duplicated,
       each_failNA,
       each_removeNA,
       each_replaceNA,
       EachCol,
       EachRow,
       failNA,
       findat,
       flipud!,
       flipud,
       Formula,
       get_groups,
       get_indices,
       gl,
       GroupApplied,
       groupby,
       GroupedDataFrame,
       head,
       Index,
       index,
       index_to_level,
       IndexedVector,
       Indexer,
       interaction_design_matrix,
       is_group,
       isna,
       letters,
       LETTERS,
       level_to_index,
       levels,
       load_df,
       matrix,
       merge,
       model_frame,
       model_matrix,
       model_response,
       ModelFrame,
       ModelMatrix,
       NA,
       NAException,
       nafilter,
       naFilter,
       NamedArray,
       names!,
       nareplace,
       naReplace,
       nas,
       NAtype,
       ncol,
       nrow,
       padNA,
       paste,
       pdatafalses,
       pdataones,
       pdatatrues,
       pdatazeros,
       percent_change,
       PooledDataVecs, # The capitalization and/or name for this is a bit inconsistent (merge_pools, maybe?). Do we want to export?
       PooledDataArray,
       PooledDataMatrix,
       PooledDataVector,
       print_table,
       range,
       rbind,
       read_minibatch,
       read_table,
       reldiff,
       removeNA,
       rename_group!,
       reorder,
       rep,
       replace!,
       rename!,
       rename,
       replaceNA,
       rowffts,
       rowmaxs,
       rowmeans,
       rowmedians,
       rowmins,
       rownorms,
       rowprods,
       rowranges,
       rowstds,
       rowsums,
       rowvars,
       save,
       set_group,
       set_groups,
       set_levels,
       set_levels!,
       SimpleIndex,
       stack,
       SubDataFrame,
       subset,
       table,
       tail,
       unique,
       unstack,
       vector,
       with,
       within!,
       within,
       without,
       write_table,
       xtab,
       xtabs,
       stack_df,
       StackedVector,
       RepeatedVector,
       EachRepeatedVector,
       melt,
       melt_df,
       pivot_table,
       RComplex,     # Vector{Complex128}
       RInteger,     # Vector{Int32} plus BitVector of NA indicators
       RLogical,     # BitVector of values and BitVector of NA indicators
       RNumeric,     # Vector{Float64}
       RList,        # Vector{Any}
       RString,      # Vector{ASCIIString}
       RSymbol,      # Symbol stored as an String b/c of embedded '.'

       class,                              # in the S3 sense of "class"
       inherits,
       read_rda,
       vecbind

##############################################################################
##
## Load files
##
##############################################################################

include("utils.jl")
include("natype.jl")
include("dataarray.jl")
include("datavector.jl")
include("datamatrix.jl")
include("pooleddataarray.jl")
include("index.jl")
include("namedarray.jl")
include("dataframe.jl")
include("merge.jl")
include("grouping.jl")
include("reshape.jl")
include("formula.jl")
include("io.jl")
include("datastream.jl")
include("linalg.jl")
include("operators.jl")
include("statistics.jl")
include("predicates.jl")
include("indexing.jl")
include("extras.jl")
include("RDA.jl")

# TODO: Remove these definitions
nafilter(x...) = error("Function removed. Please use removeNA")
nareplace(x...) = error("Function removed. Please use replaceNA")
naFilter(x...) = error("Function removed. Please use each_removeNA")
naReplace(x...) = error("Function removed. Please use each_replaceNA")

end # module DataFrames


text:
function GB.Mesh(df::GT.SignedDistanceField{3,ST,FT},
              method::AbstractMeshingAlgorithm;
              pointtype=GB.Point{3,Float32},
              facetype=GB.GLTriangleFace) where {ST, FT}

    h = df.bounds
    vts, fcs = isosurface(df.data, method, pointtype, facetype, origin=pointtype(GB.origin(h)), widths=pointtype(GB.widths(h)))
    return GB.Mesh(vts, fcs)
end

function GB.Mesh(f::Function,
                 h::GB.Rect,
                 samples::NTuple{3,T},
                 method::AbstractMeshingAlgorithm;
                 pointtype=GB.Point{3,Float32},
                 facetype=GB.GLTriangleFace) where {T <: Integer}
    vts, fcs = isosurface(f, method, pointtype, facetype; samples=samples,
                          origin=pointtype(GB.origin(h)), widths=pointtype(GB.widths(h)))

    return GB.Mesh(vts, fcs)
end

function GB.Mesh(f::Function, h::GB.Rect, method::AbstractMeshingAlgorithm;
              samples::NTuple{3,T}=_DEFAULT_SAMPLES,
              pointtype=GB.Point{3,Float32},
              facetype=GB.GLTriangleFace) where {T <: Integer}

    return GB.Mesh(f, h, samples, method; pointtype=pointtype, facetype=facetype)
end

function GB.Mesh(volume::AbstractArray{T, 3},
                 method::AbstractMeshingAlgorithm;
                 pointtype=GB.Point{3,Float32},
                 facetype=GB.GLTriangleFace,
                 kwargs...) where {T}

    vts, fcs = isosurface(volume, method, pointtype, facetype; kwargs...)
    return GB.Mesh(vts, fcs)
end


text:
using Plots
pyplot()

println("Ag) Generate and normalize continuum orbitals in a local potential.")
#
wa = Atomic.Computation("xx",  Nuclear.Model(3.); grid=JAC.Radial.Grid("grid: by given parameters"; rnt = 2.0e-5, h = 5.0e-2, hp = 2.0e-2, NoPoints = 800),
                        configs=[Configuration("1s^2 2s"), Configuration("1s^2 2p")],
                        asfSettings=AsfSettings(true, false, "meanDFS", "hydrogenic", "xxx", [1], 40, 1.0e-6, JAC.Subshell[],
                                                true, false, NoneQed(), "no-file", false, [1,2,3,4], false, JAC.LevelSymmetry[] ) )

wb = perform(wa, output=true)
#
multiplet = wb["multiplet:"]
grid      = wb["grid:"]
level     = multiplet.levels[1]
basis     = multiplet.levels[1].basis
#
println("\nCompute the nuclear and Hartree potentials.")
w1 = JAC.Nuclear.pointNucleus(5.0, grid)
## w2 = compute("radial potential: Dirac-Fock-Slater", grid, level)
## w2 = JAC.add(w1,w2)
w2 = JAC.Nuclear.nuclearPotentialDH(Nuclear.Model(5.), grid, 0.9)
w3 = JAC.Nuclear.nuclearPotentialDH(Nuclear.Model(5.), grid, 2.5)
plot("radial potentials", [w1, w2, w3], grid; N=500)

settings = Continuum.Settings(false, grid.nr-100)
setDefaults("method: continuum, Galerkin")
## setDefaults("method: normalization, pure sine")
setDefaults("method: normalization, pure Coulomb")
energy   = 30.;   sh = Subshell("3f_5/2")
wc       = JAC.Continuum.generateOrbitalLocalPotential(energy, sh, w1, settings)
wd       = JAC.Continuum.generateOrbitalLocalPotential(energy, sh, w2, settings)
we       = JAC.Continuum.generateOrbitalLocalPotential(energy, sh, w3, settings)
wavenb   = sqrt( 2energy + energy * JAC.give("alpha")^2 )
wavelgth = 2pi / wavenb
println("wavelength = $wavelgth")
plot("radial orbitals: large", [wc[1], wd[1], we[1]], grid; N=900)

error("stop here")

settings = Continuum.Settings(false, grid.nr)
setDefaults("method: continuum, pure sine")
setDefaults("method: normalization, pure sine")
wc       = JAC.Continuum.generateOrbitalLocalPotential(100.0, Subshell("3d_3/2"), w1, settings)

setDefaults("method: continuum, asymptotic Coulomb")
setDefaults("method: normalization, pure Coulomb")
wc       = JAC.Continuum.generateOrbitalLocalPotential(100.0, Subshell("3d_3/2"), w1, settings)

settings = Continuum.Settings(false, grid.nr)
setDefaults("method: continuum, nonrelativistic Coulomb")
setDefaults("method: normalization, pure sine")
wc       = JAC.Continuum.generateOrbitalNonrelativisticCoulomb(100.0, Subshell("3s_1/2"), 1.0, w1.grid, settings)

settings = Continuum.Settings(false, grid.nr)
setDefaults("method: continuum, nonrelativistic Coulomb")
setDefaults("method: normalization, pure sine")
wc       = JAC.Continuum.generateOrbitalNonrelativisticCoulomb(100.0, Subshell("3s_1/2"), w1, settings)


text:
# This is a Julia adaptation of the Matlab code from O. Ledoit and M. Wolf
# solely used to test our own implementation of their estimator.
# Their code was released under the BSD-2 license with the copyright notice
# below.
# https://www.econ.uzh.ch/en/people/faculty/wolf/publications.html#9

###########################################################################
# This file is released under the BSD 2-clause license.

# Copyright (c) 2014, Olivier Ledoit and Michael Wolf
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are
# met:
#
# 1. Redistributions of source code must retain the above copyright notice,
# this list of conditions and the following disclaimer.
#
# 2. Redistributions in binary form must reproduce the above copyright
# notice, this list of conditions and the following disclaimer in the
# documentation and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
# IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
# THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
# PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
# CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
# EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
# PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
# NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
###########################################################################

using LinearAlgebra
using Statistics

# comments marked with a `#` come from their code
# comments marked with a `#<#>` are ours

function matlab_ledoitwolf_covcor(z)
    x = copy(z)

    # de-mean returns
    t, n  = size(x)
    meanx = mean(x, dims=1)
    x -= repeat(meanx, outer=(t, 1))

    # compute sample covariance matrix
    sample = (1/t) * (x'*x)

    # compute prior
    diagsample = diag(sample)
    sqrtvar = sqrt.(diagsample)

    sqrtvar_mat = repeat(sqrtvar, outer=(1, n))

    rBar = (sum(sample./(sqrtvar_mat .* sqrtvar_mat'))-n)/(n*(n-1))
    prior = rBar * sqrtvar_mat .* sqrtvar_mat'

    #<#> replace the diagonal
    prior -= Diagonal(diag(prior))
    prior += Diagonal(diagsample)

    #<#> compute optimal shrinkage

    # what we call pi-hat
    y = x.^2;
    phiMat = y'*y/t - 2*(x'*x) .* sample/t + sample.^2
    phi = sum(phiMat)

    # what we call rho-hat
    term1    = ((x.^3)'*x)/t
    help     = x'*x/t
    helpDiag = diag(help)
    term2    = repeat(helpDiag, outer=(1, n)) .* sample
    term3    = help .* repeat(diagsample, outer=(1, n))
    term4    = repeat(diagsample, outer=(1, n)) .* sample
    thetaMat = term1 - term2 - term3 + term4
    #<#> remove the diagonal
    thetaMat -= Diagonal(diag(thetaMat))
    rho       = sum(diag(phiMat)) + rBar*sum(sqrtvar'./sqrtvar.*thetaMat)

    # what we call gamma hat
    gamma = sum((sample - prior).^2)
    kappa = (phi - rho) / gamma
    shrinkage = max(0, min(1, kappa/t))

    part_results = Dict(
        "r̄" => rBar,
        "F" => prior,
        "shrinkage" => shrinkage,
        "kappa" => kappa,
        "gamma" => gamma,
        "lwcov" => shrinkage*prior + (1-shrinkage)*sample)

    return part_results
end


text:
@doc raw"""
`PODE`: Partitioned Ordinary Differential Equation

Defines a partitioned initial value problem
```math
\begin{aligned}
\dot{q} (t) &= v(t, q(t), p(t)) , &
q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), p(t)) , &
p(t_{0}) &= p_{0} ,
\end{aligned}
```
with vector fields ``v`` and ``f``, initial conditions ``(q_{0}, p_{0})`` and the solution
``(q,p)`` taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}``.

### Parameters

* `vType <: Function`: type of `v`
* `fType <: Function`: type of `f`
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `v`: function computing the vector field ``v``
* `f`: function computing the vector field ``f``
* `invariants`: functions for the computation of invariants, either a `NamedTuple` containing the equation's invariants or `NullInvariants`
* `parameters`: type constraints for parameters, either a `NamedTuple` containing the equation's parameters or `NullParameters`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions, either a `AbstractArray` or `NullPeriodicity`

The functions `v` and `f` must have the interface
```julia
    function v(t, q, p, v, params)
        v[1] = ...
        v[2] = ...
        ...
    end
```
and
```julia
    function f(t, q, p, f, params)
        f[1] = ...
        f[2] = ...
        ...
    end
```
where `t` is the current time, `q` and `p` are the current solution vectors,
`v` and `f` are the vectors which hold the result of evaluating the vector
fields ``v`` and ``f`` on `t`, `q` and `p`, and params are additional parameters.

### Constructors

```julia
PODE(v, f, invariants, parameters, periodicity)
PODE(v, f; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity())
```

"""
struct PODE{vType <: Callable, fType <: Callable,
            invType <: OptionalInvariants,
            parType <: OptionalParameters,
            perType <: OptionalPeriodicity} <: AbstractEquationPODE{invType,parType,perType}

    v::vType
    f::fType

    invariants::invType
    parameters::parType
    periodicity::perType

    function PODE(v, f, invariants, parameters, periodicity)
        new{typeof(v), typeof(f), typeof(invariants), typeof(parameters), typeof(periodicity)}(
                v, f, invariants, parameters, periodicity)
    end
end

PODE(v, f; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = PODE(v, f, invariants, parameters, periodicity)

GeometricBase.invariants(equation::PODE) = equation.invariants
GeometricBase.parameters(equation::PODE) = equation.parameters
GeometricBase.periodicity(equation::PODE) = equation.periodicity

hasvectorfield(::PODE) = true

function check_initial_conditions(::PODE, ics::NamedTuple)
    haskey(ics, :q) || return false
    haskey(ics, :p) || return false
    eltype(ics.q) == eltype(ics.p) || return false
    typeof(ics.q) == typeof(ics.p) || return false
    axes(ics.q) == axes(ics.p) || return false
    return true
end

function check_methods(equ::PODE, tspan, ics, params)
    applicable(equ.v, tspan[begin], ics.q, ics.p, zero(ics.q), params) || return false
    applicable(equ.f, tspan[begin], ics.q, ics.p, zero(ics.p), params) || return false
    return true
end

function datatype(equ::PODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return eltype(ics.q)
end

function arrtype(equ::PODE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    typeof(ics.q)
end

_get_v(equ::PODE, params) = (t,q,p,v) -> equ.v(t, q, p, v, params)
_get_f(equ::PODE, params) = (t,q,p,f) -> equ.f(t, q, p, f, params)
_get_v̄(equ::PODE, params) = _get_v(equ, params)
_get_f̄(equ::PODE, params) = _get_f(equ, params)

_functions(equ::PODE) = (v = equ.v, f = equ.f)
_functions(equ::PODE, params::OptionalParameters) = (v = _get_v(equ, params), f = _get_f(equ, params))


text:
function ABCDQR_kfilter_s(y,u,A,B,C,D,Q,R,x10)
	#
	# STATIONARY Kalman filter for model
	#
	# x_{t+1} = A*x_{t} + B*u_{t} + w_{t}
	# y_{t}   = C*x_{t} + D*u_{t} + v_{t}
	#
	# cov(w_{t},v_{t}) = [Q 0;0 R]
	#
	# javier.cara@upm.es, 2015-10
	#

	(ny,nt) = size(y)
	nx = size(A,1)

	# allocation
	xtt = zeros(nx,nt)
	xtt1 = zeros(nx,nt+1)
	et = zeros(ny,nt)
	loglik = 0.0

	# Discrete-time Algebraic Riccati Equation (DARE) solution
	S = zeros(nx,ny)
	Ptt1 = dare(A',C',Q,R,S)

	# variance of innovations
	St = C*Ptt1*C' + R
	Stinv = inv(St)

	# Kalman gain matrix
	Kt = Ptt1*C'*Stinv

	# Var[x_t|y_t]
	Ptt = (Matrix(1.0I,nx,nx)-Kt*C)*Ptt1

	# Filter
	xtt1[:,1] = x10
	for t in 1:nt
		#  innovations
		et[:,t] = y[:,t] - C*xtt1[:,t] - D*u[:,t]

		# filtered values
		xtt[:,t] = xtt1[:,t] + Kt*et[:,t]

		# one-step ahead prediction
		xtt1[:,t+1] = A*xtt[:,t] + B*u[:,t]

		# likelihood
		l0 = et[:,t]'*Stinv*et[:,t] # typeof(l0) = Array{Float64,1}
		loglik = loglik + l0[1]
	end

	loglik = - ny*nt/2*log(2*pi) - nt/2*log(det(St)) - 1/2*loglik

	return xtt,Ptt,xtt1,Ptt1,et,St,Kt,loglik

end


text:
"""
Selectors indicate that index values are not indices, but points to
be selected from the dimension values, such as DateTime objects on a Time dimension.
"""
abstract type Selector{T} end

val(m::Selector) = m.val

"""
    At(x)

Selector that exactly matches the value on the passed-in dimensions, or throws an error.
For ranges and arrays, every value must match an existing value - not just the end points.
"""
struct At{T,A,R} <: Selector{T}
    val::T
    atol::A
    rtol::R
end
At(val::Union{Number,AbstractArray{<:Number},Tuple{<:Number,Vararg}};
   atol=zero(eltype(val)),
   rtol=(atol > zero(eltype(val)) ? zero(rtol) : Base.rtoldefault(eltype(val)))
  ) = At{typeof.((val, atol, rtol))...}(val, atol, rtol)
At(val; atol=nothing, rtol=nothing) =
    At{typeof.((val, atol, rtol))...}(val, atol, rtol)

atol(sel::At) = sel.atol
rtol(sel::At) = sel.rtol

"""
    Near(x)

Selector that selects the nearest index to its contained value(s)
"""
struct Near{T} <: Selector{T}
    val::T
end

"""
    Between(a, b)

Selector that retreive all indices located between 2 values.
"""
struct Between{T<:Union{Tuple{Any,Any},Nothing}} <: Selector{T}
    val::T
end
Between(args...) = Between{typeof(args)}(args)
Between(x::Tuple) = Between{typeof(x)}(x)

# Get the dims in the same order as the grid
# This would be called after RegularGrid and/or CategoricalGrid
# dimensions are removed
dims2indices(grid::TransformedGrid, dims::Tuple, lookups::Tuple, emptyval) =
    sel2indices(grid, dims, map(val, permutedims(dimz, dims(grid))))

sel2indices(A::AbstractArray, lookup) = sel2indices(dims(A), lookup)
sel2indices(dims::Tuple, lookup) = sel2indices(dims, (lookup,))
sel2indices(dims::Tuple, lookup::Tuple) = sel2indices(map(grid, dims), dims, lookup)
sel2indices(grids, dims::Tuple, lookup::Tuple) =
    (sel2indices(grids[1], dims[1], lookup[1]),
     sel2indices(tail(grids), tail(dims), tail(lookup))...)
    sel2indices(grids::Tuple{}, dims::Tuple{}, lookup::Tuple{}) = ()

# At selector
sel2indices(grid, dim::AbDim, sel::At) = at(dim, sel, val(sel))
sel2indices(grid, dim::AbDim, sel::At{<:Tuple}) =
    [at.(Ref(dim), Ref(sel), val(sel))...]
sel2indices(grid, dim::AbDim, sel::At{<:AbstractVector}) =
    at.(Ref(dim), Ref(sel), val(sel))

# Near selector
sel2indices(grid::T, dim::AbDim, sel::Near) where T<:Union{CategoricalGrid,UnknownGrid} =
    throw(ArgumentError("`Near` has no meaning in a `$T`. Use `At`"))
sel2indices(grid::AbstractAlignedGrid, dim::AbDim, sel::Near) =
    near(dim, val(sel))
sel2indices(grid::AbstractAlignedGrid, dim::AbDim, sel::Near{<:Tuple}) =
    [near.(Ref(dim), val(sel))...]
sel2indices(grid::AbstractAlignedGrid, dim::AbDim, sel::Near{<:AbstractVector}) =
    near.(Ref(dim), val(sel))

# Between selector
sel2indices(grid, dim::AbDim, sel::Between{<:Tuple}) = between(dim, val(sel))

# Transformed grid

# We use the transformation from the first TransformedGrid dim.
# In practice the others could be empty.
sel2indices(grids::Tuple{Vararg{<:TransformedGrid}}, dims::AbDimTuple,
            sel::Tuple{Vararg{<:Selector}}) =
    map(to_int, sel, val(dims[1])([map(val, sel)...]))

to_int(::At, x) = convert(Int, x)
to_int(::Near, x) = round(Int, x)

# Do the input values need some kind of scalar conversion?
# what is the scale of these lookup matrices?
sel2indices(grid::LookupGrid, sel::Tuple{Vararg{At}}) =
    lookup(grid)[map(val, sel)...]

at(dim::AbDim, sel::At, val) =
    relate(relationorder(dim), dim, at(dim, val, atol(sel), rtol(sel)))
at(dim::AbDim, selval, atol::Nothing, rtol::Nothing) = begin
    ind = findfirst(x -> x == selval, val(dim))
    ind == nothing ? throw(ArgumentError("$selval not found in $dim")) : ind
end
at(dim::AbDim, selval, atol, rtol) = begin
    # This is not particularly efficient.
    # It should be separated out for unordered
    # dims and otherwise treated as an ordered list.
    ind = findfirst(x -> isapprox(x, selval; atol=atol, rtol=rtol), val(dim))
    ind == nothing ? throw(ArgumentError("$selval not found in $dim")) : ind
end

near(dim::AbDim, selval) =
    relate(relationorder(dim), dim, near(indexorder(dim), dim::AbDim, selval))
near(indexorder::Unordered, dim, selval) =
    throw(ArgumentError("`Near` has no meaning in an `Unordered` grid"))
near(indexorder, dim, selval) = begin
    index = val(dim)
    i = searchsortedfirst(index, selval; rev=isrev(indexorder))
    if i > lastindex(index)
        lastindex(index)
    elseif i <= firstindex(index)
        firstindex(index)
    elseif abs(index[i] - selval) < abs(index[i-1] - selval)
        i
    else
        i - 1
    end
end

between(dim::AbDim, sel) =
    relate(relationorder(dim), dim, between(indexorder(dim), dim, sel))
between(::Unordered, dim::AbDim, sel) =
    throw(ArgumentError("Cannot use `Between` on an unordered dimension"))
between(indexorder::Forward, dim::AbDim, sel) =
    rangeorder(indexorder, searchsortedfirst(val(dim), first(sel)), searchsortedlast(val(dim), last(sel)))
between(indexorder::Reverse, dim::AbDim, sel) =
    rangeorder(indexorder, searchsortedfirst(val(dim), last(sel); rev=true),
                    searchsortedlast(val(dim), first(sel); rev=true))

rangeorder(::Forward, lower, upper) = lower:upper
rangeorder(::Reverse, lower, upper) = upper:-1:lower

relate(::Forward, dim, val) = val
relate(::Reverse, dim, val) = lastindex(dim) .- val .+ 1

# Selector indexing without dim wrappers. Must be in the right order!
Base.@propagate_inbounds Base.getindex(a::AbstractArray, I::Vararg{Selector}) =
    getindex(a, sel2indices(a, I)...)
Base.@propagate_inbounds Base.setindex!(a::AbstractArray, x, I::Vararg{Selector}) =
    setindex!(a, x, sel2indices(a, I)...)
Base.view(a::AbstractArray, I::Vararg{Selector}) =
    view(a, sel2indices(a, I)...)


text:
include("Files.jl")
include("Base.jl")
include("Edit.jl")
include("Flux.jl")


text:
@testset "init" begin
    withenv("OCTAVIAN_NUM_TASKS" => "") do
        @test Octavian._read_environment_num_tasks() == min(Threads.nthreads(), VectorizationBase.num_cores())
    end
    withenv("OCTAVIAN_NUM_TASKS" => "1") do
        @test Octavian._read_environment_num_tasks() == 1
    end
    withenv("OCTAVIAN_NUM_TASKS" => "99") do
        @test Octavian._read_environment_num_tasks() == min(Threads.nthreads(), VectorizationBase.num_cores())
    end
end


text:
using InteractiveDynamics, GLMakie, DynamicalSystems

# the second range is a convenience for intermittency example of logistic
rrange = 1:0.001:4.0
# rrange = (rc = 1 + sqrt(8); [rc, rc - 1e-5, rc - 1e-3])

lo = Systems.logistic(0.4; r=rrange[1])

interactive_cobweb(lo, rrange, 5)


text:
module TestFieldDocTables

using FieldDocTables, FieldMetadata, PrettyTables, Test, Markdown

import FieldMetadata: @default, default, @bounds, bounds, @description, description

const FIELDDOCTABLE = FieldDocTable((Description=description, Default=default, Bounds=bounds))

"""
This type tests if FIELDMETADATA is printed as a markdown table.

$(FIELDDOCTABLE)
"""
@description @bounds @default mutable struct TestStruct
   " A field doc"
   a::Int     | 2   | (1, 10)     | "an Int"
   b::Float64 | 4.0 | (2.0, 20.0) | "a Float"
end

@test Markdown.plain(@doc TestStruct) == "This type tests if FIELDMETADATA is printed as a markdown table.\n\n| Field | Description | Default |      Bounds |        Docs |\n| -----:| -----------:| -------:| -----------:| -----------:|\n|     a |      an Int |       2 |     (1, 10) | A field doc |\n|     b |     a Float |     4.0 | (2.0, 20.0) |             |\n"

end