text:
import Distributions, StatsBase, Random
@everywhere using Distributions, StatsBase, Random
@everywhere function datagen(n::Int, Ctrue::Int, N::Int, J::Int, n_ij::Vector{Int64}; b::Int=round(Int,rand()*1000), adjustmean::Bool=false, showpower::Bool=false, Xstd::Real=1.0)

    Random.seed!(200*b)
    groupindex = StatsBase.inverse_rle(1:n, n_ij)
    X = rand(Distributions.Normal(), N, 2)
    betas_true=ones(J)
    if adjustmean
        mu_all = -log(1/0.81 - 1)
        if Ctrue == 1
            mu_true = [mu_all]
            wi_true = [1.0]
            sigmas_true = [0.3493]
        elseif Ctrue == 2
            mu_true = [mu_all - 1.0, mu_all + 1.2]
            wi_true = [.3, .7]
            sigmas_true = [.6, .5]
        elseif Ctrue == 3
            mu_true = [mu_all - 1.0, mu_all, mu_all + 2.]
            wi_true = [.2, .5, .3]
            sigmas_true = [.6, .4, .6]
        end
    elseif showpower
        mu_all = log(1/0.779 - 1)
        if Ctrue == 2
            mu_true = [mu_all - 1.0, mu_all + 0.8]
            wi_true =  [.6, .4]
            sigmas_true = [1.2, .8]
        elseif Ctrue == 3
            mu_true = [mu_all - 2.0, mu_all + 1.0, mu_all + 3.5]
            wi_true = [.3, .4, .3]
            sigmas_true = [1.2, .8, .9]
        end
    else
        mu_all = log(1/0.779 - 1)
        if Ctrue == 1
            mu_true = [mu_all]
            wi_true = [1.0]
            sigmas_true = [0.5]
        elseif Ctrue == 2
            mu_true = [mu_all - 2.0, mu_all + 2.0]
            wi_true =  [.5, .5]
            sigmas_true = [1.2, .8]
        elseif Ctrue == 3
            mu_true = [mu_all - 4.0, mu_all + 1.0, mu_all + 4.0]
            wi_true = [.3, .4, .3]
            sigmas_true = [1.2, .8, .9]
        end
    end

    #Randomly data generation based on the setting on datagen.jl
    Random.seed!(b * 100)
    m = Distributions.MixtureModel(map((u, v) -> Distributions.Normal(u, v), mu_true, sigmas_true), wi_true)
    gamma_true = rand(m, n)

    prob = exp.(gamma_true[groupindex] .+ X*betas_true)
    prob= prob ./ (1 .+ prob)
    Y = Bool[Distributions.rand(Distributions.Binomial(1, prob[i])) == 1 for i in 1:N];
    X = X .- mean(X, dims=1);
    return(X, Y, groupindex, wi_true, mu_true, sigmas_true, betas_true, m, gamma_true)
end


text:
function Base.sortperm(ds::Dataset, cols; alg = HeapSortAlg(), rev = false, mapformats::Bool = true, stable = true)
    isempty(ds) && return []
    _check_consistency(ds)
    colsidx = index(ds)[cols]
    _sortperm(ds, cols, rev, a = alg, mapformats = mapformats, stable = stable)[2]
end

function Base.sort!(ds::Dataset, cols::MultiColumnIndex; alg = HeapSortAlg(), rev = false, mapformats::Bool = true, stable = true)
    isempty(ds) && return ds
    _check_consistency(ds)
    colsidx = index(ds)[cols]
    if length(rev) == 1
        revs = repeat([rev], length(colsidx))
    else
        revs = rev
    end

    @assert length(colsidx) == length(revs) "the reverse argument must be the same length as the length of selected columns"
    _check_for_fast_sort(ds, colsidx, revs, mapformats) == 0 && return ds
    _use_ds_perm = false
    if _check_for_fast_sort(ds, colsidx, revs, mapformats) == 1
        _use_ds_perm = true
    end
    skipcol = 0
    if _check_for_fast_sort(ds, colsidx, revs, mapformats) == 2
        skipcol = length(index(ds).sortedcols)
    end
    starts, perm, ngroups = _sortperm(ds, cols, revs; a = alg, mapformats = mapformats, stable = stable, skipcol = skipcol, skipcol_mkcopy = false)
    if _use_ds_perm || skipcol>0 #index(ds).perm and index(ds).starts already been updated
        empty!(index(ds).sortedcols)
        empty!(index(ds).rev)
        append!(index(ds).sortedcols, collect(colsidx))
        append!(index(ds).rev, revs)
        index(ds).ngroups[] = ngroups
        index(ds).fmt[] = mapformats
        if _use_ds_perm # current implementation needs to reset starts and use the computed one/ it should be optimised
            empty!(index(ds).starts)
            append!(index(ds).starts, starts)
        end
    elseif skipcol==0
        _reset_grouping_info!(ds)
        append!(index(ds).sortedcols, collect(colsidx))
        append!(index(ds).rev, revs)
        append!(index(ds).perm, perm)
        append!(index(ds).starts, starts)
        index(ds).ngroups[] = ngroups
        index(ds).fmt[] = mapformats
    end

    _modified(_attributes(ds))
    _permute_ds_after_sort!(ds, perm)
    ds
end

Base.sort!(ds::Dataset, col::ColumnIndex; alg = HeapSortAlg(), rev::Bool = false, mapformats::Bool = true, stable =true) = sort!(ds, [col], rev = rev, alg = alg, mapformats = mapformats, stable = stable)

function Base.sort(ds::Dataset, cols::MultiColumnIndex; alg = HeapSortAlg(), rev = false, mapformats::Bool = true, stable = true)
    isempty(ds) && return copy(ds)
    _check_consistency(ds)
    colsidx = index(ds)[cols]
    if length(rev) == 1
        revs = repeat([rev], length(colsidx))
    else
        revs = rev
    end

    @assert length(colsidx) == length(revs) "the reverse argument must be the same length as the length of selected columns"
    _check_for_fast_sort(ds, colsidx, revs, mapformats) == 0 && return copy(ds)
    skipcol = 0
    if _check_for_fast_sort(ds, colsidx, revs, mapformats) == 2
        skipcol = length(index(ds).sortedcols)
    end
    starts, perm, ngroups = _sortperm(ds, cols, revs; a = alg, mapformats = mapformats, stable = stable, skipcol = skipcol, skipcol_mkcopy = true)
    newds = ds[perm, :]
    _reset_grouping_info!(newds)
    append!(index(newds).sortedcols, collect(colsidx))
    append!(index(newds).rev, revs)
    append!(index(newds).perm, collect(perm))
    append!(index(newds).starts, starts)
    index(newds).ngroups[] = ngroups
    index(newds).fmt[] = mapformats
    newds
end

Base.sort(ds::Dataset, col::ColumnIndex; alg = HeapSortAlg(), rev::Bool = false, mapformats::Bool = true, stable =true) = sort(ds, [col], rev = rev, alg = alg, mapformats = mapformats, stable = stable)

function Base.sort(ds::SubDataset, cols::MultiColumnIndex; alg = HeapSortAlg(), rev = false, mapformats::Bool = true, stable = true, view = false)
    _check_consistency(ds)
    colsidx = index(ds)[cols]
    if rev isa AbstractVector
        @assert length(rev) == length(colsidx) "length of rev and the number of selected columns must match"
        revs = rev
    else
        revs = repeat([rev], length(colsidx))
    end
    starts, idx, last_valid_range =  _sortperm_v(ds, cols, revs, stable = stable, a = alg, mapformats = mapformats)
    if view
        Base.view(ds, idx, :)
    else
        newds = ds[idx, :]
        append!(index(newds).sortedcols, collect(colsidx))
        append!(index(newds).rev, revs)
        append!(index(newds).perm, idx)
        append!(index(newds).starts, starts)
        index(newds).ngroups[] = last_valid_range
        index(newds).fmt[] = mapformats
        newds
    end
end

Base.sort(ds::SubDataset, col::ColumnIndex; alg = HeapSortAlg(), rev::Bool = false, mapformats::Bool = true, stable =true) = sort(ds, [col], rev = rev, alg = alg, mapformats = mapformats, stable = stable)

function Base.sortperm(ds::SubDataset, cols; alg = HeapSortAlg(), rev = false, mapformats::Bool = true, stable = true)
    isempty(ds) && return []
    _check_consistency(ds)
    colsidx = index(ds)[cols]
    _sortperm_v(ds, cols, rev, a = alg, mapformats = mapformats, stable = stable)[2]
end

function unsort!(ds::Dataset)
    isempty(ds) && return ds
    if isempty(index(ds).perm) #if perm is empty everything else should be empty, here we just make sure
        empty!(index(ds).sortedcols)
        empty!(index(ds).rev)
        empty!(index(ds).starts)
        index(ds).ngroups[] = 1
        index(ds).grouped[] = false
        return ds
    else
        _permute_ds_after_sort!(ds, invperm(index(ds).perm))
        # TODO we may don't need to reset grouping info
        _reset_grouping_info!(ds)
        ds
    end
end

function Base.issorted(ds::AbstractDataset, cols::MultiColumnIndex; rev = false, mapformats = true)
    _issorted(ds, cols, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64), rev = rev, mapformats = mapformats)[1]
end
Base.issorted(ds::AbstractDataset, col::ColumnIndex; rev = false, mapformats = true) = issorted(ds, [col], rev = rev, mapformats = mapformats)

function issorted!(ds::Dataset, cols::MultiColumnIndex; rev = false, mapformats = true)
    res, starts, lastvalid, colsidx, revs, mapformats = _issorted(ds, cols, nrow(ds) < typemax(Int32) ? Val(Int32) : Val(Int64), rev = rev, mapformats = mapformats)
    if res
        _reset_grouping_info!(ds)
        append!(index(ds).sortedcols, collect(colsidx))
        append!(index(ds).rev, revs)
        append!(index(ds).perm, collect(1:nrow(ds)))
        append!(index(ds).starts, starts)
        index(ds).ngroups[] = lastvalid
        index(ds).fmt[] = mapformats
    end
    res
end

issorted!(ds::Dataset, col::ColumnIndex; rev = false, mapformats = true) = issorted!(ds, [col], rev = rev, mapformats = mapformats)

function _issorted(ds, cols::MultiColumnIndex, ::Val{T}; rev = false, mapformats = true) where T
    colsidx = index(ds)[cols]
    if rev isa AbstractVector
        @assert length(rev) == length(colsidx) "length of rev and the number of selected columns must match"
        revs = rev
    else
        revs = repeat([rev], length(colsidx))
    end
    by = Function[]

    if mapformats
        for j in 1:length(colsidx)
            push!(by, getformat(parent(ds), colsidx[j]))
        end
    else
        for j in 1:length(colsidx)
            push!(by, identity)
        end
    end
    res = true
    starts = Vector{T}(undef, nrow(ds))
    starts[1] = 1
    lastvalid = 1
    inbits = zeros(Bool, nrow(ds))
    inbits[1] = true
    for j in 1:length(colsidx)
        v = _columns(ds)[colsidx[j]]
        _ord =  Base.Order.ord(isless, by[j], revs[j])
        part_res = _issorted_check_for_each_range(v, starts, lastvalid, _ord, nrow(ds))
        !part_res && return false, starts, lastvalid, colsidx, revs, mapformats
        _find_starts_of_groups!(_columns(ds)[colsidx[j]], 1:nrow(ds), by[j], inbits)
        lastvalid = _fill_starts_from_inbits!(starts, inbits)
        lastvalid == nrow(ds) && return true, starts, lastvalid, colsidx, revs, mapformats
        # lastvalid = _fill_starts_v2!(starts, inbits, _columns(ds)[colsidx[j]], lastvalid, Base.Order.ord(isless, by[j], revs[j]), Val(T))
    end
    res, starts, lastvalid, colsidx, revs, mapformats
end

function _issorted_check_for_each_range(v, starts, lastvalid, _ord, nrows)
    part_res = ones(Bool, Threads.nthreads())
    Threads.@threads for rng in 1:lastvalid
        lo = starts[rng]
        rng == lastvalid ? hi = nrows : hi = starts[rng+1] - 1
        part_res[Threads.threadid()] = _issorted_barrier(v, _ord, lo, hi)
        !part_res[Threads.threadid()] &&  break
    end
    all(part_res)
end

function _fill_starts_from_inbits!(starts, inbits)
    lastvalid = 1
    @inbounds for i in 1:length(inbits)
        if inbits[i]
            starts[lastvalid] = i
            lastvalid += 1
        end
    end
    lastvalid - 1
end

function _issorted_barrier(v, _ord, lo, hi)
    lo >= hi && return true
    @inbounds for i in lo+1:hi
        Base.Order.lt(_ord, v[i], v[i-1]) && return false
    end
    true
end

ext: jl


text:
const TableauHSPARKprimary = AbstractTableauSPARK{:hspark_primary}
const ParametersHSPARKprimary = AbstractParametersSPARK{:hspark_primary}

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for Hamiltonian systems *EXPERIMENTAL*.

