text:
export set_mesh_field_using_materials

function set_mesh_field_using_materials(materials::Vector{Material})
    @info "Setting mesh size field using materials"
    # Setup a field with a constant mesh size material.mesh_size inside the
    # entities that make up that material
    field_ids = Vector{Int32}(undef, length(materials))
    for (i, material) in enumerate(materials)
        fid = gmsh.model.mesh.field.add("Constant")
        gmsh.model.mesh.field.set_number(fid, "VIn", material.mesh_size)
        field_ids[i] = fid
    end
    # Populate each of the fields with the highest dimensional entities in the material
    # physical group
    material_dict = Dict{String, Vector{Tuple{Int32, Int32}}}()
    material_names = ["Material: "*mat.name for mat in materials]
    for name in material_names
        material_dict[name] = Tuple{Int32, Int32}[]
    end
    for group in gmsh.model.get_physical_groups()
        gdim, gnum = group
        name = gmsh.model.get_physical_name(gdim, gnum)
        if startswith(name, "Material: ")
            tags = gmsh.model.get_entities_for_physical_group(gdim, gnum)
            dtags = [(gdim, tag) for tag in tags]
            # Empty vector, or these ents are higher dim
            if length(material_dict[name]) == 0 || gdim > material_dict[name][1][1]
                material_dict[name] = dtags
            end
        end
    end
    # Assign the entities to the fields
    for (mat_id, fid) in enumerate(field_ids)
        name = material_names[mat_id]
        dtags = material_dict[name]
        edim = dtags[1][1]
        tags = [dt[2] for dt in dtags]
        if edim == 0
            gmsh.model.mesh.field.set_numbers(fid, "PointsList", tags)
        elseif edim == 1
            gmsh.model.mesh.field.set_numbers(fid, "CurvesList", tags)
        elseif edim == 2
            gmsh.model.mesh.field.set_numbers(fid, "SurfacesList", tags)
        else # edim == 3
            gmsh.model.mesh.field.set_numbers(fid, "VolumesList", tags)
        end
    end
    # Create a field that takes the min of each and set as background mesh
    fid = gmsh.model.mesh.field.add("Min")
    gmsh.model.mesh.field.set_numbers(fid, "FieldsList", field_ids)
    gmsh.model.mesh.field.set_as_background_mesh(fid)
    gmsh.option.set_number("Mesh.MeshSizeExtendFromBoundary", 0)
    gmsh.option.set_number("Mesh.MeshSizeFromPoints", 0)
    gmsh.option.set_number("Mesh.MeshSizeFromCurvature", 0)
    return field_ids
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#### Function to define parameters

##### FUNCTION to create tuple of parameters needed for dBdt

# # assign a value for Z, I_K and T as it would be in a loop, just to check if ModelParameters works to generate parameters
# Z = 10
# I_K = 1
# T = 293.15

function ModelParameters(param, T, I_K, Z)
    ### Calculate body masses
    M = BodyMasses(Z=Z)

    ### Calculate carrying capacity
    K = carryingcapacity_BA(I_K, param, M, T)
    #@assert K[2] == 0 # defensive programming
    #@assert K[3] == 0

    ### Calculate producer growth rate
    r = growth_BA(param, M, T)
    @assert r[2] == 0.0 # defensive programming
    @assert r[3] == 0.0

    ### Calculate metabolism
    x = metabolism_BA(param, M, T)
    @assert x[1] == 0.0

    ### Calculate maximum ingestion
    y = max_ingestion_BA(param, M, Z, T)

    ### Calculate half saturation density
    B0 = half_saturation_BA(param, M, Z, T)

    # Combine parameters into Tuple
    ModelParameters = (M = M, K = K, r = r, x = x, y = y, B0 = B0)
    return ModelParameters
end

# p = ModelParameters(param, T, I_K, Z)
# typeof(p)
# p.B0

ext: jl


text:
using NewPkgEval
using Test

# determine the version to use
const version = get(ENV, "JULIA_VERSION", string(VERSION))
const julia = try
    # maybe it already refers to a version in Versions.toml
    v = VersionNumber(version)
    NewPkgEval.prepare_julia(v)
    v
catch
    # maybe it points to a build in Builds.jl
    try
        NewPkgEval.download_julia(version)
    catch
        # assume it points to something in our Git repository
        NewPkgEval.build_julia(version)
    end
end
NewPkgEval.prepare_julia(julia::VersionNumber)

@testset "sandbox" begin
    NewPkgEval.prepare_runner()
    mktemp() do path, io
        NewPkgEval.run_sandboxed_julia(julia, `-e 'print(1337)'`; stdout=io)
        close(io)
        @test read(path, String) == "1337"
    end

    # print versioninfo so we can verify in CI logs that the correct version is used
    NewPkgEval.run_sandboxed_julia(julia, `-e 'using InteractiveUtils; versioninfo()'`)
end

const pkgnames = ["TimerOutputs", "Crayons", "Example"]

@testset "low-level interface" begin
    NewPkgEval.prepare_registry()

    pkgs = NewPkgEval.read_pkgs(pkgnames)

    # timeouts
    results = NewPkgEval.run([julia], pkgs; time_limit = 0.1)
    for pkg in pkgnames
        @test results[(julia=julia,pkg=pkg)] == :kill
    end
end

@testset "main entrypoint" begin
    results = NewPkgEval.run([julia], pkgnames)
    for pkg in pkgnames
        @test results[(julia=julia,pkg=pkg)] == :ok
        output = read(joinpath(NewPkgEval.log_path(julia), "$pkg.log"), String)
        @test occursin("Testing $pkg tests passed", output)
    end
end


text:
using VLMinorityGameKit

# setup call -
memory = 3

# make call -
strategy_object = build_basic_agent_strategy(memory)


text:
module MuladdMacro

"""
    postwalk(f, expr)

Apply the function `f` to each subexpression within expr.

Example:
```
julia> double(x::Int)=2x;  # Int Literals

julia> double(x::Any)=x;   # Everything else (e.g. Expressions and Symbols)

julia> MuladdMacro.postwalk(double, :(1+2*3));

julia> MuladdMacro.postwalk(double, :(1+2*3))
:(2 + 4 * 6)
```
"""
postwalk(f, x) = f(x)
postwalk(f, expr::Expr) = f(Expr(expr.head, (postwalk(f, arg) for arg in expr.args)...))

"""
    @muladd

Convert every combination of addition/subtraction and multiplication to a call of `muladd`.

If both of the involved operators are dotted, `muladd` is applied as a dot call.
The order of summation might be changed.
"""
macro muladd(ex)
    esc(to_muladd(ex))
end

"""
    to_muladd(ex)

Convert every combination of addition/subtraction and multiplication in expression `ex` to a call of `muladd`.

If both of the involved operators are dotted, `muladd` is applied as a dot call
The order of summation might be changed.
"""
function to_muladd(ex)
    postwalk(ex) do x
        # Modify summations
        issum(x) && return sum_to_muladd(x)

        # Modify subtractions
        issub(x) && return sub_to_muladd(x)

        return x
    end
end

"""
    sum_to_muladd(ex)

Replace sum `ex` by sequence of `muladd` if possible. Hereby the order of summation might be changed.
"""
function sum_to_muladd(ex)
    # Retrieve summands
    summands = args(ex)

    # Skip further calculations if no summand is a multiplication
    dotcall = isdotcall(ex)
    any(x -> ismul(x, dotcall), summands) || return ex

    # Split summands into two groups, one with expressions
    # of multiplications and one with other expressions
    mulsummands = filter(x -> ismul(x, dotcall), summands)
    oddsummands = filter(x -> !ismul(x, dotcall), summands)

    # If all summands are multiplications the first one is not reduced
    isempty(oddsummands) && push!(oddsummands, popfirst!(mulsummands))

    # Reduce sum to a composition of muladd
    foldl(mulsummands; init = newargs(ex, oddsummands...)) do s₁, s₂
        newmuladd(splitargs(s₂)..., s₁, dotcall)
    end
end

"""
    sub_to_muladd(ex)

Replace subtraction `ex` by `muladd` if possible.
"""
function sub_to_muladd(ex)
    # Retrieve operands
    x, y = args(ex)

    # Modify subtraction if possible
    dotcall = isdotcall(ex)
    if ismul(y, dotcall)
        y₁, y₂ = splitargs(y)
        return newmuladd(:(-$y₁), y₂, x, dotcall)
    elseif ismul(x, dotcall)
        return newmuladd(splitargs(x)..., :(-$y), dotcall)
    end

    return ex
end

"""
    iscall(ex, op)

Determine whether `ex` is a call of operation `op` with at least two arguments.
"""
iscall(ex::Expr, op) =
    ex.head == :call && length(ex.args) > 2 && ex.args[1] == op
iscall(ex, op) = false

"""
    isdotcall(ex)

Determine whether `ex` is a dot call.
"""
isdotcall(ex::Expr) =
    (ex.head == :. && length(ex.args) == 2 && Meta.isexpr(ex.args[2], :tuple)) ||
    (ex.head == :call && !isempty(ex.args) && startswith(string(ex.args[1]), '.'))
isdotcall(ex) = false

"""
    isdotcall(ex, op)

Determine whether `ex` is a dot call of operation `op` with at least two arguments.
"""
isdotcall(ex::Expr, op) =
    (ex.head == :. && length(ex.args) == 2 && ex.args[1] == op && Meta.isexpr(ex.args[2], :tuple) && length(ex.args[2].args) > 1) ||
    (ex.head == :call && length(ex.args) > 2 && ex.args[1] == Symbol('.', op))
isdotcall(ex, op) = false

"""
    issum(ex)

Determine whether `ex` is a sum.
"""
issum(ex) = iscall(ex, :+) || isdotcall(ex, :+)

"""
    issub(ex)

Determine whether `ex` is a subtraction.
"""
issub(ex) = iscall(ex, :-) || isdotcall(ex, :-)

"""
    ismul(ex, dot::Bool)

Determine whether expression `ex` is a multiplication that is dotted if
`dot` is `true` and not dotted otherwise.
"""
ismul(ex, dot::Bool) = dot ? isdotcall(ex, :*) : iscall(ex, :*)

"""
    newmuladd(x, y, z, dot::Bool)

Return expression `(muladd).(x, y, z)` if `dot` is `true` and
`muladd(x, y, z)` otherwise.
"""
function newmuladd(x, y, z, dot::Bool)
    # Quoting seems to be required for the @. macro to work
    if dot
        :(($(Meta.quot(Base.muladd))).($x, $y, $z))
    else
        :($(Meta.quot(Base.muladd))($x, $y, $z))
    end
end

"""
    args(ex)

Return arguments of function call in `ex`.
"""
function args(ex::Expr)
    if ex.head == :call && length(ex.args) > 1
        return ex.args[2:end]
    end

    if ex.head == :. && length(ex.args) == 2 && Meta.isexpr(ex.args[2], :tuple) && !isempty(ex.args[2].args)
        return ex.args[2].args
    end

    error("expression is not a function call with arguments")
end

"""
    splitargs(ex)

Split arguments of function call `ex` before last argument and combine first
arguments to one expression if possible.
"""
function splitargs(ex)
    if ex.head == :call && length(ex.args) > 2
        x = ex.args[2:end-1]
        y = ex.args[end]
    elseif ex.head == :. && length(ex.args) == 2 && Meta.isexpr(ex.args[2], :tuple) && length(ex.args[2].args) > 1
        x = ex.args[2].args[1:end-1]
        y = ex.args[2].args[end]
    else
        error("cannot split arguments")
    end

    return newargs(ex, x...), y
end

"""
    newargs(ex, args)

Create new expression of function call `ex` with arguments `args`.

Unary function calls are not considered, i.e. if only one function
argument is provided it is returned.
"""
function newargs(ex::Expr, args...)
    # Return single argument
    length(args) == 1 && return args[1]

    # Create function calls with new arguments
    if ex.head == :call && !isempty(ex.args)
        return Expr(:call, ex.args[1], args...)
    end
    if ex.head == :. && length(ex.args) == 2 && Meta.isexpr(ex.args[2], :tuple)
        return Expr(:., ex.args[1], Expr(:tuple, args...))
    end

    error("expression is not a function call")
