text:
struct IsopycnalSkewSymmetricDiffusivity{K, S, M, L} <: AbstractTurbulenceClosure{ExplicitTimeDiscretization}
                    κ_skew :: K
               κ_symmetric :: S
          isopycnal_tensor :: M
             slope_limiter :: L

    function IsopycnalSkewSymmetricDiffusivity(κ_skew::K, κ_symmetric::S, isopycnal_tensor::I, slope_limiter::L) where {K, S, I, L}

        isopycnal_tensor isa SmallSlopeIsopycnalTensor ||
            error("Only isopycnal_tensor=SmallSlopeIsopycnalTensor() is currently supported.")

        return new{K, S, I, L}(κ_skew, κ_symmetric, isopycnal_tensor, slope_limiter)
    end
end

const ISSD = IsopycnalSkewSymmetricDiffusivity
const ISSDVector = AbstractVector{<:ISSD}
const FlavorOfISSD = Union{ISSD, ISSDVector}
const issd_coefficient_loc = (Center, Center, Center)

"""
    IsopycnalSkewSymmetricDiffusivity([FT=Float64;]
                                      κ_skew = 0,
                                      κ_symmetric = 0,
                                      isopycnal_tensor = SmallSlopeIsopycnalTensor(),
                                      slope_limiter = nothing)

Return parameters for an isopycnal skew-symmetric tracer diffusivity with skew diffusivity
`κ_skew` and symmetric diffusivity `κ_symmetric` that uses an `isopycnal_tensor` model for
for calculating the isopycnal slopes, and (optionally) applying a `slope_limiter` to the
calculated isopycnal slope values.

Both `κ_skew` and `κ_symmetric` may be constants, arrays, fields, or functions of `(x, y, z, t)`.
"""
IsopycnalSkewSymmetricDiffusivity(FT=Float64; κ_skew=0, κ_symmetric=0, isopycnal_tensor=SmallSlopeIsopycnalTensor(), slope_limiter=nothing) =
    IsopycnalSkewSymmetricDiffusivity(convert_diffusivity(FT, κ_skew), convert_diffusivity(FT, κ_symmetric), isopycnal_tensor, slope_limiter)

function with_tracers(tracers, closure::ISSD)
    κ_skew = !isa(closure.κ_skew, NamedTuple) ? closure.κ_skew : tracer_diffusivities(tracers, closure.κ_skew)
    κ_symmetric = !isa(closure.κ_symmetric, NamedTuple) ? closure.κ_symmetric : tracer_diffusivities(tracers, closure.κ_symmetric)
    return IsopycnalSkewSymmetricDiffusivity(κ_skew, κ_symmetric, closure.isopycnal_tensor, closure.slope_limiter)
end

# For ensembles of closures
function with_tracers(tracers, closure_vector::ISSDVector)
    arch = architecture(closure_vector)

    if arch isa Architectures.GPU
        closure_vector = Vector(closure_vector)
    end

    Ex = length(closure_vector)
    closure_vector = [with_tracers(tracers, closure_vector[i]) for i=1:Ex]

    return arch_array(arch, closure_vector)
end

#####
##### Tapering
#####

struct FluxTapering{FT}
    max_slope :: FT
end

"""
    taper_factor_ccc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, tapering::FluxTapering)

Return the tapering factor `min(1, Sₘₐₓ² / slope²)`, where `slope² = slope_x² + slope_y²`
that multiplies all components of the isopycnal slope tensor. All slopes involved in the
tapering factor are computed at the cell centers.

References
==========
R. Gerdes, C. Koberle, and J. Willebrand. (1991), "The influence of numerical advection schemes
    on the results of ocean general circulation models", Clim. Dynamics, 5 (4), 211–226.
"""
@inline function taper_factor_ccc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, tapering::FluxTapering) where FT
    # TODO: handle boundaries!
    bx = ℑxᶜᵃᵃ(i, j, k, grid, ∂x_b, buoyancy, tracers)
    by = ℑyᵃᶜᵃ(i, j, k, grid, ∂y_b, buoyancy, tracers)
    bz = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers)

    slope_x = - bx / bz
    slope_y = - by / bz
    slope² = ifelse(bz <= 0, zero(FT), slope_x^2 + slope_y^2)

    return min(one(FT), tapering.max_slope^2 / slope²)
end

"""
    taper_factor_ccc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, ::Nothing) where FT

Returns 1 for the  isopycnal slope tapering factor, that is, no tapering is done.
"""
taper_factor_ccc(i, j, k, grid::AbstractGrid{FT}, buoyancy, tracers, ::Nothing) where FT = one(FT)

# Diffusive fluxes

@inline get_tracer_κ(κ::NamedTuple, tracer_index) = @inbounds κ[tracer_index]
@inline get_tracer_κ(κ, tracer_index) = κ

# defined at fcc
@inline function diffusive_flux_x(i, j, k, grid,
                                  closure::Union{ISSD, ISSDVector}, diffusivity_fields, ::Val{tracer_index},
                                  velocities, tracers, clock, buoyancy) where tracer_index

    c = tracers[tracer_index]
    closure = getclosure(i, j, closure)

    κ_skew = get_tracer_κ(closure.κ_skew, tracer_index)
    κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index)

    κ_skewᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, clock, issd_coefficient_loc, κ_skew)
    κ_symmetricᶠᶜᶜ = κᶠᶜᶜ(i, j, k, grid, clock, issd_coefficient_loc, κ_symmetric)

    ∂x_c = ∂xᶠᶜᶜ(i, j, k, grid, c)

    # Average... of... the gradient!
    ∂y_c = ℑxyᶠᶜᵃ(i, j, k, grid, ∂yᶜᶠᶜ, c)
    ∂z_c = ℑxzᶠᵃᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c)

    R₁₁ = one(eltype(grid))
    R₁₂ = zero(eltype(grid))
    R₁₃ = isopycnal_rotation_tensor_xz_fcc(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor)

    ϵ = taper_factor_ccc(i, j, k, grid, buoyancy, tracers, closure.slope_limiter)

    return - ϵ * (              κ_symmetricᶠᶜᶜ * R₁₁ * ∂x_c +
                                κ_symmetricᶠᶜᶜ * R₁₂ * ∂y_c +
                  (κ_symmetricᶠᶜᶜ - κ_skewᶠᶜᶜ) * R₁₃ * ∂z_c)
end

# defined at cfc
@inline function diffusive_flux_y(i, j, k, grid,
                                  closure::Union{ISSD, ISSDVector}, diffusivity_fields, ::Val{tracer_index},
                                  velocities, tracers, clock, buoyancy) where tracer_index

    c = tracers[tracer_index]
    closure = getclosure(i, j, closure)

    κ_skew = get_tracer_κ(closure.κ_skew, tracer_index)
    κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index)

    κ_skewᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, clock, issd_coefficient_loc, κ_skew)
    κ_symmetricᶜᶠᶜ = κᶜᶠᶜ(i, j, k, grid, clock, issd_coefficient_loc, κ_symmetric)

    ∂y_c = ∂yᶜᶠᶜ(i, j, k, grid, c)

    # Average... of... the gradient!
    ∂x_c = ℑxyᶜᶠᵃ(i, j, k, grid, ∂xᶠᶜᶜ, c)
    ∂z_c = ℑyzᵃᶠᶜ(i, j, k, grid, ∂zᶜᶜᶠ, c)

    R₂₁ = zero(eltype(grid))
    R₂₂ = one(eltype(grid))
    R₂₃ = isopycnal_rotation_tensor_yz_cfc(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor)

    ϵ = taper_factor_ccc(i, j, k, grid, buoyancy, tracers, closure.slope_limiter)

    return - ϵ * (              κ_symmetricᶜᶠᶜ * R₂₁ * ∂x_c +
                                κ_symmetricᶜᶠᶜ * R₂₂ * ∂y_c +
                  (κ_symmetricᶜᶠᶜ - κ_skewᶜᶠᶜ) * R₂₃ * ∂z_c)
end