This integrator solves the following system of equations for the internal stages,
```math
\begin{aligned}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
0 &= \sum \limits_{j=1}^{r} \omega_{ij} \tilde{\Phi}_{n,j} , & i &= 1, ..., r-1 ,
\end{aligned}
```
with definitions
```math
\begin{aligned}
V_{n,i} &= \hphantom{-} \frac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &=           -  \frac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{aligned}
```
and update rule
```math
\begin{aligned}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) .
\end{aligned}
```
"""
struct IntegratorHSPARKprimary{DT, TT, D, S, R, PT <: ParametersHSPARKprimary{DT,TT,D,S,R},
                                                ST <: NonlinearSolver{DT},
                                                IT <: InitialGuessPODE{TT}} <: AbstractIntegratorHSPARK{DT,TT,D,S,R}
    params::PT
    solver::ST
    iguess::IT
    caches::CacheDict{PT}

    function IntegratorHSPARKprimary(params::ParametersHSPARKprimary{DT,TT,D,S,R}, solver::ST, iguess::IT, caches) where {DT,TT,D,S,R,ST,IT}
        new{DT, TT, D, S, R, typeof(params), ST, IT}(params, solver, iguess, caches)
    end

    function IntegratorHSPARKprimary{DT,D}(equations::NamedTuple, tableau::TableauHSPARKprimary{TT}, Δt::TT) where {DT,TT,D,ST}
        @assert tableau.ρ == tableau.r-1

        # get number of stages
        S = tableau.s
        R = tableau.r
        P = tableau.ρ

        N = 2*D*S + 3*D*R

        # create params
        params = ParametersHSPARKprimary{DT,D}(equations, tableau, Δt)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, N, params, caches)

        # create initial guess
        iguess = InitialGuessPODE(get_config(:ig_extrapolation), equations[:v̄], equations[:f̄], Δt)

        # create integrator
        IntegratorHSPARKprimary(params, solver, iguess, caches)
    end

    function IntegratorHSPARKprimary(equation::Union{PDAE{DT}, HDAE{DT}}, tableau::TableauHSPARKprimary, Δt; kwargs...) where {DT}
        IntegratorHSPARKprimary{DT, ndims(equation)}(get_functions(equation), tableau, Δt; kwargs...)
    end
end

GeometricBase.nconstraints(::IntegratorHSPARKprimary{DT,TT,D}) where {DT,TT,D} = D

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,D,S,R},
                                        params::ParametersHSPARKprimary{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tqᵢ::TT
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[2*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[2*(D*(i-1)+k-1)+2]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tqᵢ = params.t + params.Δt * params.tab.q.c[i]
        tpᵢ = params.t + params.Δt * params.tab.p.c[i]
        params.equs[:v](tqᵢ, cache.Qi[i], cache.Pi[i], cache.Vi[i])
        params.equs[:f](tpᵢ, cache.Qi[i], cache.Pi[i], cache.Fi[i])
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[2*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.tab.λ.c[i]
        params.equs[:u](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.equs[:g](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.equs[:ϕ](tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.q̃ .+= params.Δt .* params.tab.q.b[i] .* cache.Vi[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[i] .* cache.Fi[i]
    end
    for i in 1:R
        cache.q̃ .+= params.Δt .* params.tab.q.β[i] .* cache.Up[i]
        cache.p̃ .+= params.Δt .* params.tab.p.β[i] .* cache.Gp[i]
    end

    # compute ϕ(q,p)
    tλᵢ = params.t + params.Δt
    params.equs[:ϕ](tλᵢ, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
function Integrators.function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersHSPARKprimary{DT,TT,D,S,R,P},
                                      caches::CacheDict) where {ST,DT,TT,D,S,R,P}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(y, cache, params)

    # compute b = - [(Y-AV-AU), (Z-AF-AG)]
    for i in 1:S
        for k in 1:D
            b[2*(D*(i-1)+k-1)+1] = - cache.Yi[i][k]
            b[2*(D*(i-1)+k-1)+2] = - cache.Zi[i][k]
            for j in 1:S
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * cache.Vi[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*(D*(i-1)+k-1)+1] += params.tab.q.α[i,j] * cache.Up[j][k]
                b[2*(D*(i-1)+k-1)+2] += params.tab.p.α[i,j] * cache.Gp[j][k]
            end
        end
    end

    # compute b = - [(Y-AV-AU), (Z-AF-AG), ωΦ]
    for i in 1:R
        for k in 1:D
            b[2*D*S+3*(D*(i-1)+k-1)+1] = - cache.Yp[i][k]
            b[2*D*S+3*(D*(i-1)+k-1)+2] = - cache.Zp[i][k]
            for j in 1:S
                b[2*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * cache.Vi[j][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * cache.Fi[j][k]
            end
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+1] += params.tab.q̃.α[i,j] * cache.Up[j][k]
                b[2*D*S+3*(D*(i-1)+k-1)+2] += params.tab.p̃.α[i,j] * cache.Gp[j][k]
            end
        end
    end
    for i in 1:R-P
        for k in 1:D
            b[2*D*S+3*(D*(i-1)+k-1)+3] = 0
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,j] * cache.Φp[j][k]
            end
            b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.ω[i,R+1] * cache.ϕ̃[k]
        end
    end

    # compute b = d_λ ⋅ Λ
    for i in R-P+1:R
        for k in 1:D
            b[2*D*S+3*(D*(R-1)+k-1)+3] = 0
            for j in 1:R
                b[2*D*S+3*(D*(i-1)+k-1)+3] -= params.tab.δ[j] * cache.Λp[j][k]
            end
        end
    end
end


text:
module DD
nothing
end


text:
abstract type ThresholdStrategy end

struct FixedThresholdStrategy <: ThresholdStrategy
    threshold::Real
end

function calculate_threshold(threshold_strat::FixedThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    return threshold_strat.threshold
end

struct NumLabelTresholdStrategy <: ThresholdStrategy
    eps::Real
    function NumLabelTresholdStrategy(eps::Real)
        eps >= 0 || throw(ArgumentError("Negative eps: $(eps)."))
        eps <= 1 || throw(ArgumentError("Eps bigger than 1: $(eps)."))
        new(eps)
    end
end

function calculate_threshold(threshold_strat::NumLabelTresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    return threshold_strat.eps * length(labels)
end

struct MaxDensityThresholdStrategy <: ThresholdStrategy
    eps::Real
    function MaxDensityThresholdStrategy(eps::Real)
        eps >= 0 || throw(ArgumentError("Negative eps: $(eps)."))
        eps <= 1 || throw(ArgumentError("Eps bigger than 1: $(eps)."))
        new(eps)
    end
end

function calculate_threshold(threshold_strat::MaxDensityThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    return threshold_strat.eps * maximum(kde(kde_cache))
end

struct OutlierPercentageThresholdStrategy <: ThresholdStrategy
    eps::Union{Real, Nothing}
    function OutlierPercentageThresholdStrategy(eps::Union{Real, Nothing}=nothing)
        if eps !== nothing
            eps >= 0 || throw(ArgumentError("Negative eps: $(eps)."))
            eps <= 1 || throw(ArgumentError("Eps bigger than 1: $(eps)."))
        end
        new(eps)
    end
end

function calculate_threshold(threshold_strat::OutlierPercentageThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    if threshold_strat.eps === nothing
        eps = count(x -> x .== :outlier, labels) / length(labels)
    else
        eps = threshold_strat.eps
    end
    d_x = sort(kde(kde_cache))
    n = length(labels)
    idx = max(1, min(n, floor(Int, eps * n)))
    return (d_x[idx] + d_x[min(n, idx+1)]) / 2
end

struct GroundTruthThresholdStrategy <: ThresholdStrategy end

function calculate_threshold(threshold_strat::GroundTruthThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    eps = count(x -> x .== :outlier, labels) / length(labels)
    d_x = sort(kde(kde_cache))
    n = length(labels)
    idx = max(1, min(n, floor(Int, eps * n)))
    return (d_x[idx] + d_x[min(n, idx+1)]) / 2
end

function calculate_threshold(threshold_strat::ThresholdStrategy, data::DataSet, labels::Vector{Symbol},
                             gamma::Union{Real, Symbol})
    return calculate_threshold(threshold_strat, KDECache(data, gamma), labels)
end

function split_masks(threshold_strat::GroundTruthThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    threshold = calculate_threshold(threshold_strat, kde_cache, labels)
    inlier_mask = labels .== :inlier
    outlier_mask = .!inlier_mask
    return threshold, inlier_mask, outlier_mask
end

function split_masks(threshold_strat::Nothing, kde_cache::KDECache, labels::Vector{Symbol})
    return 0.0, trues(length(labels)), falses(length(labels))
end

function split_masks(threshold_strat::ThresholdStrategy, kde_cache::KDECache, labels::Vector{Symbol})
    threshold = calculate_threshold(threshold_strat, kde_cache, labels)
    d_x = kde(kde_cache)
    inlier_mask = d_x .>= threshold
    outlier_mask = .!inlier_mask
    return threshold, inlier_mask, outlier_mask
end

function split_masks(threshold_strat::ThresholdStrategy, data::DataSet, labels::Vector{Symbol}, gamma::Union{Nothing, Real, Symbol}=nothing)
    return split_data(threshold_strat, KDECache(data, gamma), labels)
end

mutable struct PWC
    kde_cache::KDECache
    threshold_strat::ThresholdStrategy
    threshold::Union{Real, Nothing}
    function PWC(kde_cache::KDECache, threshold_strat::ThresholdStrategy=OutlierPercentageThresholdStrategy(),
                 threshold::Union{Real, Nothing}=nothing)
        return new(kde_cache, threshold_strat, threshold)
    end
end

function PWC(data, gamma::Union{Real, Symbol}, threshold_strat::ThresholdStrategy=OutlierPercentageThresholdStrategy())
    return PWC(KDECache(data, gamma), threshold_strat)
end

function calculate_threshold!(pwc::PWC, labels::Vector{Symbol})
    pwc.threshold = calculate_threshold(pwc.threshold_strat, pwc.kde_cache, labels)
    return nothing
end

function predict(pwc::PWC, data::DataSet)
    if pwc.threshold === nothing
        throw(ArgumentError("PWC threshold not initialized."))
    end
    d_x = kde(pwc.kde_cache, data)
    pred = fill(:outlier, size(data, 2))
    pred[d_x .> pwc.threshold] .= :inlier
    return pred
end


text:
using DataFrames
using Random: randstring
using JDF
using Test

df = DataFrame([collect(1:100) for i = 1:3000]);
df[!, :strs] = [randstring(8) for i = 1:DataFrames.nrow(df)]
df[!, :substrs] = [SubString(x, 1, 3) for x in df[!, :strs]]

io = IOBuffer()

JDF.save("tmp-substring.jdf", df)

df2 = JDF.load("tmp-substring.jdf")

@test df.substrs == df2.substrs

rm("tmp-substring.jdf", force=true, recursive=true)


text:
# ### Surface (3)

#md # ![FILE_NAME.png](images/FILE_NAME.png)

f(x,y) = x^2 + y^2
x = -10:10
y = x
plot(x, y, f, linetype=:wireframe)


text:
export HRepElement, HalfSpace, HyperPlane
export VRepElement, Ray, Line
export islin, isray, ispoint, coord, lift, simplify

_vec(::Type{T}, a::AbstractVector) where {T} = AbstractArray{T}(a)
_vec(::Type{T}, a::AbstractVector{T}) where {T} = a
function _vec(::Type{T}, a::StaticArrays.SVector{N}) where {N, T}
    StaticArrays.SVector{N, T}(a)
end
_vec(::Type{T}, a::StaticArrays.SVector{N, T}) where {N, T} = a

abstract type HRepElement{T, AT} end

"""
    struct HalfSpace{T, AT} <: HRepElement{T, AT}
        a::AT
        β::T
    end

An halfspace defined by the set of points ``x`` such that ``\\langle a, x \\rangle \\le \\beta``.
"""
struct HalfSpace{T, AT <: AbstractVector{T}} <: HRepElement{T, AT}
    a::AT
    β::T
    function HalfSpace{T, AT}(a::AT, β::T) where {T, AT<:AbstractVector{T}}
        new{T, AT}(a, β)
    end
end

HalfSpace{T}(a::AT, β::T) where {T, AT <: AbstractVector{T}} = HalfSpace{T, AT}(a, β)
HalfSpace{T}(a::AbstractVector, β) where {T} = HalfSpace{T}(_vec(T, a), T(β))

"""
    struct HyperPlane{T, AT} <: HRepElement{T, AT}
        a::AT
        β::T
    end

