text:
using Neighborhood
using Distances

export beta_statistic
export mdop_embedding
export mdop_maximum_delay

"""
    mdop_embedding(s::Vector; kwargs...) → Y, τ_vals, ts_vals, FNNs, βS
MDOP (for "maximizing derivatives on projection") is a unified approach to properly
embed a timeseries or a set of timeseries (`Dataset`)
based on the paper of Chetan Nichkawde [^Nichkawde2013].

## Keyword arguments

* `τs= 0:50`: Possible delay values `τs`. For each of
  the `τs`'s the β-statistic gets computed.
* `w::Int = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `fnn_thres::Real= 0.05`: A threshold value defining a sufficiently small fraction
  of false nearest neighbors, in order to the let algorithm terminate and stop
  the embedding procedure (`0 ≤ fnn_thres < 1).
* `r::Real = 2`: The threshold for the tolerable relative increase of the distance
  between the nearest neighbors, when increasing the embedding dimension.
* `max_num_of_cycles = 50`: The algorithm will stop after that many cycles no matter what.

## Description
The method works iteratively and gradually builds the final embedding `Y`.
Based on the [`beta_statistic`](@ref) the algorithm picks an optimal delay
value `τ` for each embedding cycle as the global maximum of `β`. In case of
multivariate embedding, i.e. when embedding a set of time series (`s::Dataset`),
the optimal delay value `τ` is chosen as the maximum from all maxima's of all
considered `β`-statistics for each possible timeseries. The range of
considered delay values is determined in `τs` and for the nearest neighbor
search we respect the Theiler window `w`.

After each embedding cycle the FNN-statistic `FNNs` [^Hegger1999][^Kennel1992] is being
checked and as soon as this statistic drops below the threshold `fnn_thres`,
the algorithm terminates. In order to increase the practability of the method the
algorithm also terminates when the FNN-statistic `FNNs` increases.

The final embedding is returned as `Y`. The chosen delay values for each embedding
cycle are stored in the `τ_vals` and the according timeseries index chosen for the
the respective according delay value in `τ_vals` is stored in `ts_vals`.
`βS, FNNs` are returned for clarity and double-checking, since they are computed anyway.
In case of multivariate embedding, `βS` will store all `β`-statistics for all
available time series in each embedding cycle. To double-check the actual used
`β`-statistics in an embedding cycle 'k', simply `βS[k][:,ts_vals[k+1]]`.

[^Nichkawde2013]: Nichkawde, Chetan (2013). [Optimal state-space reconstruction using derivatives on projected manifold. Physical Review E 87, 022905](https://doi.org/10.1103/PhysRevE.87.022905).
[^Hegger1999]: Hegger, Rainer and Kantz, Holger (1999). [Improved false nearest neighbor method to detect determinism in time series data. Physical Review E 60, 4970](https://doi.org/10.1103/PhysRevE.60.4970).
[^Kennel1992]: Kennel, M. B., Brown, R., Abarbanel, H. D. I. (1992). [Determining embedding dimension for state-space reconstruction using a geometrical construction. Phys. Rev. A 45, 3403] (https://doi.org/10.1103/PhysRevA.45.3403).
"""
function mdop_embedding(s::Vector{T};
        τs = 0:50 , w::Int = 1, fnn_thres::Real = 0.05,
        max_num_of_cycles = 50,
        r::Real = 2) where {T<:Real}

    @assert 0 ≤ fnn_thres < 1 "Please select a valid breaking criterion, `fnn_thres` ∈ [0 1)"
    @assert all(x -> x ≥ 0, τs)
    metric = Euclidean()
    s_orig = s
    s = standardize(s) # especially important for fnn-computation
    # define actual phase space trajectory
    Y_act = Dataset(s)

    # compute nearest neighbor distances in the first embedding dimension for
    # FNN statistic
    vtree = KDTree(Y_act, metric)
    allNNidxs, NNdist_old = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    τ_vals = Int64[0]
    ts_vals = Int64[1]
    FNNs = Float64[]
    βS = Array{T}(undef, length(τs), max_num_of_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act, NNdist_new = mdop_embedding_cycle!(
            Y_act, flag, s, τs, w, counter, βS, τ_vals, metric, r,
            FNNs, fnn_thres, ts_vals, NNdist_old
        )
        flag = mdop_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
        counter += 1
        # for the next embedding cycle save the distances of NN for this dimension
        NNdist_old = NNdist_new
    end
    # construct final reconstruction vector
    Y_final = s_orig
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,s_orig,τ_vals[i])
    end
    return Y_final, τ_vals[1:counter-1], ts_vals[1:counter-1], FNNs, βS[:,1:counter-1]
end

function mdop_embedding(Y::AbstractDataset{D, T};
        τs = 0:50 , w::Int = 1, fnn_thres::Real = 0.05,
        max_num_of_cycles = 50,
        r::Real = 2) where {D, T<:Real}

    @assert 0 ≤ fnn_thres < 1 "Please select a valid breaking criterion, `fnn_thres` ∈ [0 1)"
    @assert all(x -> x ≥ 0, τs)
    metric = Euclidean()
    Y_orig = Y
    Y = standardize(Y) # especially important for fnn-computation
    # define actual phase space trajectory and NN-distances
    Y_act = []
    NNdist_old = []

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    τ_vals = Int64[0]
    ts_vals = Int64[]
    FNNs = Float64[]
    βS = fill(zeros(T, length(τs), size(Y,2)), 1, max_num_of_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        Y_act, NNdist_new = mdop_multivariate_embedding_cycle!(
            Y_act, flag, Y, τs, w, counter, βS, τ_vals, metric, r,
            FNNs, fnn_thres, ts_vals, NNdist_old
        )
        flag = mdop_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
        counter += 1
        # for the next embedding cycle save the distances of NN for this dimension
        NNdist_old = NNdist_new
    end
    # construct final reconstruction vector
    Y_final = Y_orig[:,ts_vals[1]]
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,Y_orig[:,ts_vals[i]],τ_vals[i])
    end
    return Y_final, τ_vals[1:counter-1], ts_vals[1:counter-1], FNNs, βS[:,1:counter-1]
end

# Here we separate the inner core loop of the mdop_embedding into another function
# not only for clarity of source code but also to introduce a function barrier
# for the type instability of `Y → Y_act`
function mdop_embedding_cycle!(
        Y, flag, s, τs, w, counter, βS, τ_vals, metric, r,
        FNNs, fnn_thres, ts_vals, NNdist_old
    )

    β = beta_statistic(Y, s, τs, w)
    βS[:, counter] = β

    # determine the optimal tau value from the β-statistic
    maxi, max_idx = findmax(β)
    # store chosen delay (and chosen time series)
    push!(τ_vals, τs[max_idx])
    push!(ts_vals, 1)

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y,s,τ_vals[counter+1])
    vtree = KDTree(Y_act, metric)
    allNNidxs, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r)
    push!(FNNs,fnns)
    return Y_act, NNdist_new
end

function mdop_multivariate_embedding_cycle!(
        Y_act, flag, Ys, τs, w, counter, βS, τ_vals, metric, r,
        FNNs, fnn_thres, ts_vals, NNdist_old
    )

    M = size(Ys,2)
    # in the 1st cycle we have to check all (size(Y,2)^2 combinations
    if counter == 1
        Y_act, NNdist_new = first_embedding_cycle_MDOP!(M, Ys, τs, w, τ_vals,
                                                ts_vals, βS, metric, FNNs, r)

    # in all other cycles we just have to check (size(Y,2)) combinations
    else

        Y_act, NNdist_new = embedding_cycle_MDOP!(Y_act, counter, NNdist_old, M, Ys, τs, w,
                                                τ_vals, ts_vals, βS, metric,
                                                FNNs, r)
    end
    return Y_act, NNdist_new
end

"""
Computes the actual reconstruction vectors `Y_act` and the according nearest
neighbour distances `NNdist_new` for the very first embedding cycle of the
multivariate MDOP-embedding procedure.
"""
function first_embedding_cycle_MDOP!(M, Ys, τs, w, τ_vals,
                                        ts_vals, βS, metric, FNNs, r)
    counter = 1
    βs = zeros(length(τs), M*M)
    beta_counter = 1
    for ts = 1:M
        for ts2 = 1:M
            βs[:,beta_counter] = beta_statistic(Dataset(Ys[:,ts]), Ys[:,ts2], τs, w)
            beta_counter += 1
        end
    end
    # determine the optimal tau value from the β-statistic
    τ_opt, ts_opt1, ts_opt2 = choose_optimal_tau1(βs, M)
    # store chosen delay, time series and according β-statistic
    push!(τ_vals, τs[τ_opt])
    push!(ts_vals, ts_opt1)
    push!(ts_vals, ts_opt2)
    βS[counter] = βs[:,(ts_vals[counter]-1)*M+1:(ts_vals[counter]-1)*M+M]

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Ys[:,ts_vals[counter]],
                                Ys[:,ts_vals[counter+1]],τ_vals[counter+1])

    # compute nearest neighbor distances in the first embedding dimension for
    # FNN statistic
    Y_old = Dataset(Ys[:,ts_vals[counter]])
    vtree = KDTree(Y_old, metric)
    _, NNdist_old = all_neighbors(vtree, Y_old, 1:length(Y_old), 1, w)

    vtree = KDTree(Y_act, metric)
    _, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r)
    push!(FNNs,fnns)

    return Y_act, NNdist_new
end

"""
Computes the actual reconstruction vectors `Y_act` and the according nearest
neighbour distances `NNdist_new` for every embedding cycle of the multivariate
MDOP-embedding procedure, but the first one.
"""
function embedding_cycle_MDOP!(Y_act, counter, NNdist_old, M, Ys, τs, w, τ_vals,
                                        ts_vals, βS, metric, FNNs, r)
    βs = zeros(length(τs), M)
    for ts = 1:M
        βs[:,ts] = beta_statistic(Y_act, Ys[:,ts], τs, w)
    end

    # determine the optimal tau value from the β-statistic
    τ_opt, ts_opt = choose_optimal_tau2(βs)
    # store chosen delay, time series and according β-statistic
    push!(τ_vals, τs[τ_opt])
    push!(ts_vals, ts_opt)
    βS[counter] = βs

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act,
                                Ys[:,ts_vals[counter+1]],τ_vals[counter+1])
    vtree = KDTree(Y_act, metric)
    _, NNdist_new = all_neighbors(vtree, Y_act, 1:length(Y_act), 1, w)

    # compute FNN-statistic and store vals, while also
    # truncating distance-vector to the "new" length of the actual embedding vector
    fnns = fnn_embedding_cycle(view(NNdist_old, 1:length(Y_act)), NNdist_new, r)
    push!(FNNs,fnns)

    return Y_act, NNdist_new
end

"""
This function chooses the optimal τ value as the maximum of the β-statistics
maxima. It returns the index of this maximum as well as the time series number
corresponding to that maximum of maxima. In this first embedding cycle it also
returns the time series number, which act as Y_act.
"""
function choose_optimal_tau1(βs::Array{T, 2}, M::Int) where {T}
    NN = size(βs,2)
    maxis = zeros(T, NN)
    max_idx = zeros(Int, NN)
    for i = 1:NN
        # determine optimal tau value and maximum vals from all β-statistic's
        maxis[i], max_idx[i] = findmax(βs[:,i])
    end
    # take the maximum of the maximums
    _, idx = findmax(maxis)
    # determine the first time series and the according second one
    a = idx/M
    ts_1 = Int(ceil(a))
    ts_2 = idx-M*(ts_1-1)
    return max_idx[idx], ts_1, ts_2
end

"""
This function chooses the optimal τ value as the maximum of the β-statistics
maxima. It returns the index of this maximum as well as the time series number
corresponding to that maximum of maxima.
"""
function choose_optimal_tau2(βs::Array{T, 2}) where {T}
    NN = size(βs,2)
    maxis = zeros(T, NN)
    max_idx = zeros(Int, NN)
    for i = 1:NN
        # determine optimal tau value and maximum vals from all β-statistic's
        maxis[i], max_idx[i] = findmax(βs[:,i])
    end
    # take the maximum of the maximums
    _, idx = findmax(maxis)

    return max_idx[idx], idx
end