end

export @muladd

end # module


text:
# QCP example
#    maximize x + y
#
#    s.t.  x, y >= 0
#          x^2 + y^2 <= 1
#
#    solution: (0.71, 0.71) objv = 1.414

using CPLEX

env = CPLEX.Env()

model = CPLEX.Model(env, "qcqp_01")
CPLEX.set_sense!(model, :Max)

CPLEX.add_vars!(model, [1., 1.], 0., Inf)

 # add_qpterms!(model, linearindices, linearcoeffs, qrowinds, qcolinds, qcoeffs, sense, rhs)
CPLEX.add_qconstr!(model, [], [], [1, 2], [1, 2], [1, 1.], '<', 1.0)

CPLEX.optimize!(model)

println("sol = $(CPLEX.get_solution(model))")
println("obj = $(CPLEX.get_objval(model))")


text:
using Test, EllZeroTrendFiltering

t = range(0, stop=4π, length=50)
g = sin

@time ℓ = compute_transition_costs(g, t);
cost_last = QuadraticPolynomial(0.0, 0.0, 0.0)

@time Λ =  pwq_dp_constrained(ℓ, cost_last, 7);

I_sols = [      [[1, 50]],
                [[1, 34, 50], [1, 17, 50]],
                [[1, 21, 30, 50]],
                [[1, 8, 19, 30, 50], [1, 21, 32, 43, 50]],
                [[1, 8, 19, 32, 43, 50]]
         ]

f_costs = [5.328255648628215, 4.8783936257642315, 1.7402065022125042, 0.9196880458290417, 0.09174442455649423]

@testset "Optimal Fit m=$m" for m = 1:5
    I, _, f = recover_solution(Λ[m, 1], ℓ, cost_last)

    @test I ∈ I_sols[m]
    @test f ≈ f_costs[m]    atol = 1e-8
end

##
@testset "Brute force m=$m" for m = 1:4
    I_bf, _, f_bf = brute_force_search(ℓ, cost_last, m);

    @test I_bf ∈  I_sols[m]
    @test f_bf ≈ f_costs[m]    atol = 1e-8
end

##

@testset "Regularize ζ=$ζ" for ζ in 0.1:0.1:2

    Λ_reg = pwq_dp_regularized(ℓ, cost_last, ζ)

    # recover_optimal_index_set returns the cost inclusive the regularization penality,
    # revober optimal solution does not do so. It is arguably more interesting
    # to test cost including regularization.
    I, _, f_reg = recover_optimal_index_set(Λ_reg[1])

    # Use the costs above to find out how many segments the solution should contain
    m_expected = argmin([f_costs[m] + m*ζ for m=1:5])

    @test m_expected == length(I) - 1
    @test I ∈ I_sols[m_expected]
    @test f_reg ≈ f_costs[m_expected] + ζ*m_expected
end


text:
abstract type Bound{T} end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, bound::Bound{T}) where {T}
    summary(io, bound)
end

### General functions
@inline getT(bound::Bound{T}) where T = T

### Constraints and derivatives
## Position level constraint wrappers
g(contact::Bound, body::Body, Δt) = g(contact, body.state, Δt)

# @inline g(bound::Bound{T}) where T = zero(T)

# @inline ∂g∂pos(bound::Bound{T}) where T = szeros(T, 6)
# @inline ∂g∂vel(bound::Bound{T}) where T = szeros(T, 6)
# @inline schurf(bound::Bound{T}) where T = szeros(T, 6)
# @inline schurD(bound::Bound{T}) where T = szeros(T, 6, 6)


text:
module Sort

using Base.Order

import
    Base.sort,
    Base.sort!,
    Base.issorted,
    Base.sortperm

export # also exported by Base
    # order-only:
    issorted,
    select,
    select!,
    searchsorted,
    searchsortedfirst,
    searchsortedlast,
    # order & algorithm:
    sort,
    sort!,
    sortperm,
    sortrows,
    sortcols,
    # algorithms:
    InsertionSort,
    QuickSort,
    MergeSort

export # not exported by Base
    Algorithm,
    DEFAULT_UNSTABLE,
    DEFAULT_STABLE,
    SMALL_ALGORITHM,
    SMALL_THRESHOLD

## functions requiring only ordering ##

function issorted(itr, order::Ordering)
    state = start(itr)
    done(itr,state) && return true
    prev, state = next(itr, state)
    while !done(itr, state)
        this, state = next(itr, state)
        lt(order, this, prev) && return false
        prev = this
    end
    return true
end
issorted(itr;
    lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) =
    issorted(itr, ord(lt,by,rev,order))

function select!(v::AbstractVector, k::Int, lo::Int, hi::Int, o::Ordering)
    lo <= k <= hi || error("select index $k is out of range $lo:$hi")
    @inbounds while lo < hi
        if hi-lo == 1
            if lt(o, v[hi], v[lo])
                v[lo], v[hi] = v[hi], v[lo]
            end
            return v[k]
        end
        pivot = v[(lo+hi)>>>1]
        i, j = lo, hi
        while true
            while lt(o, v[i], pivot); i += 1; end
            while lt(o, pivot, v[j]); j -= 1; end
            i <= j || break
            v[i], v[j] = v[j], v[i]
            i += 1; j -= 1
        end
        if k <= j
            hi = j
        elseif i <= k
            lo = i
        else
            return pivot
        end
    end
    return v[lo]
end

function select!(v::AbstractVector, r::OrdinalRange, lo::Int, hi::Int, o::Ordering)
    isempty(r) && (return v[r])
    a, b = extrema(r)
    lo <= a <= b <= hi || error("selection $r is out of range $lo:$hi")
    @inbounds while true
        if lo == a && hi == b
            sort!(v, lo, hi, DEFAULT_UNSTABLE, o)
            return v[r]
        end
        pivot = v[(lo+hi)>>>1]
        i, j = lo, hi
        while true
            while lt(o, v[i], pivot); i += 1; end
            while lt(o, pivot, v[j]); j -= 1; end
            i <= j || break
            v[i], v[j] = v[j], v[i]
            i += 1; j -= 1
        end
        if b <= j
            hi = j
        elseif i <= a
            lo = i
        else
            a <= j && select!(v, a, lo,  j, o)
            b >= i && select!(v, b,  i, hi, o)
            sort!(v, a, b, DEFAULT_UNSTABLE, o)
            return v[r]
        end
    end
end

select!(v::AbstractVector, k::Union(Int,OrdinalRange), o::Ordering) = select!(v,k,1,length(v),o)
select!(v::AbstractVector, k::Union(Int,OrdinalRange);
    lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) =
    select!(v, k, ord(lt,by,rev,order))

select(v::AbstractVector, k::Union(Int,OrdinalRange); kws...) = select!(copy(v), k; kws...)

# reference on sorted binary search:
#   http://www.tbray.org/ongoing/When/200x/2003/03/22/Binary

# index of the first value of vector a that is greater than or equal to x;
# returns length(v)+1 if x is greater than all values in v.
function searchsortedfirst(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering)
    lo = lo-1
    hi = hi+1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        if lt(o, v[m], x)
            lo = m
        else
            hi = m
        end
    end
    return hi
end

# index of the last value of vector a that is less than or equal to x;
# returns 0 if x is less than all values of v.
function searchsortedlast(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering)
    lo = lo-1
    hi = hi+1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        if lt(o, x, v[m])
            hi = m
        else
            lo = m
        end
    end
    return lo
end

# returns the range of indices of v equal to x
# if v does not contain x, returns a 0-length range
# indicating the insertion point of x
function searchsorted(v::AbstractVector, x, lo::Int, hi::Int, o::Ordering)
    lo = lo-1
    hi = hi+1
    @inbounds while lo < hi-1
        m = (lo+hi)>>>1
        if lt(o, v[m], x)
            lo = m
        elseif lt(o, x, v[m])
            hi = m
        else
            a = searchsortedfirst(v, x, max(lo,1), m, o)
            b = searchsortedlast(v, x, m, min(hi,length(v)), o)
            return a:b
        end
    end
    return lo+1:hi-1
end

function searchsortedlast{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering)
    if step(a) == 0
        lt(o, x, first(a)) ? 0 : length(a)
    else
        n = max(min(iround((x-first(a))/step(a))+1,length(a)),1)
        lt(o, x, a[n]) ? n-1 : n
    end
end

function searchsortedfirst{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering)
    if step(a) == 0
        lt(o, first(a), x) ? length(a)+1 : 1
    else
        n = max(min(iround((x-first(a))/step(a))+1,length(a)),1)
        lt(o, a[n] ,x) ? n+1 : n
    end
end

function searchsortedlast{T<:Integer}(a::Range{T}, x::Real, o::DirectOrdering)
    if step(a) == 0
        lt(o, x, first(a)) ? 0 : length(a)
    else
        max(min(fld(ifloor(x)-first(a),step(a))+1,length(a)),0)
    end
end

function searchsortedfirst{T<:Integer}(a::Range{T}, x::Real, o::DirectOrdering)
    if step(a) == 0
        lt(o, first(a), x) ? length(a)+1 : 1
    else
        max(min(-fld(ifloor(-x)+first(a),step(a))+1,length(a)+1),1)
    end
end

searchsorted{T<:Real}(a::Range{T}, x::Real, o::DirectOrdering) =
    searchsortedfirst(a,x,o):searchsortedlast(a,x,o)

for s in {:searchsortedfirst, :searchsortedlast, :searchsorted}
    @eval begin
        $s(v::AbstractVector, x, o::Ordering) = $s(v,x,1,length(v),o)
        $s(v::AbstractVector, x;
           lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) =
            $s(v,x,ord(lt,by,rev,order))
    end
end

## sorting algorithms ##

abstract Algorithm

immutable InsertionSortAlg <: Algorithm end
immutable QuickSortAlg     <: Algorithm end
immutable MergeSortAlg     <: Algorithm end

const InsertionSort = InsertionSortAlg()
const QuickSort     = QuickSortAlg()
const MergeSort     = MergeSortAlg()

const DEFAULT_UNSTABLE = QuickSort
const DEFAULT_STABLE   = MergeSort
const SMALL_ALGORITHM  = InsertionSort
const SMALL_THRESHOLD  = 20

function sort!(v::AbstractVector, lo::Int, hi::Int, ::InsertionSortAlg, o::Ordering)
    @inbounds for i = lo+1:hi
        j = i
        x = v[i]
        while j > lo
            if lt(o, x, v[j-1])
                v[j] = v[j-1]
                j -= 1
                continue
            end
            break
        end
        v[j] = x
    end
    return v
end

function sort!(v::AbstractVector, lo::Int, hi::Int, a::QuickSortAlg, o::Ordering)
    @inbounds while lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)
        mi = (lo+hi)>>>1
        if lt(o, v[mi], v[lo])
            v[mi], v[lo] = v[lo], v[mi]
        end
        if lt(o, v[hi], v[mi])
            v[hi], v[mi] = v[mi], v[hi]
        end
        if lt(o, v[mi], v[lo])
            v[mi], v[lo] = v[lo], v[mi]
        end
        v[mi], v[lo] = v[lo], v[mi]
        i, j = lo, hi
        pivot = v[lo]
        while true
            i += 1; j -= 1;
            while lt(o, v[i], pivot); i += 1; end;
            while lt(o, pivot, v[j]); j -= 1; end;
            i >= j && break
            v[i], v[j] = v[j], v[i]
        end
        v[j], v[lo] = v[lo], v[j]
        lo < (j-1) && sort!(v, lo, j-1, a, o)
        lo = j+1
    end
    return v
end