An hyperplane defined by the set of points ``x`` such that ``\\langle a, x \\rangle = \\beta``.
"""
struct HyperPlane{T, AT<:AbstractVector{T}} <: HRepElement{T, AT}
    a::AT
    β::T
    function HyperPlane{T, AT}(a::AT, β::T) where {T, AT<:AbstractVector{T}}
        new{T, AT}(a, β)
    end
end

HyperPlane{T}(a::AT, β::T) where {T, AT <: AbstractVector{T}} = HyperPlane{T, AT}(a, β)
HyperPlane{T}(a::AbstractVector, β) where {T} = HyperPlane{T}(_vec(T, a), T(β))

HalfSpace(a::AbstractVector{T}, β::T) where {T} = HalfSpace{T}(a, β)
HalfSpace(a::AbstractVector{S}, β::T) where {S, T} = HalfSpace{promote_type(S, T)}(a, β)
HyperPlane(a::AbstractVector{T}, β::T) where {T} = HyperPlane{T}(a, β)
HyperPlane(a::AbstractVector{S}, β::T) where {S, T} = HyperPlane{promote_type(S, T)}(a, β)

function Base.convert(::Type{HalfSpace{T, AT}}, h::HalfSpace) where {T, AT}
    return HalfSpace{T, AT}(convert(AT, h.a), convert(T, h.β))
end
function Base.convert(::Type{HyperPlane{T, AT}}, h::HyperPlane) where {T, AT}
    return HyperPlane{T, AT}(convert(AT, h.a), convert(T, h.β))
end
Base.convert(::Type{HalfSpace{T, AT}}, h::HalfSpace{T, AT}) where {T, AT} = h
Base.convert(::Type{HyperPlane{T, AT}}, h::HyperPlane{T, AT}) where {T, AT} = h

islin(::Union{HalfSpace, Type{<:HalfSpace}}) = false
islin(::Union{HyperPlane, Type{<:HyperPlane}}) = true

Base.:-(h1::HRepElement, h2::HRepElement) = HalfSpace(h1.a - h2.a, h1.β - h2.β)
Base.:-(h1::HyperPlane, h2::HyperPlane) = HyperPlane(h1.a - h2.a, h1.β - h2.β)

Base.:(*)(h::HyperPlane, α::Real) = HyperPlane(h.a * α, h.β * α)
Base.:(*)(α::Real, h::HyperPlane) = HyperPlane(α * h.a, α * h.β)
Base.:(*)(h::HalfSpace, α::Real) = HalfSpace(h.a * α, h.β * α)
Base.:(*)(α::Real, h::HalfSpace) = HalfSpace(α * h.a, α * h.β)

function Base.:(/)(h::ElemT, P::UniformScaling) where {T, ElemT<:HRepElement{T}}
    Tout = _promote_type(T, eltype(P))
    ElemTout = similar_type(ElemT, FullDim(h), Tout)
    ElemTout(P * _vec(Tout, h.a), Tout(h.β))
end
function Base.:(/)(h::ElemT, P::AbstractMatrix) where {T, ElemT<:HRepElement{T}}
    Tout = _promote_type(T, eltype(P))
    ElemTout = similar_type(ElemT, size(P, 2), Tout)
    ElemTout(AbstractMatrix{Tout}(P) * _vec(Tout, h.a), Tout(h.β))
end

# Point: -> A same Rep should always return the same of the two types so that when points and sympoints will have different accessors it will be type stable
# AbstractVector{T}
# Ray:
# Ray{T}
# Linear Ray:
# Line{T}

origin(::Type{<:SparseVector{T}}, N::Int) where {T} = spzeros(T, N)
origin(::Type{Vector{T}}, N::Int) where {T} = zeros(T, N)
origin(VT::Type{<:AbstractVector}, ::Int) = zeros(VT)
# Canonical basis vector
basis(::Type{<:SparseVector{T}}, N::Int, i::Int) where {T} = sparsevec([i], [one(T)], N)
function basis(::Type{Vector{T}}, N::Int, i::Int) where {T}
    v = zeros(T, N)
    v[i] = one(T)
    v
end
function basis(::Type{StaticArrays.SVector{N, T}}, ::StaticArrays.Size, i::Int) where {N, T}
    StaticArrays.SVector{N, T}(ntuple(j -> j == i ? one(T) : zero(T), Val(N)))
end

"""
    struct Ray{T, AT <: AbstractVector{T}}
        a::AT
    end

The conic hull of `a`, i.e. the set of points `λa` where `λ` is any nonnegative real number.
"""
struct Ray{T, AT <: AbstractVector{T}}
    a::AT
    function Ray{T, AT}(a::AT) where {T, AT<:AbstractVector{T}}
        new{T, AT}(a)
    end
end
Ray{T, AT}(ray::Ray) where {T, AT} = Ray{T, AT}(AT(ray.a))
Base.convert(::Type{Ray{T, AT}}, ray::Ray) where {T, AT} = Ray{T, AT}(ray)
Ray{T}(a::AT) where {T, AT<:AbstractVector{T}} = Ray{T, AT}(a)
Ray(a::AbstractVector) = Ray{eltype(a)}(a)

"""
    struct Line{T, AT <: AbstractVector{T}}
        a::AT
    end