function mdop_break_criterion(FNNs, counter, fnn_thres, max_num_of_cycles)
    flag = true
    if FNNs[counter] ≤ fnn_thres
        flag = false
        println("Algorithm stopped due to sufficiently small FNNs. "*
                "Valid embedding achieved ✓.")
    end
    if counter > 1 && FNNs[counter] > FNNs[counter-1]
        flag = false
        println("Algorithm stopped due to rising FNNs. "*
                "Valid embedding achieved ✓.")
    end
    if max_num_of_cycles == counter
        println("Algorithm stopped due to hitting max cycle number. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    return flag
end

"""
    beta_statistic(Y::Dataset, s::Vector) [, τs, w]) → β
Compute the β-statistic `β` for input state space trajectory `Y` and a timeseries `s`
according to Nichkawde [^Nichkawde2013], based on estimating derivatives on a
projected manifold. For a range of delay values `τs`, `β` gets computed and its
maximum over all considered `τs` serves as the optimal delay considered in this
embedding cycle.

Arguments `τs, w` as in [`mdop_embedding`](@ref).

## Description
The `β`-statistic is based on the geometrical idea of maximal unfolding of the
reconstructed attractor and is tightly related to the False Nearest Neighbor
method ([^Kennel1992]). In fact the method eliminates the maximum amount of
false nearest neighbors in each embedding cycle. The idea is to estimate the
absolute value of the directional derivative with respect to a possible new
dimension in the reconstruction process, and with respect to the nearest
neighbor, for all points of the state space trajectory:

ϕ'(τ) = Δϕ_d(τ) / Δx_d

Δx_d is simply the Euclidean nearest neighbor distance for a reference point
with respect to the given Theiler window `w`. Δϕ_d(τ) is the distance of the
reference point to its nearest neighbor in the one dimensional time series `s`,
for the specific τ. Δϕ_d(τ) = |s(i+τ)-s(j+τ)|, with i being the index of the
considered reference point and j the index of its nearest neighbor.

Finally,

`β` = log β(τ) = ⟨log₁₀ ϕ'(τ)⟩ ,

with ⟨.⟩ being the mean over all reference points. When one chooses the maximum
of `β` over all considered τ's, one obtains the optimal delay value for this
embedding cycle. Note that in the first embedding cycle, the input state space
trajectory `Y` can also be just a univariate time series.

[^Nichkawde2013]: Nichkawde, Chetan (2013). [Optimal state-space reconstruction using derivatives on projected manifold. Physical Review E 87, 022905](https://doi.org/10.1103/PhysRevE.87.022905).
[^Kennel1992]: Kennel, M. B., Brown, R., Abarbanel, H. D. I. (1992). [Determining embedding dimension for state-space reconstruction using a geometrical construction. Phys. Rev. A 45, 3403] (https://doi.org/10.1103/PhysRevA.45.3403).
"""
function beta_statistic(Y::AbstractDataset{D,T}, s::Vector{T}, τs = 0:50, w::Int = 1) where {D, T<:Real}
    @assert length(s) ≥ length(Y) "The length of the input time series `s` must be at least the length of the input trajectory `Y` "
    @assert all(x -> x ≥ 0, τs)
    τ_max = maximum(τs)
    metric = Euclidean()    # consider only Euclidean norm
    K = 1                   # consider only first nearest neighbor
    N = length(Y)           # length of the state space trajectory
    NN = N - τ_max          # allowed length of the trajectory w.r.t. τ_max
    vtree = KDTree(Y[1:NN], metric)
    allNNidxs, Δx = all_neighbors(vtree, Y[1:NN], 1:NN, K, w)   # Eq. 12

    Δϕ = zeros(NN, length(τs))
    @inbounds for j = 1:NN # loop over all state space points
        for (i,τ) in enumerate(τs) # loop over all considered τ's
            Δϕ[j,i] = abs(s[j+τ][1]-s[allNNidxs[j][1]+τ][1]) / Δx[j][1] # Eq. 14 & 15
        end
    end

    β = mean(log10.(Δϕ); dims=1) # β statistic is the average of the logarithms, see Eq. 16
    return vec(β)
end

"""
    mdop_maximum_delay(s, tw = 1:50, samplesize = 1.0)) -> τ_max, L
Compute an upper bound for the search of optimal delays, when using `mdop_embedding`
[`mdop_embedding`](@ref) or `beta_statistic` [`beta_statistic`](@ref).

## Description
The input time series `s` gets embedded with unit lag and increasing dimension,
for dimensions (or time windows) `tw` (`RangeObject`).
For each of such a time window the `L`-statistic from Uzal et al. [^Uzal2011]
will be computed. `samplesize` determines the fraction of points to be
considered in the computation of `L` (see [`uzal_cost`](@ref)). When this
statistic reaches its global minimum the maximum delay value `τ_max` gets
returned. When `s` is a multivariate `Dataset`, `τ_max` will becomputed for all
timeseries of that Dataset and the maximum value will be returned. The returned
`L`-statistic has size `(length(tw), size(s,2))`.

[^Nichkawde2013]: Nichkawde, Chetan (2013). [Optimal state-space reconstruction using derivatives on projected manifold. Physical Review E 87, 022905](https://doi.org/10.1103/PhysRevE.87.022905).
[^Uzal2011]: Uzal, L. C., Grinblat, G. L., Verdes, P. F. (2011). [Optimal reconstruction of dynamical systems: A noise amplification approach. Physical Review E 84, 016223](https://doi.org/10.1103/PhysRevE.84.016223).
"""
function mdop_maximum_delay(s::Vector{T}, tw=1:50, samplesize::Real=1) where {T<:Real}
    @assert all(x -> x ≥ 0, tw)
    L = zeros(T, length(tw))
    counter = 1
    for i in tw
        i==1 ? Y_act = Dataset(s) : Y_act = embed(s,i,1)
        L[counter] = uzal_cost(Y_act; samplesize = samplesize)
        counter +=1
    end
    _, τ_max = findmin(L)
    return tw[τ_max], L
end

function mdop_maximum_delay(s::AbstractDataset{D,T}, tw=1:50, samplesize::Real=1) where {D, T<:Real}
    τs_max = zeros(Int,D)
    Ls = zeros(T, length(tw), D)
    @inbounds for i = 1:D
        τs_max[i], Ls[:, i] = mdop_maximum_delay(vec(s[:,i]), tw, samplesize)
    end
    return maximum(τs_max), Ls
end


text:
function evalNaive(diag::Diagrams{W}, evalPropagator, varK, varT, root = nothing, phase = nothing, para = nothing) where {W}
    momenta = diag.momenta
    propagators = diag.propagators
    tree = diag.tree

    # calculate new momentum
    for mom in momenta
        mom.curr = varK[1] * mom.basis[1]
        for i = 2:length(mom.basis)
            if mom.basis[i] != 0
                mom.curr += varK[i] * mom.basis[i]
            end
        end
    end

    #calculate propagators
    for p in propagators
        K = momenta[p.Kidx].curr
        p.curr = evalPropagator(p.type, K, p.Tidx, varT, p.factor, para)
    end

    #calculate diagram tree
    for node in tree
        if node.operation == 1 #multiply
            node.curr = 1.0
            for pidx in node.propagators
                node.curr *= propagators[pidx].curr
            end
            for nidx in node.nodes
                node.curr *= tree[nidx].curr
            end

        elseif node.operation == 2 #sum
            node.curr = 0.0
            for pidx in node.propagators
                node.curr += propagators[pidx].curr
            end
            for nidx in node.nodes
                node.curr += tree[nidx].curr
            end
        else
            error("not implemented")
        end

        node.curr *= node.factor
        if isnothing(phase) == false && (isempty(node.extK) == false || isempty(node.extT) == false)
            node.curr *= phase(varK, varT, node.extK, node.extT)
            # println(node.id, ": ", node.curr)
        end
    end

    if isnothing(root) == false
        return [r == -1 ? W(0) : tree[r].curr for r in root]
    else
        return tree[end].curr
    end
end

ext: jl


text:
function test_construction(AT)
    @testset "construction" begin
        constructors(AT)
        conversion(AT)
        value_constructor(AT)
        iterator_constructors(AT)
    end
end

function constructors(AT)
    @testset "constructors + similar" begin
        for T in supported_eltypes()
            B = AT{T}(10)
            @test B isa AT{T,1}
            @test size(B) == (10,)
            @test eltype(B) == T

            B = AT{T}(10, 10)
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = AT{T}((10, 10))
            @test B isa AT{T,2}
            @test size(B) == (10, 10)
            @test eltype(B) == T

            B = similar(B, Int32, (11, 15))
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(B, T)
            @test B isa AT{T,2}
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(B, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(B, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T

            B = similar(AT{Int32}, (11, 15))
            @test B isa AT{Int32,2}
            @test size(B) == (11, 15)
            @test eltype(B) == Int32

            B = similar(AT{Int32, 2}, T, (11, 15))
            @test B isa AT{T,2}
            @test size(B) == (11, 15)
            @test eltype(B) == T

            B = similar(AT{T}, (5,))
            @test B isa AT{T,1}
            @test size(B) == (5,)
            @test eltype(B) == T

            B = similar(AT{T}, 7)
            @test B isa AT{T,1}
            @test size(B) == (7,)
            @test eltype(B) == T
        end
    end
end

function conversion(AT)
    @testset "conversion" begin
        for T in supported_eltypes()
            Bc = round.(rand(10, 10) .* 10.0)
            B = AT{T}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10)
            B = AT(Bc)
            @test size(B) == (10,)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(T, 10, 10)
            B = AT{T, 2}(Bc)
            @test size(B) == (10, 10)
            @test eltype(B) == T
            @test Array(B) ≈ Bc

            Bc = rand(Int32, 3, 3, 3)
            B = convert(AT{T, 3}, Bc)
            @test size(B) == (3, 3, 3)
            @test eltype(B) == T
            @test Array(B) ≈ Bc
        end
    end
end

function value_constructor(AT)
    @testset "value constructors" begin
        for T in supported_eltypes()
            x = fill(zero(T), (2, 2))
            x1 = fill(AT{T}, T(0), (2, 2))
            x2 = fill(AT{T}, T(0), (2, 2))
            x3 = fill(AT{T, 2}, T(0), (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            x = fill(T(1), (2, 2))
            x1 = fill(AT{T}, T(1), (2, 2))
            x2 = fill(AT{T}, T(1), (2, 2))
            x3 = fill(AT{T, 2}, T(1), (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            fill!(x, 0)

            x1 = fill(AT{T}, 0f0, (2, 2))
            x2 = fill(AT{T}, T(0), (2, 2))
            x3 = fill(AT{T, 2}, T(0), (2, 2))
            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x

            fill!(x1, 2f0)
            x2 = fill!(AT{Int32}((4, 4, 4)), 77f0)
            @test all(x-> x == 2f0, Array(x1))
            @test all(x-> x == Int32(77), Array(x2))

            x = Matrix{T}(I, 2, 2)

            x1 = AT{T, 2}(I, 2, 2)
            x2 = AT{T}(I, (2, 2))
            x3 = AT{T, 2}(I, (2, 2))

            @test Array(x1) ≈ x
            @test Array(x2) ≈ x
            @test Array(x3) ≈ x
        end
    end
end
function iterator_constructors(AT)
    @testset "iterator constructors" begin
        for T in supported_eltypes()
            @test AT(Fill(T(0), (10,))) == fill(AT{T}, T(0), (10,))
            @test AT(Fill(T(0), (10, 10))) == fill(AT{T}, T(0), (10, 10))
            if T <: Real
                x = AT{Float32}(Fill(T(0), (10, 10)))
                @test eltype(x) == Float32
                @test AT(Eye{T}((10, 10))) == AT{T}(I, 10, 10)
                x = AT{Float32}(Eye{T}((10, 10)))
                @test eltype(x) == Float32
            end
        end
    end
end


text:
export getfirstnotes, allowedpitches, separatepitches, firstnotes

"""
    firstnotes(notes, grid)

Return the notes that first appear in each grid point, *without quantizing them*.

This function *does not* consider the notes modulo the quarter note! Different
quarter notes have different grid points.
"""
function firstnotes(notes::Notes{N}, grid) where {N<:AbstractNote}
    isgrid(grid)

    clas = classify(notes, grid)

    # Notes with dif == 0 belong to the same grid bin
    dif = clas[2:end] - clas[1:end-1]

    posdif = [Int(notes[j].position - notes[j-1].position) for j in 2:length(notes)]
    # minpos is the minimum position difference notes can have while being in the
    # same quarter note
    tpq = notes.tpq
    minpos = tpq - tpq*(grid[end] - grid[end-1] + grid[2] - grid[1])/2

    # The first notes are notes that have dif == 0 OR posdif > minpos
    toadd = Vector{Note}()
    for j in 1:length(notes)-1
        if (dif[j] == 0 && posdif[j] > minpos) || dif[j] != 0
            push!(toadd, notes[j+1])
        end
    end

    # Take care of first note in notes:
    if clas[1] != clas[2] || notes[2].position - notes[1].position > minpos
        pushfirst!(toadd, notes[1])
    elseif notes[2].position < notes[1].position
        toadd[1] = notes[1]
    end
    return Notes(toadd, notes.tpq)
end

"""
    allowedpitches(notes::Notes, allowed) -> newnotes

Only keep the notes that have a pitch specified in `allowed` (one or many pitches).
"""
function allowedpitches(notes::Notes{N}, allowed) where {N<:AbstractNote}
    n = N[]
    for i ∈ 1:length(notes)
        notes[i].pitch ∈ allowed && push!(n, deepcopy(notes[i]))
    end
    return Notes(n, notes.tpq)
end

"""
    separatepitches(notes::Notes [, allowed])

Get a dictionary \"pitch\"=>\"notes of that pitch\".
Optionally only keep pitches that are contained in `allowed`.
"""
function separatepitches(notes::Notes{N}) where {N}
    separated = Dict{UInt8, Notes{N}}()
    for note in notes
        _add_note_to_dict!(separated, note, notes.tpq)
    end
    return separated
end

function separatepitches(notes::Notes{N}, pitches) where {N}
    separated = Dict{UInt8, Notes{N}}()
    for note in notes
        note.pitch ∈ pitches && _add_note_to_dict!(separated, note, notes.tpq)
    end
    return separated
end

function _add_note_to_dict!(separated, note::N, tpq) where {N<:AbstractNote}
    if haskey(separated, note.pitch)
        push!(separated[note.pitch], deepcopy(note))
    else
        separated[note.pitch] = Notes{N}(Vector{N}[], tpq)
        push!(separated[note.pitch], deepcopy(note))
    end
end


text:
module QuadricMeshSimplification

# Mesh Simplification Tutorial
#
# (C) by Sven Forstmann in 2014
# (C) by Steve Kelly in 2020
#
# License : MIT
# http://opensource.org/licenses/MIT
#
# Derived from:
# https://github.com/sp4cerat/Fast-Quadric-Mesh-Simplification
#
# 5/2016: Chris Rorden created minimal version for OSX/Linux/Windows compile
# 4/2020: Steve Kelly Julia port

using StaticArrays
using LinearAlgebra

export simplify_mesh!, simplify_mesh_lossless!

function SymetricMatrix(a, b, c, d)
    SVector{10,Float64}(a*a, a*b, a*c, a*d, b*b, b*c, b*d, c*c, c*d, d*d)
end

function det(m::SVector, a11, a12, a13, a21, a22, a23, a31, a32, a33)
    m[a11]*m[a22]*m[a33] + m[a13]*m[a21]*m[a32] + m[a12]*m[a23]*m[a31]
                - m[a13]*m[a22]*m[a31] - m[a11]*m[a23]*m[a32]- m[a12]*m[a21]*m[a33];
end

mutable struct Triangle
    v::SVector{3,Int}
    err::SVector{4,Float64}
    deleted::Bool
    dirty::Bool
    n::SVector{3,Float64}
end

mutable struct Vertex{PT}
    p::PT
    tstart::Int
    tcount::Int
    q::SVector{10,Float64}
    border::Bool
end

struct TRef
    tid::Int
    tvertex::Int
end

# Helper functions

"""
Main simplification function

 target_count  : target nr. of triangles
 agressiveness : sharpness to increase the threshold.
                 5..8 are good numbers
                 more iterations yield higher quality

"""
function simplify_mesh!(invts, infcs, target_count::Int, agressiveness=7, verbose=false)
    # init
    triangles = Vector{Triangle}(undef,length(infcs))
    vertices = Vector{Vertex{eltype(invts)}}(undef,length(invts))

    for i in eachindex(infcs)
        triangles[i]=Triangle(infcs[i],zero(SVector{4,Float64}),false,false,zero(SVector{3,Float64}))
    end
    for i in eachindex(invts)
        vertices[i]=Vertex(invts[i],0,0,zero(SVector{10,Float64}),false)
    end
    # main iteration loop
    deleted_triangles=0;
    deleted0 = Bool[]
    deleted1 = Bool[]
    refs = TRef[]

    triangle_count=length(triangles)

    for iteration = 0:100
        triangle_count-deleted_triangles<=target_count && break

        # update mesh once in a while
        @show iteration
        iszero(iteration%5) && update_mesh!(triangles, vertices, refs, iteration)

        # clear dirty flag
        for t in triangles
            t.dirty = false
        end

        #
        # All triangles with edges below the threshold will be removed
        #
        # The following numbers works well for most models.
        # If it does not, try to adjust the 3 parameters
        #
        threshold = 0.000000001*(iteration+3)^agressiveness

        # remove vertices & mark deleted triangles
        for i in eachindex(triangles)
            t = triangles[i]
            t.err[4]>threshold && continue
            t.deleted && continue
            t.dirty && continue

            for j = 1:3
                if t.err[j] < threshold

                    i0 = t.v[j]
                    v0 = vertices[i0]
                    i1 = t.v[(j%3)+1]
                    v1 = vertices[i1]
                    # Border check
                    v0.border != v1.border && continue

                    # Compute vertex to collapse to
                    e, p = calculate_error(vertices,i0,i1)
                    resize!(deleted0,v0.tcount) # normals temporarily
                    resize!(deleted1,v1.tcount) # normals temporarily
                    # don't remove if flipped
                    flipped(vertices,triangles,refs,p,i0,i1,v0,v1,deleted0) && continue

                    flipped(vertices,triangles,refs,p,i1,i0,v1,v0,deleted1) && continue

                    # not flipped, so remove edge
                    v0.p = p
                    v0.q = v1.q .+ v0.q
                    tstart = length(refs)

                    deleted_triangles = update_triangles(triangles,vertices,refs,i0,v0,deleted0,deleted_triangles)
                    deleted_triangles = update_triangles(triangles,vertices,refs,i0,v1,deleted1,deleted_triangles)

                    tcount = length(refs) - tstart

                    if tcount <= v0.tcount
                        # save ram
                        s = v0.tstart
                        if !iszero(tcount)
                            for ne in 0:(tcount-1)
                                refs[s+ne] = refs[tstart+ne]
                            end
                        end
                    else
                        # append
                        v0.tstart = tstart
                    end
                    v0.tcount = tcount
                    break
                end
            end
            # done?
            triangle_count-deleted_triangles <= target_count && break
        end
    end
    # clean up mesh
    compact_mesh(vertices,triangles)
    resize!(invts, length(vertices))
    resize!(infcs, length(triangles))
    for i in eachindex(vertices)
        invts[i] = vertices[i].p
    end
    for i in eachindex(triangles)
        infcs[i] = triangles[i].v
    end

    return invts, infcs
end #simplify_mesh()

"""
Check if a triangle flips when this edge is removed
"""
function flipped(vertices, triangles, refs, p, i0, i1, v0::Vertex, v1::Vertex, deleted::Vector{Bool})
    for k in 1:v0.tcount
        t=triangles[refs[v0.tstart+k].tid]

        t.deleted && continue

        s=refs[v0.tstart+k].tvertex
        id1=t.v[(s%3)+1]
        id2=t.v[((s+1)%3)+1]

        if id1==i1 || id2==i1 # delete ?
            deleted[k] = true
            continue
        end
        d1 = normalize(vertices[id1].p .- p)
        d2 = normalize(vertices[id2].p .- p)
        abs(dot(d1, d2))>0.999 && return true
        n = normalize(cross(d1,d2))
        deleted[k]=false
        dot(n,t.n) < 0.2 && return true
    end
    return false
end

# Update triangle connections and edge error after a edge is collapsed
function update_triangles(triangles, vertices, refs, i0, v::Vertex, deleted, deleted_triangles)
    dt = deleted_triangles
    for k in 1:v.tcount
        r = refs[v.tstart+k-1]
        t = triangles[r.tid]
        t.deleted && continue
        if deleted[k]
            t.deleted=true
            dt += 1
            continue
        end
        ind = r.tvertex
        tv = t.v
        if ind == 1
            triangles[r.tid].v = SVector(i0,tv[2],tv[3])
        elseif ind == 2
            triangles[r.tid].v = SVector(tv[1],i0,tv[3])
        else
            triangles[r.tid].v = SVector(tv[1],tv[2],i0)
        end
        triangles[r.tid].dirty = true
        err0, p = calculate_error(vertices, tv[1], tv[2])
        err1, p = calculate_error(vertices, tv[2], tv[3])
        err2, p = calculate_error(vertices, tv[3], tv[1])
        err3 = min(err0, err1, err2)
        t.err = SVector{4,Float64}(err0,err1,err2,err3)
        push!(refs, refs[v.tstart+k])
    end
    dt
end

function update_mesh!(triangles::Vector{Triangle}, vertices::Vector, refs::Vector{TRef}, iteration)
    if iteration > 0 # compact triangles
        dst = 1
        for i in eachindex(triangles)
            if !(triangles[i].deleted)
                triangles[dst]=triangles[i]
                dst += 1
            end
        end
        resize!(triangles,dst)
    end
    #
    # Init Quadrics by Plane & Edge Errors
    #
    # required at the beginning ( iteration == 0 )
    # recomputing during the simplification is not required,
    # but mostly improves the result for closed meshes
    #
    if iszero(iteration)
        for i in eachindex(vertices)
            vertices[i].q = zero(SVector{10,Float64})
        end

        for i in eachindex(triangles)
            t = triangles[i]
            p1 = vertices[t.v[1]].p
            p2 = vertices[t.v[2]].p
            p3 = vertices[t.v[3]].p
            n = normalize(cross(p2 .- p1, p3 .- p1))
            triangles[i].n = n
            vertices[t.v[1]].q += SymetricMatrix(n[1],n[2],n[3],dot(-n,p1))
            vertices[t.v[2]].q += SymetricMatrix(n[1],n[2],n[3],dot(-n,p1))
            vertices[t.v[3]].q += SymetricMatrix(n[1],n[2],n[3],dot(-n,p1))
        end
        for i in eachindex(triangles)
            # Calc Edge Error
            t = triangles[i]
            err1, _ = calculate_error(vertices,t.v[1],t.v[2])
            err2, _ = calculate_error(vertices,t.v[2],t.v[3])
            err3, _ = calculate_error(vertices,t.v[3],t.v[1])
            err4 = min(err1,err2,err3)
            t.err = SVector{4,Float64}(err1,err2,err3,err4)
        end
    end

    # Init Reference ID list
    for i in eachindex(vertices)
        vertices[i].tstart = 0
        vertices[i].tcount = 0
    end
    @show "hello"
    for i in eachindex(triangles)
        for j in 1:3
            #@show triangles[i].v[j], length(vertices), i, length(triangles)
            vertices[triangles[i].v[j]].tcount += 1
        end
    end
    @show "exited"
    tstart = 1
    for i in eachindex(vertices)
        vertices[i].tstart = tstart
        tstart += vertices[i].tcount
        vertices[i].tcount = 1
    end

    # Write References
    resize!(refs, length(triangles)*3)
    for i in eachindex(triangles)
        t=triangles[i]
        for j in 1:3
            v=vertices[t.v[j]]
            refs[v.tstart+v.tcount-1] = TRef(i, j)
            vertices[t.v[j]].tcount += 1
        end
    end

    # Identify boundary : vertices[].border=0,1
    if iszero(iteration)
        vcount = Int[]
        vids = Int[]

        for i in eachindex(vertices)
            vertices[i].border = false
        end

        for i in eachindex(vertices)
            v=vertices[i]
            empty!(vcount)
            empty!(vids)
            for j in 1:(v.tcount-1)
                k = refs[v.tstart+j-1].tid
                t = triangles[k]
                for k = 1:3
                    ofs = 0
                    id = t.v[k]
                    while ofs < length(vcount)
                        vids[ofs+1] == id && break
                        ofs += 1
                    end
                    if ofs == length(vcount)
                        push!(vcount,1)
                        push!(vids, id)
                    else
                        vcount[ofs+1] += 1
                    end
                end
            end
            for j in eachindex(vcount)
                if vcount[j] == 1
                    vertices[vids[j]].border = true
                end
            end
        end
    end
end

# Finally compact mesh before exiting
function compact_mesh(vertices::Vector,triangles::Vector{Triangle})
    dst = 1
    for i in eachindex(vertices)
        vertices[i].tcount=0
    end
    for i in eachindex(triangles)
        if !(triangles[i].deleted)
            t = triangles[i]
            triangles[dst]=t
            dst += 1
            for j in 1:3
                vertices[t.v[j]].tcount = 1
            end
        end
    end
    resize!(triangles, dst)
    dst = 1
    for i in eachindex(vertices)
        if !iszero(vertices[i].tcount)
            vertices[i].tstart=dst;
            vertices[dst].p = vertices[i].p
            dst += 1
        end
    end
    for i in eachindex(triangles)
        #@show triangles[i]
        #@show length(vertices), triangles[i].v[1],triangles[i].v[2],triangles[i].v[3]
        triangles[i].v = SVector(vertices[triangles[i].v[1]].tstart,
                                 vertices[triangles[i].v[2]].tstart,
                                 vertices[triangles[i].v[3]].tstart)
    end
    resize!(vertices,dst-1)
end

"""
Error between vertex and Quadric
"""
function vertex_error(q::SVector{10,Float64}, x::Float64, y::Float64, z::Float64)
    q[1]*x*x + 2*q[2]*x*y + 2*q[3]*x*z + 2*q[4]*x + q[5]*y*y + 2*q[6]*y*z + 2*q[7]*y + q[8]*z*z + 2*q[9]*z + q[10]
end

"""
Error for one edge
"""
function calculate_error(vertices::Vector, id_v1::Int, id_v2::Int)
    # compute interpolated vertex

    q = vertices[id_v1].q .+ vertices[id_v2].q
    border = vertices[id_v1].border & vertices[id_v2].border
    d = det(q,1, 2, 3, 2, 5, 6, 3, 6, 8)

    if !iszero(d) && !border
        # q_delta is invertible
        x = -1/d*(det(q,2, 3, 4, 5, 6, 7, 6, 8, 9)) # vx = A41/det(q_delta)
        y =  1/d*(det(q,1, 3, 4, 2, 6, 7, 3, 8, 9)) # vy = A42/det(q_delta)
        z = -1/d*(det(q,1, 2, 4, 2, 5, 7, 3, 6,  9)) # vz = A43/det(q_delta)

        return vertex_error(q, x, y, z), SVector{3,Float64}(x,y,z)
    else
        # det = 0 -> try to find best result
        p1 = vertices[id_v1].p
        p2 = vertices[id_v2].p
        p3 = (p1.+p2)/2
        error1 = vertex_error(q, p1[1],p1[2],p1[3])
        error2 = vertex_error(q, p2[1],p2[2],p2[3])
        error3 = vertex_error(q, p3[1],p3[2],p3[3])
        error4 = min(error1, error2, error3)
        if error1 == error4
            return error4, p1
        elseif error2 == error4
            return error4, p2
        elseif error3 == error4
            return error4, p3
        end
    end
    return 0.0, zero(SVector{3,Float64})
end

end # module


text:
# This file is used for comparison of fractal interpolation and spline intepolation

using Makie
using GeometryBasics
using FractalTools
import Interpolations
import Interpolations: Gridded, Linear

# Define test functions
f(x, y) = x^2 + y^2 + 1

# Construct gridded data to construct and test interpolant
ndata = 10
ntest = 20
data_mesh_2d = GeometryBasics.mesh(Tesselation(Rect2D(BigFloat.([-1., -1, 2., 2.])...), (ndata, ndata)))
data_points_3d = [Point(pnt[1], pnt[2], f(pnt[1], pnt[2])) for pnt in data_mesh_2d.position]     # points to construct interpolants
test_mesh_2d = GeometryBasics.mesh(Tesselation(Rect2D(BigFloat.([-1., -1, 2., 2.])...), (ntest, ntest)))
test_points_3d = [Point(pnt[1], pnt[2], f(pnt[1], pnt[2])) for pnt in test_mesh_2d.position]    # points to test interpolants

# Construct interpolants
fractal_interpolant = FractalTools.interpolate(data_points_3d)
spline_interpolant = Interpolations.interpolate(
    (unique(getindex.(data_points_3d, 1)), unique(getindex.(data_points_3d, 2))),       # x-y gridded data
    collect(reshape(getindex.(data_points_3d, 3), ndata, ndata)),                       # z values
    Gridded(Linear())                                                                   # Interpolation type
)

# Evaluations
true_data_points                        = map(pnt -> Point(f(pnt[1], pnt[2])), data_points_3d)
true_test_points                        = map(pnt -> Point(f(pnt[1], pnt[2])), test_points_3d)

fractal_interpolant_data_points         = map(pnt -> Point(fractal_interpolant(pnt[1], pnt[2])), data_points_3d)
fractal_interpolant_test_points         = map(pnt -> Point(fractal_interpolant(pnt[1], pnt[2])), test_points_3d)
fractal_interpolant_data_points_error   = Point.(
    abs.(last.(true_data_points - fractal_interpolant_data_points)) ./ abs.(last.(true_data_points)) * 100
    )
fractal_interpolant_test_points_error   = Point.(
    abs.(last.(true_test_points - fractal_interpolant_test_points)) ./ abs.(last.(true_test_points)) * 100
    )

spline_interpolant_data_points          = map(pnt -> Point(spline_interpolant(pnt[1], pnt[2])), data_points_3d)
spline_interpolant_test_points          = map(pnt -> Point(fractal_interpolant(pnt[1], pnt[2])), test_points_3d)
spline_interpolant_data_points_error   = Point.(
    abs.(last.(true_data_points - spline_interpolant_data_points)) ./ abs.(last.(true_data_points)) * 100
    )
spline_interpolant_test_points_error    = Point.(
    abs.(last.(true_test_points - spline_interpolant_test_points)) ./ abs.(last.(true_test_points)) * 100
    )

# Plots
fig = Figure()
ls11 = LScene(fig[1, 1])
ls12 = LScene(fig[1, 2])
ls13 = LScene(fig[1, 3])
ls21 = LScene(fig[2, 1])
ls22 = LScene(fig[2, 2])
ls23 = LScene(fig[2, 3])

trisurf!(ls11, data_mesh_2d.position, true_data_points, meshcolor3=last.(true_data_points))
trisurf!(ls11, data_mesh_2d.position, fractal_interpolant_data_points, meshcolor3=last.(fractal_interpolant_data_points), colormap=:heat)
wireframe!(ls11, data_mesh_2d)

trisurf!(ls12, data_mesh_2d.position, fractal_interpolant_data_points_error, meshcolor3=last.(fractal_interpolant_data_points_error))
wireframe!(ls12, data_mesh_2d)

trisurf!(ls13, test_mesh_2d.position, fractal_interpolant_test_points_error, meshcolor3=last.(fractal_interpolant_test_points_error))
wireframe!(ls13, test_mesh_2d)

trisurf!(ls21, data_mesh_2d.position, true_data_points, meshcolor3=last.(true_data_points))
trisurf!(ls21, data_mesh_2d.position, spline_interpolant_data_points, meshcolor3=last.(spline_interpolant_data_points), colormap=:heat)
wireframe!(ls21, data_mesh_2d)
trisurf!(ls22, data_mesh_2d.position, spline_interpolant_data_points_error, meshcolor3=last.(spline_interpolant_data_points_error))
wireframe!(ls22, data_mesh_2d)
trisurf!(ls23, test_mesh_2d.position, spline_interpolant_test_points_error, meshcolor3=last.(spline_interpolant_test_points_error))
wireframe!(ls23, test_mesh_2d)

display(fig)


text:
"""
    struct SimulationFree <: SimulationSpace

Simulation space for N-body simulations with all-with-all interactions.
The method implemented so far uses a brute force approach and only parallelizes over the number of agents,
hence having a minimum cost of computation of ``O(N)`` at best.
In the future we will implement parallelization over the second loop to further reduce the computational cost in favourable systems.
"""
struct SimulationFree <: SimulationSpace

    box::Array{<:FlatBoundary,1}
    medium::Array{<:Medium,1}

end

function SimulationFree(abm::Agent; box::Array{<:Any,1}=Array{FlatBoundary,1}(), medium::Array{<:Medium,1}=Array{Medium,1}())

    #Check dimensionality
    if length(box) != abm.dims && !isempty(abm.declaredSymbols["Medium"])
        error("Box has to be the same length as dimensions.")
    end
    #Make consistent box format adding Open to tuples
    box2 = Array{FlatBoundary,1}()
    for i in 1:length(box)
        if typeof(box[i])<:Tuple{Symbol,<:Real,<:Real}
            push!(box2, Bound(box[i]...))
        elseif typeof(box[i])<:FlatBoundary
            push!(box2,box[i])
        else
            error("Dimension has to be defined as a tupple with (Symbol, Real, Real) or a FlatBoundary Type.")
        end
    end

    #Check limits are correct
    vars = [i.s for i in box2]
    if !isempty(vars)
        checkIsDeclared_(abm,vars)
    end
    for i in box2
        if i.max <= i.min
            error("Superior limit is equal or smaller than inferior limit for ", i.s, ". The entry should be of the form (BoundaryType)(symbol,min,max,vargs...).")
        end
    end

    #Check medium has the same dimensions
    if abm.dims != length(medium) && length(abm.declaredSymbols["Medium"]) > 0
        error("Medium has to be specified with the same dimensions as the model. For that, it is necessary to also define a box.")
    end

    #Check symbols
    for b in box2
        for i in keys(b.addSymbols)
            for j in b.addSymbols[i]
                checkIsDeclared_(abm,j)
            end
        end
    end

    return SimulationFree(box2,medium)

end

function arguments_!(program::Program_, abm::Agent, a::SimulationFree, platform::String)

    return Nothing
end

function loop_(program::Program_, abm::Agent, a::SimulationFree, code::Expr, platform::String)

    code = vectorize_(abm, code, program)
    loop = :(
        for ic2_ in 1:N
            $code
        end)
    loop = simpleFirstLoop_(platform, loop)
    loop = subs_(loop,:nnic2_,:ic2_)

    return loop
end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
r5(x) = round(x,5)
const b = Beta(1,5)

import Base.Collections: heapify, heapify!, heappush!, heappop!
import Base: isless, zero, push!, show, display, getindex, setindex!, sort!, length, start, next, done

export Transition, ReplayMemory,SequentialReplayMemory,SortedReplayMemory, sample_greedy!, sample_beta!, sample_uniform!, push!, delete_and_push!, sort!, setindex!

type Transition
    s::Vector{Float64}
    s1::Vector{Float64}
    a::Vector{Float64}
    r::Float64
    δ::Float64
    t::Union{Int,Float64}
    i::Int
end

display(t::Transition) = println("s: ", string(r5(t.s)), " s1: ", string(r5(t.s1)), " a: ", string(r5(t.a)), " r: ", string(r5(t.r)), " δ: ", string(r5(t.δ)))
show(t::Transition) = display(t)

isless(t1::Transition,t2::Transition) = isless(abs(t1.δ),abs(t2.δ))
zero(Transition) = Transition([0.],[0.],[0.],0.,0,0,0)

abstract ReplayMemory

type SortedReplayMemory <: ReplayMemory
    mem::Vector{Transition}
    s::Int
    rs::Float64
    SortedReplayMemory(s::Int) = new(Vector{Transition}(),s, sum(1./(1:s)))
end

type SequentialReplayMemory <: ReplayMemory
    mem::Vector{Transition}
    s::Int
    _i::Int
    _filled::Bool
    _last_sampled::Vector{Int}
    SequentialReplayMemory(s::Int) = new(Vector{Transition}(s),s,0, false,Int[])
end

start(mem::ReplayMemory) = start(mem.mem)
next(mem::ReplayMemory,i) = next(mem.mem,i)
done(mem::SortedReplayMemory,i) = done(mem.mem,i)
done(mem::SequentialReplayMemory,i) = (mem._filled ? mem.s : mem._i) < i

display(mem::SortedReplayMemory) = map(display,mem.mem)
display(mem::SequentialReplayMemory) = mem._filled ? map(display,mem.mem) : [display(mem[i]) for i = 1:mem._i]
show(t::Transition) = map(display,mem.mem)

getindex(mem::ReplayMemory, ind) = mem.mem[ind]
setindex!(mem::ReplayMemory, trans::Transition, ind) = setindex!(mem.mem, trans, ind)
sort!(mem::SortedReplayMemory) = sort!(mem.mem,rev=true)
sort!(mem::SequentialReplayMemory) = nothing

function display(mem::ReplayMemory)
    map(display,mem.mem)
    nothing
end

length(mem::SortedReplayMemory) = length(mem.mem)
length(mem::SequentialReplayMemory) = mem._filled ? mem.s : mem._i

function push!(mem::SortedReplayMemory, t::Transition)
    if length(mem) >= mem.s
        return delete_and_push!(mem,t)
    end
    heappush!(mem.mem, t, Base.Order.Reverse)
    nothing
end

function push!(mem::ReplayMemory, t::Vector{Transition})
    for ti in t
        push!(mem,ti)
    end
end

function push!(mem::SequentialReplayMemory, t::Transition)
    if length(mem._last_sampled) > 0
        ind = pop!(mem._last_sampled)
        mem[ind] = t
        return
    end
    mem._i = (mem._i % mem.s) + 1
    mem[mem._i] = t
end

function delete_and_push!(mem::SortedReplayMemory, t::Transition)
    l2 = length(mem) ÷ 2
    deleteat!(mem.mem,l2+rand(1:(l2)))
    heappush!(mem.mem, t, Base.Order.Reverse)
    nothing
end

function sample_greedy!(mem::SortedReplayMemory)
    heappop!(mem.mem, Base.Order.Reverse)
    # mem[1]
end

function sample_beta!(mem::SortedReplayMemory)
    ind = ceil(Int,length(mem)*rand(b))
    retval = mem[ind]
    deleteat!(mem.mem,ind)
    retval
end

function sample_uniform!(mem::SortedReplayMemory)
    ind = rand(1:length(mem))
    retval = mem[ind]
    deleteat!(mem.mem,ind)
    retval
end

function sample_uniform!(mem::SequentialReplayMemory)
        lastind = mem._filled ? mem.s : mem._i
        ind = rand(1:lastind)
        push!(mem._last_sampled, ind)
        retval = mem[ind]
end

function sample_many!(mem::ReplayMemory,n,sample_fun)
    retvec = Vector{Transition}(n)
    for i = 1:n
        retvec[i] = sample_fun(mem)
    end
    retvec
end

sample_greedy!(mem::SortedReplayMemory,n) = sample_many!(mem::SortedReplayMemory,n,sample_greedy!)
sample_beta!(mem::SortedReplayMemory,n) = sample_many!(mem::SortedReplayMemory,n,sample_beta!)
sample_uniform!(mem::ReplayMemory,n) = sample_many!(mem::ReplayMemory,n,sample_uniform!)


text:
function pass_variables(m1, m2, f::MOI.AbstractFunction, map, tol::Number)
    vis = non_zero_variables(f, tol)
    pass_variables(m1, m2, vis, map)
    return
end

function pass_variables(m1, m2, vis, map)
    for vi in vis
        if !haskey(map, vi)
            vname = MOI.get(m1, MOI.VariableName(), vi)
            x = MOI.add_variable(m2)
            MOI.set(m2, MOI.VariableName(), x, vname)
            map[vi] = x
        end
    end
    return
end

is_approx_zero(lhs::MOI.EqualTo, tol = 1e-8) = abs(lhs.value) < tol
is_approx_zero(lhs::MOI.GreaterThan, tol = 1e-8) = abs(lhs.lower) < tol
is_approx_zero(lhs::MOI.LessThan, tol = 1e-8) = abs(lhs.upper) < tol
is_approx_zero(lhs::S, tol = 1e-8) where S<:MOI.AbstractSet = true

set_difference(lhs::S, rhs::S) where S = lhs
set_sum(lhs::S, rhs::S) where S = lhs

set_sum(lhs::MOI.EqualTo, rhs::MOI.EqualTo) = MOI.EqualTo(rhs.value + lhs.value)
set_difference(lhs::MOI.EqualTo, rhs::MOI.EqualTo) = MOI.EqualTo(lhs.value - rhs.value)

set_sum(lhs::MOI.GreaterThan, rhs::MOI.GreaterThan) = MOI.GreaterThan(rhs.lower + lhs.lower)
set_difference(lhs::MOI.GreaterThan, rhs::MOI.GreaterThan) = MOI.GreaterThan(lhs.lower - rhs.lower)

set_sum(lhs::MOI.LessThan, rhs::MOI.LessThan) = MOI.LessThan(rhs.upper + lhs.upper)
set_difference(lhs::MOI.LessThan, rhs::MOI.LessThan) = MOI.LessThan(lhs.upper - rhs.upper)

function non_zero_variables(f::MOI.ScalarAffineFunction, tol)
    vis = Set{VI}()
    to_stay = Int[]
    for i in 1:length(f.terms)
        if abs(MOI.coefficient(f.terms[i])) > tol
            push!(vis, f.terms[i].variable_index)
            push!(to_stay, i)
        end
    end
    f.terms = f.terms[to_stay]
    return vis
end
function non_zero_variables(f::MOI.SingleVariable, tol)
    vis = Set{VI}(f.variable)
    return vis
end

# TODO - Other functions
function unsafe_move_variables(set::Set, f1, f2)
    for i in set
        unsafe_move_variable(i, f1, f2)
    end
end
function unsafe_move_variable(vi, f1::MOI.ScalarAffineFunction, f2::MOI.ScalarAffineFunction) # suppose canonical
    for i in 1:length(f1.terms)
        if MOIU._hasvar(f1.terms[i], vi)
            push!(f2.terms, f1.terms[i])
            deleteat!(f1.terms, i) # suppose canonical
            return true
        end
    end
    return false
end


text:
# Doppler shift wavelengths
"""Calculate lambda in solar system barycenter frame given observed lambda and barycentric correction)."""
function lambda_ssb_from_obs(lambda_obs, z_bc)
   lambda_obs * (1+z_bc)
end

"""Calculate log lambda in solar system barycenter frame given observed log lambda and barycentric correction)."""
function log_lambda_ssb_from_obs(log_lambda_obs, z_bc)
   log_lambda_obs + log1p(z_bc)
end

"""Calculate lambda in source frame given observed lambda, barycentric correction and source z relative to solar system barycenter)."""
function lambda_source(lambda_obs, z_bc; z_source = zero(z_bc) )
   lambda_obs * (1+z_bc) / (1+z_source) # TODO: Check EQNs
end

"""Calculate lambda in source frame given observed lambda, barycentric correction and source z relative to solar system barycenter)."""
function log_lambda_source(log_lambda_obs, z_bc; z_source = zero(z_bc) )
   log_lambda_obs + log1p(z_bc) - log1p(z_source) # TODO: Check EQNs
end
export lambda_ssb_from_obs, log_lambda_ssb_from_obs, lambda_source, log_lambda_source

#----

using Statistics     # For computing sum/mean/median over pixels

# Size of data
function num_pixels_per_order(data::EOT) where {EOT <: AbstractEchelleObservation}
   size(data.flux,1)
end
function num_orders(data::EOT) where {EOT <: AbstractEchelleObservation}
   size(data.flux,2)
end

function num_pixels_per_order(data::EOST) where {EOST <: AbstractEchelleObservationSet}
   size(data.flux,1)
end
function num_orders(data::EOST) where {EOST <: AbstractEchelleObservationSet}
   size(data.flux,2)
end

function num_obs(data::EOST) where {EOST <: AbstractEchelleObservationSet}
   size(data.flux,3)
end
export num_pixels_per_order, num_orders, num_obs

#---

# Utility functions, so can work with either an EcheelObservation or an EcheelObservationSet
function get_min_log_lambda(data::T, pixel_range::AbstractRange = 1:num_pixels_per_order(data)
      ) where T<:AbstractEchelleObservationSet
   pixel_lo = max(1,pixel_range[1] - 1)
   log_lambda_lo = log(minimum(data.lambda[pixel_lo,order,:]))
end
function get_max_log_lambda(data::T, pixel_range::AbstractRange = 1:num_pixels_per_order(data)
      ) where T<:AbstractEchelleObservationSet
   pixel_hi = min(num_pixels_per_order(data), pixel_range[end]+1)
   log_lambda_hi = log(maximum(data.lambda[pixel_hi,order,:]))
end

function get_min_log_lambda(data::T, pixel_range::AbstractRange = 1:num_pixels_per_order(data)
      ) where T<:AbstractEchelleObservation
   pixel_lo = max(1,pixel_range[1] - 1)
   log_lambda_lo = log(minimum(data.lambda[pixel_lo,order]))
end
function get_max_log_lambda(data::T, pixel_range::AbstractRange = 1:num_pixels_per_order(data)
      ) where T<:AbstractEchelleObservation
   pixel_hi = min(num_pixels_per_order(data), pixel_range[end]+1)
   log_lambda_hi = log(maximum(data.lambda[pixel_hi,order]))
end
export get_min_log_lambda, get_max_log_lambda

# summary statistics for data
function flux_mean_over_pixels(data::AbstractEchelleObservationSet)
   reshape( sum(data.flux./data.var,dims=3)./sum(one(eltype(data.var))./data.var,dims=3),
            (size(data.flux,1),size(data.flux,2)) )
end

function flux_median_over_pixels(data::AbstractEchelleObservationSet)
   reshape(median(data.flux,dims=3),(size(data.flux,1),size(data.flux,2)))
end

function lambda_mean_over_pixels(data::AbstractEchelleObservationSet)
   reshape(mean(data.lambda,dims=3),(size(data.lambda,1),size(data.lambda,2)))
end

function lambda_median_over_pixels(data::AbstractEchelleObservationSet)
   reshape(median(data.lambda,dims=3),(size(data.lambda,1),size(data.lambda,2)))
end
export flux_mean_over_pixels, lambda_mean_over_pixels
export flux_median_over_pixels, lambda_median_over_pixels


text:
using Example
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@test hello("Julia") == "Hello, Julia"
@test domath(2.0) ≈ 7.0


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

"""
    nrbrevolve(crv::NURBS1D{I,F},vec::Vector{F},
               pnt::Vector{F}=[0.0;0.0;0.0],theta::F=2.0*π
              )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}

Construct a NURBS surface by revolving the profile NURBS curve around an axis
defined by a point and vector.

# Arguments:
- `crv`: NURBS curve to revolve, see [`nrbmak`](@ref).
- `vec`: Vector defining the direction of the rotation axis.
- `pnt`: Coordinate of the point used to define the axis of rotation.
- `theta`: Angle to revolve the curve, default `2*π`.

# Output:
- `srf`: Constructed surface.

# Examples:
```julia
julia> srf = nrbrevolve(crv,vec)
julia> srf = nrbrevolve(crv,vec,pnt)
julia> srf = nrbrevolve(crv,vec,pnt,ang)
```

Construct a sphere by rotating a semicircle around a x-axis.
```julia
julia> crv = nrbcirc(1.0,[0.0;0.0;0.0],0.0,pi);
julia> srf = nrbrevolve(crv,[1.0;0.0;0.0],[0.0;0.0;0.0]);
julia> nrbplot(srf,[20;20]);
```

# The algorithm:

1. vectrans the point to the origin `[0.0;0.0;0.0]`
2. rotate the vector into alignment with the z-axis
3. for each control point along the curve
   1. determine the radius and angle of control point to the z-axis
   2. construct a circular arc in the x-y plane with this radius and start
      angle and sweep angle theta
   3. combine the arc and profile, coefs and weights.
   4. next control point
4. rotate and vectrans the surface back into position by reversing 1 and 2.
"""
function nrbrevolve(crv::NURBS1D{I,F},vec::Vector{F},
                    pnt::Vector{F}=[0.0;0.0;0.0],theta::F=2.0*π
                   )::NURBS2D{I,F} where {I<:Integer,F<:AbstractFloat}
# TODO: optional arguments in julia are always at the end i switched vec and
#       pnt argument > examine if that did change the behaviour in this toolbox

if length(vec) != 3 || length(pnt) !=3
  throw(ArgumentError("All point and vector coordinates must be 3D!"));
end # if

# Translate and rotate the curve into alignment with the z-axis
T  = vectrans(-pnt);
angx = vecangle(vec[1],vec[3]);
RY = vecroty(-angx);
vectmp = RY * [vecnorm_toolbox(vec);1.0];
angy = vecangle(vectmp[2],vectmp[3]);
RX = vecrotx(angy);
curve = nrbtform(crv,RX*RY*T);

# Construct an arc
arc = nrbcirc(1.0,[0.0;0.0;0.0],0.0,theta);

# Construct the surface
coefs = zeros(F,(4,arc.number,curve.number));
angle = vecangle(curve.coefs[2:2,:],curve.coefs[1:1,:]);
radius = vecmag(curve.coefs[1:2,:]);
for i in 1:curve.number
  coefs[:,:,i] = vecrotz(angle[i]) *
                 vectrans([0.0;0.0;curve.coefs[3,i]]) *
                 vecscale([radius[i];radius[i]]) * arc.coefs;
  coefs[4,:,i] = coefs[4,:,i] * curve.coefs[4,i];
end # for i
surf = nrbmak(coefs,[arc.knots,curve.knots]);

# Rotate and vectrans the surface back into position
T = vectrans(pnt);
RX = vecrotx(-angy);
RY = vecroty(angx);
return nrbtform(surf,T*RY*RX);
end # nrbrevolve


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Tetrahedron(p1, p2, p3, p4)

A tetrahedron with points `p1`, `p2`, `p3`, `p4`.
"""
struct Tetrahedron{Dim,T,V<:AbstractVector{Point{Dim,T}}} <: Polyhedron{Dim,T}
  vertices::V
end

function measure(t::Tetrahedron)
  a, b, c, d = t.vertices
  abs((a - d) ⋅ ((b - d) × (c - d))) / 6
end


text:
using Distributed, Pkg
addprocs(8)
@everywhere using Distributed, Pkg
Pkg.activate(".")
@everywhere Pkg.activate(".")
using LTMSim, DataFrames, SimpleHypergraphs, Statistics, Plots, PyPlot, Random, Serialization, LaTeXStrings
@everywhere using LTMSim, DataFrames, SimpleHypergraphs, Statistics, Random, Serialization

nvalues = [10, 20, 40, 80]
runs = 10
data = Dict{String, Vector{Vector{Int}}}()

data["BinarySearch(H)"]=Vector{Vector{Int}}()
data["Greedy(H)"]=Vector{Vector{Int}}()
data["Greedy([H]₂)"]=Vector{Vector{Int}}()
data["SubTSS(H)"]=Vector{Vector{Int}}()

for n=nvalues
	println("n=$n")
    results1 = Vector{Int}()
    results2 = Vector{Int}()
    results3 = Vector{Int}()

    results = @distributed (append!) for run=1:runs
        h = randomHkuniform(500, 500, n)
        metaV = randMetaV(h)
        metaE = proportionalMetaE(h,0.5)

        r1 = greedy_tss_2section(h,metaV,metaE)
        r2 = bisect(h,metaV,metaE)
		r3 = greedy_tss(h,metaV,metaE)
		r4 = sub_tss_opt2(h,metaV,metaE)
		[(r1,r2,r3,r4)]

    end
    push!(data["Greedy([H]₂)"], [r[1] for r in results])
    push!(data["BinarySearch(H)"], [r[2] for r in results])
    push!(data["Greedy(H)"], [r[3] for r in results])
	push!(data["SubTSS(H)"], [r[4][1] for r in results])
    println("end ",n)

end

serialize("res/paper/exp1/random-k.data", data)
data = deserialize("res/paper/exp1/random-k.data")

labels = Dict{String, String}(
    "BinarySearch(H)" => "StaticGreedy",
    "Greedy(H)" => "DynamicGreedy",
    "Greedy([H]₂)" => L"DynamicGreedy_{[H]_2}",
	"SubTSS(H)" => "SubTSS"
)

ticks = nvalues

function set_box_color(bp, color)
    plt.setp(bp["boxes"], color=color)
    plt.setp(bp["whiskers"], color=color)
    plt.setp(bp["caps"], color=color)
    plt.setp(bp["medians"], color=color)
end

clf()

plt.figure(figsize=(7,5))

val = -0.4
c = 1

colorz=["#2C7BB6", "#D7191C", "#FF8900", "#33CC33"]

for algo in ["BinarySearch(H)", "Greedy([H]₂)", "Greedy(H)", "SubTSS(H)"]#keys(data)
    global val, c
    b = plt.boxplot(
        data[algo],
        positions=collect(range(0, stop=length(data[algo])-1)).*2.0.+val,
        sym="",
        widths=0.2
    )
    set_box_color(b, colorz[c])
    c+=1
    val+=0.4
end

# draw temporary red and blue lines and use them to create a legend
plt.plot([], c="#2C7BB6", label=labels["BinarySearch(H)"])#label=collect(keys(data))[1])
plt.plot([], c="#D7191C", label=labels["Greedy([H]₂)"])#label=collect(keys(data))[2])
plt.plot([], c="#FF8900", label=labels["Greedy(H)"])#label=collect(keys(data))[3])
plt.plot([], c="#33CC33", label=labels["SubTSS(H)"])#label=collect(keys(data))[3])
plt.legend(fontsize="x-large", loc = "lower right")

plt.xticks(range(0, length(ticks) * 2, step=2), ticks, fontsize="x-large")
plt.yticks(fontsize="x-large")

#plt.xlim(-2, length(ticks)*2)
#plt.ylim(0, 8)
plt.ylim(0, 230)

ylabel("Influence set size", fontstyle = "italic", fontsize="xx-large", labelpad=10) #fontweight="semibold",
xlabel(L"$k$", fontsize="xx-large", fontweight="semibold", labelpad=10)

plt.tight_layout()

gcf()

PyPlot.savefig("res/paper/exp1/random-k.pdf")


text:
export AsapScheduler

import Base: show, similar

"""
    AsapScheduler

`AsapScheduler` executes scheduled actions as soon as possible and does not introduce any additional logic.
`AsapScheduler` is a default scheduler for almost all observables.
"""
struct AsapScheduler <: AbstractScheduler end

Base.show(io::IO, ::AsapScheduler) = print(io, "AsapScheduler()")

Base.similar(::AsapScheduler) = AsapScheduler()

makeinstance(::Type, ::AsapScheduler) = AsapScheduler()

instancetype(::Type, ::Type{<:AsapScheduler}) = AsapScheduler

scheduled_subscription!(source, actor, instance::AsapScheduler) = on_subscribe!(source, actor, instance)

scheduled_next!(actor, value, ::AsapScheduler) = on_next!(actor, value)
scheduled_error!(actor, err, ::AsapScheduler)  = on_error!(actor, err)
scheduled_complete!(actor, ::AsapScheduler)    = on_complete!(actor)


text:
### BioSequences.jl
###
### A julia package for the representation and manipulation of biological sequences.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

module BioSequences

export
    ###
    ### Symbols
    ###

    # Types & aliases
    NucleicAcid,
    DNA,
    RNA,
    DNA_A,
    DNA_C,
    DNA_G,
    DNA_T,
    DNA_M,
    DNA_R,
    DNA_W,
    DNA_S,
    DNA_Y,
    DNA_K,
    DNA_V,
    DNA_H,
    DNA_D,
    DNA_B,
    DNA_N,
    DNA_Gap,
    ACGT,
    ACGTN,
    RNA_A,
    RNA_C,
    RNA_G,
    RNA_U,
    RNA_M,
    RNA_R,
    RNA_W,
    RNA_S,
    RNA_Y,
    RNA_K,
    RNA_V,
    RNA_H,
    RNA_D,
    RNA_B,
    RNA_N,
    RNA_Gap,
    ACGU,
    ACGUN,
    AminoAcid,
    AA_A,
    AA_R,
    AA_N,
    AA_D,
    AA_C,
    AA_Q,
    AA_E,
    AA_G,
    AA_H,
    AA_I,
    AA_L,
    AA_K,
    AA_M,
    AA_F,
    AA_P,
    AA_S,
    AA_T,
    AA_W,
    AA_Y,
    AA_V,
    AA_O,
    AA_U,
    AA_B,
    AA_J,
    AA_Z,
    AA_X,
    AA_Term,
    AA_Gap,

    # Predicates
    isGC,
    iscompatible,
    isambiguous,
    iscertain,
    isgap,
    ispurine,
    ispyrimidine,

    ###
    ### Alphabets
    ###

    # Types & aliases
    Alphabet,
    NucleicAcidAlphabet,
    DNAAlphabet,
    RNAAlphabet,
    AminoAcidAlphabet,
    CharAlphabet,

    ###
    ### BioSequences
    ###

    join!,

    # Type & aliases
    BioSequence,
    NucleotideSeq,
    AminoAcidSeq,

    # Indexing
    unsafe_setindex!,

    # Predicates
    ispalindromic,
    hasambiguity,
    isrepetitive,
    iscanonical,

    # Transformations
    canonical,
    canonical!,
    complement,
    complement!,
    reverse_complement,
    reverse_complement!,
    ungap,
    ungap!,

    # Iteration
    each,
    fwmer,
    bwmer,

    ###
    ### LongSequence
    ###

    # Type & aliases
    LongSequence,
    LongDNASeq,
    LongRNASeq,
    LongAminoAcidSeq,
    LongCharSeq,
    LongSubSeq,

    # Random
    SamplerUniform,
    SamplerWeighted,
    randseq,
    randdnaseq,
    randrnaseq,
    randaaseq,

    ###
    ### Mers
    ###

    # Type & aliases
    AbstractMer,
    Mer,
    DNAMer,
    RNAMer,
    DNAKmer,
    RNAKmer,

    BigMer,
    BigDNAMer,
    BigRNAMer,
    BigDNAKmer,
    BigRNAKmer,

    DNACodon,
    RNACodon,

    # Iteration
    neighbors,

    ###
    ### Sequence literals
    ###

    @dna_str,
    @rna_str,
    @aa_str,
    @char_str,
    @biore_str,
    @prosite_str,
    @mer_str,
    @bigmer_str,

    matched,
    captured,
    alphabet, # TODO: Resolve the use of alphabet - it's from BioSymbols.jl
    symbols,
    gap,
    mismatches,
    matches,
    n_ambiguous,
    n_gaps,
    n_certain,

    gc_content,

    eachcanonical,

    ###
    ### Composition
    ###
    Composition,
    composition,
    NucleicAcidCounts,

    translate!,
    translate,
    ncbi_trans_table,

    # Search
    ExactSearchQuery,
    ApproximateSearchQuery,
    approxsearch,
    approxsearchindex,
    approxrsearch,
    approxrsearchindex,
    PFM,
    PWM,
    maxscore,
    scoreat,

    ReferenceSequence,

    ###
    ### Demultiplexing
    ###
    Demultiplexer,
    demultiplex,

    seqmatrix,
    majorityvote,
    MinHashSketch,
    minhash,
    Site,
    Certain,
    Ambiguous,
    Gap,
    Match,
    Mismatch,
    count_pairwise

using BioGenerics
using BioSymbols
import Combinatorics
import IndexableBitVectors
import Twiddle: enumerate_nibbles,
    nibble_mask,
    count_0000_nibbles,
    count_1111_nibbles,
    count_nonzero_nibbles,
    count_00_bitpairs,
    count_01_bitpairs,
    count_10_bitpairs,
    count_11_bitpairs,
    count_nonzero_bitpairs,
    repeatpattern
using Random
using StableRNGs

BioSymbols.gap(::Type{Char}) = '-'

include("alphabet.jl")

# Load the bit-twiddling internals that optimised BioSequences methods depend on.
include("bit-manipulation/bit-manipulation.jl")

# The generic, abstract BioSequence type
include("biosequence/biosequence.jl")

# The definition of the LongSequence concrete type, and its method overloads...
include("longsequences/longsequence.jl")
include("longsequences/hash.jl")
include("longsequences/randseq.jl")

# The definition of the Skipmer concrete type, and its method overloads...
include("mers/mer.jl")

# The definition of the ReferenceSequence concrete type, and its method overloads...
include("nmask.jl")
include("refseq/refseq.jl")

# The generic iterators for any BioSequence...
include("iterators/condition.jl")
include("iterators/eachmer.jl")
include("iterators/skipmerfactory.jl")

include("composition.jl")

include("geneticcode.jl")
include("demultiplexer.jl")

# Pattern searching in sequences...
include("search/exact.jl")
include("search/approx.jl")
include("search/re.jl")
include("search/pwm.jl")

include("minhash.jl")

end  # module BioSequences


text:
# function parse_newick(io)
#     cont = read(io, String)
#     cont = replace(cont, r":\d+\.\d+" => "")
#     cont = replace(cont, r"_\d+_\d+" => "")
#     cont = replace(cont, r"\)\d+" => ")")
#     return parse_tuple(eval(Meta.parse(cont)))
# end

function parse_newicks(io)
    res = Tree[]
    skipws() = skipchars(isspace, io)
    skipws()
    while !eof(io)
        push!(res, parse_newick(io))
        accept(io, ';')
        skipws()
    end
    res
end

# function takenumber(s :: IO)
#     buffer = IOBuffer()
#     while isdigit(peek(s, Char))
#         print(buffer, peek(s, Char))
#         skip(s, 1)
#     end
#     parse(Int, take!(buffer))
# end

# function ps(s)
#     stack = Tree[]
#     while peek(s, Char) != ';'
#         if peek(s, Char) == '('
#             skip(s, 1)
#         end
#         if peek(s, Char) == ','
#             skip(s, 1)
#         end
#         if isdigit(peek(s, Char))
#             push!(stack, Leaf(takenumber(s), Multiset([tup]), nothing))
#         end
#         if peek(s, Char) == ')'
#             a = pop!(stack)
#             b = pop!(stack)
#             res = Branch(a, b, l._clades + r._clades, nothing)
#             a._parent = res
#             b._parent = res
#             push!(stack, res)
#             skip(s, 1)
#         end
#     end
#     stack
# end

# function parse_newicks(io)
#     res = Tree[]
#     tree_src = readuntil(io, ";"; keep=true)
#     skipchars(isspace, io)
#     while !isempty(tree_src)
#         push!(res, parse_newick(IOBuffer(tree_src)))
#         tree_src = readuntil(io, ";"; keep=true)
#         skipchars(isspace, io)
#     end
#     res
# end

# recursive descent parser for newicks... yeah apprantly

function takenumber(io :: IO)
    buf = IOBuffer()
    x = peek(io, Char)
    while isdigit(x)
        print(buf, x)
        skip(io, 1)
        x = peek(io, Char)
    end
    return parse(Int, String(take!(buf)))
end

function skip_suffix(io :: IO)
    while peek(io, Char) == '_'
        skip(io, 1)
        takenumber(io)
    end
end

function skiplength(io :: IO)
    if peek(io, Char) == ':' # oh boy, those pesky lengths...
        skipchars(io) do ch
            ch == ':' || ch == '.' || isdigit(ch)
        end
    end
end

function accept(io :: IO, ch :: Char)
    if peek(io, Char) != ch
        throw("$(peek(io, Char)) should not occur, it should be $ch")
    end
    skip(io, 1)
end

function buildbranch(l :: Tree, r :: Tree)
    res = Branch(l, r, +(l._clades, r._clades), nothing)
    l._parent = res
    r._parent = res
    res
end

function parse_newick(io :: IO) # return a newick value
    if isdigit(peek(io, Char))
        # it is a number value
        n = takenumber(io)
        skip_suffix(io)
        skiplength(io)
        Leaf(n, Multiset([n]), nothing)
    else
        # it is a value type
        accept(io, '(')
        l = parse_newick(io)
        accept(io, ',')
        r = parse_newick(io)
        accept(io, ')')
        if isdigit(peek(io, Char))
            takenumber(io)
            skiplength(io)
        end
        skiplength(io)
        buildbranch(l, r)
    end
end

function parse_tuple(tup)
    if isa(tup, Int)
        Leaf(tup, Multiset([tup]), nothing)
    else
        l = parse_tuple(tup[1])
        r = parse_tuple(tup[2])
        res = Branch(l, r, +(l._clades, r._clades), nothing)
        l._parent = res
        r._parent = res
        res
    end
end

newick = parse_tuple

function writetree(io::IO, tree :: Tree)
    buf = IOBuffer()
    notatetree(io, tree)
    println(io, ";")
end

function notatetree(buf :: IO, tree :: Tree)
    function f(t :: Tree)
        if isleaf(t)
            print(buf, t.value)
        else
            print(buf, '(')
            f(t.left)
            print(buf, ',')
            f(t.right)
            print(buf, ')')
        end
    end
    f(tree)
end

function notatetree(tree :: Tree)
    buf = IOBuffer()
    notatetree(buf, tree)
    take!(buf) |> String
end

function writetrees(io :: IO, trees)
    for t in trees
        writetree(io, t)
    end
end

function bintree2newick(t)
    if isa(t, Branch)
        (bintree2newick(t.left), bintree2newick(t.right))
    else
        t.value
    end
end

max_stars_repo_name: RuneBlaze/clade-phylo

max_issues_repo_name: RuneBlaze/clade-phylo

max_forks_repo_name: RuneBlaze/clade-phylo


text:
# ---
# title: 894. All Possible Full Binary Trees
# id: problem894
# author: Indigo
# date: 2021-02-01
# difficulty: Medium
# categories: Tree, Recursion
# link: <https://leetcode.com/problems/all-possible-full-binary-trees/description/>
# hidden: true
# ---
#
# A _full binary tree_  is a binary tree where each node has exactly 0 or 2
# children.
#
# Return a list of all possible full binary trees with `N` nodes.  Each element
# of the answer is the root node of one possible tree.
#
# Each `node` of each tree in the answer **must** have `node.val = 0`.
#
# You may return the final list of trees in any order.
#
#
#
# **Example 1:**
#
#
#
#     Input: 7
#     Output: [[0,0,0,null,null,0,0,null,null,0,0],[0,0,0,null,null,0,0,0,0],[0,0,0,0,0,0,0],[0,0,0,0,0,null,null,null,null,0,0],[0,0,0,0,0,null,null,0,0]]
#     Explanation:
#     ![](https://s3-lc-upload.s3.amazonaws.com/uploads/2018/08/22/fivetrees.png)
#
#
#
#
# **Note:**
#
#   * `1 <= N <= 20`
#
#
## @lc code=start
using LeetCode

function all_possible_fbt(N::Int)
    cache_res = [TreeNode{Int}[] for i in 1:N]
    cache_res[1] = [TreeNode(0)]
    if N % 2 == 0
        return TreeNode[]
    end
    for n in 3:2:N
        for i in 1:2:n-1
            lres = cache_res[i]
            rres = cache_res[n-i-1]
            append!(cache_res[n], TreeNode(0, lr, rr) for lr in lres for rr in rres)
        end
    end
    return cache_res[N]
end
## @lc code=end


text:
using FuzzySets, Test

step_size = 0.005
levels = collect(0:step_size:1)

a₁ = 2
a₂ = 5
w₁ = 0.25
w₂ = 0.75
A₁ = FuzzyNumber(levels, number=a₁)
A₂ = FuzzyNumber(levels, number=a₂)
W₁ = FuzzyNumber(levels, number=w₁)
W₂ = FuzzyNumber(levels, number=w₂)
μ = (w₁*a₁ + w₂*a₂) / (w₁ + w₂)

X⃗ = FuzzyVector([A₁, A₂])
W⃗ = FuzzyVector([W₁, W₂])
AVG = fuzzy_weighted_average(X⃗, W⃗)
@test core(AVG) == Interval(μ)
draw(AVG)

x₁ = Interval(1.8, 2.2)
x₂ = Interval(2.8, 3.2)
w₁ = Interval(0.55, 0.95)
w₂ = Interval(0.05, 0.45)
combinations = FuzzySets.getcombinations([x₁, x₂], [w₁, w₂])
@test combinations == [
    1.8  2.8  0.55  0.05;
    2.2  2.8  0.55  0.05;
    1.8  3.2  0.55  0.05;
    2.2  3.2  0.55  0.05;
    1.8  2.8  0.95  0.05;
    2.2  2.8  0.95  0.05;
    1.8  3.2  0.95  0.05;
    2.2  3.2  0.95  0.05;
    1.8  2.8  0.55  0.45;
    2.2  2.8  0.55  0.45;
    1.8  3.2  0.55  0.45;
    2.2  3.2  0.55  0.45;
    1.8  2.8  0.95  0.45;
    2.2  2.8  0.95  0.45;
    1.8  3.2  0.95  0.45;
    2.2  3.2  0.95  0.45;
]
grades = FuzzySets.fwa(combinations)
@test grades ≈ Interval(1.85, 2.65)


text:
###############
## bootstrap ##
###############

function resample(td::AbstractTimedata)
    # resample rows of Timedata object
    nRows = size(td, 1)
    inds = rand(1:nRows, nRows)

    return td[inds, :]
end

function resample(df::DataFrame)
    # resample rows of Timedata object
    nRows = size(df, 1)
    inds = rand(1:nRows, nRows)

    return df[inds, :]
end

function resample{T}(obj::Array{T, 1})
    # resample rows of Timedata object
    nRows = size(obj, 1)
    inds = rand(1:nRows, nRows)

    return obj[inds]
end

function resample{T}(obj::Array{T, 2})
    # resample rows of Timedata object
    nRows = size(obj, 1)
    inds = rand(1:nRows, nRows)

    return obj[inds, :]
end

macro bootstrap(nTimes, expr)
    # resample different types: Timematr, Arrays, DataFrames
    quote
        # evaluate real expression and stop time
        tt1 = time()
        $(esc(expr))
        tt2 = time() # in seconds

        # get type of real result
        typ = typeof($(esc(expr.args[1])))

        # estimate required resampling time
        t1 = time()
        expTime = round((tt2 - tt1) * $nTimes, 2)
        println("")
        println("Expected resampling time: $expTime seconds (rough guess)")

        # get quarters to print message
        intermediateSteps = [4 2 4/3]
        quarts = int(floor($nTimes ./ intermediateSteps))

        # get function to resample
        func = $(esc(expr.args[2].args[1]))

        # get data as first argument to function
        data = $(esc(expr.args[2].args[2]))

        bootstrVals = Array(typ, $nTimes)
        for ii=1:$nTimes
            # get sample
            samp = resample(data)

            # print message about remaining time
            if issubset(ii, quarts)
                nReps = $nTimes
                pct = floor(ii/nReps*100)
                ttSoFar = time() - t1

                remainPct = 100 - pct
                remainInUnitsOfDone = remainPct / pct
                remainTime = round(remainInUnitsOfDone*ttSoFar , 2)
                println("$pct percent done, $remainTime remaining")
            end

            # apply function to sample
            bootstrVals[ii] = func(samp)
        end

        # actually required time
        t2 = time();
        t = round(t2-t1, 2)
        println("Actual resampling time: $t seconds")
        println("")

        # try fast version of vcat(bootstrVals...) for optimized
        # output structure
        res = []
        try
            res = fastVcat(bootstrVals)
        catch
            res = bootstrVals
        end
    end
end

# fast version to convert Array{Array{Float64,2},1} to Array{Float64,
# 2} if individual entries have same length
#
# much faster than vcat(x...)
function fastVcat(x::Array{Array{Float64,2},1})
    # preallocation
    (nReps, nVals) = (size(x, 1), size(x[1], 2))
    res = Array(Float64, nReps, nVals)
    for ii=1:nReps
        res[ii, :] = x[ii]
    end
    return res
end

function fastVcat(x::Array{Array{Float64,1},1})
    # preallocation
    nReps = size(x, 1)
    nVals = 1
    res = Array(Float64, nReps, nVals)
    for ii=1:nReps
        res[ii, :] = x[ii]
    end
    return res
end

function fastVcat(x::Array{DataFrame, 1})
    # preallocation
    (nReps, nVals) = (size(x, 1), size(x[1], 2))
    res = Array(Float64, nReps, nVals)
    for ii=1:nReps
        res[ii, :] = array(x[ii])
    end
    return composeDataFrame(res, names(x[1]))
end

function fastVcat{T}(x::Array{Timematr{T}, 1})
    # preallocation
    (nReps, nVals) = (size(x, 1), size(x[1], 2))
    res = Array(Float64, nReps, nVals)
    indices = Array(T, nReps)
    for ii=1:nReps
        res[ii, :] = core(x[ii])
        indices[ii] = idx(x[ii])[1]
    end
    return Timematr(res, names(x[1]), indices)
end

######################################
## aggregating to lower frequencies ##
######################################

function aggrRets(tm::Timematr; freq = "monthly",
                  logRet = true,
                  percent = true)
    ## aggregate returns to lower frequency

    (nObs, nAss) = size(tm)
    valsArr = core(tm)

    ## define aggregation function
    aggrFun = x -> x
    if logRet == true
        if percent == true
            aggrFun = x -> sum(x, 1)
        end
    elseif logRet == false
        if percent == true
            aggrFun = x -> (prod(1 + x./100, 1) - 1)*100
        else
            aggrFun = x -> (prod(1 + x, 1) - 1)
        end
    end

    # assign equal aggregation ID for days within same period
    aggrId = Array(Float64, nObs)
    nUnique = 1
    nAggrPeriods = 1 # current number of aggregation periods
    lastDayOfPeriod = []

    dayAsAggregationPeriod = []

    if freq == "monthly"
        ## get year and month for each day
        dayAsAggregationPeriod =
            [(year(dateEntry), month(dateEntry)) for dateEntry in
             idx(tm)]
    elseif freq == "yearly"
        ## get year for each day
        dayAsAggregationPeriod =
            [year(dateEntry) for dateEntry in idx(tm)]
    end

    for ii=1:nObs
        aggrId[ii] = nUnique
        if ii < nObs
            if dayAsAggregationPeriod[ii] != dayAsAggregationPeriod[ii+1]
                # next day is of next aggregation period
                nUnique = nUnique + 1
                lastDayOfPeriod = [lastDayOfPeriod; ii]
            end
        end
    end
    lastDayOfPeriod = [0; lastDayOfPeriod; nObs]
    nPeriods = length(lastDayOfPeriod)-1

   # for each period, get aggregated values
    aggrVals = Array(Float64, nPeriods, nAss)
    origDates = idx(tm)
    newDates = origDates[lastDayOfPeriod[2:end]]

    for ii=1:nPeriods
        # get values
        currPeriod =
            valsArr[(lastDayOfPeriod[ii]+1):lastDayOfPeriod[ii+1], :]
        # apply aggregation function
        aggrVals[ii, :] = aggrFun(currPeriod)
    end

    # put everything together in timematr
    df = composeDataFrame(aggrVals, names(tm))
    tmNew = Timematr(df, newDates)
end


text:
# ---
# title: 658. Find K Closest Elements
# id: problem658
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Binary Search
# link: <https://leetcode.com/problems/find-k-closest-elements/description/>
# hidden: true
# ---
#
# Given a **sorted** integer array `arr`, two integers `k` and `x`, return the
# `k` closest integers to `x` in the array. The result should also be sorted in
# ascending order.
#
# An integer `a` is closer to `x` than an integer `b` if:
#
#   * `|a - x| < |b - x|`, or
#   * `|a - x| == |b - x|` and `a < b`
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [1,2,3,4,5], k = 4, x = 3
#     Output: [1,2,3,4]
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,2,3,4,5], k = 4, x = -1
#     Output: [1,2,3,4]
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= arr.length`
#   * `1 <= arr.length <= 104`
#   * Absolute value of elements in the array and `x` will not exceed `104`
#
#
## @lc code=start
using LeetCode

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


text:
using KnuthMorrisPratt
using Test

@testset "KMP_search" begin
    test_sample1 = KMP_search("ABC ABCDAB ABCDABCDABDE", "ABCDABD")
    test_sample2 = KMP_search("babcbabcabcaabcabcabcacabc", "abcabcacab")
    test_sample3 = KMP_search("ABAB CDCDAB ABDE", "AB")
    @test test_sample1 == 16
    @test test_sample2 == 16
    @test test_sample3 == 1
end

@testset "generate_KMP_table" begin
    test_sample1 = KnuthMorrisPratt.generate_KMP_table("ABCDABD")
    test_sample2 = KnuthMorrisPratt.generate_KMP_table("AAAAAAA")
    test_sample3 = KnuthMorrisPratt.generate_KMP_table("abcabcacab")
    test_len1 = KnuthMorrisPratt.generate_KMP_table("a")
    test_len2 = KnuthMorrisPratt.generate_KMP_table("ab")
    @test test_sample1 == [0, 1, 1, 1, 1, 2, 3]
    @test test_sample2 == [0, 1, 2, 3, 4, 5, 6]
    @test test_sample3 == [0, 1, 1, 1, 2, 3, 4, 5, 1, 2]
    @test test_len1 == [0]
    @test test_len2 == [0, 1]
end


text:
function rock_temperature(S, a, e)
    T = 0

    # Your code goes here!

    return T
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
# Nonlinear Projection
## Helper 1 - construct model around the units
function embedUnits!(ena, col, knn, weight1, seed)

    ### Seed
    Random.seed!(seed)

    ### Prepare data
    colVals = Vector{Float64}(ena.metadata[!, col])
    colVals = colVals .- minimum(colVals)
    colVals /= maximum(colVals)
    accumX = transpose(Matrix{Float64}(hcat(colVals, ena.accumModel[!, ena.networkModel[!, :relationship]])))

    ### Model the units
    weightn = (1-weight1) / nrow(ena.networkModel)
    weights = [weight1, [weightn for row in eachrow(ena.networkModel)]...]
    metric = WeightedEuclidean(weights/sum(weights))
    model = UMAP_(accumX; n_neighbors=knn, min_dist=0.0000000001, metric=metric)

    ### Add to accumModel
    ena.accumModel[!, :umap_x] = model.embedding[1, :]
    ena.accumModel[!, :umap_y] = model.embedding[2, :]

    ### Unseed and return model
    Random.seed!(Dates.value(now()))
    return model
end

## Helper 2 - add the networks to the model
function embedNetwork!(ena, model, seed)

    ### Seed
    Random.seed!(seed)

    ### Prepare data
    lineX = nothing
    for (i, networkRow) in enumerate(eachrow(ena.networkModel))
        j, k = ena.relationshipMap[networkRow[:relationship]]
        pointA = Vector{Float64}([-9999, ena.codeModel[j, ena.networkModel[!, :relationship]]...])
        pointB = Vector{Float64}([-9999, ena.codeModel[k, ena.networkModel[!, :relationship]]...])
        point = 0.5*pointA + 0.5*pointB
        if isnothing(lineX)
            lineX = point
        else
            lineX = hcat(lineX, point)
        end
    end

    ### Project network components into the embedding
    weightn = 1 / nrow(ena.networkModel)
    weights = [0, [weightn for row in eachrow(ena.networkModel)]...]
    metric = WeightedEuclidean(weights/sum(weights))
    networkEmbedding = UMAP.transform(model, lineX; n_neighbors=1, min_dist=0.0000000001, metric=metric)

    ### Add to networkModel
    ena.networkModel[!, :umap_x] = networkEmbedding[1, :]
    ena.networkModel[!, :umap_y] = networkEmbedding[2, :]

    ### Add to codeModel
    #### Only care about non-empty networks
    nonZeroRows = map(eachrow(ena.accumModel)) do row
        if sum(row[ena.networkModel[!, :relationship]]) == 0
            return false
        else
            return true
        end
    end

    #### Regression model for placing the code dots
    X = Matrix{Float64}(zeros(nrow(ena.accumModel[nonZeroRows, :]), nrow(ena.codeModel)))
    for (i, unitRow) in enumerate(eachrow(ena.accumModel[nonZeroRows, :]))
        for r in keys(ena.relationshipMap)
            a, b = ena.relationshipMap[r]
            X[i, a] += unitRow[r] / 2
            X[i, b] += unitRow[r] / 2
        end
    end

    X = (transpose(X) * X)^-1 * transpose(X)

    #### Running regressions
    meanX = mean(ena.accumModel[nonZeroRows, :umap_x])
    meanY = mean(ena.accumModel[nonZeroRows, :umap_y])
    ena.codeModel[!, :umap_x] = X * Vector{Float64}(ena.accumModel[nonZeroRows, :umap_x] .- meanX)
    ena.codeModel[!, :umap_y] = X * Vector{Float64}(ena.accumModel[nonZeroRows, :umap_y] .- meanY)

    ### Unseed
    Random.seed!(Dates.value(now()))
end

## Helper 3 - plot UMAP+ENA projections
function plotUMAP(ena, colorMap, col; group=nothing, colormode=:spectral, lineSize=8, codeSize=8, unitSize=4, xlims=(-1.5, 1.5), ylims=(-1.5, 1.5))

    ### Empty plot
    p = plot(; size=(800,800))
    xticks!(p, [-1, 1])
    yticks!(p, [-1, 1])
    xlims!(p, xlims)
    ylims!(p, ylims)

    ### Hide empty networks, and optionally focus on just one group
    nonZeroRows = map(eachrow(ena.accumModel)) do row
        if sum(row[ena.networkModel[!, :relationship]]) == 0
            return false
        else
            return true
        end
    end

    displayRows = nonZeroRows
    if !isnothing(group)
        displayRows = map(enumerate(eachrow(ena.metadata))) do (i, row)
            if nonZeroRows[i] && row[:LABEL] == group
                return true
            else
                return false
            end
        end
    end

    displayAccums = ena.accumModel[displayRows, :]

    ### Center and scale points
    meanX = mean(ena.accumModel[nonZeroRows, :umap_x])
    meanY = mean(ena.accumModel[nonZeroRows, :umap_y])
    scaleX = 1 / maximum(abs.(ena.accumModel[nonZeroRows, :umap_x] .- meanX))
    scaleY = 1 / maximum(abs.(ena.accumModel[nonZeroRows, :umap_y] .- meanY))
    xs = Vector{Float64}(displayAccums[!, :umap_x])
    xs = xs .- meanX
    xs *= scaleX
    ys = Vector{Float64}(displayAccums[!, :umap_y])
    ys = ys .- meanY
    ys *= scaleY

    ### Finding network component sizes
    allLineWidths = map(eachrow(ena.networkModel)) do networkRow
        return sum(ena.accumModel[!, networkRow[:relationship]])
    end

    allCodeWidths = zeros(nrow(ena.codeModel))
    for (i, networkRow) in enumerate(eachrow(ena.networkModel))
        j, k = ena.relationshipMap[networkRow[:relationship]]
        allCodeWidths[j] += allLineWidths[i]
        allCodeWidths[k] += allLineWidths[i]
    end

    lineWidths = map(eachrow(ena.networkModel)) do networkRow
        return sum(displayAccums[!, networkRow[:relationship]])
    end

    codeWidths = zeros(nrow(ena.codeModel))
    for (i, networkRow) in enumerate(eachrow(ena.networkModel))
        j, k = ena.relationshipMap[networkRow[:relationship]]
        codeWidths[j] += lineWidths[i]
        codeWidths[k] += lineWidths[i]
    end

    ### Finding node positions
    codeXs = ena.codeModel[!, :umap_x] * scaleX
    codeYs = ena.codeModel[!, :umap_y] * scaleY

    ### Draw the network
    lineWidths *= lineSize / maximum(allLineWidths)
    codeWidths *= codeSize / maximum(codeWidths)
    for (i, networkRow) in enumerate(eachrow(ena.networkModel))
        j, k = ena.relationshipMap[networkRow[:relationship]]
        pointA = [codeXs[j], codeYs[j]]
        pointB = [codeXs[k], codeYs[k]]
        pointT = [(networkRow[:umap_x] - meanX) * scaleX, (networkRow[:umap_y] - meanY) * scaleY]
        points = hcat(pointA, pointT, pointT, pointT, pointB)
        plot!(p,
            points[1, :],
            points[2, :],
            label=nothing,
            seriestype=:curves,
            linewidth=lineWidths[i],
            linecolor=:grey)
    end

    labels = map(label->text(label, :top, default(:xtickfontsize)), ena.codeModel[!, :code])
    plot!(p,
        codeXs,
        codeYs,
        label=nothing,
        seriestype=:scatter,
        series_annotations=labels,
        markershape=:circle,
        markersize=codeWidths,
        markercolor=:white,
        markerstrokecolor=:black
    )

    ### Draw the units
    colVals = Vector{Float64}(ena.metadata[displayRows, col])
    colVals = colVals .- minimum(colVals)
    colVals /= maximum(colVals)
    colors = map(zip(eachrow(ena.metadata[displayRows, :]), colVals)) do (row, colVal)
        if colormode == :spectral
            return HSL(colVal*240, 1, 0.5)
        elseif colormode == :label
            return colorMap[row[:LABEL]]
        else
            error("Unrecognized color mode $(colormode)")
        end
    end

    plot!(p,
        (displayAccums[!, :umap_x] .- meanX) * scaleX,
        (displayAccums[!, :umap_y] .- meanY) * scaleY,
        label=nothing,
        seriestype=:scatter,
        markersize=unitSize,
        markercolor=colors,
        markerstrokecolor=colors
    )

    ### Draw the legend
    if colormode == :label
        if isnothing(group)
            for mylabel in sort(unique(ena.metadata[!, :LABEL]))
                color = colorMap[mylabel]
                plot!(p, [-9999], [-9999],
                    label=mylabel,
                    seriestype=:scatter,
                    markersize=4,
                    markercolor=color,
                    markerstrokecolor=color
                )
            end
        else
            color = colorMap[group]
            plot!(p, [-9999], [-9999],
                label=group,
                seriestype=:scatter,
                markersize=4,
                markercolor=color,
                markerstrokecolor=color
            )
        end
    end

    ### Done, return
    return p
end

# Clustering
## Helper 4 - detect clusters from umap_x and umap_y positions
function autocluster!(ena, data, colorMap, epsval, min_cluster_size, min_neighbors)

    ### Prepare data
    X = Matrix{Float64}(ena.accumModel[!, [:umap_x, :umap_y]])
    X = Matrix{Float64}(transpose(X))

    ### Using DBSCAN to detect labels
    results = dbscan(X, epsval; min_cluster_size=min_cluster_size, min_neighbors=min_neighbors)

    ### Applying the labels
    ena.metadata[!, :LABEL] .= "No Label"
    for (i, cluster) in enumerate(results)
        mylabel = "#$(i)"
        # mylabel = "Auto Cluster #$(i)"
        ena.metadata[cluster.core_indices, :LABEL] .= mylabel
        ena.metadata[cluster.boundary_indices, :LABEL] .= mylabel
        if !haskey(colorMap, mylabel)
            if i <= length(EpistemicNetworkAnalysis.DEFAULT_EXTRA_COLORS)
                colorMap[mylabel] = EpistemicNetworkAnalysis.DEFAULT_EXTRA_COLORS[i]
            else
                colorMap[mylabel] = RGB(rand(0.4:0.01:0.9), rand(0.4:0.01:0.9), rand(0.4:0.01:0.9))
            end
        end
    end

    ### Saving labels back into original dataframe
    data[!, :LABEL] .= "No Label"
    for row in eachrow(ena.metadata)
        dataRows = [join(dataRow[ena.units], ".") == row[:ENA_UNIT] for dataRow in eachrow(data)]
        data[dataRows, :LABEL] .= row[:LABEL]
    end
end

# Descriptive Statistics
## Helper 5 - create df containing upper and lower bound for each labeled group
function findLabelBounds(ena, col)
    labels = String[]
    sizes = Any[]
    starts = Any[]
    ends = Any[]
    for label in sort(unique(ena.metadata[!, :LABEL]))
        if label != "No Label"
            labelRows = ena.metadata[!, :LABEL] .== label
            push!(labels, label)
            push!(sizes, length(ena.metadata[labelRows, col]))
            push!(starts, first(ena.metadata[labelRows, col]))
            push!(ends, last(ena.metadata[labelRows, col]))
        end
    end

    return DataFrame(
        :LABEL => labels,
        :N => sizes,
        :LowerBound => starts,
        :UpperBound => ends
    )
end

## Helper 6 - plot cdfs so we can see how the groups overlap
function plotCDFs(ena, col, colorMap)
    dayLabelMap = Dict(row[col] => row[:LABEL] for row in eachrow(ena.metadata))
    labelCounts = Dict(label => 0 for label in unique(ena.metadata[!, :LABEL]))
    labelPrevX = Dict(label => 0 for label in unique(ena.metadata[!, :LABEL]))
    p = plot(; size=(800,800))
    for (x, day) in enumerate(sort(ena.metadata[!, col]))
        label = dayLabelMap[day]
        if label != "No Label"
            plot!(p,
                [labelPrevX[label], x],
                [labelCounts[label], labelCounts[label]],
                label=nothing,
                seriestype=:line,
                linecolor=colorMap[label])

            plot!(p,
                [x, x],
                [labelCounts[label], labelCounts[label] + 1],
                label=nothing,
                seriestype=:line,
                linecolor=colorMap[label])

            labelCounts[label] += 1
            labelPrevX[label] = x
        end
    end

    return p
end

# Preprocessing
## Helper 7 -
function derivedCode!(data, newCol, oldCols...)
    data[!, newCol] = ones(nrow(data))
    for col in oldCols
        data[!, newCol] = data[!, newCol] .* (1 .- data[!, col])
    end

    data[!, newCol] = 1 .- data[!, newCol]
end