# defined at ccf
@inline function diffusive_flux_z(i, j, k, grid,
                                  closure::Union{ISSD, ISSDVector}, diffusivity_fields, ::Val{tracer_index},
                                  velocities, tracers, clock, buoyancy) where tracer_index

    c = tracers[tracer_index]
    closure = getclosure(i, j, closure)

    κ_skew = get_tracer_κ(closure.κ_skew, tracer_index)
    κ_symmetric = get_tracer_κ(closure.κ_symmetric, tracer_index)

    κ_skewᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, clock, issd_coefficient_loc,κ_skew)
    κ_symmetricᶜᶜᶠ = κᶜᶜᶠ(i, j, k, grid, clock, issd_coefficient_loc, κ_symmetric)

    ∂z_c = ∂zᶜᶜᶠ(i, j, k, grid, c)

    # Average... of... the gradient!
    ∂x_c = ℑxzᶜᵃᶠ(i, j, k, grid, ∂xᶠᶜᶜ, c)
    ∂y_c = ℑyzᵃᶜᶠ(i, j, k, grid, ∂yᶜᶠᶜ, c)

    R₃₁ = isopycnal_rotation_tensor_xz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor)
    R₃₂ = isopycnal_rotation_tensor_yz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor)
    R₃₃ = isopycnal_rotation_tensor_zz_ccf(i, j, k, grid, buoyancy, tracers, closure.isopycnal_tensor)

    ϵ = taper_factor_ccc(i, j, k, grid, buoyancy, tracers, closure.slope_limiter)

    return - ϵ * ((κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₁ * ∂x_c +
                  (κ_symmetricᶜᶜᶠ + κ_skewᶜᶜᶠ) * R₃₂ * ∂y_c +
                                κ_symmetricᶜᶜᶠ * R₃₃ * ∂z_c)
end

@inline viscous_flux_ux(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_uy(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_uz(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))

@inline viscous_flux_vx(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_vy(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_vz(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))

@inline viscous_flux_wx(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_wy(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))
@inline viscous_flux_wz(i, j, k, grid, closure::Union{ISSD, ISSDVector}, args...) = zero(eltype(grid))

calculate_diffusivities!(diffusivity_fields, closure::Union{ISSD, ISSDVector}, model) = nothing

#####
##### Show
#####

Base.show(io::IO, closure::ISSD) =
    print(io, "IsopycnalSkewSymmetricDiffusivity: " *
              "(κ_symmetric=$(closure.κ_symmetric), κ_skew=$(closure.κ_skew), " *
              "(isopycnal_tensor=$(closure.isopycnal_tensor), slope_limiter=$(closure.slope_limiter))")


text:
using Catlab
using Catlab.Theories
using Catlab.Syntax
using Metatheory, Metatheory.EGraphs

# ============================================================

# GATExpr => normal Expr in MT
function gat_to_expr(ex::ObExpr{:generator})
    @assert length(ex.args) == 1
    return ex.args[1]
end
function gat_to_expr(ex::ObExpr{H}) where {H}
    return Expr(:call, head(ex), map(gat_to_expr, ex.args)...)
end
function gat_to_expr(ex::HomExpr{H}) where {H}
    @assert length(ex.type_args) == 2
    expr = Expr(:call, head(ex), map(gat_to_expr, ex.args)...)
    type_ex = Expr(:call, :Hom, map(gat_to_expr, ex.type_args)...)
    return Expr(:call, :~, expr, type_ex)
end
function gat_to_expr(ex::HomExpr{:generator})
    f = ex.args[1]
    type_ex = Expr(:call, :Hom, map(gat_to_expr, ex.type_args)...)
    return Expr(:call, :~, f, type_ex)
end

const Code = Union{Symbol, Expr}
const TTags = Dict{Code, Tuple{Code, Symbol}}

# ============================================================

# infer type of morphisms and objects
# a morphism f: A → B will be typed as f ~ Hom(A, B)
# an object A will be typed as Ob(A)
function get_concrete_type_expr(theory, x::Symbol, ctx, type_tags = TTags())
    t = ctx[x]
    @show(t)
    # t === :Ob && (t = Expr(:call, :Ob, x))
    if t === :Ob
        type_tags[x] = (x, t)
        return (x, t)
    else
        @assert t.args[1] == :Hom
        type_tags[x] = (t, t.args[1])
        return (t, t.args[1])
    end
end

function get_concrete_type_expr(theory, x::Expr, ctx, type_tags = TTags())
    @assert exprhead(x) == :call
    f = x.args[1]
    rest = x.args[2:end]
    # recursion case - inductive step (?)
    for a in rest
        (t, sort) = get_concrete_type_expr(theory, a, ctx, type_tags)
        type_tags[a] = (t, sort)
        println("$a ~ $t")
    end
    # get the corresponding TermConstructor from theory.terms
    # for each arg in `rest`, instantiate the term.params with term.context
    # instantiate term.typ

    (t, sort) = gat_type_inference(theory, f, [type_tags[a] for a in rest])
    type_tags[x] = (t, sort)
    # println("$x ~ $(type_tags[x])")
    return (t, sort)
end

function is_context_match(t, head, args)
    # t isa TermConstructor
    # println(repeat("=", 30))
    # println("is_context_match")
    # @show t
    # @show head
    # @show args
    # println(repeat("=", 30))

    # TODO fixme!

    t.name !== head && return false
    n = length(t.params)
    n != length(args) && return false
    for i = 1:n
        arg, sort = args[i]

        if t.context[t.params[i]] === :Ob
            if sort !== :Ob
                return false
            end
        else
            if sort === :Ob
                return false
            end
        end
    end
    return true
end

function gat_type_inference(theory, head, args)
    for t in theory.terms
        if is_context_match(t, head, args)
            # @show t, head, args
            return gat_type_inference(t, head, args)
        end
    end
    # @show theory, head, args
    @error "can not find $(Expr(:call, head, args...)) in the theory"
end

function gat_type_inference(t::GAT.TermConstructor, head, args)
    @assert length(t.params) == length(args) && t.name === head
    bindings = Dict()

    println(args)
    texprs = map(first, args)
    sorts = map(last, args)

    for i = 1:length(args)
        template = t.context[t.params[i]]
        template === :Ob && (template = t.params[i])
        # @show template
        update_bindings!(bindings, template, texprs[i])
    end
    # @show bindings
    r = GAT.replace_types(bindings, t)
    if r.typ == :Ob
        return Expr(:call, head, texprs...), r.typ
        #     # return Expr(:call, :Ob, Expr(:call, head, args...))
    #     Expr(:call, head, args...)
    else
        @show(r.typ)
        return r.typ, r.typ.args[1]
    end
    # end
end
function update_bindings!(bindings, template::Expr, target::Expr)
    for i = 1:length(template.args)
        update_bindings!(bindings, template.args[i], target.args[i])
    end
end
function update_bindings!(bindings, template, target)
    bindings[template] = target
end

function tag_expr(x::Expr, axiom, theory)
    texpr, sort = get_concrete_type_expr(theory, x, axiom.context)
    start = exprhead(x) == :call ? 2 : 1

    nargs = Any[tag_expr(y, axiom, theory) for y in x.args[start:end]]

    if start == 2
        pushfirst!(nargs, x.args[1])
    end

    z = Expr(exprhead(x), nargs...)

    (sort === :Ob) && (return z)
    :($z ~ $texpr)
end

function tag_expr(x::Symbol, axiom, theory)
    (texpr, sort) = get_concrete_type_expr(theory, x, axiom.context)
    (sort === :Ob) && (return x)
    # return (t == x ? x : :($x ~ $t))
    return :($x ~ $texpr)
end

# ============================================================
# Convert Catlab Axioms to rules
# ============================================================

function axiom_to_rule(theory, axiom)
    op = axiom.name
    @assert op == :(==)
    lhs = tag_expr(axiom.left, axiom, tt) |> Pattern
    rhs = tag_expr(axiom.right, axiom, tt) |> Pattern

    pvars = patvars(lhs) ∪ patvars(rhs)
    extravars = setdiff(pvars, patvars(lhs) ∩ patvars(rhs))
    if !isempty(extravars)
        if extravars ⊆ patvars(lhs)
            println(lhs)
            println(rhs)
            return RewriteRule(lhs, rhs)
        else
            return RewriteRule(rhs, lhs)
        end
    end
    # println("$lhs == $rhs")
    EqualityRule(lhs, rhs)
end

function gen_theory(t::Catlab.GAT.Theory)
    [axiom_to_rule(t, ax) for ax in t.axioms]
end

# =========================================================
# Utility Functions
# =========================================================

function Base.show(io::IO, a::Catlab.GAT.AxiomConstructor)
    print(io, a.left)
    print(io, ' ', a.name, ' ' )
    print(io, a.right)
    print(io, " where ")
    n = length(a.context)
    ctx = collect(a.context)
    for i in 1:n
        (k,v) = ctx[i]
        print(io, "$k => $v")
        if i !== n
            print(io, ", ")
        end
    end
end

ax

ext: jl


text:
using SnoopCompile, Pkg
project_path = dirname(dirname(@__FILE__))
rm(joinpath(project_path, "src", "precompile.jl"); force=true)
open(joinpath(project_path, "src", "precompile.jl"), "w+") do io
    println(io, "_precompile_() = nothing")
end

### Log the compiles
# This only needs to be run once (to generate "/tmp/comonicon_compiles.log")

SnoopCompile.@snoopc ["--project=$(project_path)"] "/tmp/ion_compiles.log" begin
    using IonBase, Comonicon

    cd(tempdir()) do
        rm("Foo"; recursive=true, force=true)
        IonBase.create("Foo"; user="xyz")
        rm("Foo"; recursive=true, force=true)
        IonBase.create("Foo"; user="xyz", template="comonicon")
        rm("Foo"; recursive=true, force=true)
    end

    include(Comonicon.PATH.project(IonBase, "test", "runtests.jl"))
    IonBase.search("Yao")
    IonBase.Doc.build(Comonicon.PATH.project(IonBase))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts

data = SnoopCompile.read("/tmp/ion_compiles.log")
pc = SnoopCompile.parcel(reverse!(data[2]))

open(joinpath(project_path, "src", "precompile.jl"), "w") do io
    if any(str->occursin("__lookup", str), pc[:IonBase])
        println(io, lookup_kwbody_str)
    end
    println(io, "function _precompile_()")
    println(io, "    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing")
    for ln in sort(pc[:IonBase])
        println(io, "    ", ln)
    end
    println(io, "end")
end


text:
const PARSED_SECTION = [
    "Intro",
    "Introduction", # long description
    "Arguments",
    "Args", # arguments
    "Options",
    "Flags", # kwargs
]

function split_docstring(f::Function)
    doc = Base.Docs.doc(f)
    return split_docstring(doc)
end

function split_docstring(m::Module)
    doc = Base.Docs.doc(m)
    has_docstring(doc) || return "", ""
    return read_description(doc), read_intro(doc)
end

function split_docstring(doc::Markdown.MD)
    has_docstring(doc) || return JLMD()
    desc = read_description(doc)
    intro = read_intro(doc)
    args = read_arguments(doc)
    flags = read_flags(doc)
    options = read_options(doc)
    return JLMD(desc, intro, args, options, flags)
end

function has_docstring(doc::Markdown.MD)
    paragraph = first(read_content(doc))
    flag = paragraph isa Markdown.Paragraph && paragraph.content == Any["No documentation found."]
    return !flag
end

function read_intro(md::Markdown.MD)
    sec = read_section(md, ["Intro", "Introduction"])
    isempty(sec) && return ""
    length(sec) == 1 || error("section Intro/Introduction can only have one paragraph")
    return strip(md_to_string(Markdown.MD(sec, md.meta)))
end

function read_arguments(md::Markdown.MD)
    args = Dict{String,String}()
    sec = read_section(md, ["Arguments", "Args"])
    isempty(sec) && return args
    length(sec) == 1 || error("section Arguments/Args can only have one paragraph")
    sec = sec[]

    for each in sec.items
        name, doc = read_item(each[1])
        args[name] = doc
    end
    return args
end

function read_flags(md::Markdown.MD)
    flags = Dict{String,JLMDFlag}()
    sec = read_section(md, "Flags")
    isempty(sec) && return flags
    length(sec) == 1 || error("section Flags can only have one paragraph")
    sec = sec[]

    for each in sec.items
        name, doc = read_item(each[1])
        name, short, hint = split_option(name)
        hint === nothing || error("flag cannot have hint")
        flags[name] = JLMDFlag(doc, short !== nothing)
    end
    return flags
end

function read_options(md::Markdown.MD)
    options = Dict{String,JLMDOption}()
    sec = read_section(md, "Options")
    isempty(sec) && return options
    length(sec) == 1 || error("section Options can only have one paragraph")
    sec = sec[]

    for each in sec.items
        code, doc = read_item(each[1])
        name, short, hint = split_option(code)
        options[name] = JLMDOption(hint, doc, !isnothing(short))
    end
    return options
end

function read_description(md::Markdown.MD)
    intro = []
    lines = read_content(md)
    # ignore julia function signature
    if lines[1] isa Markdown.Code
        lines = lines[2:end]
    end

    for line in lines
        if line isa Markdown.Header{1} && line.text[1] in PARSED_SECTION
            break
        else
            push!(intro, line)
        end
    end
    isempty(intro) && return ""
    return strip(md_to_string(Markdown.MD(intro, md.meta)))
end

function read_section(md::Markdown.MD, title::Vector{String})
    for each in title
        sec = read_section(md, each)
        isempty(sec) || return sec
    end
    return []
end

function read_section(md::Markdown.MD, title::String)::Vector{Any}
    ct = read_content(md)
    nlines = length(ct)
    content = []
    for k in 1:nlines
        line = ct[k]
        if line isa Markdown.Header{1} && line.text[1] == title
            k + 1 ≤ nlines || return content # return on last line
            for idx in k+1:nlines
                ct[idx] isa Markdown.Header{1} && return content # another title
                push!(content, ct[idx])
            end
        end
    end
    return content
end

read_content(x) = x
function read_content(md::Markdown.MD)
    if md.content[1] isa Markdown.MD
        return read_content(md.content[1])
    else
        return md.content
    end
end

function read_item(raw::Markdown.Paragraph)
    raw.content[1] isa Markdown.Code ||
        throw(Meta.ParseError("command argument name should be marked by inline code"))
    name = raw.content[1].code

    raw_doc = md_to_string(Markdown.MD(Markdown.Paragraph(raw.content[2:end])))
    doc = read_docstring(raw_doc)
    return name, doc
end

function read_docstring(doc::String)
    m = match(r"^: *(.*)", strip(doc))
    m === nothing && throw(Meta.ParseError("invalid docstring format: $doc"))
    return String(m[1])
end

"""
    rm_format(md)

Remove Markdown DOM and flatten to strings.
"""
function rm_format end
rm_format(x::Markdown.Paragraph) = join(map(rm_format, x.content))
rm_format(x::Markdown.Code) = x.code
rm_format(x::String) = x
rm_format(x::Markdown.MD) = rm_format(x.content[1])

function docstring(x)
    return sprint(Base.Docs.doc(x); context = :color => true) do io, x
        show(io, MIME"text/plain"(), x)
    end
end

function md_to_string(md::Markdown.MD) # don't print newlines for one paragraph
    return sprint(md; context = (:color => true, :displaysize => (typemax(Int), 1000))) do io, x
        show(io, MIME"text/plain"(), x)
    end
end

function split_option(content::String)
    content = strip(content)
    startswith(content, "-") || throw(Meta.ParseError("expect --option[,-o], got $content"))
    names = map(strip, split(content, ",")) # rm space

    if length(names) == 1 # long option
        long, hint = split_hint(lstrip(names[1], '-'))
        short = nothing
    elseif length(names) == 2 # short option
        if startswith(names[1], "--") && !startswith(names[2], "--") && startswith(names[2], '-') # --option, -o
            long, (short, hint) = lstrip(names[1], '-'), split_hint(lstrip(names[2], '-'))
        elseif startswith(names[2], "--") && !startswith(names[1], "--") && startswith(names[1], '-') # -o, --option
            (long, hint), short = split_hint(lstrip(names[2], '-')), lstrip(names[1], '-')
        else
            error("expect --option[, -o], got $content")
        end

        length(short) == 1 ||
            throw(Meta.ParseError("short option can only use one letter, got $content"))
        first(short) == first(long) ||
            throw(Meta.ParseError("short option must use the same first letter, got $content"))
    else
        throw(Meta.ParseError("too much inputs, expect --option[,-o], got $content"))
    end
    long = replace(long, '_' => '-')
    return long, short, hint
end

function split_hint(content::AbstractString)
    msg = "expect --option[,-o]=<hint> or --option[,-o] <hint>, got $content"
    content = strip(content)
    m = match(r"([^\s]+)(?:\s+|=)(<.+>)", content)

    # no hint, just -o, --option
    # shouldn't contain inner space
    if m === nothing
        any(x -> isspace(x) || isequal(x, '='), content) && throw(Meta.ParseError(msg))
        return content, nothing
    end

    return m[1], strip(m[2], ['<', '>'])
end


text:
export integral,
    integrate

"""
The suggested quadrature strategy based on the `integral` arguments.

By default, adaptive quadrature is chosen.
"""
suggestedstrategy(args...) = QuadAdaptive()

returntype(integrand, ::Type{S}) where {S} = Base.Core.Compiler.return_type(integrand, (S,))

zero_error(x) = zero_error(typeof(x))
zero_error(::Type{T}) where {T} = zero(prectype(T))
unknown_error(T) = -one(prectype(T))
unknown_error(z::Number) = -one(z)

tolerance(d::Domain{T}) where {T} = tolerance(T)
tolerance(::Type{T}) where {T<:AbstractFloat} = 10eps(T)
tolerance(::Type{T}) where {T} = tolerance(prectype(T))

"""
Compute an integral of the given integrand on the given domain,
or with the given measure.

Example:
```
integral(cos, 0.0..1.0)
```
"""
function integral end

"Like integral, but also returns an error estimate (if applicable)."
function integrate end

function integral(args...)
    I, E = integrate(args...)
    I
end

# associate a domain with a measure: use its support
associated_domain(μ::Measure) = support(μ)

# associate a measure with a domain
# we try to avoid memory allocations, hence LebesgeeSpace{T} is the default
associated_measure(domain::Domain{T}) where {T} = Lebesgue{T}()
# these cases have a known allocation-free lebesgue measure
associated_measure(domain::ChebyshevInterval) = lebesguemeasure(domain)
associated_measure(domain::UnitInterval) = lebesguemeasure(domain)
# sometimes the T of an interval is an integer (e.g. in 0..1)
associated_measure(domain::AbstractInterval{T}) where {T} = Lebesgue{float(T)}()

# Process the arguments until there is a domain, a measure and zero or more property objects.
process_arguments(measure::Measure, properties::Property...) =
    process_arguments(associated_domain(measure), measure, properties...)
process_arguments(domain::Domain, properties::Property...) =
    process_arguments(domain, associated_measure(domain), properties...)
process_arguments(domain::Domain, measure::LebesgueDomain{T}, properties::Property...) where {T} =
    process_arguments(domain ∩ support(measure), Lebesgue{T}(), properties...)

process_arguments(domain::Domain{S}, measure::Measure{T}, properties::Property...) where {S,T} =
    process_arguments(convert(Domain{promote_type(S,T)}, domain),
        convert(Measure{promote_type(S,T)}, measure), properties...)

# all good now
process_arguments(domain::Domain{T}, measure::Measure{T}, properties::Property...) where {T} =
    (domain, measure, properties...)

process_arguments(args...) = error("Arguments to integral or integrate functions not understood.")

integrate(qs::ChebyshevIntervalRule{T}, integrand) where {T} =
    integrate(qs, integrand, ChebyshevInterval{T}())
integrate(qs::HalfLineRule{T}, integrand) where {T} =
    integrate(qs, integrand, HalfLine{T}())
integrate(qs::RealLineRule{T}, integrand) where {T} =
    integrate(qs, integrand, FullSpace{T}())

## The use of generators
integrate(gen::Base.Generator, args...) =
    integrate(process_generator(gen)..., args...)
# - catch something like f(x) for x in domain
process_generator(gen::Base.Generator{<:Domain}) = (gen.f, gen.iter)
# - catch something like f(x,y) for x in domain1, y in domain2
process_generator(gen::Base.Generator{<:Base.Iterators.ProductIterator}) =
    process_generator(gen, gen.iter.iterators)
function process_generator(gen, iterators::Tuple{Vararg{Domain}})
    domain = productdomain(iterators...)
    dims = map(dimension, iterators)
    (gen.f, domain)
end

integrate(integrand, args...) =
    integrate(integrand, process_arguments(args...)...)

integrate(qs::QuadratureStrategy, integrand, args...) =
    integrate_start(qs, integrand, process_arguments(args...)...)

integrate(integrand, domain::Domain, measure::Measure, properties::Property...) =
    integrate(suggestedstrategy(domain, measure, properties...), integrand, domain, measure, properties...)

# The process is as follows:
# - The argument list is completed (missing domain or measure are added)
#   -> integrate_property is invoked
# - integrate_property processes properties using dispatch
#   -> integrate_measure is invoked
# - integrate_measure processes measures using dispatch
#   -> integrate_domain is invoked
# - integrate_domain processes domains using dispatch
#   -> apply_quad is invoked
# - If apply_quad is not intercepted, a fallback routine is invoked.

integrate_start(qs, integrand, domain, measure, properties...) =
    integrate_property(qs, Integrand{promote_type(numtype(domain),codomaintype(measure))}(integrand), domain, measure, properties...)

integrate_start(qs, integrand::AbstractIntegrand, domain, measure, properties...) =
    integrate_property(qs, integrand, domain, measure, properties...)

function sum_integrals(qs, integrand::AbstractIntegrand, domains, measure, properties...)
    Itot, Etot = zero_result(integrand)
    for domain in domains
        I,E = integrate_start(qs, integrand, domain, measure, filter_singularities(properties...)...)
        Itot += I
        Etot += I
    end
    Itot, Etot
end

# Step 1: process any singularities in the list of properties
function integrate_property(qs::AdaptiveStrategy, integrand, domain, measure, property, properties...)
    if property_splits_domain(property, properties...)
        domains = property_split_domain(domain, property, properties...)
        sum_integrals(qs, integrand, domains, measure, filter_singularities(property, properties...)...)
    else
        transformed_arguments = process_property(qs, integrand, domain, measure, property, properties...)
        integrate_measure(qs, transformed_arguments...)
    end
end

integrate_property(qs, integrand, domain, measure, properties...) =
    integrate_measure(qs, integrand, domain, measure, properties...)

# STEP M: dispatch on the measure
function integrate_measure(qs, integrand, domain, measure, properties...)
    transformed_arguments = process_measure(qs, integrand, domain, measure, properties...)
    integrate_domain(qs, transformed_arguments...)
end

# STEP D: dispatch on the domain
function integrate_domain(qs, integrand, domain, measure, properties...)
    if domain_splits(domain)
        domains = domain_split(domain)
        sum_integrals(qs, integrand, domains, measure, properties...)
    else
        transformed_arguments = process_domain(qs, integrand, domain, measure, properties...)
        integrate_done(qs, transformed_arguments...)
    end
end

integrate_done(qs, integrand::Integrand, domain, measure, properties...) =
    integrate_done(qs, integrand.fun, domain, measure, properties...)
integrate_done(qs, integrand::Function, domain, measure, properties...) =
    select_quad(qs, integrand, domain, measure, properties...)

# Selection: use a suitable quadrature strategy

# Generic adaptive quadrature in 1D invokes QuadGK, everywhere else it uses hcubature
select_quad(qs::QuadAdaptive, integrand, domain::Domain{T}, measure, properties...) where {T<:Number} =
    apply_quad(Q_quadgk(qs), integrand, domain, measure, properties...)
select_quad(qs::QuadAdaptive, integrand, domain::Domain{T}, measure, properties...) where {T} =
    apply_quad(Q_hcubature(qs), integrand, domain, measure, properties...)
select_quad(qs, integrand, domain, measure, properties...) =
    apply_quad(qs, integrand, domain, measure, properties...)

# FINAL STEP: invoke apply_quad, with a fallback if necessary
apply_quad(qs, integrand, domain, measure, properties...) =
    fallback_integrate(qs, integrand, domain, measure, properties...)

fallback_integrate(qs, integrand, domain, measure, properties...) =
    error("Don't know how to integrate on $domain with strategy $(qs).")

# For quadgk, we only know how to compute intervals
function apply_quad(qs::Q_quadgk, integrand, domain::AbstractInterval, measure::LebesgueMeasure, properties...)
    if isempty(domain)
        zero_result(integrand, prectype(domain))
    else
        quadgk(integrand, extrema(domain)...; atol = qs.atol, rtol = qs.rtol, maxevals = qs.maxevals)
    end
end

function apply_quad(qs::Q_hcubature, integrand, domain, measure, properties...)
    if islebesguemeasure(measure)
        apply_hcubature(qs, integrand, domain)
    else
        error("HCubature is invoked, but not with a Lebesgue measure: $(measure)")
    end
end

function apply_hcubature(qs::Q_hcubature, integrand, domain)
    f = x -> integrand(x) * DomainSets.indomain(x, domain)
    apply_hcubature(qs, f, boundingbox(domain))
end

function apply_hcubature(qs::Q_hcubature, integrand, domain::ProductDomain)
    if all(DomainSets.isinterval, components(domain))
        a = map(float, map(infimum, components(domain)))
        b = map(float, map(supremum, components(domain)))
        hcubature(integrand, a, b; rtol=qs.rtol, atol=qs.atol, maxevals = qs.maxevals)
    else
        apply_hcubature(qs, x -> integrand(x)*DomainSets.indomain(x, domain), boundingbox(domain))
    end
end

# Given a rule on [-1,1] and a different interval, we map the rule
function apply_quad(qs::ChebyshevIntervalRule, integrand, domain::AbstractInterval, measure::LebesgueMeasure, properties...)
    a = infimum(domain)
    b = supremum(domain)
    D = (b-a)/2
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(a + (x[i]+1)*D)*D for i in 1:length(x))
    z, unknown_error(z)
end

function apply_quad(qs::UnitIntervalRule, integrand, domain::AbstractInterval, measure::LebesgueMeasure, properties...)
    a = leftendpoint(domain)
    b = rightendpoint(domain)
    D = b-a
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(a + x[i]*D)*D for i in 1:length(x))
    z, unknown_error(z)
end

function apply_quad(qs::IntervalRule, integrand, interval::AbstractInterval, measure::LebesgueMeasure, properties...)
    a, b = extrema(interval)
    A, B = extrema(domain(qs))
    m = mapto(A..B, a..b)
    J = (b-a)/(B-A)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(m(x[i]))*J for i in 1:length(x))
    z, unknown_error(z)