The conic hull of `a` and `-a`, i.e. the set of points `λa` where `λ` is any real number.
"""
struct Line{T, AT<:AbstractVector{T}}
    a::AT
    function Line{T, AT}(a::AT) where {T, AT<:AbstractVector{T}}
        new{T, AT}(a)
    end
end
Line{T, AT}(line::Line) where {T, AT} = Line{T, AT}(AT(line.a))
Base.convert(::Type{Line{T, AT}}, line::Line) where {T, AT} = Line{T, AT}(line)
Line{T}(a::AT) where {T, AT<:AbstractVector{T}} = Line{T, AT}(a)
Line(a::AbstractVector) = Line{eltype(a)}(a)

Base.:(==)(a::Ray, b::Ray) = coord(a) == coord(b)
Base.:(==)(a::Line, b::Line) = coord(a) == coord(b)

const VStruct{T, AT} = Union{Line{T, AT}, Ray{T, AT}}
Base.getindex(x::VStruct, i) = x.a[i]
Base.vec(x::VStruct) = vec(x.a)

Base.:-(h::ElemT) where {ElemT<:Union{HyperPlane, HalfSpace}} = ElemT(-h.a, -h.β)
Base.:-(elem::ElemT) where ElemT<:VStruct = ElemT(-coord(elem))
# Used in remproj
Base.:-(p::AbstractVector, l::Line) = p - coord(l)
# `Ray - Line` and `Line - Line` are done in `remproj`
Base.:-(r::Ray, s::Union{Ray, Line}) = Ray(r.a - s.a)
Base.:-(r::Line, s::Line) = Line(r.a - s.a)
Base.:+(r::Ray, s::Ray) = Ray(r.a + s.a)
Base.:+(p::AbstractVector, r::Ray) = p + coord(r)

for op in [:dot, :cross]
    @eval begin
        LinearAlgebra.$op(x::VStruct, y) = $op(x.a, y)
        LinearAlgebra.$op(x, y::VStruct) = $op(x, y.a)
        LinearAlgebra.$op(x::VStruct, y::VStruct) = $op(x.a, y.a)
    end
end

for ElT in (:HyperPlane, :HalfSpace, :Line, :Ray)
    @eval begin
        Base.promote_rule(::Type{$ElT{T, AT}}, ::Type{$ElT{T, AT}}) where {T, AT} = $ElT{T, AT}
        # Allowing mixing e.g. sparse vector with vector would not be helpful as code meant
        # to use sparse polyhedra would lose sparsity silently. Same thing for StaticArrays.SVector
        Base.promote_rule(::Type{$ElT{T, VT}}, ::Type{$ElT{T, WT}}) where {T, VT, WT} = error("Cannot mix Polyhedra elements of vector type $VT and $WT")
        function Base.promote_rule(::Type{$ElT{S, AS}}, ::Type{$ElT{T, AT}}) where {S, T, AS, AT}
            U = promote_type(S, T)
            VU = similar_type(AS, U)
            WU = similar_type(AT, U)
            # Using promote_type instead of promote_rule here will promote to `Vector`
            # when `VU` is `Vector` and `WU` is `SVector` in Julia v1.2 instead of throwing
            # the error above.
            return promote_rule($ElT{U, VU}, $ElT{U, WU})
        end
    end
end

Base.:(*)(α, r::T) where T<:VStruct = T(α * r.a)
Base.:(*)(r::T, α) where T<:VStruct = T(r.a * α)
Base.:(/)(r::T, α) where T<:VStruct = T(r.a / α)

const VRepElement{T} = Union{VStruct{T}, AbstractVector{T}}
const RepElement{T} = Union{HRepElement{T}, VRepElement{T}}
const StructElement{T, AT} = Union{VStruct{T, AT}, HRepElement{T, AT}}

vectortype(::Type{<:StructElement{T, AT}}) where {T, AT} = AT
vectortype(AT::Type{<:AbstractVector}) = AT

FullDim(::Type{<:StructElement{T, AT}}) where {T, AT} = FullDim(AT)
FullDim(el::StructElement) = FullDim(coord(el))
coefficient_type(::Union{RepElement{T}, Type{<:RepElement{T}}}) where {T} = T

islin(::Union{Line, Type{<:Line}}) = true
islin(::Union{AbstractVector, Ray, Type{<:Union{AbstractVector, Ray}}}) = false
ispoint(::Union{AbstractVector, Type{<:AbstractVector}}) = true
ispoint(::Union{Line, Ray, Type{<:Union{Line, Ray}}}) = false
isray(v) = !ispoint(v)

coord(v::AbstractVector) = v
coord(v::Union{HRepElement, VStruct}) = v.a

function Base.:*(P::AbstractMatrix, v::ElemT) where {T, ElemT<:VStruct{T}}
      Tout = _promote_type(T, eltype(P))
      ElemTout = similar_type(ElemT, size(P, 1), Tout)
      return ElemTout(P * v.a)
end

function zeropad(a::AbstractSparseVector{T}, n::Integer, start::Integer) where T
    if iszero(n)
        return a
    elseif iszero(start)
        return zeropad(a, -n)
    elseif start == length(a)
        return zeropad(a, n)
    else
        return [a[1:start]; spzeros(T, n); a[start+1:end]]::AbstractSparseVector{T}
    end
end

function zeropad(a::AbstractSparseVector{T}, n::Integer) where T
    if iszero(n)
        return a
    elseif n < 0
        # Add type assert to check that vcat does not make it dense
        return [spzeros(T, -n); a]::AbstractSparseVector{T}
    else
        return [a; spzeros(T, n)]::AbstractSparseVector{T}
    end
end
function (zeropad(a::StaticArrays.SVector{N1, T}, ::StaticArrays.Size{N2})::StaticArrays.SVector{N1+abs(N2[1]), T}) where {T, N1, N2}
    # TODO, should probably de a generated function to make it type stable
    if iszero(N2[1])
        return a
    else
        z = StaticArrays.@SVector zeros(T, abs(N2[1]))
        if N2[1] < 0
            return vcat(z, a)::StaticArrays.SVector{N1-N2[1], T}
        else
            return vcat(a, z)::StaticArrays.SVector{N1+N2[1], T}
        end
    end
end
function zeropad(a::Vector{T}, n::Integer) where T
    if iszero(n)
        return a
    elseif n < 0
        return [zeros(T, -n); a]
    else
        return [a; zeros(T, n)]
    end
end
zeropad(h::HRepElement, d, start) = constructor(h)(zeropad(h.a, d, start), h.β)
zeropad(h::HRepElement, d::FullDim) = constructor(h)(zeropad(h.a, d), h.β)
zeropad(v::VStruct, d::FullDim)     = constructor(v)(zeropad(v.a, d))
# Called in cartesian product of two polyhedra, only the second using static arrays.
zeropad(a::Vector, ::StaticArrays.Size{N}) where N = zeropad(a, N[1])
zeropad(a::StaticArrays.SVector, n::Int) = zeropad(collect(a), n)

for ElemT in [:HalfSpace, :HyperPlane, :Ray, :Line]
    @eval begin
        function similar_type(::Type{$ElemT{T, AT}}, dout::FullDim, ::Type{Tout}) where {T, AT, Tout}
            return $ElemT{Tout, similar_type(AT, dout, Tout)}
        end
    end
end

ininterior(r::Ray, h::HalfSpace) = _neg(h.a ⋅ r)
ininterior(l::Line, h::HalfSpace) = _neg(h.a ⋅ l)
ininterior(p::AbstractVector, h::HalfSpace) = _lt(h.a ⋅ p, h.β)

inrelativeinterior(p::VRepElement, h::HalfSpace) = ininterior(p, h)

Base.in(r::Ray, h::HalfSpace) = _nonpos(h.a ⋅ r)
Base.in(l::Line, h::HalfSpace) = _nonpos(h.a ⋅ l)
Base.in(p::AbstractVector, h::HalfSpace) = _leq(h.a ⋅ p, h.β)

ininterior(p::VRepElement, h::HyperPlane) = false
inrelativeinterior(p::VRepElement, h::HyperPlane) = p in h

Base.in(r::Ray, h::HyperPlane) = isapproxzero(h.a ⋅ r)
Base.in(l::Line, h::HyperPlane) = isapproxzero(h.a ⋅ l)
Base.in(p::AbstractVector, h::HyperPlane) = _isapprox(h.a ⋅ p, h.β)

#function Base.vec(x::FixedVector{T}) where {T}
#    y = Vector{T}(N)
#    for i in 1:N
#        y[i] = x[i]
#    end
#    y
#end
#Base.vec{ElemT::VStruct}(x::ElemT) = ElemT(vec(x.a))

function pushbefore(a::AbstractSparseVector{T}, β::T) where T
    b = spzeros(T, length(a)+1)
    b[1] = β
    b[2:end] = a
    b
end
pushbefore(a::AbstractVector, β) = [β; a]
function pushbefore(a::StaticArrays.SVector, β)
    StaticArrays.SVector(β, a...)
end

function push_after(a::AbstractSparseVector{T}, β::T) where T
    b = spzeros(T, length(a)+1)
    b[end] = β
    b[1:end-1] = a
    return b
end
push_after(a::AbstractVector, β) = [a; β]
function push_after(a::StaticArrays.SVector, β)
    StaticArrays.SVector(a..., β)
end

constructor(::HyperPlane) = HyperPlane
constructor(::HalfSpace) = HalfSpace
constructor(::Ray) = Ray
constructor(::Line) = Line

function lift(h::HRepElement{T}, pre::Bool=true) where T
    a = pre ? pushbefore(h.a, -h.β) : push_after(h.a, -h.β)
    constructor(h)(a, zero(T))
end
function complement_lift(h::HRepElement{T}, pre::Bool=true) where T
    a = pre ? pushbefore(h.a, h.β) : push_after(h.a, h.β)
    constructor(h)(a, h.β)
end

lift(v::VStruct{T}) where {T} = constructor(v)(pushbefore(v.a, zero(T)))
lift(v::AbstractVector{T}) where {T} = pushbefore(v, one(T))

translate(p::AbstractVector, v) = p + v
translate(r::VStruct, v) = r

translate(h::ElemT, p) where {ElemT<:HRepElement} = ElemT(h.a, h.β + h.a ⋅ p)

_simplify(a::AbstractVector) = a
_simplify(a::AbstractVector, β) = a, β

iszo(g) = iszero(g) || g == 1 # TODO isone in v0.7

function _simplify(a::AbstractVector{<:Integer})
    g = gcd(a)
    if !iszo(g)
        div.(a, g)
    else
        a
    end
end
function _simplify(a::AbstractVector{<:Integer}, β::Integer)
    g = gcd(gcd(a), β)
    if !iszo(g)
        div.(a, g), div(β, g)
    else
        a, β
    end
end
function _simplify(a::AbstractVector{<:Rational})
    g = gcd(numerator.(a))
    if !iszo(g)
        a = a ./ g
    end
    g = gcd(denominator.(a))
    if !iszo(g)
        a .* g
    else
        a
    end
end
function _simplify(a::AbstractVector{<:Rational}, β::Rational)
    g = gcd(gcd(numerator.(a)), numerator(β))
    if !iszo(g)
        a = a ./ g
        β = β / g
    end
    g = gcd(gcd(denominator.(a)), denominator(β))
    if !iszo(g)
        a .* g, β * g
    else
        a, β
    end
end

simplify(h::HalfSpace{T}) where {T} = HalfSpace{T}(_simplify(h.a, h.β)...)
simplify(h::HyperPlane{T}) where {T} = HyperPlane{T}(_simplify(h.a, h.β)...)
simplify(r::VStruct) = constructor(r)(_simplify(coord(r)))
# Cannot scale points
simplify(p::AbstractVector) = p


text:
# export countt
# export collision
# export bordercollisions
# export hastop
# export hasleft
# export hasright
# export hasbot
# export mapzeroes

function countt(x)	#Counting function, counts the number of particles in a cell be performing integer
			#Division of it's number and the particle direction primes.
    count=0
    if x==1		#1 Corresponds to an empty cell.
        return 0
    end
    for k in [2,3,5,7]
        if x%k==0
            count+=1
        end
    end
    return count
end
function collision(x)	#This function computes the result of a particle-particle collision.
    if (x==1) | (x==0)	#x=1 corresponds to an empty cell, x=0 corresponds to an inert object.
        return x
    end
    if x%10==0 	 	#10=2 x 5
        return 21
    elseif x%21==0		#21=7 x 3
        return 10
    else
        return x
    end
end

function bordercollisions(x)		#This function computes the result of a particle-wall collision.
    if ((x%11!=0) & (x%13!=0) & (x%17!=0) & (x%19!=0)) | (x==0)
        return x
    end
    R=x
    if x%22==0			#22=11 x 2	55=11 x 5
        R=(R/22)*55
    end
    if x%39==0			#39=13 x 3	91=13 x 7
        R=(R/39)*91
    end
    if x%85==0             	#85=17 x 5	34=17 x 2
        R=R/(85)*34
    end
    if x%133==0            	#133=19 x 7	57= 19 x 3
        #print("BOOM")
        R=R/(133)*57

    end
    return Int(R)		#R is turned into an integer to accomodate it within the integer array.
end

#=========== Checking  functions ===================#

function hasleft(x)	#This function checks the right cell, for an incoming particle, or for a wall.
    if x==0		#If the right cell contains a 0 ( i.e it's a wall ) the current cell is multiplied by 19.
			#Which means a right border. This is useful for the bordercollisions function.
        return 19
    end
    X=x%3		#The value is divided by 3. which is an incoming particle. If the remainder is 0, the cell is 			      # multiplied by 3, as there's a particle heading towards it.
    ret=1
    if X==0
        return ret*3
    else
        return 1	#1 Is used as the neutral element for multiplication and division.
    end
end
function hastop(x)	#This function does the same for the bottom cell, with 2 as the incoming cell ( reversed ),
			# and 17 as a border value.
    X=x%2
    ret=1
    if x==0
        return 17
    end

    if X==0
        return ret*2
    else
        return 1
    end
end
function hasright(x)	#Right cell, 7 for incoming particle, 13 for border.
    X=x%7
    ret=1
    if x==0
        return 13
    end
    if X==0
        return ret*7
    else
        return 1
    end
end
function hasbot(x)	#Top cell, 5 for incoming particle, 11 for border.
    #X=x%5
    return (1 + (x%5)==0 * 4 + x==0 * 10)
    # return ret
    # if x==0
    #     return 11
    # end
    # if X==0
    #     return ret*5
    # else
    #     return 1
    # end
end

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

function mapzeroes(A1,A2)	#This function maps the zeroes of an array into another array.
				#Necessary for the functioning of intert objects. Otherwise the particles
				#Will just ignore the object.

    for i in enumerate(A1)

        if A1[i[1]]==0
            A2[i[1]]=0

        end
    end
    return A2
end


text:
# Test name methods
@testset "Basics" begin
    # initialize model and variable
    m = InfiniteModel()
    num = Float64(0)
    info = VariableInfo(false, num, false, num, false, num, false, num, false, false)
    new_info = VariableInfo(true, 0., true, 0., true, 0., true, 0., true, false)
    var = ScalarVariable(info)
    object = VariableData(var, "test")
    idx = FiniteVariableIndex(1)
    vref = FiniteVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, FiniteVariableIndex)
    bad_vref = FiniteVariableRef(m, FiniteVariableIndex(-1))
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(vref) === m
        @test owner_model(gvref) === m
    end
    # JuMP.index
    @testset "JuMP.index" begin
        @test index(vref) == idx
        @test index(gvref) == idx
    end
    # dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test dispatch_variable_ref(m, idx) == vref
        @test dispatch_variable_ref(gvref) == vref
    end
    # _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object) == idx
    end
    # _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(m, FiniteVariable) === m.finite_vars
        @test InfiniteOpt._data_dictionary(vref) === m.finite_vars
        @test InfiniteOpt._data_dictionary(gvref) === m.finite_vars
    end
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, vref)
        @test is_valid(m, gvref)
    end
    # _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(vref) === object
        @test InfiniteOpt._data_object(gvref) === object
        @test_throws ErrorException InfiniteOpt._data_object(bad_vref)
    end
    # _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(vref) === var
        @test InfiniteOpt._core_variable_object(gvref) === var
    end
    @testset "_variable_info" begin
        @test InfiniteOpt._variable_info(vref) == info
    end
    # _update_variable_info
    @testset "_update_variable_info" begin
        @test isa(InfiniteOpt._update_variable_info(vref, new_info), Nothing)
        @test InfiniteOpt._variable_info(vref) == new_info
    end
    # _set_core_variable_object
    @testset "_set_core_variable_object" begin
        @test InfiniteOpt._set_core_variable_object(vref, var) isa Nothing
    end
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(vref) == "test"
        @test name(gvref) == "test"
        @test name(bad_vref) == ""
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(vref, "new"), Nothing)
        @test name(vref) == "new"
        @test isa(set_name(gvref, "new2"), Nothing)
        @test name(vref) == "new2"
        @test_throws ErrorException set_name(bad_vref, "")
    end
    # test parameter_refs
    @testset "parameter_refs" begin
        @test parameter_refs(vref) == ()
        @test parameter_refs(gvref) == ()
    end
    # _make_variable_ref
    @testset "_make_variable_ref" begin
        @test InfiniteOpt._make_variable_ref(m, idx) == gvref
    end
    # _var_name_dict
    @testset "_var_name_dict" begin
        @test InfiniteOpt._var_name_dict(m) isa Nothing
    end
    # _update_var_name_dict
    @testset "_update_var_name_dict" begin
        m.name_to_var = Dict{String, ObjectIndex}()
        dict = InfiniteOpt._data_dictionary(vref)
        @test InfiniteOpt._update_var_name_dict(m, dict) isa Nothing
        @test InfiniteOpt._var_name_dict(m) == Dict(name(vref) => idx)
        m.name_to_var = nothing
    end
    # parameter_by_name
    @testset "JuMP.variable_by_name" begin
        # test normal
        @test variable_by_name(m, "new2") == gvref
        @test variable_by_name(m, "bob") isa Nothing
        # prepare variable with same name
        idx2 = FiniteVariableIndex(2)
        @test InfiniteOpt._add_data_object(m, object) == idx2
        vref2 = FiniteVariableRef(m, idx2)
        @test set_name(vref2, "new2") isa Nothing
        # test multiple name error
        @test_throws ErrorException variable_by_name(m, "new2")
    end
    # _delete_data_object
    @testset "_delete_data_object" begin
        @test InfiniteOpt._delete_data_object(vref) isa Nothing
        @test length(InfiniteOpt._data_dictionary(vref)) == 1
        @test !is_valid(m, vref)
    end
end

# Test variable definition methods
@testset "Definition" begin
    # initialize model and info
    m = InfiniteModel()
    num = Float64(0)
    info = VariableInfo(false, num, false, num, false, num, false, num, false, false)
    info2 = VariableInfo(true, num, true, num, true, num, true, num, true, false)
    info3 = VariableInfo(true, num, true, num, true, num, true, num, false, true)
    # build_variable
    @testset "JuMP.build_variable" begin
        # test normal
        @test build_variable(error, info).info == info
    end
    # _check_and_make_variable_ref
    @testset "_check_and_make_variable_ref" begin
        # test normal
        v = build_variable(error, info)
        idx = FiniteVariableIndex(1)
        vref = FiniteVariableRef(m, idx)
        @test InfiniteOpt._check_and_make_variable_ref(m, v, "") == vref
    end
    # add_variable
    @testset "JuMP.add_variable" begin
        idx = FiniteVariableIndex(2)
        vref = FiniteVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        v = build_variable(error, info)
        @test add_variable(m, v, "name") == gvref
        @test haskey(InfiniteOpt._data_dictionary(vref), idx)
        @test name(vref) == "name"
        # prepare variable with all the possible info additions
        v = build_variable(error, info2)
        # test info addition functions
        idx = FiniteVariableIndex(3)
        vref = FiniteVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test add_variable(m, v, "name") == gvref
        @test !optimizer_model_ready(m)
        # lower bound
        cindex = InfOptConstraintIndex(1)
        cref = InfOptConstraintRef(m, cindex)
        @test has_lower_bound(vref)
        @test InfiniteOpt._lower_bound_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.GreaterThan{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # upper bound
        cindex = InfOptConstraintIndex(2)
        cref = InfOptConstraintRef(m, cindex)
        @test has_upper_bound(vref)
        @test InfiniteOpt._upper_bound_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.LessThan{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # fix
        cindex = InfOptConstraintIndex(3)
        cref = InfOptConstraintRef(m, cindex)
        @test is_fixed(vref)
        @test InfiniteOpt._fix_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.EqualTo{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # binary
        cindex = InfOptConstraintIndex(4)
        cref = InfOptConstraintRef(m, cindex)
        @test is_binary(vref)
        @test InfiniteOpt._binary_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.ZeroOne}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        @test InfiniteOpt._constraint_dependencies(vref) == [InfOptConstraintIndex(i)
                                                             for i = 1:4]
    end
end

# Test the finite variable macro if no bounds are provided
@testset "Macro" begin
    # initialize model
    m = InfiniteModel()
    # test the deprecations
    @testset "@hold_variable" begin
        @test_macro_throws ErrorException @hold_variable(m, w)
    end
    @testset "@variable" begin
        # test regular
        idx = FiniteVariableIndex(1)
        vref = FiniteVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @variable(m, x >= 1, Bin) == gvref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test anan
        idx = FiniteVariableIndex(2)
        vref = FiniteVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @variable(m, binary = true, lower_bound = 1,
                        base_name = "x") == gvref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test array
        idxs = [FiniteVariableIndex(3), FiniteVariableIndex(4)]
        vrefs = [FiniteVariableRef(m, idx) for idx in idxs]
        gvrefs = [InfiniteOpt._make_variable_ref(m, idx) for idx in idxs]
        @test @variable(m, y[1:2] == 2, Int) == gvrefs
        @test name(vrefs[1]) == "y[1]"
        @test fix_value(vrefs[2]) == 2
        @test is_integer(vrefs[1])
        # test errors
        @test_macro_throws ErrorException @variable(m, x >= 1, Bin)
    end
end

# Test variable(s) constrained on creation
@testset "Creation Constraints" begin
    # initialize model and stuff
    m = InfiniteModel()
    info = VariableInfo(false, NaN, false, NaN, false, NaN, false, NaN, false, false)
    v1 = build_variable(error, info)
    # test add_variable with VariableConstrainedOnCreation
    @testset "Adding VariableConstrainedOnCreation" begin
        v = VariableConstrainedOnCreation(v1, MOI.Integer())
        vref = GeneralVariableRef(m, 1, FiniteVariableIndex)
        @test add_variable(m, v, "test") == vref
        @test name(vref) == "test"
        @test num_constraints(m, MOI.Integer) == 1
    end
    # test add_variable with VariablesConstrainedOnCreation
    @testset "Adding VariablesConstrainedOnCreation" begin
        v = VariablesConstrainedOnCreation([v1, v1, v1], MOI.SecondOrderCone(3),
                                           VectorShape())
        vrefs = [GeneralVariableRef(m, i, FiniteVariableIndex) for i in 2:4]
        names = ["test[$i]" for i in 1:3]
        @test add_variable(m, v, names) == vrefs
        @test name.(vrefs) == names
        @test num_constraints(m, MOI.SecondOrderCone) == 1
    end
end

# test usage methods
@testset "Usage" begin
    # initialize model and stuff
    m = InfiniteModel()
    @variable(m, y0)
    vref = dispatch_variable_ref(y0)
    # test used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(vref)
        push!(InfiniteOpt._measure_dependencies(vref), MeasureIndex(1))
        @test used_by_measure(y0)
        @test used_by_measure(vref)
        empty!(InfiniteOpt._measure_dependencies(vref))
    end
    # test used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(vref)
        push!(InfiniteOpt._constraint_dependencies(vref), InfOptConstraintIndex(1))
        @test used_by_constraint(y0)
        @test used_by_constraint(vref)
        empty!(InfiniteOpt._constraint_dependencies(vref))
    end
    # test used_by_objective
    @testset "used_by_objective" begin
        @test !used_by_objective(y0)
        @test !used_by_objective(vref)
        InfiniteOpt._data_object(vref).in_objective = true
        @test used_by_objective(vref)
        InfiniteOpt._data_object(vref).in_objective = false
    end
    # test is_used
    @testset "is_used" begin
        # test not used
        @test !is_used(vref)
        # test used by constraint and/or measure
        push!(InfiniteOpt._constraint_dependencies(vref), InfOptConstraintIndex(1))
        @test is_used(y0)
        empty!(InfiniteOpt._constraint_dependencies(vref))
        # test used by objective
        InfiniteOpt._data_object(vref).in_objective = true
        @test is_used(vref)
    end
end

# Test Deletion helpers
@testset "Deletion Helpers" begin
    m = InfiniteModel()
    @variable(m, 1 <= z <= 2)
    @test InfiniteOpt._delete_variable_dependencies(dispatch_variable_ref(z)) isa Nothing
    @test !has_lower_bound(z)
    @test !has_upper_bound(z)
end

# test deprecations
@testset "Parameter Bound Deprecations" begin
    # setup
    m = InfiniteModel()
    @infinite_parameter(m, t in [0, 1])
    @variable(m, z)
    rs = DomainRestrictions(t => 0)
    @test_throws ErrorException has_parameter_bounds(z)
    @test_throws ErrorException parameter_bounds(z)
    @test_throws ErrorException set_parameter_bounds(z, rs)
    @test_throws ErrorException add_parameter_bounds(z, rs)
    @test_throws ErrorException delete_parameter_bounds(z)
end


text:
##interlace block operators
function isboundaryrow(A,k)
    for j=1:size(A,2)
        if isafunctional(A[k,j])
            return true
        end
    end

    return false
end

domainscompatible(A::AbstractMatrix{T}) where {T<:Operator} = domainscompatible(map(domainspace,A))

function spacescompatible(A::AbstractMatrix{T}) where T<:Operator
    for k=1:size(A,1)
        if !spacescompatible(map(rangespace,A[k,:]))
            return false
        end
    end
    for k=1:size(A,2)
        if !spacescompatible(map(domainspace,A[:,k]))
            return false
        end
    end
    true
end

spacescompatible(A::AbstractVector{T}) where {T<:Operator} = spacescompatible(map(domainspace,A))

function domainspace(A::AbstractMatrix{T}) where T<:Operator
    if !spacescompatible(A)
        error("Cannot construct domainspace for $A as spaces are not compatible")
    end

    spl=map(domainspace,A[1,:])
    Space(spl)
end

function rangespace(A::AbstractVector{T}) where T<:Operator
    if !spacescompatible(A)
        error("Cannot construct rangespace for $A as domain spaces are not compatible")
    end

    spl=map(rangespace,A)
    Space(spl)
end

promotespaces(A::AbstractMatrix{T}) where {T<:Operator} = promotespaces(Matrix(A))

function promotespaces(A::Matrix{T}) where T<:Operator
    isempty(A) && return A
    ret = similar(A) #TODO: promote might have different Array type
    for j=1:size(A,2)
        ret[:,j] = promotedomainspace(A[:,j])
    end
    for k=1:size(A,1)
        ret[k,:] = promoterangespace(ret[k,:])
    end

    # do a second loop as spaces might have been inferred
    # during range space
    for j=1:size(A,2)
        ret[:,j] = promotedomainspace(ret[:,j])
    end
    ret
end

## Interlace operator

struct InterlaceOperator{T,p,DS,RS,DI,RI,BI} <: Operator{T}
    ops::Array{Operator{T},p}
    domainspace::DS
    rangespace::RS
    domaininterlacer::DI
    rangeinterlacer::RI
    bandwidths::BI
end

const VectorInterlaceOperator = InterlaceOperator{T,1,DS,RS} where {T,DS,RS<:Space{D,R}} where {D,R<:AbstractVector}
const MatrixInterlaceOperator = InterlaceOperator{T,2,DS,RS} where {T,DS,RS<:Space{D,R}} where {D,R<:AbstractVector}

InterlaceOperator(ops::Array{T,p},ds,rs,di,ri,bi) where {T,p} =
    InterlaceOperator{T,p,typeof(ds),typeof(rs),
                        typeof(di),typeof(ri),typeof(bi)}(ops,ds,rs,di,ri,bi)

function InterlaceOperator(ops::AbstractMatrix{Operator{T}},ds::Space,rs::Space) where T
    # calculate bandwidths TODO: generalize
    p=size(ops,1)
    dsi = interlacer(ds)
    rsi = interlacer(rs)

    if size(ops,2) == p && all(isbanded,ops) &&# only support nblocks 1 for now
            all(i->isa(i,AbstractFill) && getindex_value(i) == 1, dsi.blocks) &&
            all(i->isa(i,AbstractFill) && getindex_value(i) == 1, rsi.blocks)
        l,u = 0,0
        for k=1:p,j=1:p
            l=max(l,p*bandwidth(ops[k,j],1)+k-j)
        end
        for k=1:p,j=1:p
            u=max(u,p*bandwidth(ops[k,j],2)+j-k)
        end
    elseif p == 1 && size(ops,2) == 2 && size(ops[1],2) == 1
        # special case for example
        l,u = max(bandwidth(ops[1],1),bandwidth(ops[2],1)-1),bandwidth(ops[2],2)+1
    else
        l,u = (1-dimension(rs),dimension(ds)-1)  # not banded
    end

    InterlaceOperator(ops,ds,rs,
                        cache(dsi),
                        cache(rsi),
                        (l,u))
end

function InterlaceOperator(ops::Vector{Operator{T}},ds::Space,rs::Space) where T
    # calculate bandwidths
    p=size(ops,1)
    if all(isbanded,ops)
        l,u = 0,0
        #TODO: this code assumes an interlace strategy that might not be right
        for k=1:p
            l=max(l,p*bandwidth(ops[k],1)+k-1)
        end
        for k=1:p
            u=max(u,p*bandwidth(ops[k],2)+1-k)
        end
    else
        l,u = (1-dimension(rs),dimension(ds)-1)  # not banded
    end

    InterlaceOperator(ops,ds,rs,
                        cache(BlockInterlacer(tuple(blocklengths(ds)))),
                        cache(interlacer(rs)),
                        (l,u))
end

function InterlaceOperator(opsin::AbstractMatrix{Operator{T}}) where {T}
    isempty(opsin) && throw(ArgumentError("Cannot create InterlaceOperator from empty Matrix"))
    ops=promotespaces(opsin)
    InterlaceOperator(ops,domainspace(ops),rangespace(ops[:,1]))
end

function InterlaceOperator(opsin::AbstractMatrix{Operator{T}},::Type{DS},::Type{RS}) where {T,DS<:Space,RS<:Space}
    isempty(opsin) && throw(ArgumentError("Cannot create InterlaceOperator from empty Matrix"))
    ops=promotespaces(opsin)
    InterlaceOperator(ops,DS(components(domainspace(ops))),RS(rangespace(ops[:,1]).spaces))
end

function InterlaceOperator(opsin::RowVector{Operator{T}}) where {T}
    isempty(opsin) && throw(ArgumentError("Cannot create InterlaceOperator from empty Matrix"))

    ops=promotespaces(opsin)
    InterlaceOperator(ops,domainspace(ops),rangespace(ops[1]))
end

function InterlaceOperator(opsin::RowVector{Operator{T}},::Type{DS},::Type{RS}) where {T,DS<:Space,RS<:Space}
    isempty(opsin) && throw(ArgumentError("Cannot create InterlaceOperator from empty Matrix"))
    ops=promotespaces(opsin)
    InterlaceOperator(ops,DS(components(domainspace(ops))),rangespace(ops[1]))
end

InterlaceOperator(opsin::AbstractMatrix{Operator{T}},::Type{DS}) where {T,DS<:Space} =
    InterlaceOperator(opsin,DS,DS)

InterlaceOperator(opsin::AbstractMatrix,S...) =
    InterlaceOperator(Matrix{Operator{mapreduce(eltype,promote_type,opsin)}}(promotespaces(opsin)),S...)

function InterlaceOperator(opsin::Vector{Operator{T}}) where T
    ops=promotedomainspace(opsin)
    InterlaceOperator(ops,domainspace(first(ops)),rangespace(ops))
end

InterlaceOperator(ops::AbstractArray{T,p}) where {T,p} =
    InterlaceOperator(Array{Operator{mapreduce(eltype,promote_type,ops)},p}(ops))

function convert(::Type{Operator{T}},S::InterlaceOperator) where T
    if T == eltype(S)
        S
    else
        ops=Array{Operator{T}}(undef, size(S.ops)...)
        for j=1:size(S.ops,2),k=1:size(S.ops,1)
            ops[k,j]=S.ops[k,j]
        end
        InterlaceOperator(ops,domainspace(S),rangespace(S),
                            S.domaininterlacer,S.rangeinterlacer,S.bandwidths)
    end
end

#TODO: More efficient to save bandwidth
bandwidths(M::InterlaceOperator) = M.bandwidths

blockbandwidths(M::InterlaceOperator) =
    (mapreduce(op->blockbandwidth(op,1),max,M.ops),
     mapreduce(op->blockbandwidth(op,2),max,M.ops))

isblockbanded(M::InterlaceOperator) = all(isblockbanded,M.ops)

function blockcolstop(M::InterlaceOperator,J::Integer)
    if isblockbandedbelow(M)
        Block(J + blockbandwidth(M,1))
    else
        mapreduce(op->blockcolstop(op,J),max,M.ops)
    end
end

function colstop(M::InterlaceOperator{T},j::Integer) where T
#    b=bandwidth(M,1)
    if isbandedbelow(M)
        min(j+bandwidth(M,1)::Int,size(M,1))::Int
    elseif isblockbandedbelow(M)
        J=block(domainspace(M), j)::Block{1}
        blockstop(rangespace(M), blockcolstop(M,J)::Block{1})::Int
    else #assume is raggedbelow
        K = 0
        (J,jj) = M.domaininterlacer[j]
        for N = 1:size(M.ops,1)
            cs = colstop(M.ops[N,J],jj)::Int
            if cs > 0
                K = max(K,findfirst(M.rangeinterlacer,(N,cs))::Int)
            end
        end
        K
    end
end

israggedbelow(M::InterlaceOperator) = all(israggedbelow,M.ops)

getindex(op::InterlaceOperator,k::Integer,j::Integer) =
    error("Higher tensor InterlaceOperators not supported")

function getindex(op::InterlaceOperator{T,2},k::Integer,j::Integer) where {T}
    M,J = op.domaininterlacer[j]
    N,K = op.rangeinterlacer[k]
    op.ops[N,M][K,J]::T
end

# the domain is not interlaced
function getindex(op::InterlaceOperator{T,1},k::Integer,j::Integer) where T
    N,K = op.rangeinterlacer[k]
    op.ops[N][K,j]::T
end

function getindex(op::InterlaceOperator{T},k::Integer) where T
    if size(op,1) == 1
        op[1,k]
    elseif size(op,2) == 1
        op[k,1]
    else
        error("Only implemented for row/column operators.")
    end
end

findsub(cr,ν) = findall(x->x[1]==ν,cr)

function getindex(L::InterlaceOperator{T},kr::UnitRange) where T
    ret=zeros(T,length(kr))

    if size(L,1) == 1
        ds=domainspace(L)
        cr=cache(interlacer(ds))[kr]
    elseif size(L,2) == 1
        rs=rangespace(L)
        cr=cache(interlacer(rs))[kr]
    else
        error("Only implemented for row/column operators.")
    end

    for ν=1:length(L.ops)
        # indicies of ret
        ret_kr=findsub(cr,ν)

        # block indices
        if !isempty(ret_kr)
            sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]

            LinearAlgebra.axpy!(1.0,L.ops[ν][sub_kr],view(ret,ret_kr))
        end
    end
    ret
end

# overwritten for functions
# this won't work in 0.4 as expected, though the user
# should call vec anyways for 0.5 compatibility
function getindex(L::InterlaceOperator,k::Integer,j)
    if k==1 && size(L,1) == 1
        L[j]
    else
        defaultgetindex(L,k,j)
    end
end

function getindex(L::InterlaceOperator,k,j::Integer)
    if j==1 && size(L,2) == 1
        L[k]
    else
        defaultgetindex(L,k,j)
    end
end

#####
# optimized copy routine for when there is a single domainspace
# and no interlacing of the columns is necessary
# this is especially important for \
######
for TYP in (:BandedMatrix, :BlockBandedMatrix, :BandedBlockBandedMatrix, :RaggedMatrix,
                :Matrix)
    @eval begin
        function $TYP(S::SubOperator{T,InterlaceOperator{T,1,SS,PS,DI,RI,BI},
                            Tuple{UnitRange{Int},UnitRange{Int}}}) where {SS,PS,DI,RI,BI,T}
            kr,jr=parentindices(S)
            L=parent(S)

            ret=$TYP(Zeros, S)

            ds=domainspace(L)
            rs=rangespace(L)
            cr=cache(interlacer(rs))[kr]
            for ν=1:length(L.ops)
                # indicies of ret
                ret_kr=findsub(cr,ν)

                # block indices
                if !isempty(ret_kr)
                    sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]

                    LinearAlgebra.axpy!(1.0,view(L.ops[ν],sub_kr,jr),view(ret,ret_kr,:))
                end
            end
            ret
        end

        function $TYP(S::SubOperator{T,InterlaceOperator{T,2,SS,PS,DI,RI,BI},
                      Tuple{UnitRange{Int},UnitRange{Int}}}) where {SS,PS,DI,RI,BI,T}
            kr,jr=parentindices(S)
            L=parent(S)

            ret=$TYP(Zeros, S)

            if isempty(kr) || isempty(jr)
                return ret
            end

            ds=domainspace(L)
            rs=rangespace(L)
            cr=L.rangeinterlacer[kr]
            cd=L.domaininterlacer[jr]
            for ν=1:size(L.ops,1),μ=1:size(L.ops,2)
                # indicies of ret
                ret_kr=findsub(cr,ν)
                ret_jr=findsub(cd,μ)

                # block indices
                if !isempty(ret_kr) && !isempty(ret_jr)
                    sub_kr=cr[ret_kr[1]][2]:cr[ret_kr[end]][2]
                    sub_jr=cd[ret_jr[1]][2]:cd[ret_jr[end]][2]

                    LinearAlgebra.axpy!(1.0,view(L.ops[ν,μ],sub_kr,sub_jr),
                                   view(ret,ret_kr,ret_jr))
                end
            end
            ret
        end
    end
end

## Build block-by-block
function blockbanded_interlace_convert!(S,ret)
    T = eltype(S)
    KR,JR = parentindices(S)
    l,u=blockbandwidths(S)::Tuple{Int,Int}

    M = map(op -> begin
                KR_size = Block.(Int(first(KR)):min(Int(last(KR)),nblocks(op,1)))
                JR_size = Block.(Int(first(JR)):min(Int(last(JR)),nblocks(op,2)))
                BlockBandedMatrix(view(op, KR_size, JR_size))
            end, parent(S).ops)

    for J=Block(1):Block(nblocks(ret,2)),K=blockcolrange(ret,Int(J))
        Bs=view(ret,K,J)
        j = 0
        for ξ=1:size(M,2)
            k = 0
            m = 0
            for κ=1:size(M,1)
                if K.n[1] ≤ nblocks(M[κ,ξ],1) && J.n[1] ≤ nblocks(M[κ,ξ],2)
                    MKJ = M[κ,ξ][K,J]::Matrix{T}
                    n,m = size(MKJ)
                    Bs[k+1:k+n,j+1:j+m] = MKJ
                    k += n
                end
            end
            j += m
        end
    end
    ret
end

for d in (:1,:2)
    @eval BlockBandedMatrix(S::SubOperator{T,InterlaceOperator{T,$d,SS,PS,DI,RI,BI},
                          Tuple{BlockRange1,BlockRange1}}) where {SS,PS,DI,RI,BI,T} =
    blockbanded_interlace_convert!(S, BlockBandedMatrix(Zeros, S))
end

domainspace(IO::InterlaceOperator) = IO.domainspace
rangespace(IO::InterlaceOperator) = IO.rangespace

#tests whether an operator can be made into a column
iscolop(op) = isconstop(op)
iscolop(::Multiplication) = true

promotedomainspace(A::InterlaceOperator{T,1},sp::Space) where {T} =
    InterlaceOperator(map(op->promotedomainspace(op,sp),A.ops))

interlace(A::AbstractArray{T}) where {T<:Operator} = InterlaceOperator(A)

const OperatorTypes = Union{Operator,Fun,Number,UniformScaling}

operators(A::InterlaceOperator) = A.ops
operators(A::Matrix{<:Operator}) = A
operators(A::Operator) = [A]

Base.vcat(A::MatrixInterlaceOperator...) =
    InterlaceOperator(vcat(map(operators,A)...))
function _vcat(A::OperatorTypes...)
    Av = Vector{Operator{mapreduce(eltype,promote_type,A)}}()
    for a in A
        if a isa VectorInterlaceOperator
            append!(Av,a.ops)
        else
            push!(Av,a)
        end
    end
    InterlaceOperator(vnocat(Av...))
end

Base.hcat(A::Union{VectorInterlaceOperator,MatrixInterlaceOperator}...) =
    InterlaceOperator(hcat(map(A->A.ops,A)...))
_hcat(A::OperatorTypes...) = InterlaceOperator(hnocat(A...))
function _hvcat(rows::Tuple{Vararg{Int}},as::OperatorTypes...)
    # Based on Base
    nbr = length(rows)  # number of block rows
    rs = Array{Any,1}(undef, nbr)
    a = 1
    for i = 1:nbr
        rs[i] = hcat(map(op -> convert(Operator,op),as[a:a-1+rows[i]])...)
        a += rows[i]
    end
    vcat(rs...)
end

Base.vcat(A::Operator, B::OperatorTypes...) = _vcat(A, B...)
Base.hcat(A::Operator, B::OperatorTypes...) = _hcat(A, B...)
Base.hvcat(rows::Tuple{Vararg{Int}}, A::Operator, B::OperatorTypes...) =
    _hvcat(rows, A, B...)

Base.vcat(C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) = _vcat(C, A, B...)
Base.hcat(C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) = _hcat(C, A, B...)
Base.hvcat(rows::Tuple{Vararg{Int}}, C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) =
    _hvcat(rows, C, A, B...)

Base.vcat(D::Union{Fun,Number,UniformScaling}, C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) = _vcat(D, C, A, B...)
Base.hcat(D::Union{Fun,Number,UniformScaling}, C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) = _hcat(D, C, A, B...)
Base.hvcat(rows::Tuple{Vararg{Int}}, D::Union{Fun,Number,UniformScaling}, C::Union{Fun,Number,UniformScaling}, A::Operator, B::OperatorTypes...) =
    _hvcat(rows, D, C, A, B...)

## Convert Matrix operator to operators

Operator(M::AbstractArray{OO}) where {OO<:Operator} = InterlaceOperator(M)

function interlace_choosedomainspace(ops,sp::UnsetSpace)
    # this ensures correct dispatch for unino
    sps = Vector{Space}(
        filter(x->!isambiguous(x),map(choosedomainspace,ops)))
    if isempty(sps)
        UnsetSpace()
    else
        union(sps...)
    end
end

function interlace_choosedomainspace(ops,rs::Space)
    # this ensures correct dispatch for unino
    sps = Vector{Space}(
        filter(x->!isambiguous(x),map((op)->choosedomainspace(op,rs),ops)))
    if isempty(sps)
        UnsetSpace()
    else
        union(sps...)
    end
end

choosedomainspace(A::InterlaceOperator{T,1},rs::Space) where {T} =
    interlace_choosedomainspace(A.ops,rs)

choosedomainspace(A::InterlaceOperator{T,2},rs::Space) where {T} =
    Space([interlace_choosedomainspace(A.ops[:,k],rs) for k=1:size(A.ops,2)])


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

mutable struct SlowActor
    t::Union{Nothing,Task}
    sticky::Bool
    act::_ACT
end

struct Mailbox
    box::Deque
    lock::ReentrantLock
    A::Union{Nothing,SlowActor}
end

"""
    newLink()