function sort!(v::AbstractVector, lo::Int, hi::Int, a::MergeSortAlg, o::Ordering, t=similar(v))
    @inbounds if lo < hi
        hi-lo <= SMALL_THRESHOLD && return sort!(v, lo, hi, SMALL_ALGORITHM, o)

        m = (lo+hi)>>>1
        sort!(v, lo,  m,  a, o, t)
        sort!(v, m+1, hi, a, o, t)

        i, j = 1, lo
        while j <= m
            t[i] = v[j]
            i += 1
            j += 1
        end

        i, k = 1, lo
        while k < j <= hi
            if lt(o, v[j], t[i])
                v[k] = v[j]
                j += 1
            else
                v[k] = t[i]
                i += 1
            end
            k += 1
        end
        while k < j
            v[k] = t[i]
            k += 1
            i += 1
        end
    end

    return v
end

## generic sorting methods ##

defalg(v::AbstractArray) = DEFAULT_STABLE
defalg{T<:Number}(v::AbstractArray{T}) = DEFAULT_UNSTABLE

sort!(v::AbstractVector, alg::Algorithm, order::Ordering) = sort!(v,1,length(v),alg,order)
sort!(v::AbstractVector; alg::Algorithm=defalg(v),
    lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) =
    sort!(v, alg, ord(lt,by,rev,order))

sort(v::AbstractVector; kws...) = sort!(copy(v); kws...)

## sortperm: the permutation to sort an array ##

sortperm(v::AbstractVector; alg::Algorithm=DEFAULT_UNSTABLE,
    lt::Function=isless, by::Function=identity, rev::Bool=false, order::Ordering=Forward) =
    sort!([1:length(v)], alg, Perm(ord(lt,by,rev,order),v))

## sorting multi-dimensional arrays ##

sort(A::AbstractArray, dim::Integer; kws...) = mapslices(a->sort(a; kws...), A, [dim])

function sortrows(A::AbstractMatrix; kws...)
    c = 1:size(A,2)
    rows = [ sub(A,i,c) for i=1:size(A,1) ]
    p = sortperm(rows; kws..., order=Lexicographic)
    A[p,:]
end

function sortcols(A::AbstractMatrix; kws...)
    r = 1:size(A,1)
    cols = [ sub(A,r,i) for i=1:size(A,2) ]
    p = sortperm(cols; kws..., order=Lexicographic)
    A[:,p]
end

## fast clever sorting for floats ##

module Float
using ..Sort
using ...Order

import Core.Intrinsics: unbox, slt_int
import ..Sort: sort!
import ...Order: lt, DirectOrdering, uint_mapping

typealias Floats Union(Float32,Float64)

immutable Left <: Ordering end
immutable Right <: Ordering end

left(::DirectOrdering) = Left()
right(::DirectOrdering) = Right()

left(o::Perm) = Perm(left(o.order), o.data)
right(o::Perm) = Perm(right(o.order), o.data)

lt{T<:Floats}(::Left, x::T, y::T) = slt_int(unbox(T,y),unbox(T,x))
lt{T<:Floats}(::Right, x::T, y::T) = slt_int(unbox(T,x),unbox(T,y))

isnan(o::DirectOrdering, x::Floats) = (x!=x)
isnan(o::Perm, i::Int) = isnan(o.order,o.data[i])

function nans2left!(v::AbstractVector, o::Ordering, lo::Int=1, hi::Int=length(v))
    i = lo
    @inbounds while i <= hi && isnan(o,v[i])
        i += 1
    end
    j = i + 1
    @inbounds while j <= hi
        if isnan(o,v[j])
            v[i], v[j] = v[j], v[i]
            i += 1
        end
        j += 1
    end
    return i, hi
end
function nans2right!(v::AbstractVector, o::Ordering, lo::Int=1, hi::Int=length(v))
    i = hi
    @inbounds while lo <= i && isnan(o,v[i])
        i -= 1
    end
    j = i - 1
    @inbounds while lo <= j
        if isnan(o,v[j])
            v[i], v[j] = v[j], v[i]
            i -= 1
        end
        j -= 1
    end
    return lo, i
end

nans2end!(v::AbstractVector, o::ForwardOrdering) = nans2right!(v,o)
nans2end!(v::AbstractVector, o::ReverseOrdering) = nans2left!(v,o)
nans2end!{O<:ForwardOrdering}(v::AbstractVector{Int}, o::Perm{O}) = nans2right!(v,o)
nans2end!{O<:ReverseOrdering}(v::AbstractVector{Int}, o::Perm{O}) = nans2left!(v,o)

issignleft(o::ForwardOrdering, x::Floats) = lt(o, x, zero(x))
issignleft(o::ReverseOrdering, x::Floats) = lt(o, x, -zero(x))
issignleft(o::Perm, i::Int) = issignleft(o.order, o.data[i])

function fpsort!(v::AbstractVector, a::Algorithm, o::Ordering)
    i, j = lo, hi = nans2end!(v,o)
    @inbounds while true
        while i <= j &&  issignleft(o,v[i]); i += 1; end
        while i <= j && !issignleft(o,v[j]); j -= 1; end
        i <= j || break
        v[i], v[j] = v[j], v[i]
        i += 1; j -= 1
    end
    sort!(v, lo, j,  a, left(o))
    sort!(v, i,  hi, a, right(o))
    return v
end

sort!{T<:Floats}(v::AbstractVector{T}, a::Algorithm, o::DirectOrdering) = fpsort!(v,a,o)
sort!{O<:DirectOrdering,T<:Floats}(v::Vector{Int}, a::Algorithm, o::Perm{O,Vector{T}}) = fpsort!(v,a,o)

end # module Sort.Float

end # module Sort

max_stars_repo_name: porterjamesj/julia

max_issues_repo_name: porterjamesj/julia

max_forks_repo_name: porterjamesj/julia


text:
### This file contains convenience functions for constructing LoopSets.

# function strip_unneeded_const_deps!(ls::LoopSet)
#     for op ∈ operations(ls)
#         if isconstant(op) && iszero(length(reducedchildren(op)))
#             op.dependencies = NODEPENDENCY
#         end
#     end
#     ls
# end

function Base.copyto!(ls::LoopSet, q::Expr)
  q.head === :for || throw("Expression must be a for loop.")
  add_loop!(ls, q, 8)
    # strip_unneeded_const_deps!(ls)
end

function add_ci_call!(q::Expr, @nospecialize(f), args, syms, i, valarg = nothing, mod = nothing)
    call = if f isa Core.SSAValue
        Expr(:call, syms[f.id])
    else
        Expr(:call, f)
    end
    for arg ∈ @view(args[2:end])
        if arg isa Core.SSAValue
            push!(call.args, syms[arg.id])
        else
            push!(call.args, arg)
        end
    end
    if mod !== nothing # indicates it's `vmaterialize(!)`
        push!(call.args, Expr(:call, Expr(:curly, :Val, QuoteNode(mod))), valarg)
    end
    push!(q.args, Expr(:(=), syms[i], call))
end

function substitute_broadcast(q::Expr, mod::Symbol, inline::Bool, u₁::Int8, u₂::Int8, v::Int8, threads::Int, warncheckarg::Int)
  ci = first(Meta.lower(LoopVectorization, q).args).code
  nargs = length(ci)-1
  ex = Expr(:block)
  syms = [gensym() for _ ∈ 1:nargs]
  configarg = (inline,u₁,u₂,v,true,threads,warncheckarg)
  unroll_param_tup = Expr(:call, lv(:avx_config_val), :(Val{$configarg}()), staticexpr(0))
  for n ∈ 1:nargs
    ciₙ = ci[n]
    ciₙargs = ciₙ.args
    f = first(ciₙargs)
    if ciₙ.head === :(=)
      push!(ex.args, Expr(:(=), f, syms[((ciₙargs[2])::Core.SSAValue).id]))
    elseif isglobalref(f, Base, :materialize!)
      add_ci_call!(ex, lv(:vmaterialize!), ciₙargs, syms, n, unroll_param_tup, mod)
    elseif isglobalref(f, Base, :materialize)
      add_ci_call!(ex, lv(:vmaterialize), ciₙargs, syms, n, unroll_param_tup, mod)
    else
      add_ci_call!(ex, f, ciₙargs, syms, n)
    end
  end
  esc(ex)
end

function LoopSet(q::Expr, mod::Symbol = :Main)
  contract_pass!(q)
  ls = LoopSet(mod)
  copyto!(ls, q)
  resize!(ls.loop_order, num_loops(ls))
  ls
end
LoopSet(q::Expr, m::Module) = LoopSet(macroexpand(m, q)::Expr, Symbol(m))

function loopset(q::Expr) # for interactive use only
    ls = LoopSet(q)
    set_hw!(ls)
    ls
end

function check_macro_kwarg(arg, inline::Bool, check_empty::Bool, u₁::Int8, u₂::Int8, v::Int8, threads::Int, warncheckarg::Int)
    ((arg.head === :(=)) && (length(arg.args) == 2)) || throw(ArgumentError("macro kwarg should be of the form `argname = value`."))
    kw = (arg.args[1])::Symbol
    value = (arg.args[2])
    if kw === :inline
        inline = value::Bool
    elseif kw === :unroll
        if value isa Integer
            u₁ = convert(Int8, value)::Int8
        elseif Meta.isexpr(value,:tuple,2)
            u₁ = convert(Int8, value.args[1])::Int8
            u₂ = convert(Int8, value.args[2])::Int8
        else
            throw(ArgumentError("Don't know how to process argument in `unroll=$value`."))
        end
    elseif kw === :vectorize
        v = convert(Int8, value)
    elseif kw === :check_empty
        check_empty = value::Bool
    elseif kw === :thread
        if value isa Bool
            threads = Core.ifelse(value::Bool, -1, 1)
        elseif value isa Integer
            threads = max(1, convert(Int,value)::Int)
        else
            throw(ArgumentError("Don't know how to process argument in `thread=$value`."))
        end
    elseif kw === :warn_check_args
        warncheckarg = convert(Int, value)::Int
    else
        throw(ArgumentError("Received unrecognized keyword argument $kw. Recognized arguments include:\n`inline`, `unroll`, `check_empty`, and `thread`."))
    end
    inline, check_empty, u₁, u₂, v, threads, warncheckarg
end
function process_args(args; inline::Bool = false, check_empty::Bool = false, u₁::Int8 = zero(Int8), u₂::Int8 = zero(Int8), v::Int8 = zero(Int8), threads::Int = 1, warncheckarg::Int = 1)
  for arg ∈ args
    inline, check_empty, u₁, u₂, v, threads, warncheckarg = check_macro_kwarg(arg, inline, check_empty, u₁, u₂, v, threads, warncheckarg)
  end
  inline, check_empty, u₁, u₂, v, threads, warncheckarg
end
function turbo_macro(mod, src, q, args...)
  q = macroexpand(mod, q)
  if q.head === :for
    ls = LoopSet(q, mod)
    inline, check_empty, u₁, u₂, v, threads, warncheckarg = process_args(args)
    esc(setup_call(ls, q, src, inline, check_empty, u₁, u₂, v, threads, warncheckarg))
  else
    inline, check_empty, u₁, u₂, v, threads, warncheckarg = process_args(args, inline=true)
    substitute_broadcast(q, Symbol(mod), inline, u₁, u₂, v, threads, warncheckarg)
  end