end

# Make sure that a half line rule is applied to a half line integral
function apply_quad(qs::HalfLineRule, integrand, domain::HalfLine, measure::LebesgueMeasure, properties...)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end

# Make sure that a real line rule is applied to a real line integral
function apply_quad(qs::RealLineRule, integrand, domain::FullSpace, measure::LebesgueMeasure, properties...)
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end

# Any other quadrature rule, we just apply it
function apply_quad(qs::GenericDomainRule, integrand, domain, measure, properties...)
    @assert domain(qs) == domain
    x = points(qs)
    w = weights(qs)
    z = sum(w[i]*integrand(x[i]) for i in 1:length(x))
    z, unknown_error(z)
end


text:
@testset "Faro" begin
    @test Faro === Weave

    @test infaro === inweave === Faro{:in}()
    @test outfaro === outweave === Faro{:out}()

    x = [1, 2, 3, 4, 5, 6]
    shuffle!(x, Faro{:in}()) == x == [4, 1, 5, 2, 6, 3]

    y = [1, 2, 3, 4, 5, 6]
    shuffle([1, 2, 3, 4, 5, 6], Faro{:out}()) == [1, 4, 2, 5, 3, 6] != y

    arr = rand(Int, 1000)
    inshuffled = shuffle(arr, Faro{:in}())
    outshuffled = shuffle(arr, Faro{:out}())

    @test inshuffled[1:2:end] == arr[501:end]
    @test inshuffled[2:2:end] == arr[1:500]

    @test outshuffled[2:2:end] == arr[501:end]
    @test outshuffled[1:2:end] == arr[1:500]

    arr = nshuffle(collect(1:52), 26, Faro{:in}())
    @test arr == collect(52:-1:1)

    revarr = collect(52:-1:1)
    nshuffle!(revarr, 26, Faro{:in}())
    @test revarr == collect(1:52)

    @test nshuffle(collect(1:52), 8, Faro{:out}()) == collect(1:52)

    OA = OffsetArray(1:6, 6:11)
    @test parent(shuffle(OA, Faro{:out}())::OffsetArray) == [1, 4, 2, 5, 3, 6]
    @test parent(shuffle(OA, Faro{:in}())::OffsetArray) == [4, 1, 5, 2, 6, 3]
end

@testset "Cut" begin
    @test shuffle([1, 2, 3, 4, 5], Cut(2)) == [3, 4, 5, 1, 2]

    @test shuffle([1, 2, 3, 4, 5], Cut(0)) == [1, 2, 3, 4, 5]

    @test shuffle([1, 2, 3, 4, 5], Cut(5)) == [1, 2, 3, 4, 5]

    @test_throws DomainError shuffle([1, 2, 3, 4, 5], Cut(7))

    @test nshuffle(collect(1:100), 7, Cut(10)) == vcat(71:100, 1:70)