Create a mailbox (without an actor) which can be used to
communicate with slow actors.
"""
newLink() = Link(Mailbox(Deque{Any}(), ReentrantLock(), nothing),
                 myid(), :local)

Base.push!(mbx::Mailbox, items...) = lock(mbx.lock) do
    push!(mbx.box, items...)
end
Base.popfirst!(mbx::Mailbox) = lock(mbx.lock) do
    popfirst!(mbx.box)
end
Base.isempty(mbx::Mailbox) = lock(mbx.lock) do
    isempty(mbx.box)
end


text:
Lynx.init("SnakeGame", 800, 600)

@setup function()
    time = 0
    points = 0

    colors = Settings.colors
    rez = Settings.resolution

    w = floor(Int, width / rez)
    h = floor(Int, height / rez)

    snake = Snake()
    food = Food()

    spawn!(food)
    framerate!(Settings.framerate)
end

onkeypress(@window) do event
    key = keyname(event) |> uppercase
    dir = nothing;
    if key in KeyMap.LEFT
        dir = Dir.LEFT
    elseif key in KeyMap.RIGHT
        dir = Dir.RIGHT
    elseif key in KeyMap.DOWN
        dir = Dir.DOWN
    elseif key in KeyMap.UP
        dir = Dir.UP
    end
    if !isnothing(dir) && dir != -snake.dir
        setdir!(snake, dir)
    end
end

function gameover!()
    background(colors.gameover)
    origin()
    fontsize(28)
    sethue("white")
    text("Game Over", halign = :center)
end

function gameinfo()
    origin(0, 0)
    fontsize(16)
    t = floor(Int, time)
    text("Time: $(t)s", 40, 40)
    text("Score: $points", 120, 40)
end

function update(dt::Real)
    background(colors.bg)
    origin(0, 0)
    gameinfo()
    scale(rez)

    update!(snake)
    draw(food)
    draw(snake)
    if endgame(snake)
        gameover!()
        gameinfo()
        @info "Game over"
        return loop!(false)
    end

    if eat(snake, food)
        grow(snake)
        spawn!(food)
        global points += 1
    end

    global time += dt
end

max_stars_repo_name: jorge-brito/LynxGallery

max_issues_repo_name: jorge-brito/LynxGallery

max_forks_repo_name: jorge-brito/LynxGallery


text:
### SELECTED DATA FILE INDEXES REQUIRED FOR TGFb MODEL ONLY

# INDEX VALUES CHANGE WITH UPDATES TO MODEL.NET

initial_cond_update_array = [];
rate_constant_update_array = [];
control_parameter_update_array = [];

# initial conditions
#push!(initial_cond_update_array,[tgfb12_initial_cond_index,0.0]) # TGFB12 initial conditions

# This gets hacky, we are using the TGFB infastructure to run the inhibitors etc.

# rates
# TGFB dosing
push!(rate_constant_update_array,[yreg_translation_rate_index,10.0]) # addition_TGFb12: [] =([])=> TGFB12
push!(rate_constant_update_array,[yreg_deg_rate_index,yreg_ox_deg_rate]) # addition_TGFb12: [] =([])=> TGFB12

# Run snail slug instead
#push!(rate_constant_update_array,[152,10.0]) # induction_gene_SNAIL_SLUG: gene_SNAIL_SLUG = mRNA_SNAIL_SLUG
#push!(rate_constant_update_array,[ss_virus_forward_rate_index,4.0])
#push!(rate_constant_update_array,[ss_virus_reverse_rate_index,0.1])

# give MAPK inhibitor

#push!(rate_constant_update_array,[U0126_foward_rate_index,4.0])
#push!(rate_constant_update_array,[U0126_reverse_rate_index,0.1])

# ECADHERINchange this to transcription rate
# push!(rate_constant_update_array,[ecad_trans_rate_index,ecad_trans_rate*5.0])

# Block TGFB3
#push!(rate_constant_update_array,[TGFB3_Ab_foward_rate_index,4.0])
#push!(rate_constant_update_array,[TGFB3_Ab_reverse_rate_index,0.1])

# TGFB3 dosing

#push!(rate_constant_update_array,[TGFB3_foward_rate_index,4.0])
#push!(rate_constant_update_array,[TGFB3_reverse_rate_index,0.1])

# DN LEF1 temoprary
#lef1_induciton_rate_index = 141
#lef1_induciton_rate = 0.01

# DN SMAD temoprary
#SMAD_induciton_rate_index = 134
#SMAD_induciton_rate = 0.01

# push modifications to RAM
#UpdateArray(data_dictionary,"INITIAL_CONDITION_ARRAY",initial_cond_update_array)
#UpdateArray(data_dictionary,"RATE_CONSTANT_ARRAY",rate_constant_update_array)

# mRNA degradation test
#tubulin_induction_index = 145
#push!(rate_constant_update_array,[tubulin_induction_index,0.0])


text:
using MPITestImages
using Documenter

DocMeta.setdocmeta!(MPITestImages, :DocTestSetup, :(using MPITestImages); recursive=true)

makedocs(;
    modules=[MPITestImages],
    authors="Jonas Schumacher <jonas.schumacher@imte.fraunhofer.de> and contributors",
    repo="https://github.com/MagneticParticleImaging/MPITestImages.jl/blob/{commit}{path}#{line}",
    sitename="MPITestImages.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://github.com/MagneticParticleImaging/MPITestImages.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(repo   = "github.com/MagneticParticleImaging/MPITestImages.jl.git",
           target = "build")


text:
# ---
# title: 515. Find Largest Value in Each Tree Row
# id: problem515
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/find-largest-value-in-each-tree-row/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _an array of the largest value in
# each row_ of the tree **(0-indexed)**.
#
#
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/21/largest_e1.jpg)
#
#
#
#     Input: root = [1,3,2,5,3,null,9]
#     Output: [1,3,9]
#
#
# **Example 2:**
#
#
#
#     Input: root = [1,2,3]
#     Output: [1,3]
#
#
# **Example 3:**
#
#
#
#     Input: root = [1]
#     Output: [1]
#
#
# **Example 4:**
#
#
#
#     Input: root = [1,null,2]
#     Output: [1,2]
#
#
# **Example 5:**
#
#
#
#     Input: root = []
#     Output: []
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree will be in the range `[0, 104]`.
#   * `-231 <= Node.val <= 231 - 1`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Point2DTools

using StaticArrays

export Point2D, convex, convexperm, centroid, bounds, scalefns

const Point2D = SVector{2}

Base.atan(p::Point2D) = atan(p.y, p.x)
Base.atan(p1::Point2D, p2::Point2D) = atan(p2-p1)

dtan(p::Point2D) = rad2deg(a)

"""
    convex(points)