end
"""
    @turbo

Annotate a `for` loop, or a set of nested `for` loops whose bounds are constant across iterations, to optimize the computation. For example:

    function AmulB!(C, A, B)
        @turbo for m ∈ indices((A,C), 1), n ∈ indices((B,C), 2) # indices((A,C),1) == axes(A,1) == axes(C,1)
            Cₘₙ = zero(eltype(C))
            for k ∈ indices((A,B), (2,1)) # indices((A,B), (2,1)) == axes(A,2) == axes(B,1)
                Cₘₙ += A[m,k] * B[k,n]
            end
            C[m,n] = Cₘₙ
        end
    end

The macro models the set of nested loops, and chooses an ordering of the three loops to minimize predicted computation time.

Current limitations:

1. It assumes that loop iterations are independent.
2. It does not perform bounds checks.
3. It assumes that each loop iterates at least once. (Use `@turbo check_empty=true` to lift this assumption.)
4. That there is only one loop at each level of the nest.

It may also apply to broadcasts:

```jldoctest
julia> using LoopVectorization

julia> a = rand(100);

julia> b = @turbo exp.(2 .* a);

julia> c = similar(b);

julia> @turbo @. c = exp(2a);

julia> b ≈ c
true
```

# Extended help

Advanced users can customize the implementation of the `@turbo`-annotated block
using keyword arguments:

```julia
@turbo inline=false unroll=2 body
```

where `body` is the code of the block (e.g., `for ... end`).

`thread` is either a Boolean, or an integer.
The integer's value indicates the number of threads to use.
It is clamped to be between `1` and `min(Threads.nthreads(),LoopVectorization.num_cores())`.
`false` is equivalent to `1`, and `true` is equivalent to `min(Threads.nthreads(),LoopVectorization.num_cores())`.

`inline` is a Boolean. When `true`, `body` will be directly inlined
into the function (via a forced-inlining call to `_turbo_!`).
When `false`, it wont force inlining of the call to `_turbo_!` instead, letting Julia's own inlining engine
determine whether the call to `_turbo_!` should be inlined. (Typically, it won't.)
Sometimes not inlining can lead to substantially worse code generation, and >40% regressions, even in very
large problems (2-d convolutions are a case where this has been observed).
One can find some circumstances where `inline=true` is faster, and other circumstances
where `inline=false` is faster, so the best setting may require experimentation. By default, the macro
tries to guess. Currently the algorithm is simple: roughly, if there are more than two dynamically sized loops
or and no convolutions, it will probably not force inlining. Otherwise, it probably will.

`check_empty` (default is `false`) determines whether or not it will check if any of the iterators are empty.
If false, you must ensure yourself that they are not empty, else the behavior of the loop is undefined and
(like with `@inbounds`) segmentation faults are likely.

`unroll` is an integer that specifies the loop unrolling factor, or a
tuple `(u₁, u₂) = (4, 2)` signaling that the generated code should unroll more than
one loop. `u₁` is the unrolling factor for the first unrolled loop and `u₂` for the next (if present),
but it applies to the loop ordering and unrolling that will be chosen by LoopVectorization,
*not* the order in `body`.
`uᵢ=0` (the default) indicates that LoopVectorization should pick its own value,
and `uᵢ=-1` disables unrolling for the correspond loop.

The `@turbo` macro also checks the array arguments using `LoopVectorization.check_args` to try and determine
if they are compatible with the macro. If `check_args` returns false, a fall back loop annotated with `@inbounds`
and `@fastmath` is generated. Note that `VectorizationBase` provides functions such as `vadd` and `vmul` that will
ignore `@fastmath`, preserving IEEE semantics both within `@turbo` and `@fastmath`.
`check_args` currently returns false for some wrapper types like `LinearAlgebra.UpperTriangular`, requiring you to
use their `parent`. Triangular loops aren't yet supported.

Setting the keyword argument `warn_check_args=true` (e.g. `@turbo warn_check_args=true for ...`) in a loop or
broadcast statement will cause it to warn once if `LoopVectorization.check_args` fails and the fallback
loop is executed instead of the LoopVectorization-optimized loop.
Setting it to an integer > 0 will warn that many times, while setting it to a negative integer will warn
an unlimited amount of times. The default is `warn_check_args = 0`.
"""
macro turbo(args...)
  turbo_macro(__module__, __source__, last(args), Base.front(args)...)
end
"""
    @tturbo

Equivalent to `@turbo`, except it adds `thread=true` as the first keyword argument.
Note that later arguments take precendence.

Meant for convenience, as `@tturbo` is shorter than `@turbo thread=true`.
"""
macro tturbo(args...)
    turbo_macro(__module__, __source__, last(args), :(thread=true), Base.front(args)...)
end

function def_outer_reduct_types!(ls::LoopSet)
  for or ∈ ls.outer_reductions
    op = operations(ls)[or]
    pushpreamble!(ls, Expr(:(=), outer_reduct_init_typename(op), typeof_expr(op)))
  end
end
"""
    @_turbo

This macro mostly exists for debugging/testing purposes. It does not support many of the use cases of [`@turbo`](@ref).
It emits loops directly, rather than punting to an `@generated` function, meaning it doesn't have access to type
information when generating code or analyzing the loops, often leading to bad performance.

This macro accepts the `inline` and `unroll` keyword arguments like `@turbo`, but ignores the `check_empty` argument.
"""
macro _turbo(q)
  q = macroexpand(__module__, q)
  ls = LoopSet(q, __module__)
  set_hw!(ls)
  def_outer_reduct_types!(ls)
  esc(Expr(:block, ls.prepreamble, lower_and_split_loops(ls, -1)))
end
macro _turbo(arg, q)
  @assert q.head === :for
  q = macroexpand(__module__, q)
  inline, check_empty, u₁, u₂, v = check_macro_kwarg(arg, false, false, zero(Int8), zero(Int8), zero(Int8), 1, 0)
  ls = LoopSet(q, __module__)
  set_hw!(ls)
  def_outer_reduct_types!(ls)
  esc(Expr(:block, ls.prepreamble, lower(ls, u₁ % Int, u₂ % Int, v % Int, -1)))
end

"""
    @turbo_debug

Returns a `LoopSet` object instead of evaluating the loops. Useful for debugging and introspection.
"""
macro turbo_debug(q)
    q = macroexpand(__module__, q)
    ls = LoopSet(q, __module__)
    esc(LoopVectorization.setup_call_debug(ls))
end

# define aliases
const var"@avx" = var"@turbo"
const var"@avxt" = var"@tturbo"
const var"@avx_debug" = var"@turbo_debug"


text:
## Interface

const Projection{T <: Real} = AbstractMatrix{T}

"""
Abstract type for dimensionality reduction methods
"""
abstract type AbstractDimensionalityReduction end

"""
    vertices(R::AbstractDimensionalityReduction)

Returns vertices of largest connected component in the model `R`.
"""
vertices(R::AbstractDimensionalityReduction) = Int[]

"""
    neighbors(R::AbstractDimensionalityReduction)

Returns the number of nearest neighbors used for aproximate local subspace
"""
neighbors(R::AbstractDimensionalityReduction) = 0

"""
    fit(AbstractDimensionalityReduction, X)

Perform model fitting given the data `X`
"""
fit(::Type{AbstractDimensionalityReduction}, X::AbstractMatrix) = throw("Model fitting is not implemented")

"""
    transfrom(R::AbstractDimensionalityReduction)

Returns a reduced space representation of the data given the model `R` inform of  the projection matrix (of size ``(d, n)``), where `d` is a dimension of the reduced space and `n` in the number of the observations. Each column of the projection matrix corresponds to an observation in projected reduced space.
"""
transform(R::AbstractDimensionalityReduction) = throw("Data transformation is not implemented")

"""
    outdim(R::AbstractDimensionalityReduction)

Returns a dimension of the reduced space for the model `R`
"""
outdim(R::AbstractDimensionalityReduction) = 0

"""
    eigvals(R::AbstractDimensionalityReduction)

Returns eignevalues of the reduced space reporesentation for the model `R`
"""
eigvals(R::AbstractDimensionalityReduction) = Float64[]

# Auxiliary functions

show(io::IO, ::MIME"text/plain", R::T) where {T<:AbstractDimensionalityReduction} = summary(io, R)
function show(io::IO, R::T) where {T<:AbstractDimensionalityReduction}
    summary(io, R)
    io = IOContext(io, :limit=>true)
    println(io)
    println(io, "connected component: ")
    Base.show_vector(io, vertices(R))
    println(io)
    println(io, "eigenvalues: ")
    Base.show_vector(io, eigvals(R))
    println(io)
    println(io, "projection:")
    Base.print_matrix(io, transform(R), "[", ",","]")
end


text:
module CameraOffsetTrigger

using ..Ahorn, Maple

const placements = Ahorn.PlacementDict(
    "Camera Offset" => Ahorn.EntityPlacement(
        Maple.CameraOffsetTrigger,
        "rectangle"
    )
)

end


text:
using Test
using RadauBig
using LinearAlgebra
using Polynomials: Poly, polyval, polyder, coeffs
using PolynomialRoots: roots

include("test_exports.jl")
include("test_T_lambda.jl")
include("test_generate_butcher_table.jl")

# TODO: add tests for calculating of Eigencalculatino of T and λ


text:
# this file is for short recipes

# XXX it would be great to have a 1-1 aspect ratio here,
# but this seems like something that should be done upstream
function Makie.convert_arguments(P::Type{<:Makie.QQNorm}, x::MixedModel, args...; qqline=:R, kwargs...)
    # this addresses a bug in the current Makie release:
    # https://github.com/JuliaPlots/Makie.jl/pull/1277
    return convert_arguments(QQPlot, Distributions.Normal(0, 1), residuals(x), args...; qqline=qqline, kwargs...)
end
function Makie.convert_arguments(
    P::Type{<:Makie.QQPlot}, d::Distributions.Distribution, x::MixedModel, args...; kwargs...
)
    return convert_arguments(P, d, residuals(x), args...; kwargs...)
end

# Makie.convert_arguments(P::PointBased, x::MixedModel) = convert_arguments(P, response(x), fitted(x))
# Makie.convert_arguments(P::Type{<:Makie.Scatter}, x::MixedModel) = convert_arguments(P, fitted(x), residuals(x))
# Makie.plottype(::MixedModel) = Scatter


text:
include(joinpath(@__DIR__, "..", "src/dynamics", "hopper_2D", "visuals.jl"))
T = Float64
vis = Visualizer()
open(vis)

# qet hopper model
model = get_model("hopper_2D")
nq = model.dim.q
nu = model.dim.u
nc = model.dim.c
nb = model.dim.b
nd = nq + nc + nb
nr = nq + nu + nc + nb + nd

H = 119
h = 0.02377
κ = 1.0e-8

# Design open-loop control trajectory
q0_ref = SVector{nq,T}([0.0, 0.50, 0.0, 0.5])
q1_ref = SVector{nq,T}([0.0, 0.50, 0.0, 0.5])

α = 0.02
u_ref = []
push!(u_ref,  [[0.0,  5.00*α*model.g*(model.mb+model.ml)/2] for k=1:8]...);
push!(u_ref,  [[0.0, -1.00*α*model.g*(model.mb+model.ml)/2] for k=1:10]...);
push!(u_ref,  [[0.0, -0.55*α*model.g*(model.mb+model.ml)/2] for k=1:15]...);
push!(u_ref,  [[0.0,  0.14*α*model.g*(model.mb+model.ml)/2] for k=1:35]...);
push!(u_ref,  [[0.0,  2.105*α*model.g*(model.mb+model.ml)/2] for k=1:H-length(u_ref)]...);

# Simulate
sim = simulator(model, q0_ref, q1_ref, h, H;
    p = open_loop_policy(u_ref),
    d = no_disturbances(model),
    r! = model.res.r!,
    rz! = model.res.rz!,
    rθ! = model.res.rθ!,
    rz = model.spa.rz_sp,
    rθ = model.spa.rθ_sp,
    ip_opts = InteriorPointOptions(r_tol=1e-9, κ_tol=2κ),
    sim_opts = SimulatorOptions())