end

@testset "RandomShuffle" begin

    @test randshuffle === RandomShuffle()

    let mt1 = MersenneTwister(132), mt2 = MersenneTwister(132), x = rand(1000)
        @test shuffle(mt1, x, RandomShuffle()) == Random.shuffle(mt2, x)
    end

    # probably not the same after shuffling
    arr = rand(1000)
    arr_copy = copy(arr)
    @test shuffle!(arr, RandomShuffle()) == arr != arr_copy

    narr = rand(1000)
    narr_copy = copy(narr)
    @test nshuffle!(narr, 5, RandomShuffle()) == narr != narr_copy

    let mt1 = MersenneTwister(42), mt2 = MersenneTwister(42), x = rand(1000)
        @test shuffle(mt1, x, RandomShuffle()) ==
              shuffle!(mt2, x, RandomShuffle())
    end

    let mt1 = MersenneTwister(12), mt2 = MersenneTwister(12), x = rand(1000)
        @test shuffle(mt1, x, RandomShuffle()) ==
              shuffle!(mt2, x, RandomShuffle())
    end

    let mt1 = MersenneTwister(12), mt2 = MersenneTwister(12),
            mt3 = MersenneTwister(12), x = rand(1000)

        y = shuffle(mt1, x, RandomShuffle())
        shuffle!(mt1, y, RandomShuffle())
        shuffle!(mt1, y, RandomShuffle())

        @test nshuffle(mt2, x, 3, RandomShuffle()) == y

        @test nshuffle!(mt3, x, 3, RandomShuffle()) == y
    end
end

@testset "GilbertShannonReeds" begin

    @test gsrshuffle === GilbertShannonReeds()

    # probably not the same after shuffling
    arr = rand(1000)
    @test shuffle(arr, GilbertShannonReeds()) != arr

    narr = rand(1000)
    narr_copy = copy(narr)
    @test nshuffle!(narr, 10, GilbertShannonReeds()) == narr != narr_copy

    let mt1 = MersenneTwister(32), mt2 = MersenneTwister(32), x = rand(1000)
        @test shuffle(mt1, x, GilbertShannonReeds()) ==
              shuffle(mt2, x, GilbertShannonReeds())
    end

    let mt1 = MersenneTwister(153), mt2 = MersenneTwister(153), x = rand(1000)
        @test nshuffle(mt1, x, 7, GilbertShannonReeds()) ==
              nshuffle!(mt2, x, 7, GilbertShannonReeds())
    end
end


text:
export
    FLowProblems

module FlowProblems

# using ..CRCBS
using LightGraphs, MetaGraphs, JuMP

"""
    `extend_G_through_time(G,T::Int;render_depth=4)`

    Returns a time-extended (to horizon of T steps) version of G.
    Extension is performed via the "gadget" introduced by Yu and LaValle.

"""
function extend_G_through_time(G,T::Int;render_depth=4)
    N = nv(G)
    Gx = MetaDiGraph(nv(G)*(T+1))
    for v in vertices(G)
        for t in 0:T
            set_props!(Gx,v+N*t,Dict(:x=>get_prop(G,v,:x),:y=>render_depth*t))
        end
    end
    for t in 1:T
        edge_list = collect(edges(G))
        for e in edge_list
            x_mid = .5*(get_prop(Gx,e.src,:x)+get_prop(Gx,e.dst,:x))
            # gadget:
            #   o     o
            #   | \ / |
            #   |  o  |
            #   |  |  |
            #   |  o  |
            #   | / \ |
            #   o     o
            #   |     |
            #   o     o
            add_vertex!(Gx,Dict(:x=>x_mid,:y=>render_depth*(t-1) + render_depth*0.25))
            v = nv(Gx)
            add_edge!(Gx,e.src+N*(t-1),v)
            add_edge!(Gx,e.dst+N*(t-1),v)
            add_vertex!(Gx,Dict(:x=>x_mid,:y=>render_depth*(t-1) + render_depth*0.5))
            vp = nv(Gx)
            add_edge!(Gx,v,vp)
            add_vertex!(Gx,Dict(:x=>get_prop(Gx,e.src,:x),:y=>render_depth*(t-1) + render_depth*0.75))
            add_vertex!(Gx,Dict(:x=>get_prop(Gx,e.dst,:x),:y=>render_depth*(t-1) + render_depth*0.75))
            add_edge!(Gx,vp,nv(Gx)-1)
            add_edge!(Gx,vp,nv(Gx))
            # remain at node
            add_edge!(Gx,nv(Gx)-1,e.src+N*t)
            add_edge!(Gx,nv(Gx),e.dst+N*t)
            # skip
            add_edge!(Gx,e.src+N*(t-1),nv(Gx)-1)
            add_edge!(Gx,e.dst+N*(t-1),nv(Gx))
        end
    end
    Gx
end

"""
    `construct_ILP_flow_model(G,T::Int)`

    Constructs a time-extended version of the graph G out to a horizon
    of T steps, then returns both the ILP optimization model and the
    time extended graph Gx
"""
function construct_ILP_flow_model(G,T::Int,starts::Vector{Int},goals::Vector{Int},optimizer;TimeLimit=50, OutputFlag=0)
    Gx = extend_G_through_time(G,T;render_depth=4)
    model = Model(with_optimizer(optimizer,TimeLimit=TimeLimit,OutputFlag=OutputFlag))
    @variable(model, x[1:ne(Gx)], binary=true)
    M = incidence_matrix(Gx)
    srcs_sinks = zeros(Int,nv(Gx),1)
    for v in starts
        srcs_sinks[v] = 1
    end
    for v in goals
        srcs_sinks[v+nv(G)*T] = -1
    end
    # flow constraints between nodes (unit capacity) + sources and - sinks
    @constraint(model, (M*x .+ srcs_sinks) .== 0)
    @objective(model, Min, sum(x))
    return Gx, model, x
end

"""
    `extract_solution_from_flow_matrix(Gx,x)`

    Returns the sequences of vertices in the graph `Gx` that match the
    edge sequences encoded in `x`, the solution matrix of the flow problem
"""
function extract_solution_from_flow_matrix(Gx,x,starts::Vector{Int},goals::Vector{Int})
    edge_list = Vector{Edge}()
    for (i,e) in enumerate(edges(Gx))
        if x[i] == 1
            push!(edge_list,e)
        end
    end
    paths = [[v] for v in starts]
    while length(edge_list) > 0
        for (i,path) in enumerate(paths)
            if length(edge_list) <= 0
                break
            end
            for v in outneighbors(Gx,path[end])
                e = Edge(path[end],v)
                if e in edge_list
                    push!(path,v)
                    setdiff!(edge_list, [e])
                end
            end
        end
    end
    paths
end

"""
    `extract_paths_from_flow_solution(G,paths)`

    Takes the original graph G and the flow solution paths, and returns
    the corresponding edge_sequence through the original graph.
"""
function extract_paths_from_flow_solution(G,paths::Vector{Vector{Int}})
    true_paths = Vector{Vector{Int}}()
    for path in paths
        true_path = Vector{Int}()
        k = 1
        for i in path
            if i <= nv(G)*k
                push!(true_path, i-nv(G)*(k-1))
                k = k + 1
            end
        end
        push!(true_paths, true_path)
    end
    true_paths
end

end


text:
using Test
using Word2Vec

include("inputprocessing_tests.jl")
include("huffman_tests.jl")
include("contextdata_tests.jl")
# include("network_tests.jl")


text:
module Day2

export day2

struct PassWdRules
    range::Tuple{Int,Int}
    ch::Char
    passwd::String
end

function count(char, str)
    counter = 0
    for c in str
        if c == char
            counter += 1
        end
    end
    counter
end

function parse_rule(str)
    v = split(str,' ')
    r = map(t -> parse(Int,t), split(v[1],'-'))
    PassWdRules(tuple(r...), v[2][1], strip(v[3]))
end

function day2(::Val{:a}, filename)
    w = map(parse_rule, readlines(filename))
    counter = 0
    for p in w
        if p.range[1] <= count(p.ch, p.passwd) <= p.range[2]
            counter += 1
        end
    end
    counter
end

function day2(::Val{:b}, filename)
    w = map(parse_rule, readlines(filename))
    counter = 0
    for p in w
        if xor(p.passwd[p.range[1]] == p.ch, p.passwd[p.range[2]] == p.ch)
            counter += 1
        end
    end
    counter
end

end


text:
using CPLEX, Base.Test

@testset "LP 03" begin
    # Test get/set objective coefficients in LP

    env = CPLEX.Env()

    # original model
    #
    #   maximize  2x + 2y
    #
    # s.t. 0.2 <= x, y <= 1
    #

    model = CPLEX.cplex_model(env;
      name="lp_03",
      sense=:Max,
      f=[2.0, 2.0],
      lb=[0.2, 0.2],
      ub=[1.0, 1.0])

    lb_ = CPLEX.get_varLB(model)
    ub_ = CPLEX.get_varUB(model)
    c_ = CPLEX.get_obj(model)

    @test lb_ == [0.2, 0.2]
    @test ub_ == [1.0, 1.0]
    @test c_ == [2.0, 2.0]

    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test sol[1] ≈ 1
    @test sol[2] ≈ 1
    @test CPLEX.get_objval(model) ≈ 4

    # change objective (warm start)
    #
    # maximize x - y
    #
    # s.t. 0.2 <= x, y <= 1
    #

    CPLEX.set_obj!(model, [1, -1])

    c_ = CPLEX.get_obj(model)
    @test c_ == [1.0, -1.0]

    CPLEX.optimize!(model)

    sol = CPLEX.get_solution(model)
    @test sol[1] ≈ 1
    @test sol[2] ≈ 0.2
    @test CPLEX.get_objval(model) ≈ 0.8
end