Return the list of Point2D points which form the convex hull of the supplied points.
"""
convex(points) = points[convexperm(points)]

×(p1::Point2D, p2::Point2D) = p1.x * p2.y - p2.x * p1.y

direction(p1, p2, p3) = ×(p3-p1, p2-p1)

right(p1, p2, p3) = direction(p1, p2, p3) > 0
left(p1, p2, p3) = direction(p1, p2, p3) < 0
collinear(p1, p2, p3) = direction(p1, p2, p3) == 0 # is 0 too large ?
"""
    convexperm(points)

 Return the list of indexes of the Point2D points which form the convex hull of the supplied points.
"""
function convexperm(points)
    ips = sortperm(points)

    Lupper = [ips[1], ips[2]]

    for i in 3:length(ips)
        push!(Lupper, ips[i])
        while length(Lupper) > 2 && !right(points[Lupper[end-2:end]]...)
            Lupper[end-1] = Lupper[end]
            pop!(Lupper)
        end
    end
    Lower = [ips[end], ips[end-1]]
    for i in length(ips)-2:-1:1
        push!(Lower, ips[i])
        while length(Lower) > 2 && !right(points[Lower[end-2:end]]...)
            Lower[end-1] = Lower[end]
            pop!(Lower)
        end
    end
    vcat(Lupper, Lower[2:end])
end

"""
    centroid(points)