simulate!(sim)
plot(hcat(Vector.(sim.traj.u)...)')
plot(hcat(Vector.(sim.traj.q)...)')
visualize_robot!(vis, model, sim.traj)

plot_surface!(vis, model.env)
for t = 1:10
    @show norm(residual(model, sim.traj.z[t], sim.traj.θ[t], [κ]))
end

# Check ~perfect loop
@show round.(sim.traj.q[end-2] - q0_ref, digits=3)
@show round.(sim.traj.q[end-1] - q1_ref, digits=3)

# Save trajectory
traj = deepcopy(sim.traj)
gait_path = joinpath(@__DIR__, "..", "dynamics", "hopper_2D", "gaits", "gait_in_place_flip.jld2")
@save gait_path traj

# Reload trajectory
res = JLD2.jldopen(gait_path)
loaded_traj = res["traj"]

traj = get_trajectory("hopper_2D", "gait_in_place_flip", load_type=:joint_traj)
plot(hcat(Vector.(traj.q)...)')


text:
export GridRegressor

const NG = 16 # number of grids

@with_kw mutable struct GridRegressor <: BaseEstimator
    W::Array{Float32} = zeros(Float32, 0)
    α::Float32 = 1e-4
    λ::Float32 = 0
    max_iter::Int = 1000
    tol::Float32 = 1e-3
    verbose::Int = 1
    η0::Float32 = 1e-1
    power_t::Float32 = 0.25
    dim::Int = 2
end

is_classifier(::GridRegressor) = false

function gridparams(m::GridRegressor)
    grid = ["α" => [0.001, 0.01, 0.1, 1, 10]]
    gridparams(grid)
end

function fit!(m::GridRegressor, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    @unpack α, λ, max_iter, tol, verbose, η0, power_t, dim = m
    x = reshape(x, size(x, 1), :)
    F, N = size(x, 1), size(x, 2)
    dims = ntuple(i -> i <= dim ? F : NG, 2dim)
    m.W = W = zeros(Float32, dims)
    ŷ, Δ, ∂W = zero(y), zero(y), zero(W)
    fs = CartesianIndices(ntuple(i -> 1:F, dim))
    loss = prevloss = Inf32
    for t in 1:max_iter
        fill!(∂W, 0)
        predict!(ŷ, m, x)
        Δ .= ŷ .- y
        loss = mean(abs2, Δ)
        # loss > prevloss - tol && break
        prevloss = loss
        if verbose == 2 || (verbose == 1 && t % 10 == 0)
            @printf("t: %d, Norm: %.2g, Loss: %.2g\n", t, norm(W), loss)
        end
        Threads.@threads for n in 1:N
            xn = view(x, :, n)
            for f in fs
                g = gridize(xn, f)
                ∂W[f, g] += Δ[n]
            end
        end
        ∂W ./= N
        α > 0 && (∂W .+= α .* W)
        λ > 0 && (∂W .+= λ .* sign.(W))
        W .-= (η0 / t^power_t) .* ∂W
    end
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

function predict!(ŷ, m::GridRegressor, x)
    fill!(ŷ, 0)
    @unpack W, dim = m
    isempty(W) && return ŷ
    F, N = size(x)
    fs = CartesianIndices(ntuple(i -> 1:F, dim))
    Threads.@threads for n in 1:N
        xn = view(x, :, n)
        for f in fs
            g = gridize(xn, f)
            ŷ[n] += W[f, g]
        end
    end
    return ŷ
end

function predict(m::GridRegressor, x::Array{T}) where T
    @unpack W, dim = m
    x = reshape(x, size(x, 1), :)
    ŷ = zeros(T, size(x, 2))
    predict!(ŷ, m, x)
end

modelhash(m::GridRegressor) = hash(m.W)

function visualize(m::GridRegressor, columns)
    @unpack W, dim = m
    columns′ = fill("", size(m.W))
    for I in CartesianIndices(columns′)
        features = [columns[I[i]] for i in 1:dim]
        bins =  [I[i + dim ] for i in 1:dim]
        columns′[I] = join(join.(zip(features, bins), ':'), 'x')
    end
    write_feaimpt(m.W, columns′)
end

gridize(x) = (n½ = NG ÷ 2 - 1; clamp(floor(Int, n½ * x / 1.96) + n½  + 2, 1, NG))

@generated function gridize(x, f::CartesianIndex{N}) where N
    :(CartesianIndex(($([:(gridize(x[f[$i]])) for i in 1:N]...),)))
end


text:
# Tokenize
# ========

# A set of lexing parameters.
struct LexerParameters
    delim::Union{UInt8,Nothing}
    quot::UInt8
    trim::Bool
    lzstring::Bool
    skip::Int
    skipblank::Bool
    comment::String
    colnames::Union{Vector{Symbol},Nothing}
    normalizenames::Bool
    hasheader::Bool
    chunkbits::Int

    function LexerParameters(
            delim::Union{Char,Nothing}, quot::Union{Char,Nothing}, trim::Bool, lzstring::Bool,
            skip::Integer, skipblank::Bool, comment::String,
            colnames::Any, normalizenames::Bool,
            hasheader::Bool, chunkbits::Integer)
        if delim ∉ ALLOWED_DELIMITERS && delim !== nothing
            throw(ArgumentError("delimiter $(repr(delim)) is not allowed"))
        elseif quot ∉ ALLOWED_QUOTECHARS && quot !== nothing
            throw(ArgumentError("quotation character $(repr(quot)) is not allowed"))
        elseif delim == quot && !(delim === quot === nothing)
            throw(ArgumentError("delimiter and quotation character cannot be the same character"))
        elseif delim == ' ' && trim
            throw(ArgumentError("delimiting with space and space trimming are exclusive"))
        elseif quot == ' ' && trim
            throw(ArgumentError("quoting with space and space trimming are exclusive"))
        elseif skip < 0
            throw(ArgumentError("skip cannot be negative"))
        elseif occursin(r"[\r\n]", comment)
            throw(ArgumentError("comment cannot contain newline characters"))
        elseif chunkbits < 0
            throw(ArgumentError("chunkbits cannot be negative"))
        elseif chunkbits != 0 && !(MINIMUM_CHUNK_BITS ≤ chunkbits ≤ MAXIMUM_CHUNK_BITS)
            throw(ArgumentError("chunkbits must be zero or between $(MINIMUM_CHUNK_BITS) and $(MAXIMUM_CHUNK_BITS) (inclusive)"))
        end
        if colnames != nothing
            colnames = Symbol.(collect(colnames))
        end
        return new(
            delim isa Char ? UInt8(delim) : delim,
            quot isa Char ? UInt8(quot) : NO_QUOTECHAR,
            trim,
            lzstring,
            skip,
            skipblank,
            comment,
            colnames,
            normalizenames,
            hasheader,
            chunkbits,
        )
    end
end

# Token type
# ----------

# token kind
const STRING  = 0b0000
const INTEGER = 0b0001
const FLOAT   = 0b0010
const BOOL    = 0b0100
const QSTRING = 0b1000  # string with quotation marks
const MISSING = 0b1111  # missing can be any data type

const MAX_TOKEN_START = 2^36 - 1
const MAX_TOKEN_LENGTH = 2^24 - 1

struct Token
    # From most significant
    #    4bit: kind (+ missing)
    #   36bit: start position (64 GiB)
    #   24bit: length (16 MiB)
    value::UInt64

    function Token(kind::UInt8, start::Int, len::Int)
        @assert start ≤ MAX_TOKEN_START
        @assert len ≤ MAX_TOKEN_LENGTH
        return new((UInt64(kind) << 60) | (UInt64(start) << 24) | UInt64(len))
    end
end

const TOKEN_NULL = Token(0x00, 0, 0)

function kind(token::Token)
    return (token.value >> 60) % UInt8
end

function ismissing(token::Token)
    return (token.value & (UInt64(1) << 63)) != 0
end

function location(token::Token)
    x = token.value & (~UInt64(0) >> 4)
    return (x >> 24) % Int, (x & (~UInt64(0) >> 40)) % Int
end

function Base.length(token::Token)
    return (token.value & (~UInt64(0) >> 40)) % Int
end

# Line scanner
# ------------

struct ReadError <: Exception
    msg::String
end

function Base.show(io::IO, error::ReadError)
    print(io, summary(error), ": ", error.msg)
end

macro state(name, ex)
    @assert name isa Symbol
    @assert ex isa Expr && ex.head == :block
    quote
        @label $(name)
        #println($(QuoteNode(name)))
        #@show quoted
        pos += 1
        #@assert 1 ≤ pos ≤ lastindex(mem)
        @inbounds c = mem[pos]
        #@show Char(c)
        #println()
        $(ex)
        # This area is unreachable because transition must be exhaustive.
        @assert false
    end |> esc
end

macro multibytestring()
    quote
        # multibyte UTF8 character
        @inbounds if 0b110_00000 ≤ c ≤ 0b110_11111
            if pos + 1 ≤ pos_end && (mem[pos+1] >> 6) ≤ 0b10  # same as: 0b10_000000 ≤ mem[pos+1] ≤ 0b10_111111
                pos += 1
                @goto STRING
            end
        elseif 0b1110_0000 ≤ c ≤ 0b1110_1111
            if pos + 2 ≤ pos_end && (max(mem[pos+1], mem[pos+2]) >> 6) ≤ 0b10
                pos += 2
                @goto STRING
            end
        elseif 0b11110_000 ≤ c ≤ 0b11110_111
            if pos + 3 ≤ pos_end && (max(mem[pos+1], mem[pos+2], mem[pos+3]) >> 6) ≤ 0b10
                pos += 3
                @goto STRING
            end
        else
            @goto ERROR
        end
    end |> esc
end

macro follows(s, opts...)
    @assert s isa String && isascii(s)
    casesensitive = :casesensitive ∈ opts
    i = 0
    foldl(s, init = :(pos + $(sizeof(s)) ≤ pos_end)) do ex, c
        i += 1
        if casesensitive
            :($(ex) && mem[pos+$(i)] == $(UInt8(c)))
        else
            up = UInt8(uppercase(c))
            lo = UInt8(lowercase(c))
            :($(ex) && (mem[pos+$(i)] == $(up) || mem[pos+$(i)] == $(lo)))
        end
    end |> esc
end

macro begintoken()
    esc(:(start = pos))
end

macro recordtoken(kind)
    quote
        token = Token($(kind), start, pos - start)
    end |> esc
end

macro endheadertoken()
    quote
        push!(tokens, token)
        quoted = false
        qstring = false
    end |> esc
end

function scanheader(mem::Memory, params::LexerParameters)
    # Check parameters.
    delim, quot, trim = params.delim, params.quot, params.trim
    @assert delim != quot
    @assert !trim || delim != SP
    @assert !trim || quot != SP
    @assert mem[end] == CR || mem[end] == LF

    # Initialize variables.
    tokens = Token[]
    token = TOKEN_NULL
    quoted = false
    qstring = false
    pos = 0
    start = 0
    pos_end = lastindex(mem)

    @state BEGIN begin
        @begintoken
        if c == quot
            if c == NO_QUOTECHAR
                # invalid UTF-8 encoding
                @goto ERROR
            end
            if quoted
                @recordtoken STRING
                @goto QUOTE_END
            end
            quoted = true
            @goto BEGIN
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endheadertoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @goto BEGIN
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING begin
        if quoted && c == quot
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @endheadertoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                if qstring
                    @recordtoken QSTRING
                else
                    @recordtoken STRING
                end
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING_SPACE begin
        if c == SP
            @goto STRING_SPACE
        elseif c == delim
            @endheadertoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state QUOTE_END begin
        if c == delim
            @endheadertoken
            @goto BEGIN
        elseif c == quot
            qstring = true
            @goto STRING
        elseif c == SP
            if trim
                @goto QUOTE_END_SPACE
            end
            @goto ERROR
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @state QUOTE_END_SPACE begin
        if c == delim
            @endheadertoken
            @goto BEGIN
        elseif c == SP
            @goto QUOTE_END_SPACE
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @label ERROR
    throw(ReadError("invalid file header format"))

    @label CR
    if quoted
        throw(ReadError("quoted multiline string is not allowed in the header"))
    elseif pos + 1 ≤ pos_end && mem[pos + 1] == LF
        pos += 1
    end

    @label LF
    if quoted
        throw(ReadError("quoted multiline string is not allowed in the header"))
    end
    @endheadertoken

    @label END
    return pos, tokens
end

macro endtoken()
    quote
        if i ≥ ncols
            msg = "unexpected number of columns at line $(line)"
            @goto ERROR
        end
        @inbounds tokens[(i+=1)+base] = token
        quoted = false
        qstring = false
    end |> esc
end

# Scan a line in mem; mem must include one or more lines.
function scanline!(
        # output info
        tokens::Matrix{Token}, row::Int,
        # input info
        mem::Memory, pos::Int, line::Int,
        # lexer parameters
        params::LexerParameters
    )

    # Check parameters.
    delim, quot, trim, lzstring = params.delim, params.quot, params.trim, params.lzstring
    @assert delim != quot
    @assert !trim || delim != SP
    @assert !trim || quot != SP
    @assert length(mem) > 0 && (mem[end] == CR || mem[end] == LF)

    # Initialize variables.
    pos_end = lastindex(mem)
    ncols = size(tokens, 1)
    quoted = false
    qstring = false
    token = TOKEN_NULL
    blank = true  # blank line?
    msg = ""  # error message
    start = 0  # the starting position of a token
    i = 0  # the column of a token
    base = ncols * (row - 1)  # the base index of tokens in the current row

    if !isempty(params.comment)
        q = 1
        while pos + q ≤ lastindex(mem) && q ≤ sizeof(params.comment) && mem[pos+q] == codeunit(params.comment, q)
            q += 1
        end
        if q > sizeof(params.comment)
            # found a line starting with a comment sequence
            while pos + q ≤ lastindex(mem) && mem[pos+q] != CR && mem[pos+q] != LF
                q += 1
            end
            if pos + q + 1 ≤ lastindex(mem) && mem[pos+q] == CR && mem[pos+q+1] == LF
                pos += q + 1
            else
                pos += q
            end
            return pos, 0, true
        end
    end

    @state BEGIN begin
        blank &= (params.trim && c == SP) || c == CR || c == LF
        @begintoken
        if c == quot
            if c == NO_QUOTECHAR
                # invalid UTF-8 encoding
                @goto ERROR
            end
            if quoted
                @recordtoken MISSING
                @goto QUOTE_END
            end
            quoted = true
            @goto BEGIN
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken MISSING
            @endtoken
            @goto BEGIN
        elseif c == UInt8('-') || c == UInt8('+')
            @goto SIGN
        elseif c == UInt8('0')
            if lzstring
                @goto ZERO
            end
            @goto INTEGER
        elseif UInt8('1') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == SP
            if trim && !quoted
                @goto BEGIN
            end
            @goto STRING
        elseif c == UInt8('.')
            @goto DOT
        elseif c == UInt8('I') || c == UInt8('i')
            # Infinity
            if @follows("NFINITY")
                pos += 7
                @goto SPECIAL_FLOAT
            elseif @follows("NF")
                pos += 2
                @goto SPECIAL_FLOAT
            end
            @goto STRING
        elseif (c == UInt8('N') || c == UInt8('n'))
            if @follows("AN")  # case-insensitive
                # NaN
                pos += 2  # for 'A' and 'N'
                @goto SPECIAL_FLOAT
            elseif c == UInt8('N') && @follows("A", casesensitive)
                # NA
                pos += 1  # for 'A'
                @goto NA
            end
            @goto STRING
        elseif c == UInt8('F') || c == UInt8('f')
            if @follows("alse")
                pos += 4
            end
            @goto BOOL
        elseif c == UInt8('T') || c == UInt8('t')
            if @follows("rue")
                pos += 3
            end
            @goto BOOL
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            if i + 1 == ncols
                @recordtoken MISSING
            end
            @goto LF
        elseif c == CR
            if i + 1 == ncols
                @recordtoken MISSING
            end
            @goto CR
        else
            @multibytestring
        end
    end

    @state SIGN begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == UInt8('.')
            @goto DOT
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif (c == UInt8('N') || c == UInt8('n')) && @follows("AN")
            pos += 2
            @goto SPECIAL_FLOAT
        elseif c == UInt8('I') || c == UInt8('i')
            # Infinity
            if @follows("NFINITY")
                pos += 7
                @goto SPECIAL_FLOAT
            elseif @follows("NF")
                pos += 2
                @goto SPECIAL_FLOAT
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state ZERO begin
        if quoted && c == quot
            @recordtoken INTEGER|FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken INTEGER|FLOAT
            @endtoken
            @goto BEGIN
        elseif c == UInt8('.')
            @goto POINT_FLOAT
        elseif c == UInt8('e') || c == UInt8('E')
            @goto EXPONENT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken INTEGER|FLOAT
                @goto INTEGER_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken INTEGER|FLOAT
            @goto LF
        elseif c == CR
            @recordtoken INTEGER|FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state INTEGER begin
        if quoted && c == quot
            @recordtoken INTEGER|FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken INTEGER|FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto INTEGER
        elseif c == UInt8('.')
            @goto POINT_FLOAT
        elseif c == SP
            if trim && !quoted
                @recordtoken INTEGER|FLOAT
                @goto INTEGER_SPACE
            end
            @goto STRING
        elseif c == UInt8('e') || c == UInt8('E')
            @goto EXPONENT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken INTEGER|FLOAT
            @goto LF
        elseif c == CR
            @recordtoken INTEGER|FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state INTEGER_SPACE begin
        if c == SP
            @goto INTEGER_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state DOT begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto POINT_FLOAT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state POINT_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto POINT_FLOAT
        elseif c == UInt8('e') || c == UInt8('E')
            @goto EXPONENT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT
            @goto LF
        elseif c == CR
            @recordtoken FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif c == UInt8('-') || c == UInt8('+')
            @goto EXPONENT_SIGN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT_SIGN begin
        if quoted && c == quot
            @recordtoken STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken STRING
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                @recordtoken STRING
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            @recordtoken STRING
            @goto LF
        elseif c == CR
            @recordtoken STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state EXPONENT_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT
            @endtoken
            @goto BEGIN
        elseif UInt8('0') ≤ c ≤ UInt8('9')
            @goto EXPONENT_FLOAT
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT
            @goto LF
        elseif c == CR
            @recordtoken FLOAT
            @goto CR
        else
            @multibytestring
        end
    end

    @state SPECIAL_FLOAT begin
        if quoted && c == quot
            @recordtoken FLOAT|STRING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken FLOAT|STRING
            @endtoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @recordtoken FLOAT|STRING
                @goto FLOAT_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken FLOAT|STRING
            @goto LF
        elseif c == CR
            @recordtoken FLOAT|STRING
            @goto CR
        else
            @multibytestring
        end
    end

    @state FLOAT_SPACE begin
        if c == SP
            @goto FLOAT_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state BOOL begin
        if quoted && c == quot
            @recordtoken BOOL
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken BOOL
            @endtoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @recordtoken BOOL
                @goto BOOL_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken BOOL
            @goto LF
        elseif c == CR
            @recordtoken BOOL
            @goto CR
        else
            @multibytestring
        end
    end

    @state BOOL_SPACE begin
        if c == SP
            @goto BOOL_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken BOOL
            @goto LF
        elseif c == CR
            @recordtoken BOOL
            @goto CR
        else
            @multibytestring
        end
    end

    @state NA begin
        if quoted && c == quot
            @recordtoken MISSING
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            @recordtoken MISSING
            @endtoken
            @goto BEGIN
        elseif c == SP
            if trim && !quoted
                @recordtoken MISSING
                @goto NA_SPACE
            end
            @goto STRING
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken MISSING
            @goto LF
        elseif c == CR
            @recordtoken MISSING
            @goto CR
        else
            @multibytestring
        end
    end

    @state NA_SPACE begin
        if c == SP
            @goto NA_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @recordtoken MISSING
            @goto LF
        elseif c == CR
            @recordtoken MISSING
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING begin
        if quoted && c == quot
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto QUOTE_END
        elseif c == delim
            if quoted
                @goto STRING
            end
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == SP
            if trim && !quoted
                if qstring
                    @recordtoken QSTRING
                else
                    @recordtoken STRING
                end
                @goto STRING_SPACE
            end
            @goto STRING
        elseif c == LF
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto LF
        elseif c == CR
            if qstring
                @recordtoken QSTRING
            else
                @recordtoken STRING
            end
            @goto CR
        else
            @multibytestring
        end
    end

    @state STRING_SPACE begin
        if c == SP
            @goto STRING_SPACE
        elseif c == delim
            @endtoken
            @goto BEGIN
        elseif UInt8('!') ≤ c ≤ UInt8('~')
            @goto STRING
        elseif c == LF
            @goto LF
        elseif c == CR
            @goto CR
        else
            @multibytestring
        end
    end

    @state QUOTE_END begin
        if c == delim
            @endtoken
            @goto BEGIN
        elseif c == quot  # e.g. xxx," foo ""bar""",xxx
            qstring = true
            @goto STRING
        elseif c == SP
            if trim
                @goto QUOTE_END_SPACE
            end
            @goto ERROR
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @state QUOTE_END_SPACE begin
        if c == delim
            @endtoken
            @goto BEGIN
        elseif c == SP
            @goto QUOTE_END_SPACE
        elseif c == LF
            quoted = false
            @goto LF
        elseif c == CR
            quoted = false
            @goto CR
        else
            @goto ERROR
        end
    end

    @label ERROR
    if isempty(msg)
        # default error message
        msg = "invalid file format at line $(line), column $(i + 1) "
        if c ≤ 0x7f  # ASCII
            msg = string(msg, "(found $(repr(Char(c))))")
        else
            msg = string(msg, "(found $(repr(c)))")
        end
    end
    throw(ReadError(msg))

    @label CR  # carriage return
    if pos + 1 ≤ pos_end && mem[pos + 1] == LF
        pos += 1
        # fall through
    elseif quoted
        # maybe a quoted multiline string
        if pos == pos_end
            # need more data
            return 0, 0, false
        end
        @goto STRING
    else
        @endtoken
        @goto END
    end

    @label LF  # line feed
    if quoted
        # maybe a quoted multiline string
        if pos == pos_end
            # need more data
            return 0, 0, false
        end
        @goto STRING
    else
        @endtoken
        # fall through
    end

    @label END
    return pos, i, params.skipblank && blank
end


text:
runall = false;

@testset "Testing Substitution Models, P and Q matrices" begin

m1 = BinaryTraitSubstitutionModel(1.0, 2.0);
show(m1)
m1 = BinaryTraitSubstitutionModel(1.0,2.0, SVector("carnivory", "non-carnivory"));
@test nStates(m1)==2
show(m1)
@test_throws ErrorException PhyloNetworks.BinaryTraitSubstitutionModel(-1.0,2.0)
m2 = EqualRatesSubstitutionModel(4, 3.0);
@test nStates(m2)==4
m2 = EqualRatesSubstitutionModel(4, 3.0, ["S1","S2","S3","S4"]);
show(m2)
m3 = TwoBinaryTraitSubstitutionModel([2.0,1.2,1.1,2.2,1.0,3.1,2.0,1.1],
["carnivory", "noncarnivory", "wet", "dry"]);
show(m3)

@test Q(m1) == SMatrix{2,2}(-1.0, 2.0, 1.0, -2.0)
@test Q(m2) == SMatrix{4,4}(-9.0, 3, 3, 3, 3, -9, 3, 3, 3, 3, -9, 3, 3, 3, 3, -9)
@test Q(m3) ≈ SMatrix{4,4}(-3.0, 3.1, 1.2, 0.0, 1.0, -4.2, 0.0, 2.2, 2.0, 0.0, -3.2, 1.1, 0.0, 1.1, 2.0, -3.3) atol=1e-4
@test P(m1, 0.5) ≈ SMatrix{2,2}(0.7410433867161432,0.5179132265677134,0.2589566132838567,0.4820867734322865) atol=1e-15
@test P(m2, 0.1) ≈ SMatrix{4,4}(0.4758956589341516,0.17470144702194956,0.17470144702194962,0.17470144702194945,0.17470144702194956,0.47589565893415153,0.17470144702194967,0.17470144702194945,0.17470144702194962,0.17470144702194967,0.4758956589341516,0.17470144702194945,0.17470144702194945,0.17470144702194945,0.17470144702194945,0.4758956589341518) atol=1e-15
@test P(m3, 0.5) ≈ SMatrix{4,4}(0.39839916380463375,0.36847565707393248,0.23055614536582461,0.22576141081414305,0.1545971371424259,0.25768553619230444,0.14816051303688715,0.24300762748855972,0.29194735222005136,0.20198250750421617,0.35349416558860958,0.20267178083856716,0.15505634683288913,0.17185629922954704,0.26778917600867863,0.32855918085873009) atol=1e-15
@test P(m1, [0.02,0.01]) ≈ StaticArrays.SArray{Tuple{2,2},Float64,2,4}[[0.980588 0.0194118; 0.0388236 0.961176], [0.990149 0.00985149; 0.019703 0.980297]] atol=1e-6
@test P(m2, [0.02,0.01]) ≈ Array{Float64,2}[[0.839971 0.053343 0.053343 0.053343; 0.053343 0.839971 0.053343 0.053343; 0.053343 0.053343 0.839971 0.053343; 0.053343 0.053343 0.053343 0.839971], [0.91519 0.0282699 0.0282699 0.0282699; 0.0282699 0.91519 0.0282699 0.0282699; 0.0282699 0.0282699 0.91519 0.0282699; 0.0282699 0.0282699 0.0282699 0.91519]] atol=1e-6

end

@testset "Testing random discrete trait simulation" begin

m1 = BinaryTraitSubstitutionModel(1.0,2.0, SVector("carnivory", "non-carnivory"));
m2 = EqualRatesSubstitutionModel(4, 3.0, ["S1","S2","S3","S4"]);

info("\ton a single branch")
srand(12345);
@test randomTrait(m1, 0.2, [1,2,1,2,2]) == [1,2,1,1,2]
srand(12345);
@test randomTrait(m2, 0.05, [1,3,4,2,1]) == [1,3,4,2,1]

info("\ton a network")
net = readTopology("(A:1.0,(B:1.0,(C:1.0,D:1.0):1.0):1.0);")
srand(12345);
a,b = randomTrait(m1, net)
@test a == [1 2 1 1 1 1 2]
@test b == ["-2", "-3", "-4", "D", "C", "B", "A"]
if runall
    for e in net.edge e.length = 10.0; end
    @time a,b = randomTrait(m1, net; ntraits=100000) # ~ 0.014 seconds
    mean(a[:,1]) # expect 1.5 at root
    mean(a[:,2]) # expect 1.333 at other nodes
    @time a,b = randomTrait(m2, net; ntraits=100000) # ~ 0.02 seconds
    length([x for x in a[:,1] if x==4])/length(a[:,1]) # expect 0.25
    length([x for x in a[:,2] if x==4])/length(a[:,2])
    length([x for x in a[:,3] if x==4])/length(a[:,3])
    length([x for x in a[:,4] if x==4])/length(a[:,4])
    length([x for x in a[:,5] if x==4])/length(a[:,5])
    length([x for x in a[:,6] if x==4])/length(a[:,6])
    length([x for x in a[:,7] if x==4])/length(a[:,7]) # expect 0.25
end

net2 = readTopology("(((A:4.0,(B:1.0)#H1:1.1::0.9):0.5,(C:0.6,#H1:1.0):1.0):3.0,D:5.0);")
srand(12345);
a,b = randomTrait(m1, net2; keepInternal=false)
@test a == [1  1  1  2]
@test b == ["D", "C", "B", "A"]
srand(12345);
a,b = randomTrait(m1, net2; keepInternal=true)
@test a == [1  2  1  1  1  1  1  1  1]
@test b == ["-2", "D", "-3", "-6", "C", "-4", "#H1", "B", "A"]
if runall
    for e in net2.edge
        if e.hybrid
            e.length = 0.0
        end
    end
    a,b = randomTrait(m1, net2; ntraits=100000)
    # plot(net2, showNodeNumber=true) shows: H1 listed 7th, parents listed 4th and 6th
    c = map( != , a[:, 4],a[:, 6] ); # traits when parents have different traits
    n1 = sum(map( ==, a[c,7],a[c,6] )) # 39644 traits: hybrid ≠ major parent
    n2 = sum(map( ==, a[c,7],a[c,4] )) #  4401 traits: hybrid ≠ minor parent
    n1/sum(c) # expected 0.9
    n2/sum(c) # expected 0.1
    for e in net2.edge
        e.length = 0.0
    end
    net2.edge[4].length = 10.0
    a,b = randomTrait(m1, net2; ntraits=100000);
    a[:, 1] == a[:, 2]  # true: root = leaf D, as expected
    a[:, 1] == a[:, 5]  # true: root = leaf C
    mean(a[:, 6]) # expected 1.3333
    a[:, 6] == a[:, 9] # true: major hybrid parent node = leaf A
end

end


text:
__precompile__(true)

module EconUtils
using Compat

using Missings, CategoricalArrays, DataFrames, StatsBase, StatsModels

include.((
	"utils.jl",
	"firstdifference.jl",
	"between.jl",
	"dropsingletons.jl",
	"within.jl",
	"EconFormula.jl",
	"modelframe.jl",
	"tsls.jl"
	))

end


text:
include("utilities.jl")

using SEAL
using Printf

"""
    example_batch_encoder()

Example to demonstrate encoding and decoding vectors of integers with the BFV scheme.
This function is based on the file `native/examples/2_encoders.cpp` of the original
SEAL library and should yield the exact same output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/2_encoders.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/2_encoders.cpp)

See also: [`example_ckks_encoder`](@ref), [`example_encoders`](@ref)
"""
function example_batch_encoder()
  print_example_banner("Example: Encoders / Batch Encoder")

  parms = EncryptionParameters(SchemeType.bfv)
  poly_modulus_degree = 8192
  set_poly_modulus_degree!(parms, poly_modulus_degree)
  set_coeff_modulus!(parms, coeff_modulus_bfv_default(poly_modulus_degree))

  set_plain_modulus!(parms, plain_modulus_batching(poly_modulus_degree, 20))
  context = SEALContext(parms)
  print_parameters(context)
  println()

  context_data = first_context_data(context)
  epq = qualifiers(context_data)
  println("Batching enabled: ", using_batching(epq))

  keygen = KeyGenerator(context)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  secret_key_ = secret_key(keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)
  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  batch_encoder = BatchEncoder(context)
  slot_count_ = slot_count(batch_encoder)
  row_size = div(slot_count_, 2)
  println("Plaintext matrix row size: ", row_size)

  pod_matrix = zeros(UInt64, slot_count_)
  pod_matrix[1] = 0
  pod_matrix[2] = 1
  pod_matrix[3] = 2
  pod_matrix[4] = 3
  pod_matrix[row_size + 1] = 4
  pod_matrix[row_size + 2] = 5
  pod_matrix[row_size + 3] = 6
  pod_matrix[row_size + 4] = 7

  println("Input plaintext matrix:")
  print_matrix(pod_matrix, row_size)

  plain_matrix = Plaintext()
  print_line(@__LINE__)
  println("Encode plaintext matrix:")
  encode!(plain_matrix, pod_matrix, batch_encoder)

  println("    + Decode plaintext matrix ...... Correct.")
  pod_result = similar(pod_matrix)
  decode!(pod_result, plain_matrix, batch_encoder)
  print_matrix(pod_result, row_size)

  encrypted_matrix = Ciphertext()
  print_line(@__LINE__)
  println("Encrypt plain_matrix to encrypted_matrix.")
  encrypt!(encrypted_matrix, plain_matrix, encryptor)
  println("    + Noise budget in encrypted_matrix: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")

  pod_matrix2 = ones(UInt64, slot_count_)
  pod_matrix2[2:2:slot_count_] .= 2
  plain_matrix2 = Plaintext()
  encode!(plain_matrix2, pod_matrix2, batch_encoder)
  println()
  println("Second input plaintext matrix:")
  print_matrix(pod_matrix2, row_size)

  print_line(@__LINE__)
  println("Sum, square, and relinearize.")
  add_plain_inplace!(encrypted_matrix, plain_matrix2, evaluator)
  square_inplace!(encrypted_matrix, evaluator)
  relinearize_inplace!(encrypted_matrix, relin_keys_, evaluator)
  println("    + Noise budget in result: ",
          invariant_noise_budget(encrypted_matrix, decryptor),
          " bits")

  plain_result = Plaintext()
  print_line(@__LINE__)
  println("Decrypt and decode result.")
  decrypt!(plain_result, encrypted_matrix, decryptor)
  decode!(pod_result, plain_result, batch_encoder)
  println("    + Result plaintext matrix ...... Correct.")
  print_matrix(pod_result, row_size)
end

"""
    example_ckks_encoder()

Example to demonstrate encoding and decoding vectors of floating point values with the CKKS scheme.
This function is based on the file `native/examples/2_encoders.cpp` of the original
SEAL library and should yield the exact same output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/2_encoders.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/2_encoders.cpp)

See also: [`example_batch_encoder`](@ref), [`example_encoders`](@ref)
"""
function example_ckks_encoder()
  print_example_banner("Example: Encoders / CKKS Encoder")

  parms = EncryptionParameters(SchemeType.ckks)

  poly_modulus_degree = 8192
  set_poly_modulus_degree!(parms, poly_modulus_degree)
  set_coeff_modulus!(parms, coeff_modulus_create(poly_modulus_degree, [40, 40, 40, 40, 40]))

  context = SEALContext(parms)
  print_parameters(context)
  println()

  keygen = KeyGenerator(context)
  public_key_ = PublicKey()
  create_public_key!(public_key_, keygen)
  secret_key_ = secret_key(keygen)
  relin_keys_ = RelinKeys()
  create_relin_keys!(relin_keys_, keygen)

  encryptor = Encryptor(context, public_key_)
  evaluator = Evaluator(context)
  decryptor = Decryptor(context, secret_key_)

  encoder = CKKSEncoder(context)

  slot_count_ = slot_count(encoder)
  println("Number of slots: ", slot_count_)

  input = Float64[0.0, 1.1, 2.2, 3.3]
  println("Input vector: ")
  print_vector(input)

  plain = Plaintext()
  initial_scale = 2.0^30
  print_line(@__LINE__)
  println("Encode input vectors.")
  encode!(plain, input, initial_scale, encoder)

  output = Vector{Float64}(undef, slot_count_)
  println("    + Decode input vector ...... Correct.")
  decode!(output, plain, encoder)
  print_vector(output)

  encrypted = Ciphertext()
  print_line(@__LINE__)
  println("Encrypt input vector, square, and relinearize.")
  encrypt!(encrypted, plain, encryptor)

  square_inplace!(encrypted, evaluator)
  relinearize_inplace!(encrypted, relin_keys_, evaluator)

  println("    + Scale in squared input: ", scale(encrypted),
          " (", log2(scale(encrypted)), " bits)")

  print_line(@__LINE__)
  println("Decrypt and decode.")
  decrypt!(plain, encrypted, decryptor)
  decode!(output, plain, encoder)
  println("    + Result vector ...... Correct.")
  print_vector(output)

end

"""
    example_encoders()

Demonstrate how vectors of integers (BFV) or floating point values (CKKS) can be encoded/decoded
with the BFV and CKKS schemes.
This function is based on the file `native/examples/2_encoders.cpp` of the original
SEAL library and should yield the exact same output.

* [SEAL](https://github.com/microsoft/SEAL)
* [native/examples/2_encoders.cpp](https://github.com/microsoft/SEAL/blob/master/native/examples/2_encoders.cpp)

See also: [`example_batch_encoder`](@ref), [`example_ckks_encoder`](@ref)
"""
function example_encoders()
  print_example_banner("Example: Encoders")

  example_batch_encoder()
  example_ckks_encoder()
end


text:
module SaveVTK

using StaticArrays
using WriteVTK
using Random
using Parameters

export write_vtp

"""
    Struct to hold simulation data
"""
@with_kw mutable struct SimData
    Points::Any#::Array{Float32,1} = Array{Float32}(undef, 0)
    Idp::Any#::Array{Int32,1} = Array{Int32,1}()
    Vel::Any#::Array{Float32,1} = Array{Float32}(undef, 0)
    Rhop::Any#::Array{Float32,1} = Array{Float32,1}()
end

"""
    write_vtp(filename,points,attribute)

Saves a polydata file (.vtp) given a filename, a 2xN point array and and 1xN attribute array.
"""
function write_vtp(filename::String, sim_arr::SimData)

    x = sim_arr.Points[1:3:end]
    y = sim_arr.Points[2:3:end]
    z = sim_arr.Points[3:3:end]

    points = hcat(x, y, z)'

    polys = empty(MeshCell{WriteVTK.PolyData.Polys,UnitRange{Int64}}[])
    verts = empty(MeshCell{WriteVTK.PolyData.Verts,UnitRange{Int64}}[])

    # Note: the order of verts, lines, polys and strips is not important.
    # One doesn't even need to pass all of them.
    all_cells = (verts, polys)

    #allocations here, not my library
    vtk = vtk_grid(filename, points, all_cells..., compress = true, append = false)

    if !isempty(sim_arr.Idp)
        vtk["Idp"] = sim_arr.Idp
    end
    if !isempty(sim_arr.Vel)
        vx = sim_arr.Vel[1:3:end]
        vy = sim_arr.Vel[2:3:end]
        vz = sim_arr.Vel[3:3:end]
        velocity = hcat(vx, vy, vz)'
        vtk["Velocity"] = velocity
    end
    if !isempty(sim_arr.Rhop)
        vtk["Rhop"] = sim_arr.Rhop
    end

    vtk_save(vtk)
end

end


text:
LarA = Lar.Arrangement

using IntervalTrees
using SparseArrays
using NearestNeighbors
using DataStructures

#-------------------------------------------------------------------------------
#   MODEL BOUNDING BOXES
#-------------------------------------------------------------------------------
"""
	getModelBoundingBoxes(
		model::Lar.Model,
		deg::Int
	)::Array{Array{Float64,2},1}

Evaluates the bounding boxes related to each `deg`-cell of the `model`

The function builds a vector of matrices containing the two distincive vertices
of the axis-coherent hyper-cube that contains each `deg`-cell.

See also: [`Lar.Arrangement.spatialIndex`](@ref)
It uses:
---
# Examples
```jldoctest
julia> Lar.Model(hcat([
    [1.0, 0.0], [0.0, 1.0], [0.0, 0.5], [0.5, 1.0], [1.0, 1.0]
]...));
julia> model.T[1] = LAR.coboundary_0([[1, 2], [2, 5], [3, 4], [4, 5]]);
```
"""

function getModelBoundingBoxes(
		model::Lar.Model,
		deg::Int
	)::Array{Array{Float64,2},1}

	function findBBox(pts::Lar.Points)::Tuple{Array{Float64,2},Array{Float64,2}}
		minimum = mapslices(x->min(x...), pts, dims=2)
		maximum = mapslices(x->max(x...), pts, dims=2)
		return minimum, maximum
	end

	cellpts = [ hcat(Lar.getModelCellVertices(model, deg, c)...)
				for c = 1 : size(model, deg, 1)
			  ]
	return [hcat(findBBox(c)...) for c in cellpts]
end

#-------------------------------------------------------------------------------
#   SPATIAL_INDEXING METHODS
#	 buildBBIntervalTree
#	 extractBoxCovering
#-------------------------------------------------------------------------------

"""
	buildBBIntervalTree(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::IntervalTrees.IntervalMap{Float64, Array}

Build the Interval Tree of a set of bounding boxes along the `dim` coordinate.
"""
function buildBBIntervalTree(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::IntervalTrees.IntervalMap{Float64, Array}

	# Build the dictionary related to the `dim`-coordinate
	boxDict = OrderedDict{Array{Float64,1}, Array{Int,1}}()
	for (h, box) in enumerate(bboxes)
		key = box[dim, :]
		if haskey(boxDict, key) == false
			boxDict[key] = [h]
		else
			push!(boxDict[key], h)
		end
	end

	# Build the related Interval Tree
	iT = IntervalTrees.IntervalMap{Float64, Array}()
	for (key, boxSet) in boxDict
		iT[tuple(key...)] = boxSet
	end

	return iT
end

"""
	extractBoxCovering(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::Array{Array{Int,1},1}

Evaluates possible `bboxes` intersections over the `dim`-th coordinate

The function evaluates the pairwise intersection between the `dim` coordinate
of each `bboxes` couple. The operation is performed via Interval Tree.
"""
function extractBoxCovering(
		bboxes::Array{Array{Float64,2},1},
		dim::Int
	)::Array{Array{Int,1},1}

	# Build the related Interval Tree
	tree = LarA.buildBBIntervalTree(bboxes, dim)

	# Build the covers
	covers = [[] for k = 1 : length(bboxes)]
	for (i, bbox) in enumerate(bboxes)
		extent = bboxes[i][dim, :]
		iterator = IntervalTrees.intersect(tree, tuple(extent...))
		for x in iterator
			append!(covers[i], x.value)
		end
	end

	return covers
end

#-------------------------------------------------------------------------------
#   SPATIAL_INDEXING
#-------------------------------------------------------------------------------

"""
	spaceIndex(model::Lar.Model)::Array{Array{Int,1},1}

Generation of *space indexes* for all ``(d-1)``-dim cell members of `model`.

*Spatial index* made by ``d`` *interval-trees* on
bounding boxes of ``sigma in S_{d−1}``. Spatial queries solved by
intersection of ``d`` queries on IntervalTrees generated by
bounding-boxes of geometric objects (LAR cells).

The return value is an array of arrays of `int`s, indexing cells whose
containment boxes are intersecting the containment box of the first cell.
According to Hoffmann, Hopcroft, and Karasick (1989) the worst-case complexity of
Boolean ops on such complexes equates the total sum of such numbers.

# Examples 2D

```
julia> V = hcat([[0.,0],[1,0],[1,1],[0,1],[2,1]]...);

julia> EV = [[1,2],[2,3],[3,4],[4,1],[1,5]];

julia> Sigma = Lar.spaceindex((V,EV))
5-element Array{Array{Int64,1},1}:
 [4, 5, 2]
 [1, 3, 5]
 [4, 5, 2]
 [1, 3, 5]
 [4, 1, 3, 2]
```

From `model2d` value, available in `?input_collection` docstring:

```julia
julia> Sigma =  spaceindex(model2d);
```

# Example 3D

```julia
model = model3d
Sigma =  spaceindex(model3d);
Sigma
```
"""
function spaceIndex(model::Lar.Model)::Array{Array{Int,1},1}

	bboxes = LarA.getModelBoundingBoxes(model, 1)
	# Build Box Coverings for each dimension and intersecting them
	covers = LarA.extractBoxCovering(bboxes, 1)
	for d = 2 : length(model)
		dcovers = LarA.extractBoxCovering(bboxes, d)
		covers  = [intersect(pair...) for pair in zip(covers, dcovers)]
	end

	# Remove each cell from its cover
	for k = 1 : length(covers)
		covers[k] = setdiff(covers[k],[k])
	end

	return covers
end


text:
# empty types for correlations

abstract type CorrelationType end

struct KendallCorrelation <: CorrelationType end
struct SpearmanCorrelation <: CorrelationType end

"""
  Debye(x, k::Int)

Returns Debye function Dₖ(x) value
"""
Debye(x, k::Int=1) = k/x^k*(quadgk(i -> i^k/(exp(i)-1), 0, x)[1])

# kendall's τ to copulas parameters

"""
  τ2θ(τ, copula::String)

Returns Float, a single parameter of Archimedean copula, given the Kenalss τ correlation
"""
function τ2θ(τ, copula)
  if copula == "gumbel"
    return 1/(1-τ)
  elseif copula == "clayton"
    return 2*τ/(1-τ)
  elseif copula == "frank"
    return frankτ2θ(τ)
  elseif copula == "amh"
    return AMHτ2θ(τ)
  else
  throw(AssertionError("$(copula) not supported"))
  end
end

"""
  frankτ2θ(τ)

Returns a Frank copula θ parameter, givem Kendall's τ
"""
function frankτ2θ(τ)
  f(θ) = 1+4*(Debye(θ)-1)/θ - τ
  if τ > 0.
    return fzero(f, τ, 100.)
  elseif τ < 0.
    return fzero(f, -100., τ)
  end
  throw(DomainError("τ = 0 not supported"))
end

"""
  AMHτ2θ(τ)

Returns Ali-Mikhail-Haq copula θ parameter, givem Kendall's τ
"""
function AMHτ2θ(τ)
  f(θ) = (1 - 2*(*(1-θ)*(1-θ)log(1-θ) + θ)/(3*θ^2))-τ
  if -0.01 < τ < 0.01
    return 0.0000000000001
  elseif τ >= 0.28
    return 0.9999
  elseif 0. < τ < 0.28
    return fzero(f, 0.005, 0.995)
  elseif -2/11 < τ < 0.
    return fzero(f, -0.995, -0.005)
  end
  -0.9999
end

# Spearman ρ to copulas parameter

"""
 ρ2θ(ρ, copula)

 Returns a Float, an Archimedean copula parameter given expected Spermann correlation
 ρ and a copula.

"""
function ρ2θ(ρ, copula)
  if copula == "gumbel"
    return gumbelρ2θ(ρ)
  elseif copula == "clayton"
    return claytonρ2θ(ρ)
  elseif copula == "frank"
    return frankρ2θ(ρ)
  elseif copula == "amh"
    return AMHρ2θ(ρ)
  end
  throw(AssertionError("$(copula) not supported"))
end

### Clayton and Gumbel copulas

function Ccl(x, θ)
  if θ > 0
    return (x[1]^(-θ)+x[2]^(-θ)-1)^(-1/θ)
  else
    return (maximum([x[1]^(-θ)+x[2]^(-θ)-1, 0]))^(-1/θ)
  end
end

Cg(x, θ) = exp(-((-log(x[1]))^θ+(-log(x[2]))^θ)^(1/θ))

dilog(x) = quadgk(t -> log(t)/(1-t), 1, x)[1]

# converts parameter to correlations and vice versa

gumbelθ2ρ(θ) = 12*hcubature(x-> Cg(x, θ), [0,0],[1,1])[1]-3

function gumbelρ2θ(ρ)
  if ρ < 0.01
    return 1.
  else
    f(θ) = gumbelθ2ρ(θ)-ρ
    return fzero(f, 1.000001, 100.)
  end
end

 claytonθ2ρ(θ) = 12*hcubature(x-> Ccl(x, θ), [0,0],[1,1])[1]-3

function claytonρ2θ(ρ)
  f(θ) = claytonθ2ρ(θ)-ρ
  if ρ > .038
    return fzero(f, .01, 10.)
  elseif ρ < -.038
    return fzero(f, -1., -0.01)
  else
    return 0.052
  end
end

function AMHρ2θ(ρ)
  f(p) = (12*(1+p)*dilog(1-p)-24*(1-p)*log(1-p))/p^2-3*(p+12)/p-ρ
  if -0.01 < ρ  < 0.01
    return 0.00001
  elseif ρ <= -0.272
    return -0.999999
  elseif 0. < ρ < 0.47
    return fzero(f, 0.005, 0.995)
  elseif -0.272 < ρ < 0.
    return fzero(f, -0.995, -0.005)
  end
  0.999999
end

function frankρ2θ(ρ)
  f(θ) = 1+12*(Debye(θ, 2)- Debye(θ))/θ-ρ
  if ρ > 0.00001
    return fzero(f, ρ, 100.)
  elseif ρ < -0.00001
    return fzero(f, -100., ρ)
  end
  throw(DomainError("ρ = 0 not supported"))
end
Real