text:
@testset "1D LDB" begin
    X, y = generateclassdata(ClassData(:tri, 5, 5, 5))
    wt = wavelet(WT.haar)

    # AsymmetricRelativeEntropy + TimeFrequency + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # SymmetricRelativeEntropy + TimeFrequency + FishersClassSeparability
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, dm=SymmetricRelativeEntropy(),
        dp=FishersClassSeparability(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # LpDistance + TimeFrequency + RobustFishersClassSeparability
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, dm=LpDistance(),
        dp=RobustFishersClassSeparability(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # HellingerDistance + ProbabilityDensity + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, dm=HellingerDistance(),
        en=ProbabilityDensity(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # EarthMoverDistance + Signatures(equal weight) + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, dm=EarthMoverDistance(),
        en=Signatures(:equal), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # EarthMoverDistance + Signatures(pdf weight) + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(wt=wt, max_dec_level=4, dm=EarthMoverDistance(),
        en=Signatures(:pdf), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (32, 15)

    # change number of features
    @test_nowarn change_nfeatures(f, Xc, 5)
    x = change_nfeatures(f, Xc, 5)
    @test size(x) == (5, 15)
    @test_logs (:warn, "Proposed n_features larger than currently saved n_features. Results will be less accurate since inverse_transform and transform is involved.") change_nfeatures(f, Xc, 10)
    @test_throws ArgumentError change_nfeatures(f, Xc, 10)
end

@testset "2D LDB" begin
    X = cat(rand(Normal(0,1), 8,8,5), rand(Normal(1,1), 8,8,5), rand(Normal(2,1), 8,8,5), dims=3)
    y = [repeat([1], 5); repeat([2],5); repeat([3],5)]

    # AsymmetricRelativeEntropy + TimeFrequency + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(max_dec_level=2, top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # SymmetricRelativeEntropy + TimeFrequency + FishersClassSeparability
    f = LocalDiscriminantBasis(max_dec_level=2, dm=SymmetricRelativeEntropy(),
        dp=FishersClassSeparability(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # LpDistance + TimeFrequency + RobustFishersClassSeparability
    f = LocalDiscriminantBasis(max_dec_level=2, dm=LpDistance(),
        dp=RobustFishersClassSeparability(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # HellingerDistance + ProbabilityDensity + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(max_dec_level=2, dm=HellingerDistance(),
        en=ProbabilityDensity(), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # EarthMoverDistance + Signatures(equal weight) + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(max_dec_level=2, dm=EarthMoverDistance(),
        en=Signatures(:equal), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # EarthMoverDistance + Signatures(pdf weight) + BasisDiscriminantMeasure
    f = LocalDiscriminantBasis(max_dec_level=2, dm=EarthMoverDistance(),
        en=Signatures(:pdf), top_k=5, n_features=5)
    @test typeof(f) == LocalDiscriminantBasis
    @test_nowarn fit_transform(f, X, y)
    @test_nowarn fit!(f, X, y)
    @test_nowarn transform(f, X)
    Xc = transform(f, X)
    @test size(Xc) == (5,15)
    @test_nowarn inverse_transform(f, Xc)
    X̂ = inverse_transform(f, Xc)
    @test size(X̂) == (8, 8, 15)

    # change number of features
    @test_nowarn change_nfeatures(f, Xc, 5)
    x = change_nfeatures(f, Xc, 5)
    @test size(x) == (5, 15)
    @test_logs (:warn, "Proposed n_features larger than currently saved n_features. Results will be less accurate since inverse_transform and transform is involved.") change_nfeatures(f, Xc, 10)
    @test_throws ArgumentError change_nfeatures(f, Xc, 10)
end


text:
"""
Start a notebook server using:

```julia
julia> Pluto.run()
```

Have a look at the FAQ:
https://github.com/fonsp/Pluto.jl/wiki
"""
module Pluto

import Pkg

project_relative_path(xs...) = normpath(joinpath(dirname(dirname(pathof(Pluto))), xs...))

const PLUTO_VERSION = VersionNumber(Pkg.TOML.parsefile(project_relative_path("Project.toml"))["version"])
const PLUTO_VERSION_STR = 'v' * string(PLUTO_VERSION)
const JULIA_VERSION_STR = 'v' * string(VERSION)

include("./Configuration.jl")

include("./evaluation/Tokens.jl")
include("./runner/PlutoRunner.jl")
# import .PlutoRunner
# @eval Main begin
#   PlutoRunner = $(PlutoRunner)
# end
include("./analysis/ExpressionExplorer.jl")
include("./analysis/ReactiveNode.jl")

include("./notebook/PathHelpers.jl")
include("./notebook/Cell.jl")
include("./notebook/Notebook.jl")
include("./webserver/Session.jl")
include("./webserver/PutUpdates.jl")

include("./analysis/Errors.jl")
include("./analysis/Parse.jl")
include("./analysis/Topology.jl")

include("./evaluation/WorkspaceManager.jl")
include("./evaluation/Update.jl")
include("./evaluation/Run.jl")

include("./webserver/MsgPack.jl")
include("./webserver/SessionActions.jl")
include("./webserver/Static.jl")
include("./webserver/Dynamic.jl")
include("./webserver/REPLTools.jl")
include("./webserver/WebServer.jl")

if get(ENV, "JULIA_PLUTO_SHOW_BANNER", "1") !== "0"
@info """\n
    Welcome to Pluto $(PLUTO_VERSION_STR) 🎈
    Start a notebook server using:

  julia> Pluto.run()

    Have a look at the FAQ:
    https://github.com/fonsp/Pluto.jl/wiki
\n"""
end

end


text:
###
### Particle Filtering and Particle MCMC Samplers.
###

####
#### Generic Sequential Monte Carlo sampler.
####

"""
$(TYPEDEF)

Sequential Monte Carlo sampler.

# Fields

$(TYPEDFIELDS)
"""
struct SMC{space, R} <: ParticleInference
    resampler::R
end

"""
    SMC(space...)
    SMC([resampler = AdvancedPS.ResampleWithESSThreshold(), space = ()])
    SMC([resampler = AdvancedPS.resample_systematic, ]threshold[, space = ()])

Create a sequential Monte Carlo sampler of type [`SMC`](@ref) for the variables in `space`.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
function SMC(resampler = AdvancedPS.ResampleWithESSThreshold(), space::Tuple = ())
    return SMC{space, typeof(resampler)}(resampler)
end

# Convenient constructors with ESS threshold
function SMC(resampler, threshold::Real, space::Tuple = ())
    return SMC(AdvancedPS.ResampleWithESSThreshold(resampler, threshold), space)
end
SMC(threshold::Real, space::Tuple = ()) = SMC(AdvancedPS.resample_systematic, threshold, space)

# If only the space is defined
SMC(space::Symbol...) = SMC(space)
SMC(space::Tuple) = SMC(AdvancedPS.ResampleWithESSThreshold(), space)

struct SMCTransition{T,F<:AbstractFloat}
    "The parameters for any given sample."
    θ::T
    "The joint log probability of the sample (NOTE: does not work, always set to zero)."
    lp::F
    "The weight of the particle the sample was retrieved from."
    weight::F
end

function SMCTransition(vi::AbstractVarInfo, weight)
    theta = tonamedtuple(vi)

    # This is pretty useless since we reset the log probability continuously in the
    # particle sweep.
    lp = getlogp(vi)

    return SMCTransition(theta, lp, weight)
end

metadata(t::SMCTransition) = (lp = t.lp, weight = t.weight)

DynamicPPL.getlogp(t::SMCTransition) = t.lp

struct SMCState{P,F<:AbstractFloat}
    particles::P
    particleindex::Int
    # The logevidence after aggregating all samples together.
    average_logevidence::F
end

function getlogevidence(samples, sampler::Sampler{<:SMC}, state::SMCState)
    return state.average_logevidence
end

function AbstractMCMC.sample(
    rng::AbstractRNG,
    model::AbstractModel,
    sampler::Sampler{<:SMC},
    N::Integer;
    chain_type=MCMCChains.Chains,
    resume_from=nothing,
    progress=PROGRESS[],
    kwargs...
)
    if resume_from === nothing
        return AbstractMCMC.mcmcsample(rng, model, sampler, N;
                                       chain_type=chain_type,
                                       progress=progress,
                                       nparticles=N,
                                       kwargs...)
    else
        return resume(resume_from, N;
                      chain_type=chain_type, progress=progress, nparticles=N, kwargs...)
    end
end

function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:SMC},
    vi::AbstractVarInfo;
    nparticles::Int,
    kwargs...
)
    # Reset the VarInfo.
    reset_num_produce!(vi)
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!!(vi)
    empty!!(vi)

    # Create a new set of particles.
    particles = AdvancedPS.ParticleContainer(
        [AdvancedPS.Trace(model, spl, vi) for _ in 1:nparticles],
    )

    # Perform particle sweep.
    logevidence = AdvancedPS.sweep!(rng, particles, spl.alg.resampler)

    # Extract the first particle and its weight.
    particle = particles.vals[1]
    weight = AdvancedPS.getweight(particles, 1)

    # Compute the first transition and the first state.
    transition = SMCTransition(particle.f.varinfo, weight)
    state = SMCState(particles, 2, logevidence)

    return transition, state
end

function AbstractMCMC.step(
    ::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:SMC},
    state::SMCState;
    kwargs...
)
    # Extract the index of the current particle.
    index = state.particleindex

    # Extract the current particle and its weight.
    particles = state.particles
    particle = particles.vals[index]
    weight = AdvancedPS.getweight(particles, index)

    # Compute the transition and the next state.
    transition = SMCTransition(particle.f.varinfo, weight)
    nextstate = SMCState(state.particles, index + 1, state.average_logevidence)

    return transition, nextstate
end

####
#### Particle Gibbs sampler.
####

"""
$(TYPEDEF)

Particle Gibbs sampler.

Note that this method is particle-based, and arrays of variables
must be stored in a [`TArray`](@ref) object.

# Fields

$(TYPEDFIELDS)
"""
struct PG{space,R} <: ParticleInference
    """Number of particles."""
    nparticles::Int
    """Resampling algorithm."""
    resampler::R
end

"""
    PG(n, space...)
    PG(n, [resampler = AdvancedPS.ResampleWithESSThreshold(), space = ()])
    PG(n, [resampler = AdvancedPS.resample_systematic, ]threshold[, space = ()])

Create a Particle Gibbs sampler of type [`PG`](@ref) with `n` particles for the variables
in `space`.

If the algorithm for the resampling step is not specified explicitly, systematic resampling
is performed if the estimated effective sample size per particle drops below 0.5.
"""
function PG(
    nparticles::Int,
    resampler = AdvancedPS.ResampleWithESSThreshold(),
    space::Tuple = (),
)
    return PG{space, typeof(resampler)}(nparticles, resampler)
end

# Convenient constructors with ESS threshold
function PG(nparticles::Int, resampler, threshold::Real, space::Tuple = ())
    return PG(nparticles, AdvancedPS.ResampleWithESSThreshold(resampler, threshold), space)
end
function PG(nparticles::Int, threshold::Real, space::Tuple = ())
    return PG(nparticles, AdvancedPS.resample_systematic, threshold, space)
end

# If only the number of particles and the space is defined
PG(nparticles::Int, space::Symbol...) = PG(nparticles, space)
function PG(nparticles::Int, space::Tuple)
    return PG(nparticles, AdvancedPS.ResampleWithESSThreshold(), space)
end

const CSMC = PG # type alias of PG as Conditional SMC

struct PGTransition{T,F<:AbstractFloat}
    "The parameters for any given sample."
    θ::T
    "The joint log probability of the sample (NOTE: does not work, always set to zero)."
    lp::F
    "The log evidence of the sample."
    logevidence::F
end

function PGTransition(vi::AbstractVarInfo, logevidence)
    theta = tonamedtuple(vi)

    # This is pretty useless since we reset the log probability continuously in the
    # particle sweep.
    lp = getlogp(vi)

    return PGTransition(theta, lp, logevidence)
end

metadata(t::PGTransition) = (lp = t.lp, logevidence = t.logevidence)

DynamicPPL.getlogp(t::PGTransition) = t.lp

function getlogevidence(samples, sampler::Sampler{<:PG}, vi::AbstractVarInfo)
    return mean(x.logevidence for x in samples)
end

function DynamicPPL.initialstep(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:PG},
    vi::AbstractVarInfo;
    kwargs...
)
    # Reset the VarInfo before new sweep
    reset_num_produce!(vi)
    set_retained_vns_del_by_spl!(vi, spl)
    resetlogp!!(vi)

    # Create a new set of particles
    num_particles = spl.alg.nparticles
    particles = AdvancedPS.ParticleContainer(
        [AdvancedPS.Trace(model, spl, vi) for _ in 1:num_particles],
    )

    # Perform a particle sweep.
    logevidence = AdvancedPS.sweep!(rng, particles, spl.alg.resampler)

    # Pick a particle to be retained.
    Ws = AdvancedPS.getweights(particles)
    indx = AdvancedPS.randcat(rng, Ws)
    reference = particles.vals[indx]

    # Compute the first transition.
    _vi = reference.f.varinfo
    transition = PGTransition(_vi, logevidence)

    return transition, _vi
end

function AbstractMCMC.step(
    rng::AbstractRNG,
    model::AbstractModel,
    spl::Sampler{<:PG},
    vi::AbstractVarInfo;
    kwargs...
)
    # Reset the VarInfo before new sweep.
    reset_num_produce!(vi)
    resetlogp!!(vi)

    # Create reference particle for which the samples will be retained.
    reference = AdvancedPS.forkr(AdvancedPS.Trace(model, spl, vi))

    # For all other particles, do not retain the variables but resample them.
    set_retained_vns_del_by_spl!(vi, spl)

    # Create a new set of particles.
    num_particles = spl.alg.nparticles
    x = map(1:num_particles) do i
        if i != num_particles
            return AdvancedPS.Trace(model, spl, vi)
        else
            return reference
        end
    end
    particles = AdvancedPS.ParticleContainer(x)

    # Perform a particle sweep.
    logevidence = AdvancedPS.sweep!(rng, particles, spl.alg.resampler, reference)

    # Pick a particle to be retained.
    Ws = AdvancedPS.getweights(particles)
    indx = AdvancedPS.randcat(rng, Ws)
    newreference = particles.vals[indx]

    # Compute the transition.
    _vi = newreference.f.varinfo
    transition = PGTransition(_vi, logevidence)

    return transition, _vi
end

DynamicPPL.use_threadsafe_eval(::SamplingContext{<:Sampler{<:Union{PG,SMC}}}, ::AbstractVarInfo) = false

function DynamicPPL.assume(
    rng,
    spl::Sampler{<:Union{PG,SMC}},
    dist::Distribution,
    vn::VarName,
    __vi__::AbstractVarInfo
)
    local vi
    try
        vi = AdvancedPS.current_trace().f.varinfo
    catch e
        # NOTE: this heuristic allows Libtask evaluating a model outside a `Trace`.
        if e == KeyError(:__trace) || current_task().storage isa Nothing
            vi = __vi__
        else
            rethrow(e)
        end
    end
    if inspace(vn, spl)
        if ~haskey(vi, vn)
            r = rand(rng, dist)
            push!!(vi, vn, r, dist, spl)
        elseif is_flagged(vi, vn, "del")
            unset_flag!(vi, vn, "del")
            r = rand(rng, dist)
            vi[vn] = vectorize(dist, r)
            DynamicPPL.setgid!(vi, spl.selector, vn)
            setorder!(vi, vn, get_num_produce(vi))
        else
            DynamicPPL.updategid!(vi, vn, spl)
            r = vi[vn]
        end
    else # vn belongs to other sampler <=> conditionning on vn
        if haskey(vi, vn)
            r = vi[vn]
        else
            r = rand(rng, dist)
            push!!(vi, vn, r, dist, DynamicPPL.Selector(:invalid))
        end
        lp = logpdf_with_trans(dist, r, istrans(vi, vn))
        acclogp!!(vi, lp)
    end
    return r, 0, vi
end

function DynamicPPL.observe(spl::Sampler{<:Union{PG,SMC}}, dist::Distribution, value, vi)
    Libtask.produce(logpdf(dist, value))
    return 0, vi
end

# Convenient constructor
function AdvancedPS.Trace(
    model::Model,
    sampler::Sampler{<:Union{SMC,PG}},
    varinfo::AbstractVarInfo,
)
    newvarinfo = deepcopy(varinfo)
    DynamicPPL.reset_num_produce!(newvarinfo)
    f = Turing.Essential.TracedModel(model, sampler, newvarinfo)
    return AdvancedPS.Trace(f)
end


text:
using Plots, JuMP, Ipopt

mutable struct SVM
    X # training features
    y # training labels
    α
    β
    C
    η
end

SVM(X,y) = SVM(X,y,0,0,1,0)

function dualfit!(S::SVM; C = 1.0)
    S.C = C
    nrows, ncols = size(S.X)
    model = Model(with_optimizer(Ipopt.Optimizer, print_level=0))
    @variable(model, η[1:nrows])
    @variable(model, α)
    @objective(model, Max, -1/2 * (η .* S.y)' * S.X * S.X' * (η .* S.y) + sum(η))
    @constraint(model, 0 .≤ η[1:nrows] .≤ S.C)
    @constraint(model, η' * S.y == 0)
    optimize!(model)
    S.η = value.(η)
    support_vectors = 1e-2 .<  S.η .< S.C - 1e-2
    S.β = S.X' * (S.η .* S.y)
    S.α = mean((S.y - X * S.β)[support_vectors])
    nothing
end

import Plots.plot

function plot(S::SVM)
    zone(t) = t > 1 ? 2 : (t > 0 ? 1 : (t > -1 ? -1 : -2))
    xmin, xmax = extrema(S.X[:,1])
    ymin, ymax = extrema(S.X[:,2])
    xmin -= (xmax - xmin)/20; xmax += (xmax - xmin)/20
    ymin -= (ymax - ymin)/20; ymax += (ymax - ymin)/20
    xgrid = range(xmin, xmax, length=512)
    ygrid = range(ymin, ymax, length=512)
    scatter(X[:,1], X[:,2], group = y, markersize = 2)
    contour!(xgrid, ygrid, (x,y) -> S.β'*[x,y] + S.α,
              levels = [-1, 0, 1])
    heatmap!(xgrid, ygrid,
             (x,y) -> zone(S.β'*[x,y] + S.α), fillopacity = 0.5,
             fillcolor = cgrad([:blue, :lightblue, :pink, :red]))
end

function plot(S::SVM, f::Function)
    P = plot(S)
    annotate!(P, [(S.X[i,1],S.X[i,2]-0.2,
                  text("$(abs(round(f(S,i),digits=2)))",8)) for i in 1:size(S.X,1)])
    plot!(P, size = (600,400), colorbar = false, legend = false)
    P
end

max_stars_repo_name: data-gymnasia/julia-binder

max_issues_repo_name: data-gymnasia/julia-binder

max_forks_repo_name: data-gymnasia/julia-binder


text:
# test functions for developing optimizers
# see https://en.wikipedia.org/wiki/Test_functions_for_optimization

# optimal = (0, 0, 0, ...)
function sphere(x)
    return sum(x^2)
end

# optimal = (1,1,1,...)
function rosenbrock(x)
    x1 = x[1:end-1]
    x2 = x[2:end]
    return sum(100 * (x2-x1^2)^2 + (1-x1)^2)
end

# optimal = (3, 0.5)
function beale(p)
    x, y = p[1], p[2]
    return (1.5 - x + x*y)^2 + (2.25 - x + x*y^2)^2 +
            (2.625 - x + x*y^3)
end

# optimal = (1, 3)
function booth(p)
    x, y = p[1], p[2]
    return (x + 2y - 7)^2 + (2x+y-5)^2
end

#optimal = (-0.2903, ...)
function stybliski_tang(x)
    return sum(x^4 - 16x^2 + 5x)/2
end

# optimal = (1,1)
function levi(p)
    x, y = p[1], p[2]
    sin(3π*x)^2 + (x-1)^2 * (1+sin(3π*y)^2) + (y-1)^2 * (1+sin(2π*y)^2)
end


text:
# Autogenerated wrapper script for SuiteSparse_GPU_jll for x86_64-linux-gnu
export libGPUQREngine, libSuiteSparse_GPURuntime, libamd, libbtf, libcamd, libccolamd, libcholmod, libcolamd, libklu, libldl, librbio, libsliplu, libspqr, libsuitesparseconfig, libumfpack

using libblastrampoline_jll
using METIS_jll
using MPFR_jll
using GMP_jll
JLLWrappers.@generate_wrapper_header("SuiteSparse_GPU")
JLLWrappers.@declare_library_product(libGPUQREngine, "libGPUQREngine.so.1")
JLLWrappers.@declare_library_product(libSuiteSparse_GPURuntime, "libSuiteSparse_GPURuntime.so.1")
JLLWrappers.@declare_library_product(libamd, "libamd.so.2")
JLLWrappers.@declare_library_product(libbtf, "libbtf.so.1")
JLLWrappers.@declare_library_product(libcamd, "libcamd.so.2")
JLLWrappers.@declare_library_product(libccolamd, "libccolamd.so.2")
JLLWrappers.@declare_library_product(libcholmod, "libcholmod.so.3")
JLLWrappers.@declare_library_product(libcolamd, "libcolamd.so.2")
JLLWrappers.@declare_library_product(libklu, "libklu.so.1")
JLLWrappers.@declare_library_product(libldl, "libldl.so.2")
JLLWrappers.@declare_library_product(librbio, "librbio.so.2")
JLLWrappers.@declare_library_product(libsliplu, "libsliplu.so.1")
JLLWrappers.@declare_library_product(libspqr, "libspqr.so.2")
JLLWrappers.@declare_library_product(libsuitesparseconfig, "libsuitesparseconfig.so.5")
JLLWrappers.@declare_library_product(libumfpack, "libumfpack.so.5")
function __init__()
    JLLWrappers.@generate_init_header(libblastrampoline_jll, METIS_jll, MPFR_jll, GMP_jll)
    JLLWrappers.@init_library_product(
        libGPUQREngine,
        "lib/libGPUQREngine.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libSuiteSparse_GPURuntime,
        "lib/libSuiteSparse_GPURuntime.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libamd,
        "lib/libamd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libbtf,
        "lib/libbtf.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcamd,
        "lib/libcamd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libccolamd,
        "lib/libccolamd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcholmod,
        "lib/libcholmod.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcolamd,
        "lib/libcolamd.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libklu,
        "lib/libklu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libldl,
        "lib/libldl.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        librbio,
        "lib/librbio.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsliplu,
        "lib/libsliplu.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libspqr,
        "lib/libspqr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsuitesparseconfig,
        "lib/libsuitesparseconfig.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libumfpack,
        "lib/libumfpack.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module JungleHelperFallingKillbox

using ..Ahorn, Maple

@mapdef Entity "JungleHelper/FallingKillbox" FallingKillbox(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, fallSpeed::Number=100.0)

const placements = Ahorn.PlacementDict(
    "Falling Killbox (Jungle Helper)" => Ahorn.EntityPlacement(
        FallingKillbox,
        "rectangle"
    ),
)

Ahorn.minimumSize(entity::FallingKillbox) = 8, 0
Ahorn.resizable(entity::FallingKillbox) = true, false

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

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

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::FallingKillbox, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = 32

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.8, 0.4, 0.4, 0.8), (0.0, 0.0, 0.0, 0.0))
end

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
fs(f, s) = map(f, s)
f1(x) = 2x
f2(x) = x^2
fsf1(s) = fs(f1, s)
fsf2(s) = fs(f2, s)

s1 = [0, 1, 2 ,3]
s2 = [2, 4, 6, 8]
println("fsf1 of s1 is $(fsf1(s1))")
println("fsf2 of s1 is $(fsf2(s1))")
println("fsf1 of s2 is $(fsf1(s2))")
println("fsf2 of s2 is $(fsf2(s2))")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# von Mises-Fisher distribution is useful for directional statistics
#
# The implementation here follows:
#
#   - Wikipedia:
#     http://en.wikipedia.org/wiki/Von_Mises–Fisher_distribution
#
#   - R's movMF package's document:
#     http://cran.r-project.org/web/packages/movMF/vignettes/movMF.pdf
#
#   - Wenzel Jakob's notes:
#     http://www.mitsuba-renderer.org/~wenzel/files/vmf.pdf
#

immutable VonMisesFisher <: ContinuousMultivariateDistribution
    μ::Vector{Float64}
    κ::Float64
    logCκ::Float64

    function VonMisesFisher(μ::Vector{Float64}, κ::Float64; checknorm::Bool=true)
        if checknorm
            isunitvec(μ) || error("μ must be a unit vector")
        end
        κ > 0 || error("κ must be positive.")
        new(μ, κ, vmflck(length(μ), κ))
    end
end

@compat VonMisesFisher{T<:Real}(μ::Vector{T}, κ::Real) = VonMisesFisher(Float64(μ), Float64(κ))

VonMisesFisher(θ::Vector{Float64}) = (κ = vecnorm(θ); VonMisesFisher(scale(θ, 1.0 / κ), κ))
VonMisesFisher{T<:Real}(θ::Vector{T}) = VonMisesFisher(Float64(θ))

show(io::IO, d::VonMisesFisher) = show(io, d, (:μ, :κ))

### Basic properties

length(d::VonMisesFisher) = length(d.μ)

meandir(d::VonMisesFisher) = d.μ
concentration(d::VonMisesFisher) = d.κ

insupport{T<:Real}(d::VonMisesFisher, x::DenseVector{T}) = isunitvec(x)

### Evaluation

function _vmflck(p, κ)
    hp = 0.5 * p
    q = hp - 1.0
    q * log(κ) - hp * log(2π) - log(besseli(q, κ))
end
_vmflck3(κ) = log(κ) - log2π - κ - log1mexp(-2.0 * κ)
vmflck(p, κ) = (p == 3 ? _vmflck3(κ) : _vmflck(p, κ))::Float64

_logpdf{T<:Real}(d::VonMisesFisher, x::DenseVector{T}) = d.logCκ + d.κ * dot(d.μ, x)

### Sampling

sampler(d::VonMisesFisher) = VonMisesFisherSampler(d.μ, d.κ)

_rand!(d::VonMisesFisher, x::DenseVector) = _rand!(sampler(d), x)
_rand!(d::VonMisesFisher, x::DenseMatrix) = _rand!(sampler(d), x)

### Estimation

function fit_mle(::Type{VonMisesFisher}, X::Matrix{Float64})
    r = vec(sum(X, 2))
    n = size(X, 2)
    r_nrm = vecnorm(r)
    μ = scale!(r, 1.0 / r_nrm)
    ρ = r_nrm / n
    κ = _vmf_estkappa(length(μ), ρ)
    VonMisesFisher(μ, κ)
end

@compat fit_mle{T<:Real}(::Type{VonMisesFisher}, X::Matrix{T}) = fit_mle(VonMisesFisher, Float64(X))

function _vmf_estkappa(p::Int, ρ::Float64)
    # Using the fixed-point iteration algorithm in the following paper:
    #
    #   Akihiro Tanabe, Kenji Fukumizu, and Shigeyuki Oba, Takashi Takenouchi, and Shin Ishii
    #   Parameter estimation for von Mises-Fisher distributions.
    #   Computational Statistics, 2007, Vol. 22:145-157.
    #

    const maxiter = 200
    half_p = 0.5 * p

    ρ2 = abs2(ρ)
    κ = ρ * (p - ρ2) / (1 - ρ2)
    i = 0
    while i < maxiter
        i += 1
        κ_prev = κ
        a = (ρ / _vmfA(half_p, κ))
        # println("i = $i, a = $a, abs(a - 1) = $(abs(a - 1))")
        κ *= a
        if abs(a - 1.0) < 1.0e-12
            break
        end
    end
    return κ
end

_vmfA(half_p::Float64, κ::Float64) = besseli(half_p, κ) / besseli(half_p - 1.0, κ)


text:
include("pass_helpers.jl")
include("pass1.jl")
include("pass2.jl")
include("pass3.jl")

function transpile_py_2_jl(project)

    # pass1(project)
    # pass2(project)
    pass3(project)

end


text:
abstract type Queryable end

QueryOperators.query(x::Queryable) = x

IteratorInterfaceExtensions.isiterable(x::Queryable) = true
TableTraits.isiterabletable(x::Queryable) = true

function IteratorInterfaceExtensions.getiterator(x::Queryable)
    return x.getiterator(x)
end


text:
# MLJ interface for Decision Trees/Random Forests models

import MLJModelInterface       # It seems that having done this in the top module is not enought
const MMI = MLJModelInterface  # We need to repoeat it here

export DecisionTreeRegressor, RandomForestRegressor, DecisionTreeClassifier, RandomForestClassifier

# ------------------------------------------------------------------------------
# Model Structure declarations..

mutable struct DecisionTreeRegressor <: MMI.Deterministic
    maxDepth::Int64
    minGain::Float64
    minRecords::Int64
    maxFeatures::Int64
    splittingCriterion::Function
    rng::AbstractRNG
end
DecisionTreeRegressor(;
   maxDepth=0, #typemax(Int)
   minGain=0.0,
   minRecords=2,
   maxFeatures=0,
   splittingCriterion=variance,
   rng = Random.GLOBAL_RNG,
   ) = DecisionTreeRegressor(maxDepth,minGain,minRecords,maxFeatures,splittingCriterion,rng)

mutable struct DecisionTreeClassifier <: MMI.Probabilistic
   maxDepth::Int64
   minGain::Float64
   minRecords::Int64
   maxFeatures::Int64
   splittingCriterion::Function
   rng::AbstractRNG
end
DecisionTreeClassifier(;
  maxDepth=0,
  minGain=0.0,
  minRecords=2,
  maxFeatures=0,
  splittingCriterion=gini,
  rng = Random.GLOBAL_RNG,
  ) = DecisionTreeClassifier(maxDepth,minGain,minRecords,maxFeatures,splittingCriterion,rng)

mutable struct RandomForestRegressor <: MMI.Deterministic
   nTrees::Int64
   maxDepth::Int64
   minGain::Float64
   minRecords::Int64
   maxFeatures::Int64
   splittingCriterion::Function
   β::Float64
   rng::AbstractRNG
end
RandomForestRegressor(;
  nTrees=30,
  maxDepth=0,
  minGain=0.0,
  minRecords=2,
  maxFeatures=0,
  splittingCriterion=variance,
  β=0.0,
  rng = Random.GLOBAL_RNG,
  ) = RandomForestRegressor(nTrees,maxDepth,minGain,minRecords,maxFeatures,splittingCriterion,β,rng)

mutable struct RandomForestClassifier <: MMI.Probabilistic
    nTrees::Int64
    maxDepth::Int64
    minGain::Float64
    minRecords::Int64
    maxFeatures::Int64
    splittingCriterion::Function
    β::Float64
    rng::AbstractRNG
end
RandomForestClassifier(;
    nTrees=30,
    maxDepth=0,
    minGain=0.0,
    minRecords=2,
    maxFeatures=0,
    splittingCriterion=gini,
    β=0.0,
    rng = Random.GLOBAL_RNG,
) = RandomForestClassifier(nTrees,maxDepth,minGain,minRecords,maxFeatures,splittingCriterion,β,rng)

#=
# skipped for now..
# ------------------------------------------------------------------------------
# Hyperparameters ranges definition (for automatic tuning)

MMI.hyperparameter_ranges(::Type{<:DecisionTreeRegressor}) = (
#    (range(Float64, :alpha, lower=0, upper=1, scale=:log),
#     range(Int, :beta, lower=1, upper=Inf, origin=100, unit=50, scale=:log),
#         nothing)
    range(Int64,:maxDepth,lower=0,upper=Inf,scale=:log),
    range(Float64,:minGain,lower=0,upper=Inf,scale=:log),
    range(Int64,:minRecords,lower=0,upper=Inf,scale=:log),
    range(Int64,:maxFeatures,lower=0,upper=Inf,scale=:log),
    nothing
)
=#

# ------------------------------------------------------------------------------
# Fit functions...

function MMI.fit(model::Union{DecisionTreeRegressor,RandomForestRegressor}, verbosity, X, y)
   x = MMI.matrix(X)                     # convert table to matrix
   maxDepth         = model.maxDepth == 0 ? size(x,1) : model.maxDepth
   if (typeof(model) == DecisionTreeRegressor)
       maxFeatures = model.maxFeatures == 0 ? size(x,2) : model.maxFeatures
       fitresult   = buildTree(x, y, maxDepth=maxDepth, minGain=model.minGain, minRecords=model.minRecords, maxFeatures=maxFeatures, splittingCriterion=model.splittingCriterion,rng=model.rng)
   else
       maxFeatures = model.maxFeatures == 0 ? Int(round(sqrt(size(x,2)))) : model.maxFeatures
       fitresult   = buildForest(x, y, model.nTrees, maxDepth=maxDepth, minGain=model.minGain, minRecords=model.minRecords, maxFeatures=maxFeatures, splittingCriterion=model.splittingCriterion, β=model.β,rng=model.rng)
   end
   cache=nothing
   report=nothing
   return fitresult, cache, report
end

function MMI.fit(model::Union{DecisionTreeClassifier,RandomForestClassifier}, verbosity, X, y)
   x                = MMI.matrix(X)                        # convert table to matrix
   a_target_element = y[1]                                 # a CategoricalValue or CategoricalString
   #y_plain          = MMI.int(y) .- 1                     # integer relabeling should start at 0
   yarray           = convert(Vector{eltype(levels(y))},y) # convert to a simple Array{T}
   maxDepth         = model.maxDepth == 0 ? size(x,1) : model.maxDepth
   if (typeof(model) == DecisionTreeClassifier)
       maxFeatures   = model.maxFeatures == 0 ? size(x,2) : model.maxFeatures
       fittedmodel   = buildTree(x, yarray, maxDepth=maxDepth, minGain=model.minGain, minRecords=model.minRecords, maxFeatures=maxFeatures, splittingCriterion=model.splittingCriterion, forceClassification=true,rng=model.rng)
   else
       maxFeatures   = model.maxFeatures == 0 ? Int(round(sqrt(size(x,2)))) : model.maxFeatures
       fittedmodel   = buildForest(x, yarray, model.nTrees, maxDepth=maxDepth, minGain=model.minGain, minRecords=model.minRecords, maxFeatures=maxFeatures, splittingCriterion=model.splittingCriterion, forceClassification=true, β=model.β,rng=model.rng)
   end
   cache            = nothing
   report           = nothing
   fitresult        = (fittedmodel,a_target_element)
   return (fitresult, cache, report)
end

# ------------------------------------------------------------------------------
# Predict functions....

MMI.predict(model::Union{DecisionTreeRegressor,RandomForestRegressor}, fitresult, Xnew) = Trees.predict(fitresult, MMI.matrix(Xnew))
function MMI.predict(model::Union{DecisionTreeClassifier,RandomForestClassifier}, fitresult, Xnew)
    fittedModel      = fitresult[1]
    a_target_element = fitresult[2]
    decode           = MMI.decoder(a_target_element)
    classes          = MMI.classes(a_target_element)
    #println(typeof(classes))
    nLevels          = length(classes)
    nRecords         = MMI.nrows(Xnew)
    treePredictions  = Trees.predict(fittedModel, MMI.matrix(Xnew),rng=model.rng)
    predMatrix       = zeros(Float64,(nRecords,nLevels))
    # Transform the predictions from a vector of dictionaries to a matrix
    # where the rows are the PMF of each record
    for n in 1:nRecords
        for (c,cl) in enumerate(classes)
            predMatrix[n,c] = get(treePredictions[n],cl,0.0)
        end
    end
    #predictions = [MMI.UnivariateFinite(classes, predMatrix[i,:])
    #               for i in 1:nRecords]
    predictions = MMI.UnivariateFinite(classes, predMatrix)
    return predictions
end

# ------------------------------------------------------------------------------
# Model metadata for registration in MLJ...

MMI.metadata_model(DecisionTreeRegressor,
    input_scitype    = MMI.Table(MMI.Missing, MMI.Known),           # also ok: MMI.Table(Union{MMI.Missing, MMI.Known}),
    target_scitype   = AbstractVector{<: MMI.Continuous},           # for a supervised model, what target?
    supports_weights = false,                                       # does the model support sample weights?
    descr            = "A simple Decision Tree for regression with support for Missing data, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Trees.DecisionTreeRegressor"
    )
MMI.metadata_model(RandomForestRegressor,
    input_scitype    = MMI.Table(MMI.Missing, MMI.Known),
    target_scitype   = AbstractVector{<: MMI.Continuous},
    supports_weights = false,
    descr            = "A simple Random Forest ensemble for regression with support for Missing data, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Trees.RandomForestRegressor"
    )
MMI.metadata_model(DecisionTreeClassifier,
    input_scitype    = MMI.Table(MMI.Missing, MMI.Known),
    target_scitype   = AbstractVector{<: Union{MMI.Missing,MMI.Finite}},
    supports_weights = false,
    descr            = "A simple Decision Tree for classification with support for Missing data, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Trees.DecisionTreeClassifier"
    )
MMI.metadata_model(RandomForestClassifier,
    input_scitype    = MMI.Table(MMI.Missing, MMI.Known),
    target_scitype   = AbstractVector{<: Union{MMI.Missing,MMI.Finite}},
    supports_weights = false,
    descr            = "A simple Random Forest ensemble for classification with support for Missing data, from the Beta Machine Learning Toolkit (BetaML).",
	load_path        = "BetaML.Trees.RandomForestClassifier"
    )


text:
# Autogenerated wrapper script for Qt6Base_jll for i686-w64-mingw32-cxx11
export libqt6concurrent, libqt6core, libqt6dbus, libqt6gui, libqt6network, libqt6opengl, libqt6printsupport, libqt6sql, libqt6test, libqt6widgets, libqt6xml

using Xorg_libXext_jll
using Xorg_libxcb_jll
using Xorg_xcb_util_wm_jll
using Xorg_xcb_util_image_jll
using Xorg_xcb_util_keysyms_jll
using Xorg_xcb_util_renderutil_jll
using Xorg_libXrender_jll
using xkbcommon_jll
using Libglvnd_jll
using Fontconfig_jll
using Glib_jll
using Zlib_jll
using CompilerSupportLibraries_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("Qt6Base")
JLLWrappers.@declare_library_product(libqt6concurrent, "Qt6Concurrent.dll")
JLLWrappers.@declare_library_product(libqt6core, "Qt6Core.dll")
JLLWrappers.@declare_library_product(libqt6dbus, "Qt6DBus.dll")
JLLWrappers.@declare_library_product(libqt6gui, "Qt6Gui.dll")
JLLWrappers.@declare_library_product(libqt6network, "Qt6Network.dll")
JLLWrappers.@declare_library_product(libqt6opengl, "Qt6OpenGL.dll")
JLLWrappers.@declare_library_product(libqt6printsupport, "Qt6PrintSupport.dll")
JLLWrappers.@declare_library_product(libqt6sql, "Qt6Sql.dll")
JLLWrappers.@declare_library_product(libqt6test, "Qt6Test.dll")
JLLWrappers.@declare_library_product(libqt6widgets, "Qt6Widgets.dll")
JLLWrappers.@declare_library_product(libqt6xml, "Qt6Xml.dll")
function __init__()
    JLLWrappers.@generate_init_header(Xorg_libXext_jll, Xorg_libxcb_jll, Xorg_xcb_util_wm_jll, Xorg_xcb_util_image_jll, Xorg_xcb_util_keysyms_jll, Xorg_xcb_util_renderutil_jll, Xorg_libXrender_jll, xkbcommon_jll, Libglvnd_jll, Fontconfig_jll, Glib_jll, Zlib_jll, CompilerSupportLibraries_jll, OpenSSL_jll)
    JLLWrappers.@init_library_product(
        libqt6concurrent,
        "bin\\Qt6Concurrent.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6core,
        "bin\\Qt6Core.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6dbus,
        "bin\\Qt6DBus.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6gui,
        "bin\\Qt6Gui.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6network,
        "bin\\Qt6Network.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6opengl,
        "bin\\Qt6OpenGL.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6printsupport,
        "bin\\Qt6PrintSupport.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6sql,
        "bin\\Qt6Sql.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6test,
        "bin\\Qt6Test.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6widgets,
        "bin\\Qt6Widgets.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqt6xml,
        "bin\\Qt6Xml.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Test the TranscriptionData datatype
@testset "TranscriptionData" begin
    @test TranscriptionData isa DataType
    @test TranscriptionData().finite_to_constr isa Dict
end

# Test basic definition and queries
@testset "Basic Definition and Queries" begin
    # initialize needed data
    mockoptimizer = () -> MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()),
                                             eval_objective_value=false)
    # test TranscriptionModel (no optimizer)
    @testset "TranscriptionModel (Default)" begin
        @test haskey(TranscriptionModel().ext, :TransData)
    end
    # test TranscriptionModel (with optimizer)
    @testset "TranscriptionModel (Optimizer)" begin
        @test haskey(TranscriptionModel(mockoptimizer).ext, :TransData)
    end
    # test is_transcription_model
    @testset "is_transcription_model" begin
        @test is_transcription_model(TranscriptionModel())
        @test !is_transcription_model(Model())
    end
    # test transcription_data
    @testset "transcription_data" begin
        @test transcription_data(TranscriptionModel()).point_to_var isa Dict
        @test_throws ErrorException transcription_data(Model())
    end
end

# Test variable mapping queries
@testset "Variable Mapping Queries" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, x(par))
    @point_variable(m, x(0), x0)
    @hold_variable(m, y)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    # test transcription_variable (HoldVariableRef)
    @testset "transcription_variable (Hold)" begin
        # test error
        @test_throws ErrorException transcription_variable(tm, y)
        # test normal
        tm.ext[:TransData].hold_to_var[y] = a
        @test transcription_variable(tm, y) == a
    end
    # test transcription_variable (InfiniteVariableRef)
    @testset "transcription_variable (Infinite)" begin
        # test error
        @test_throws ErrorException transcription_variable(tm, x)
        # test normal
        tm.ext[:TransData].infinite_to_vars[x] = [b, c]
        @test transcription_variable(tm, x) == [b, c]
    end
    # test transcription_variable (PointVariableRef)
    @testset "transcription_variable (Point)" begin
        # test error
        @test_throws ErrorException transcription_variable(tm, x0)
        # test normal
        tm.ext[:TransData].point_to_var[x0] = b
        @test transcription_variable(tm, x0) == b
    end
    # test transcription_variable (Single argument)
    @testset "transcription_variable (Single)" begin
        @test transcription_variable(y) == a
        @test transcription_variable(x) == [b, c]
        @test transcription_variable(x0) == b
    end
    # test optimizer_model_variable extension
    @testset "optimizer_model_variable" begin
        @test optimizer_model_variable(y, Val(:TransData)) == a
        @test optimizer_model_variable(x, Val(:TransData)) == [b, c]
        @test optimizer_model_variable(x0, Val(:TransData)) == b
    end
    # test variable_supports for infinite variable with 2 inputs
    @testset "variable_supports (Model, Infinite Variable)" begin
        # test error
        @test_throws ErrorException variable_supports(tm, x)
        # test normal
        tm.ext[:TransData].infvar_to_supports[x] = [(0.,), (1.,)]
        @test variable_supports(tm, x) == [(0.,), (1.,)]
        # undo changes
        delete!(tm.ext[:TransData].infvar_to_supports, x)
    end
    # test supports for infinite variable
    @testset "supports (Infinite Variable)" begin
        # test error
        @test_throws ErrorException supports(x)
        # test normal
        tm.ext[:TransData].infvar_to_supports[x] = [(0.,), (1.,)]
        @test supports(x) == [(0.,), (1.,)]
    end
end

# Test constraint mapping queries
@testset "Constraint Mapping Queries" begin
    # initialize models
    m = InfiniteModel()
    @infinite_parameter(m, 0 <= par <= 1, supports = [0, 1])
    @infinite_variable(m, x(par))
    @point_variable(m, x(0), x0)
    @hold_variable(m, y)
    data = DiscreteMeasureData(par, [0.5, 0.5], [0, 1])
    @constraint(m, c1, x + y - 2 <= 0)
    @constraint(m, c2, measure(x, data) == 0)
    @constraint(m, c3, x0 + y == 5)
    tm = optimizer_model(m)
    @variable(tm, a)
    @variable(tm, b)
    @variable(tm, c)
    @constraint(tm, tc1, b + a <= 2)
    @constraint(tm, tc2, c + a <= 2)
    @constraint(tm, tc3, 0.5b + 0.5c == 0)
    @constraint(tm, tc4, b + a == 5)
    # test transcription_constraint (InfiniteConstraintRef)
    @testset "transcription_constraint (Infinite)" begin
        # test error
        @test_throws ErrorException transcription_constraint(tm, c1)
        # test normal
        tm.ext[:TransData].infinite_to_constrs[c1] = [tc1, tc2]
        @test transcription_constraint(tm, c1) == [tc1, tc2]
    end
    # test transcription_constraint (MeasureConstraintRef)
    @testset "transcription_constraint (Measure)" begin
        # test error
        @test_throws ErrorException transcription_constraint(tm, c2)
        # test normal
        tm.ext[:TransData].measure_to_constrs[c2] = [tc3]
        @test transcription_constraint(tm, c2) == [tc3]
    end
    # test transcription_constraint (FiniteConstraintRef)
    @testset "transcription_constraint (Finite)" begin
        # test error
        @test_throws ErrorException transcription_constraint(tm, c3)
        # test normal
        tm.ext[:TransData].finite_to_constr[c3] = tc4
        @test transcription_constraint(tm, c3) == tc4
    end
    # test transcription_constraint (Single argument)
    @testset "transcription_constraint (Single)" begin
        @test transcription_constraint(c1) == [tc1, tc2]
        @test transcription_constraint(c2) == [tc3]
        @test transcription_constraint(c3) == tc4
    end
    # test optimizer_model_constraint extension
    @testset "optimizer_model_constraint" begin
        @test optimizer_model_constraint(c1, Val(:TransData)) == [tc1, tc2]
        @test optimizer_model_constraint(c2, Val(:TransData)) == [tc3]
        @test optimizer_model_constraint(c3, Val(:TransData)) == tc4
    end
    # test constraint_supports for infinite constraint with 2 inputs
    @testset "constraint_supports (Model, Infinite)" begin
        # test error
        @test_throws ErrorException constraint_supports(tm, c1)
        # test normal
        tm.ext[:TransData].infconstr_to_supports[c1] = [(0.,), (1.,)]
        @test constraint_supports(tm, c1) == [(0.,), (1.,)]
        # undo changes
        delete!(tm.ext[:TransData].infconstr_to_supports, c1)
    end
    # test supports for infinite constraint with 1 input
    @testset "supports (Infinite)" begin
        # test error
        @test_throws ErrorException supports(c1)
        # test normal
        tm.ext[:TransData].infconstr_to_supports[c1] = [(0.,), (1.,)]
        @test supports(c1) == [(0.,), (1.,)]
    end
    # test supports for measure constraint with 2 inputs
    @testset "constraint_supports (Model, Measure)" begin
        # test error
        @test_throws ErrorException constraint_supports(tm, c2)
        # test normal
        tm.ext[:TransData].measconstr_to_supports[c2] = [(0.,), (1.,)]
        @test constraint_supports(tm, c2) == [(0.,), (1.,)]
        # undo changes
        delete!(tm.ext[:TransData].measconstr_to_supports, c2)
    end
    # test supports for measure constraint with 1 input
    @testset "supports (Measure)" begin
        # test error
        @test_throws ErrorException supports(c2)
        # test normal
        tm.ext[:TransData].measconstr_to_supports[c2] = [(0.,), (1.,)]
        @test supports(c2) == [(0.,), (1.,)]
    end
    # test constraint_parameter_refs for infinite constraint with 2 inputs
    @testset "constraint_parameter_refs (Model, Infinite)" begin
        # test error
        @test_throws ErrorException constraint_parameter_refs(tm, c1)
        # test normal
        tm.ext[:TransData].infconstr_to_params[c1] = (par,)
        @test constraint_parameter_refs(tm, c1) == (par,)
        # undo changes
        delete!(tm.ext[:TransData].infconstr_to_params, c1)
    end
    # test parameter_refs for infinite constraint with 1 input
    @testset "parameter_refs (Infinite)" begin
        # test error
        @test_throws ErrorException parameter_refs(c1)
        # test normal
        tm.ext[:TransData].infconstr_to_params[c1] = (par,)
        @test parameter_refs(c1) == (par,)
    end
    # test constraint_parameter_refs for measure constraint with 2 inputs
    @testset "constraint_parameter_refs (Model, Measure)" begin
        # test error
        @test_throws ErrorException constraint_parameter_refs(tm, c2)
        # test normal
        tm.ext[:TransData].measconstr_to_params[c2] = (par,)
        @test constraint_parameter_refs(tm, c2) == (par,)
        # undo changes
        delete!(tm.ext[:TransData].measconstr_to_params, c2)
    end
    # test parameter_refs for measure constraint with 1 input
    @testset "parameter_refs (Measure)" begin
        # test error
        @test_throws ErrorException parameter_refs(c2)
        # test normal
        tm.ext[:TransData].measconstr_to_params[c2] = (par,)
        @test parameter_refs(c2) == (par,)
    end
end


text:
module CoreRenderer

import DataFrames

import Dance.HTMLRenderer
import Dance.JSONRenderer
import Dance.StaticRenderer
import Dance.Router

"""
    render(;headers::Dict{String, String}=Dict("" => ""), status_code::Int16, endpoint::String, data::Union{DataFrames.DataFrame, Dict, String}, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}

Generic web output rendering function

Depending on `endpoint` field, render JSON or HTML string

Cannot set `headers type to `Dict{String, String} here, as can be blank (status_code !=200)
"""
function render(;headers::Dict{String, String}=Dict("" => ""), status_code::Int16, endpoint::String, data::Union{DataFrames.DataFrame, Dict, String}, html_file::String) :: Dict{Symbol, Union{Dict{String, String}, Int16, String}}
    if endpoint==Router.EP_HTML
        HTMLRenderer.render(;headers=headers, status_code=status_code, data=data, html_file=html_file)
    elseif endpoint==Router.EP_JSON
        if isa(data, String)
            data::Dict = Dict(:error => data)
        end
        JSONRenderer.render(;headers=headers, status_code=status_code, data=data)
    else
        StaticRenderer.render(;headers=headers, status_code=status_code, data=data)
    end
end

end


text:
import Base: getindex, first, last, size, show

struct ContinuousRange{T<:AbstractFloat} <: AbstractRange{T}
    a::T
    b::T

    function ContinuousRange(a::T, b::T) where {T<:AbstractFloat}
        a = clamp(a, 0., Inf)
        b = clamp(b, 0., Inf)
        new{T}(a, b)
    end
end

(::Colon)(a::T, b::T) where {T<:AbstractFloat} = ContinuousRange(a,b)

ϵ = 1e-15;
start(c::ContinuousRange) = c.a
stop(c::ContinuousRange) = c.b
firstidx(c::ContinuousRange) = Int(ceil(start(c)+ϵ))
lastidx(c::ContinuousRange) = Int(ceil(stop(c)))
indexify(c::ContinuousRange) = firstidx(c):lastidx(c)
size(c::ContinuousRange) = c.b - c.a
show(io::IO, c::ContinuousRange) = print(io, repr(start(c)), ':', repr(stop(c)))

function getindex(v::AbstractArray{T}, c::ContinuousRange, rest...) where {T<:AbstractFloat}
    ai,bi = firstidx(c), lastidx(c)
    ai = (ai <= 0) ? 1 : ai
    bi = (bi == 0) ? -1 : bi

    if ai > bi
        if (length(rest) > 0) && (typeof(rest[1]) == EllipsisNotation.Ellipsis)
            return zeros(1, size(v)[2:end]...)
        else
            u = [i[1] for i in size.(rest) if length(i)!=0]
            return zeros(1, u...)
        end
    end

    rv = getindex(v, ai:bi, rest...)
    if ai == bi
        rv[1,..] .*= size(c)
    else
        rv[1,..] .*= (firstidx(c) - start(c))
        rv[end,..] .*= (stop(c) - lastidx(c) + 1)
    end
    return rv
end

getindex(a::AbstractArray, c::ContinuousRange, rest...) = getindex(convert(Array{Float64}, a), c, rest...)
getindex(a::AbstractArray, i::Int, c::ContinuousRange, rest...) = getindex(a[i,..], c, rest...)