Calculate the unweighted centre of the closed polygon defined by `points`.
# Arguments
- `points` a Vector of points, assumed to include the "start" and "end" point twice, so points[1] is ignored
"""
centroid(points) = length(points) < 2 ? points[1] : sum(points[2:end]) / length(points[2:end])

"""
    bounds(points)
Return the [Minimum::Point2D, Maximum::Point2D] (so you can vcat with a list of lists)
# Arguments
- `points` Vector of points with .x and .y properties
# Examples
julia> b1 = bounds([Point2D(-10,0), Point2D(1,2), Point2D(2,1)])
2-element Vector{StaticArrays.SVector{2, Int64}}:
 [-10, 0]
 [2, 2]
julia> b2 = bounds([Point2D(0,0), Point2D(1,2), Point2D(2,11)])
2-element Vector{StaticArrays.SVector{2, Int64}}:
 [0, 0]
 [2, 11]
 julia> bounds(vcat(b1, b2))
 2-element Vector{StaticArrays.SVector{2, Int64}}:
  [-10, 0]
  [2, 11]
"""
function bounds(points)
    minx, miny, maxx, maxy = Inf, Inf, -Inf, -Inf
    for p in points
        if p[1] < minx
            minx = p[1]
        elseif p[1] > maxx
            maxx = p[1]
        end
        if p[2] < miny
            miny = p[2]
        elseif p.y > maxy
            maxy = p[2]
        end
    end
    [Point2D(minx, miny), Point2D(maxx, maxy)]
end

"""
    scalefns(points, width, height)
    scalefns(min_max, width, height)
return functions to apply to fx and fy of the points or bounds to scale to specified `width` and `height` - preserving aspect
# Arguments
- `points` Vector of points with .x and .y properties
"""
scalefns(points, width, height) = scalefns(bounds(points)..., width, height)
function scalefns(minxy, maxxy, width, height)
    scale = min(width, height) / min(maxxy.x - minxy.x, maxxy.y - minxy.y)
    x -> scale * (x - minxy.x), y -> scale * (y - minxy.y)
end

###
end


text:
# Autogenerated wrapper script for cif_api_jll for armv7l-linux-musleabihf-cxx11
export libcif

using ICU_jll
using SQLite_jll
JLLWrappers.@generate_wrapper_header("cif_api")
JLLWrappers.@declare_library_product(libcif, "libcif.so.1")
function __init__()
    JLLWrappers.@generate_init_header(ICU_jll, SQLite_jll)
    JLLWrappers.@init_library_product(
        libcif,
        "lib/libcif.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using VUMPS
using VUMPS: fixedpoint
using Test

struct StopFunction{T,S}
    o::Ref{T}
    tol::S
end

mutable struct StopFunction2
    counter::Int
    limit::Int
end

@testset "fixedpoint" begin
    # example squareroot
    (st::StopFunction)(v) = abs(v - st.o[]) < st.tol ? true : (st.o[] = v; false)

    init = 9
    next(guess, n) = 1/2*(guess + n/guess)
    stopfun = StopFunction(Ref(Inf), 1e-9)

    #do nothing
    @test fixedpoint(x->next(x, 9), 9, x -> true) ≈ 9
    #evaluate once
    (st::StopFunction2)(v) = (st.counter += 1; st.counter == st.limit)
    @test fixedpoint(x->next(x, 9), 9, StopFunction2(0,2)) ≈ 1/2*(9 + 1)
end


text:
using Pkg
Pkg.activate(".")

using BenchmarkTools, LinearAlgebra

function dot1(x, y)
    s = 0.0
    for i in 1:length(x)
        s += x[i]*y[i]
    end
    s
end

function dot2(x, y)
    s = 0.0
    @simd for i in 1:length(x)
        @inbounds s += x[i]*y[i]
    end
    s
end

x = 1000*rand(10000);
y = 1000*rand(10000);

res1 = @btime dot1($x, $y)
res2 = @btime dot2($x, $y)
res3 = @btime dot($x, $y)

println([res2, res3] .- res1)


text:
function encodeMTF(str::AbstractString, symtable::Vector{Char}=collect('a':'z'))
    function encode(ch::Char)
        r = findfirst(symtable, ch)
        deleteat!(symtable, r)
        prepend!(symtable, ch)
        return r
    end
    collect(encode(ch) for ch in str)
end

function decodeMTF(arr::Vector{Int}, symtable::Vector{Char}=collect('a':'z'))
    function decode(i::Int)
        r = symtable[i]
        deleteat!(symtable, i)
        prepend!(symtable, r)
        return r
    end
    join(decode(i) for i in arr)
end

testset = ["broood", "bananaaa", "hiphophiphop"]
encoded = encodeMTF.(testset)
decoded = decodeMTF.(encoded)
for (str, enc, dec) in zip(testset, encoded, decoded)
    println("Test string: $str\n -> Encoded: $enc\n -> Decoded: $dec")
end

using Base.Test
@testset "Decoded string equal to original" begin
    for (str, dec) in zip(testset, decoded)
        @test str == dec
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Distributions

m = constructdice()
run(m)

mcs = @defsim begin
    # The parameters of the log‐normal distribution fit to Olsen et al.
    # are μ = 1.10704 and σ = 0.264 (Gillingham et al. 2015)
    t2xco2 = LogNormal(1.10704, 0.264)

    save(damages.DAMAGES)
end

res = run(mcs, m, 1000; trials_output_filename="/tmp/dice-2013/trialdata.csv", results_output_dir="/tmp/dice-2013")


text:
@inline function compose_axes(::Tuple{}, x::AbstractArray{<:Any,N}) where {N}
    if N === 0
        return ()
    elseif N === 1 && can_change_size(x)
        return (compose_axis(OneToMRange(length(x))),)
    else
        return map(compose_axis, axes(x))
    end
end
function compose_axes(ks::Tuple{Vararg{<:Any,N}}, x::AbstractArray{<:Any,N}) where {N}
    if N === 0
        return ()
    elseif N === 1 && can_change_size(x)
        return compose_axes(ks, (OneToMRange(length(x)),))
    else
        return compose_axes(ks, axes(x))
    end
end
function compose_axes(ks::Tuple, x::AbstractArray{<:Any,N}) where {N}
    throw(DimensionMismatch("Number of axis arguments provided ($(length(ks))) does " *
                            "not match number of parent axes ($N)."))
end
@inline function compose_axes(ks::Tuple{Vararg{<:Any,N}}, inds::Tuple{Vararg{<:Any,N}}) where {N}
    return (
        compose_axis(first(ks), first(inds)),
        compose_axes(tail(ks), tail(inds))...
    )
end
compose_axes(::Tuple{}, ::Tuple{}) = ()
compose_axes(::Tuple{}, inds::Tuple) = map(compose_axis, inds)
compose_axes(axs::Tuple, ::Tuple{}) = map(compose_axis, axs)

###
### compose_axis
###
compose_axis(x::Integer) = SimpleAxis(x)
compose_axis(x) = Axis(x)
compose_axis(x::AbstractAxis) = x
function compose_axis(x::AbstractUnitRange{I}) where {I<:Integer}
    if known_first(x) === one(eltype(x))
        return SimpleAxis(x)
    else
        return OffsetAxis(x)
    end
end
compose_axis(x::IdentityUnitRange) = compose_axis(x.indices)

# 3-args
compose_axis(::Nothing, inds) = compose_axis(inds)
compose_axis(ks::Function, inds) = ks(inds)
function compose_axis(ks::Integer, inds)
    if ks isa StaticInt
        return SimpleAxis(known_first(inds):ks)
    else
        return SimpleAxis(inds)
    end
end
function compose_axis(ks, inds)
    check_axis_length(ks, inds)
    return _compose_axis(ks, inds)
end
function _compose_axis(ks::AbstractAxis, inds)
    # if the indices are the same then don't reconstruct
    if first(parent(ks)) == first(inds)
        return copy(ks)
    else
        return unsafe_reconstruct(ks, inds)
    end
end
@inline function _compose_axis(ks, inds)
    start = known_first(ks)
    if known_step(ks) === 1
        if known_first(ks) === nothing
            return OffsetAxis(first(ks) - static_first(inds), inds)
        elseif known_first(ks) === known_first(inds)
            # if we don't know the length of `inds` but we know the length of `ks` then we
            # should reconstruct `inds` so that it has a static length
            if known_last(inds) === nothing && known_last(ks) !== nothing
                return set_length(inds, static_length(ks))
            else
                return copy(inds)
            end
        else
            return OffsetAxis(static_first(ks) - static_first(inds), inds)
        end
    else
        check_unique_keys(ks)
        return initialize_axis(ks, inds)
    end
end

"""
    AxisArray{T,N,P,AI}

An array struct that wraps any parent array and assigns it an `AbstractAxis` for
each dimension. The first argument is the parent array and the second argument is
a tuple of subtypes to `AbstractAxis` or keys that will be converted to subtypes
of `AbstractAxis` with the provided keys.
"""
struct AxisArray{T,N,D,Axs<:Tuple{Vararg{<:Any,N}}} <: AbstractArray{T,N}
    data::D
    axes::Axs

    # TODO robust checking of indices should happen at this level
    # FIXME this needs to check that all axs are AbstractAxis
    function AxisArray{T,N,P,A}(p::P, axs::A) where {T,N,P,A}
        for i in OneTo(N)
            check_axis_length(axs[i], axes(p, i))
        end
        return new{T,N,P,A}(p, axs)
    end

    ###
    ### AxisArray{T,N,P}
    ###
    function AxisArray{T,N,P}(x::P, axs::Tuple) where {T,N,P<:AbstractArray{T,N}}
        axs = compose_axes(axs, x)
        return new{T,N,P,typeof(axs)}(x, axs)
    end

    AxisArray{T,N,P}(A::AxisArray{T,N,P}) where {T,N,P} = A

    AxisArray{T,N,P}(A::AbstractArray, args...) where {T,N,P} = AxisArray{T,N,P}(A, args)

    function AxisArray{T,N,P}(A::AxisArray; kwargs...) where {T,N,P}
        return AxisArray{T,N,P}(convert(P, parent(A)), axes(A); kwargs...)
    end

    function AxisArray{T,N,P}(x::AbstractArray, axs::Tuple) where {T,N,P}
        return AxisArray{T,N,P}(convert(P, x), axs)
    end

    # TODO fix/clean up these docs
    """
        AxisArray{T,N}(undef, dims::NTuple{N,Integer})
        AxisArray{T,N}(undef, keys::NTuple{N,AbstractVector})

    Construct an uninitialized `N`-dimensional array containing elements of type `T` were
    the size of each dimension is equal to the corresponding integer in `dims`.

    Construct an uninitialized `N`-dimensional array containing elements of type `T` were
    the size of each dimension is determined by the length of the corresponding collection
    in `keys.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> size(AxisArray{Int,2}(undef, (2,2)))
    (2, 2)

    julia> size(AxisArray{Int,2}(undef, (["a", "b"], [:one, :two])))
    (2, 2)

    """
    function AxisArray{T,N}(A::AbstractArray, ks::Tuple) where {T,N}
        if eltype(A) <: T
            axs = compose_axes(ks, A)
            return new{T,N,typeof(A),typeof(axs)}(p, axs)
        else
            p = AbstractArray{T}(A)
            axs = compose_axes(ks, p)
            return new{T,N,typeof(p),typeof(axs)}(p, axs)
        end
    end
    function AxisArray{T,N}(x::AbstractArray{T,N}, axs::Tuple) where {T,N}
        axs = compose_axes(axs, x)
        return new{T,N,typeof(x),typeof(axs)}(x, axs)
    end
    function AxisArray{T,N}(A::AxisArray, ks::Tuple) where {T,N}
        if eltype(A) <: T
            axs = compose_axes(ks, A)
            return new{T,N,parent_type(A),typeof(axs)}(p, axs)
        else
            p = AbstractArray{T}(parent(A))
            axs = compose_axes(ks, A)
            return new{T,N,typeof(p),typeof(axs)}(p, axs)
        end
    end
    function AxisArray{T,N}(init::ArrayInitializer, args...; kwargs...) where {T,N}
        return AxisArray{T,N}(init, args; kwargs...)
    end
    AxisArray{T,N}(x::AbstractArray, args...) where {T,N} = AxisArray{T,N}(x, args)
    function AxisArray{T,N}(init::ArrayInitializer, ks::Tuple{Vararg{<:Any,N}}) where {T,N}
        axs = map(compose_axis, ks)
        p = init_array(T, init, axs)
        return new{T,N,typeof(p),typeof(axs)}(p, axs)
    end

    ### AxisArray{T}
    """
        AxisArray{T}(undef, keys::NTuple{N,AbstractVector})

    Construct an uninitialized `N`-dimensional array containing elements of type `T` were
    the size of each dimension is determined by the length of the corresponding collection
    in `keys.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> size(AxisArray{Int}(undef, (["a", "b"], [:one, :two])))
    (2, 2)
    ```
    """
    function AxisArray{T}(x::AbstractArray, axs::Tuple) where {T}
        return AxisArray{T,ndims(x)}(x, axs)
    end
    function AxisArray{T}(x::AbstractArray, axs::Vararg) where {T}
        return AxisArray{T,ndims(x)}(x, axs)
    end
    function AxisArray{T}(init::ArrayInitializer, axs::Tuple) where {T}
        return AxisArray{T,length(axs)}(init, axs)
    end
    function AxisArray{T}(init::ArrayInitializer, axs::Vararg) where {T}
        return AxisArray{T,length(axs)}(init, axs)
    end

    """
        AxisArray(parent::AbstractArray, axes::Tuple)

    Construct an `AxisArray` using `parent` and explicit subtypes of `AbstractAxis`.
    If `check_length` is `true` then each dimension of parent's length is checked to match
    the length of the corresponding axis (e.g., `size(parent 1) == length(axes[1])`.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> AxisArray(ones(2,2), (SimpleAxis(2), SimpleAxis(2)))
    2×2 AxisArray(::Array{Float64,2}
      • axes:
         1 = 1:2
         2 = 1:2
    )
         1    2
      1  1.0  1.0
      2  1.0  1.0

    julia> AxisArray(ones(2,2), (["a", "b"], ["one", "two"]))
    2×2 AxisArray(::Array{Float64,2}
      • axes:
         1 = ["a", "b"]
         2 = ["one", "two"]
    )
           "one"   "two"
      "a"  1.0     1.0
      "b"  1.0     1.0

    ```
    """
    function AxisArray(x::AbstractArray{T,N}, ks::Tuple) where {T,N}
        axs = compose_axes(ks, x)
        return new{T,N,typeof(x),typeof(axs)}(x, axs)
    end

    """
        AxisArray(parent::AbstractArray, args...) -> AxisArray(parent, tuple(args))

    Passes `args` to a tuple for constructing an `AxisArray`.

    ## Examples
    ```jldoctest
    julia> using AxisIndices

    julia> A = AxisArray(reshape(1:9, 3,3), 2:4, 3.0:5.0);

    julia> A[1, 1]
    1

    julia> A[==(2), ==(3.0)]
    1

    julia> A[1:2, 1:2] == [1 4; 2 5]
    true

    julia> A[<(4), <(5.0)] == [1 4; 2 5]
    true
    ```
    """
    AxisArray(x::AbstractArray, args...) = AxisArray(x, args)

    #= TODO delete this?
    function AxisArray(x::AbstractVector{T}; kwargs...) where {T}
        if can_change_size(x)
            axs = (SimpleAxis(OneToMRange(length(x))),)
        else
            axs = (SimpleAxis(axes(x, 1)),)
        end
        return new{T,1,typeof(x),typeof(axs)}(x, axs)
    end

    function AxisArray(x::AbstractVector{T}, ks::AbstractAxis) where {T}
        check_axis_length(ks, axes(x, 1))
        return new{T,1,typeof(x),Tuple{typeof(ks)}}(x, (ks,))
    end

    function AxisArray(x::AbstractVector{T}, ks::AbstractVector) where {T}
        if can_change_size(x)
            axs = (Axis(ks, OneToMRange(length(x))),)
        else
            axs = (Axis(ks, axes(x, 1)),)
        end
        return new{T,1,typeof(x),typeof(axs)}(x, axs)
    end

    AxisArray(x::AbstractVector{T}, ks::Tuple{}; kwargs...) where {T} = AxisArray(x)

    AxisArray(x::AbstractVector{T}, ks::Tuple; kwargs...) where {T} = AxisArray{T}(x, ks)

    function AxisArray(x::AbstractArray{T,0}, axs::Tuple{}=()) where {T}
        return new{T,0,typeof(x),Tuple{}}(x, ())
    end

    function AxisArray(x::AbstractVector{T}, ks::Tuple) where {T}
        return new{T,1,typeof(x),typeof(axs)}(x, axs)
    end
    =#
end

function initialize_axis_array(data, axs)
    return unsafe_initialize(
        AxisArray{eltype(data),ndims(data),typeof(data),typeof(axs)},
        (data, axs)
    )
end

Base.axes(x::AxisArray) = getfield(x, :axes)

Base.parent(x::AxisArray) = getfield(x, :data)

ArrayInterface.parent_type(::Type{T}) where {P,T<:AxisArray{<:Any,<:Any,P}} = P
@inline function ArrayInterface.can_change_size(::Type{T}) where {D,Axs,T<:AxisArray{<:Any,<:Any,D,Axs}}
    if can_change_size(D)
        return _can_change_axes_size(Axs)
    else
        return false
    end
end
ArrayInterface.can_setindex(::Type{T}) where {T<:AxisArray} = can_setindex(parent_type(T))

@generated function _can_change_axes_size(::Type{T}) where {T<:Tuple}
    for i in T.parameters
        can_change_size(i) && return true
    end
    return false
end

Base.parentindices(x::AxisArray) = parentindices(parent(x))

Base.length(x::AxisArray) = prod(size(x))
@generated function ArrayInterface.known_length(::Type{T}) where {Axs,T<:AxisArray{<:Any,<:Any,<:Any,Axs}}
    out = 1
    for axis in Axs.parameters
        known_length(axis) === nothing && return nothing
        out = out * known_length(axis)
    end
    return out
end

Base.size(x::AxisArray) = map(static_length, axes(x))

@inline function Base.axes(x::AxisArray, i::Integer)
    if i < 1
        error("BoundsError: attempt to access $(typeof(x)) at dimension $i")
    else
        return unsafe_axes(axes(x), i)
    end
end

@inline unsafe_axes(axs::Tuple, ::StaticInt{I}) where {I} = unsafe_axes(axs, I)
@inline function unsafe_axes(axs::Tuple, i)
    if i > length(axs)
        return SimpleAxis(1)
    else
        return getfield(axs, i)
    end
end

Base.eachindex(A::AxisArray) = eachindex(IndexStyle(A), A)

Base.eachindex(::IndexCartesian, A::AxisArray{T,N}) where {T,N} = CartesianIndices(axes(A))
function Base.eachindex(S::IndexLinear, A::AxisArray{T,N}) where {T,N}
    if N === 1
        return axes(A, 1)
    else
        return compose_axis(eachindex(S, parent(A)))
    end
end

Base.eachindex(A::AxisArray{T,1}) where {T} = axes(A, 1)

function ArrayInterface.unsafe_reconstruct(A::AxisArray, data; axes=nothing, kwargs...)
    return _unsafe_reconstruct(A, data, axes)
end

# TODO function _unsafe_reconstruct(A, data, ::Nothing) end
_unsafe_reconstruct(A, data, axs) = initialize_axis_array(data, axs)

###
### getindex
###
@inline function ArrayInterface.unsafe_get_element(A::AxisArray, inds)
    if is_dense_wrapper(A)
        return @inbounds(parent(A)[apply_offsets(A, inds)...])
    else
        return _unsafe_get_element(A, apply_offsets(A, inds))
    end
end
# other methods in padded_axis.jl
@inline _unsafe_get_element(A, inds::Tuple{Vararg{Integer}}) = @inbounds(parent(A)[inds...])

function ArrayInterface.unsafe_get_collection(A::AxisArray, inds)
    axs = to_axes(A, inds)
    dest = AxisArray(similar(parent(A), length.(axs)), axs)
    if map(length, axes(dest)) == map(length, axs)
        Base._unsafe_getindex!(dest, A, inds...) # usually a generated function, don't allow it to impact inference result
    else
        Base.throw_checksize_error(dest, axs)
    end
    return dest
end

function ArrayInterface.unsafe_set_element!(A::AxisArray, value, inds)
    return @inbounds(setindex!(parent(A), value, apply_offsets(A, inds)...))
end

# * apply_offsets helps apply offsets for offset axes and if there is a padded region
#   applies the padding.
# * If padding returns a function instead of redirecting to a another index, then we
#   filter out the function and apply it to the array to grab the value that represents
#   padded regions
@inline function apply_offsets(A::AxisArray, inds::Tuple)
    if length(inds) === 1 && ndims(A) > 1 # linear indexing
        return (_sub_offset(eachindex(A), first(inds)),)
    else
        return apply_offsets(axes(A), inds)
    end

    #=
    if new_inds isa Tuple{Vararg{<:Integer}}
        return new_inds
    else
        return _filter_pad(new_inds)
    end
    =#
end

_filter_pad(inds::Tuple{I,Vararg}) where {I<:Function} = first(inds)
@inline _filter_pad(inds::Tuple{<:Integer,Vararg}) = _filter_pad(tail(inds))

###
### SubArray
###
#=
@propagate_inbounds function Base.getindex(A::SubArray{T,N,<:AxisArray{T,N}}, args...) where {T,N}
    return ArrayInterface.getindex(A, args...)
end
@propagate_inbounds function Base.getindex(A::SubArray{T,N,<:AxisArray{T,N}}, args::Vararg{Int,N}) where {T,N}
    return ArrayInterface.getindex(A, args...)
end

function ArrayInterface.unsafe_set_element!(A::AxisArray, value, inds)
    return @inbounds(setindex!(parent(A), value, apply_offsets(A, inds)...))
end

=#
@inline function ArrayInterface.to_axes(A::AxisArray, inds::Tuple)
    if ndims(A) === 1
        return (to_axis(axes(A, 1), first(inds)),)
    elseif ArrayInterface.is_linear_indexing(A, inds)
        return (to_axis(eachindex(IndexLinear(), A), first(inds)),)
    else
        return to_axes(A, axes(A), inds)
    end
end

"""
    is_dense_wrapper(::Type{T}) where {T} -> Bool

Do all the indices of `T` map to a unique indice of the parent data that is wrapped?
This is not true for padded axes.
"""
is_dense_wrapper(x) = is_dense_wrapper(typeof(x))
is_dense_wrapper(::Type{T}) where {T} = true
@generated function is_dense_wrapper(::Type{T}) where {Axs,T<:AxisArray{<:Any,<:Any,<:Any,Axs}}
    for i in Axs.parameters
        is_dense_wrapper(i) || return false
    end
    return true
end

@inline function Base.IndexStyle(::Type{A}) where {A<:AxisArray}
    if is_dense_wrapper(A)
        return IndexStyle(parent_type(A))
    else
        return IndexCartesian()
    end
end


text:
abstract type MagnusPropagator end

export MagnusPropagator