text:
# default labels and starting values for simulation
const labels = [L"V" L"N" L"Ca" L"Ca_{er}" L"Ca_m" "ADP" "F6P" "FBP"];
const y0 = [-60; 0; 0.1; 185; 100; 780; 60; 40];
# y0_stat produces a stationary solution
const y0_stat = [ -60.486843436763024
                    0.0001367295835481462
                    0.06376844044216665
                  127.60067923018624
                   25.57114498828348
                  809.2663464858219
                    9.141692725854208
                    4.432631767186038e-6 ];

# defines a struct with default values
# A new system can be instancated with keyword arguments
#
# Example
# System(t = 120) -> Struct with default parameters but 120 min length
#
@with_kw struct System
    ode_func :: Function = sys
    y0 :: AbstractVecOrMat{Float64} = y0
    time :: Integer = 60
    plot_args :: Dict = Dict()
    params :: Dict{String,Any} = deepcopy(params)
    plot_params :: AbstractVector{Integer} = [1,2,3,4,5,6,7,8]
    change_params :: Dict{String,Number} = Dict()
    meds :: Vector{Meds} = []
end

# This function attempts to solve the ode system.
# If an DomainError exception occures it retries with
# increased accuracy. If an unexpted exception is thrown
# it is printed to stdout
# Refere to the DifferentialEquations.jl docs
# https://diffeq.sciml.ai/stable/tutorials/ode_example/
#trysolve(system :: System, callback, iteration) = trysolve(to_named_tuple(trysolve), callback, iteration)
function trysolve(system, callback, iteration)
    tol = 1e-1^(iteration*6)
    backup = deepcopy(system) # Create copy for further iterations.
			      # This is necessary because the Med and ExpMed objects are destroyed
    problem = ODEProblem(
        system.ode_func,
        system.y0,
        (0.0, (
            system.time < 1000)
                ? system.time * 6000
                : system.time),
        system.params)
    if tol <= 1e-19 # maximum tolerance
        return nothing
    end
    try
	    # the following callback is called every 100 timesteps and checks for medications to apply
        cb = (callback) ? PeriodicCallback(parse_medications, 100) : nothing
	    # https://diffeq.sciml.ai/stable/basics/common_solver_opts/
        return solve(problem,
                Tsit5(),
                saveat   = 1,
                maxiters = 1e6,
                reltol   = tol,
                abstol   = tol,
                callback = cb
        )
    catch err
        if isa(err, DomainError)
            return trysolve(backup, callback, iteration+1)
        else
            println(err)
        end
    end
end

function add_meds_to_plot(p, meds :: Vector{Union{Activa,ExpMed,Med}}, matrix)
    plot_max = findmax(matrix)[1][1]
    p
    println("ok: " * string(plot_max))
    for (i, m) ∈ enumerate(meds)
            plot!(
                [m.time, m.time + m.duration],
                [plot_max + 20 * i, plot_max + 20 * i],
                label = m.name,
                linestyle = :dash
                )
            annotate!(m.time + m.duration / 2, plot_max + 20 * (i + 0.4), string(m.dose))
            #println("med annotations cannot be added with legacy Med or ExpMed structs")
    end
    return p
end

function gen_solution_plot(solution :: ODESolution, system :: System)
    gr()
    # apply scaling factors to fit all parameters in one plot,
    # since they are numerically unreliable anyways
    matrix = scale_solution_columns(solution)
    p = plot(
        solution.t/6000,
        matrix[system.plot_params];
        label = hcat(labels[system.plot_params]...),
        format=:png,
        fontfamily = "Computer Modern",
        titlefontsize  = 14,
        guidefontsize  = 14,
        tickfontsize   = 14,
        legendfontsize = 13,
        framestyle = :default,
        width=1.2,
        size = (900,600),
        palette = :tab10,
        grid = false,
        legend = (1.07,0.8),
        rightmargin = 42mm,
        system.plot_args...
        )

    print(p)
    p = add_meds_to_plot(p, system.meds, matrix)
    return p
end

# Takes a settings object (NamedTuple) and changes to parameters for the simulation
# accordingly. It returns DiffEq Solution object, the time series for the 8 values
# and a plot fo the solution.

function simulate(system :: System; iteration = 1)
    for (key,val) ∈ system.change_params
        system.params[key] = val;
    end
    system.params["meds"] = system.meds

    @time solution = trysolve(system, system.meds |> length > 0, iteration)

    (solution === nothing) ? nothing : (
        solution,
        rowtocolumn(solution.u),
        gen_solution_plot(solution,system)
    )
end

# maintina legacy compatablity with old named tuple args
simulate(system) = simulate(System(system...))

# can be called to put multiple vales for one parameter
# and produce one solution for each
# If multiple threads are available the simulations
# will be run in parallel

function loopvals(key::String, vals::Array, system)
    sols = Dict()
    sols_mat = Dict()
    plots = Dict()

    @threads for val ∈ vals
        local s = deepcopy(system)
        s.params[key] = val
        s.plot_args[:title] = key * " = " * string(val)
        try
            sol, mat, p = simulate(s)
            sols[val] = sol
            sols_mat[val] = mat
            plots[val] = p
        catch
            println("simulate " * key * " = " * string(val) * " failed")
        end
    end
    acc_plot = plot(
        sortedvalues(plots)...,
        dpi=300,
        label = hcat(labels[system.plot_params]...),
    )
    return (sols, sols_mat, plots, acc_plot)
end

# same as above but for two chaning variables
function loopvals(key1::String, vals1::Vector, key2::String, vals2::Vector, system)
    local sols::Dict{Any, Dict}     = Dict()
    local sols_mat::Dict{Any, Dict} = Dict()
    local plots::Dict{Any, Dict}    = Dict()

    @threads for val1 ∈ vals1
    #for val1 ∈ vals1
            sols[val1] = Dict()
            sols_mat[val1] = Dict()
            plots[val1] = Dict()
        @threads for val2 ∈ vals2
        #for val2 ∈ vals2
            local s = deepcopy(system)
            s.params[key1] = val1
            s.params[key2] = val2
            s.plot_args[:title] = key1 * " = " * string(val1) * ", " * key2 * " = " * string(val2)
            try
                sol, mat, p = simulate(s)
                sols[val1][val2] = sol
                sols_mat[val1][val2] = mat
                plots[val1][val2] = p
            catch
                println("simulate " * key1 * " = " * string(val1) * " " * key2 * " = " * string(val2) *" failed")
            end
        end
    end
    acc_plot = plot(
        sortedvalues(plots)...,
        dpi=300,
        label = hcat(labels[system.plot_params]...),
    )

    return (sols, sols_mat, plots, acc_plot)
end

# use multiple dispatch to avoid calls to loopvals
simulate(system :: System, key1::String, vals1::Vector, key2::String, vals2::Vector) = loopvals(key1, vals1, key2, vals2, system)
simulate(system :: System, key1::String, vals1::Vector) = loopvals(key1, vals1, system)


text:
using Test
using NSGAII

include("crowdingDistanceAssignement.jl")
include("evolution.jl")
include("fastNonDominatedSort.jl")

ext: jl


text:
function win_installer(builddir; name = "nothing",
				license = "$JULIA_HOME/../License.md")

	# check = success(`makensis`)
	# !check && throw(ErrorException("NSIS not found in path. Exiting."))

	nsis_commands = """
	# set the name of the installer
	Outfile "$(name)_Installer.exe"

	# Default install directory
	InstallDir "\$LOCALAPPDATA"

	Page license
	Page directory
	Page instfiles

	LicenseData "$license"

	# create a default section.
	Section "Install"

		SetOutPath "$(joinpath("\$INSTDIR", name))"
		File /nonfatal /a /r "$builddir"

		CreateShortcut "$(joinpath("\$INSTDIR", name, "$name.lnk"))" "$(joinpath(builddir, "core", "blink.exe"))"

	SectionEnd
	"""

	@info "Creating installer at $builddir"
	nsis_file = joinpath(builddir, "..", "$name.nsi")
	open(nsis_file, "w") do f
		write(f, nsis_commands)
	end
	run(`makensis $nsis_file`)

	@info "Created installer successfully."

end


text:
abstract type ParallelSortAlgorithm <: Base.Sort.Algorithm end
(alg::ParallelSortAlgorithm)(; kw...) = setproperties(alg; kw...)

Base.@kwdef struct ParallelMergeSortAlg{Alg,SmallSize,BaseSize} <: ParallelSortAlgorithm
    smallsort::Alg = Base.Sort.DEFAULT_STABLE
    smallsize::SmallSize = nothing  # lazily determined
    basesize::BaseSize = nothing  # lazily determined
end

function Base.sort!(
    v::AbstractVector,
    lo::Integer,
    hi::Integer,
    a::ParallelMergeSortAlg,
    o::Ordering,
)
    if a.basesize === nothing
        a = @set a.basesize = default_basesize(hi - lo + 1)
    end
    if a.smallsize === nothing
        a = @set a.smallsize = a.basesize
    end
    _mergesort!(view(v, lo:hi), a, o)
    return v
end

function mergesorted!(dest, left, right, order, basesize)
    @assert length(dest) == length(left) + length(right)
    if length(left) <= basesize || length(right) <= basesize
        return mergesorted_basecase!(dest, left, right, order)
    end
    if length(left) < length(right)
        c, d = halve(right)
        a, b = halveat(left, searchsortedlast(left, last(c), order) + 1)
    else
        a, b = halve(left)
        c, d = halveat(right, searchsortedfirst(right, first(b), order))
    end
    # length(c) > 0 && length(b) > 0 && @assert Base.lt(order, last(c), first(b))  # stable sort
    ac, bd = halve(dest, length(a) + length(c))
    task = let ac = ac, c = c, a = a
        @spawn mergesorted!(ac, a, c, order, basesize)
    end
    mergesorted!(bd, b, d, order, basesize)
    wait(task)
    return dest
end

@inline function _copyto!(ys, xs)
    for i in eachindex(ys, xs)
        @inbounds ys[i] = xs[i]
    end
    return ys
end

function mergesorted_basecase!(dest::D, left, right, order) where {D}
    # @assert issorted(left; order = order)
    # @assert issorted(right; order = order)
    if isempty(left)
        _copyto!(dest, right)
        return dest
    elseif isempty(right)
        _copyto!(dest, left)
        return dest
    end
    @assert length(dest) == length(left) + length(right)
    i = firstindex(left)
    j = firstindex(right)
    k = firstindex(dest)
    a = @inbounds left[i]
    b = @inbounds right[j]
    while true
        if Base.lt(order, b, a)
            @inbounds dest[k] = b
            j += 1
            k += 1
            if j > lastindex(right)
                _copyto!((@view dest[k:end]), (@view left[i:end]))
                break
            end
            b = @inbounds right[j]
        else
            @inbounds dest[k] = a
            i += 1
            k += 1
            if i > lastindex(left)
                _copyto!((@view dest[k:end]), (@view right[j:end]))
                break
            end
            a = @inbounds left[i]
        end
    end
    # @assert issorted(dest; order = order)
    return dest
end

halve(arr::AbstractArray, mid = length(arr) ÷ 2) = halveat(arr, firstindex(arr) + mid)

function halveat(arr::AbstractArray, i)
    left = @view arr[firstindex(arr):i-1]
    right = @view arr[i:end]
    return (left, right)
end

function _mergesort!(xs, alg, order, tmp = nothing)
    if length(xs) <= alg.smallsize
        sort!(xs, alg.smallsort, order)
        return xs
    end
    left, right = halve(xs)
    left_tmp, right_tmp = halve(tmp === nothing ? similar(xs) : tmp)
    task = @spawn _mergesort!(left, alg, order, left_tmp)
    _mergesort!(right, alg, order, right_tmp)
    wait(task)
    mergesorted!(
        xs,
        _copyto!(left_tmp, left),
        _copyto!(right_tmp, right),
        order,
        alg.basesize,
    )
    return xs
end

ParallelSortAlgorithm(alg::ParallelSortAlgorithm) = alg
ParallelSortAlgorithm(::typeof(MergeSort)) = ParallelMergeSortAlg()
ParallelSortAlgorithm(::typeof(QuickSort)) = ParallelQuickSortAlg()

# bypass `fpsort!`
_sort!(v, a, o) = sort!(v, first(axes(v, 1)), last(axes(v, 1)), a, o)
# TODO: Implement fast `NaN` support as done in `fpsort!`.
# TODO: Use `_sort!` instead of `sort!` for small case (once the `NaN`
# support is implemented).

defalg(v::AbstractArray) = ParallelMergeSortAlg()
defalg(v::AbstractArray{<:Union{Number, Missing}}) = ParallelQuickSortAlg()

"""
    ThreadsX.sort(xs; [smallsort, smallsize, basesize, alg, lt, by, rev, order])

See also [`ThreadsX.sort!`](@ref).
"""
ThreadsX.sort(xs; kwargs...) = ThreadsX.sort!(Base.copymutable(xs); kwargs...)

function ThreadsX.sort!(
    xs;
    smallsort = nothing,
    smallsize = nothing,
    basesize = nothing,
    alg::Base.Sort.Algorithm = defalg(xs),
    lt = isless,
    by = identity,
    rev::Union{Bool,Nothing} = nothing,
    order::Base.Ordering = Base.Forward,
)
    alg = ParallelSortAlgorithm(alg)
    if basesize !== nothing
        alg = @set alg.basesize = basesize
    end
    if smallsort !== nothing
        alg = @set alg.smallsort = smallsort
    end
    if smallsize !== nothing
        alg = @set alg.smallsize = smallsize
    end
    ordr = Base.ord(lt, by, rev, order)
    if maybe_counting_sort!(xs, ordr)
        return xs
    end
    return _sort!(xs, alg, ordr)
end


text:
module Scat

using FFTW

include("./filters/filter.jl")
include("./filters/filterBank.jl")
include("./transform/scatteringTransform.jl")
include("./transform/cwt.jl")
include("./transform/scattered.jl")
include("./transform/scatter1d.jl")
end


text:
# p39.m - eigenmodes of biharmonic on a square with clamped BCs
#         (compare p38.m)

# Construct spectral approximation to biharmonic operator:
N = 17; (D,x) = cheb(N); D2 = D^2; D2 = D2[2:N,2:N];
S = diagm([0; 1 ./(1-x[2:N].^2); 0]);
D4 = (diagm(1-x.^2)*D^4 - 8*diagm(x)*D^3 - 12*D^2)*S;
D4 = D4[2:N,2:N]; I = eye(N-1);
L = kron(I,D4) + kron(D4,I) + 2*kron(D2,I)*kron(I,D2);

# Find and plot 25 eigenmodes:
(Lam,V) = eig(-L); Lam = -real(Lam);
ii = sortperm(Lam)[1:25]; Lam = Lam[ii]; V = real(V[:,ii]);
Lam = sqrt.(Lam/Lam[1]);
y = x; xxx = yyy = -1:.01:1;
sq = [1+1im,-1+1im,-1-1im,1-1im,1+1im]; clf();
for i = 1:25
    uu = zeros(N+1,N+1); uu[2:N,2:N] = reshape(V[:,i],N-1,N-1);
    subplot(5,5,i); plot(real(sq),imag(sq));
    s = interpolate((x[end:-1:1],y[end:-1:1]),reduce(flipdim,uu,1:2),Gridded(Linear()));
    contour(xxx,yyy,s[xxx,yyy],levels=[0],color="k"); axis("square");
    axis(1.25*[-1,1,-1,1]); axis("off");
    text(-.3,1.15,"$(signif(Lam[i],5))",fontsize=7);
end


text:
######
# Output TODO
######

### Summaries

# "Summary statistics of a run."
# function summarize_run(thetas::Matrix; theta_true=similar(thetas,0), names=["$i" for i=1:size(thetas,1)],
#                        eff_samples=nothing, mode=nothing)
#     nt = size(thetas,1)
#     ns = size(thetas,2)
#     if length(theta_true)>0
#         cols = eff_samples==nothing ? Any[[], [],[],[],[],[]] : Any[[], [],[],[],[],[],[]]
#         header = eff_samples==nothing ? [:var, :err, :median, :mean, :mode, :std] : [:var, :err, :median, :mean, :mode, :std, :eff_samples]

#         for i=1:nt
#             push!(cols[1], Symbol(names[i]))
#             t = theta_true[i]
#             m = median(view(thetas,i,1:ns))
#             push!(cols[2], abs(t-m))
#             push!(cols[3], median(view(thetas,i,1:ns)))
#             push!(cols[4], mean(view(thetas,i,1:ns)))
#             push!(cols[5], mod==nothing ? nothing : mode[i])
#             push!(cols[6], std(view(thetas,i,1:ns)))
#             eff_samples!=nothing && push!(cols[7], eff_samples[i])
#         end
#     else
#         cols = eff_samples==nothing ? Any[[], [],[],[],[]] : Any[[], [],[],[],[], []]
#         header = eff_samples==nothing ? [:var, :median, :mean, :mode, :std] : [:var, :median, :mean, :mode, :std, :eff_samples]
#         for i=1:nt
#             push!(cols[1], Symbol(names[i]))
#             push!(cols[2], median(view(thetas,i,1:ns)))
#             push!(cols[3], mean(view(thetas,i,1:ns)))
#             push!(cols[4], mode==nothing ? nothing : mode[i])
#             push!(cols[5], std(view(thetas,i,1:ns)))
#             eff_samples!=nothing && push!(cols[6], eff_samples[i])
#         end

#     end
#     #DataFrame(cols, header)
# end

# "Print result summary"
# function print_results(thetas::Matrix, accept_ratio, eff_samples=nothing, mode=nothing; theta_true=similar(thetas,0),
#                        names=["$i" for i=1:size(thetas,1)],
#                        maxvar=45,
#                        title="")
#     println(title)
#     println("Ratio of accepted/total steps: $accept_ratio\n")
#     out = summarize_run(thetas, theta_true=theta_true, names=names, eff_samples=eff_samples)
#     show(out[1:min(size(out,1),maxvar),:])
#     nothing
# end

### Evaluation

# """
#     evaluate_convergence(thetas1, thetas2; indices=1:size(thetas)[1], walkernr=1)

# Evaluates:
# - R̂ (potential scale reduction): should be <1.1
# - total effective sample size of all chains combined, and
# - the average thinning factor

# See Gelman etal 2014, page 281-287

# Note that https://dfm.io/posts/autocorr/ says "you should not compute
# the G–R statistic using multiple chains in the same emcee ensemble because
# the chains are not independent!"  Thus this needs input from two separate emcee runs.

# Example

#     thetas1, _, accept_ratio, logposteriors = emcee(x->sum(-x.^2), ([1.0,2, 5], 0.1))
#     thetas2, _, accept_ratio, logposteriors = emcee(x->sum(-x.^2), ([1.0,2, 5], 0.1))
#     Rhat, sample_size, nthin = evaluate_convergence(thetas1, thetas2)
# """
# function evaluate_convergence(thetas... ; indices=1:size(thetas[1])[1],
#                               walkernr=1)
#     # @show [size(t) for t in thetas]
#     # @assert all([size(t)==size(thetas[1]) for t in thetas])
#     Rs = Float64[]
#     sample_size = Float64[]
#     split = size(thetas[1])[2]÷2
#     for i in indices
#         #chains = [thetas1[i,1:split,walkernr], thetas2[i,split+1:2*split,walkernr]]
#         chains = [t[i,1:split,walkernr] for t in thetas]
#         push!(Rs, MCMCDiagnostics.potential_scale_reduction(chains...))
#         push!(sample_size, sum(MCMCDiagnostics.effective_sample_size.(chains)))
#     end
#     nwalkers = length(size(thetas[1]))==3 ? size(thetas[1])[3] : 1
#     nthin = size(thetas[1])[2]*nwalkers/mean(sample_size)
#     return Rs, sample_size, isnan(nthin) ? -1 : round(Int, nthin)
# end

# """
#     int_acorr(thetas::Array{<:Any,3}; c=5, warn=true, warnat=50)

# Returns:
# - Estimated integrated autocorrelation time τ (in number of steps) for each theta
# - Indication whether the estimate has converged as ratio between chain-length and
#   τ.  Should probably be larger than 50 or 100.  (again for each theta)

# τ gives the factor by which the variance of the calculated expectation of y
# is larger than for independent samples.

# In other words:
# - τ is the number of steps that are needed before the chain "forgets" where it started.
# - N/τ is the effective number of samples (N is total number of samples)

# This means that, if you can estimate τ, then you can estimate the number
# of samples that you need to generate to reduce the relative error on your
# target integral to (say) a few percent.

# NOTE: the employed algorithm is not that great at estimating τ
# with short chain lengths.  To get good estimates a chain length greater than
# 50τ is recommended (which may well be longer than what you actually need).
# By default this function will warn if the estimate is likely not converged.

# Optionals:
# - c -- window width to search when calculating autocor (5)
# - warn -- warn if chain length is too small for an accurate τ estimate (true)

# Notes:
# - It is important to remember that f depends on the specific
#   function f(θ). This means that there isn't just one integrated
#   autocorrelation time for a given Markov chain. Instead, you must
#   compute a different τ for any integral you estimate using the samples.
# - works equally with thinned an non-thinned chains
# - whilst it is not recommended by the referenced work, it is probably ok
#   to apply this also to a squashed chain.  Then run it as:
#   `int_acorr(reshape(thetas, (size(thetas,1), size(thetas,2), 1))`

# Ref:
# https://emcee.readthedocs.io/en/latest/tutorials/autocorr/
# Adapted from various bits of the code there.
# """
# function int_acorr(thetas::Array{<:Any,3}; c=5, warn=true, warnat=50)
#     @assert c>1
#     sz = size(thetas)
#     ntheta, nsamples, nchains = size(thetas)
#     out = Float64[]
#     for n = 1:ntheta
#         # mean autocorr of all chains:
#         rho = zeros(nsamples÷2)
#         for cc = 1:nchains
#             rho .+= acor1d(thetas[n,:,cc])
#         end
#         rho ./= nchains
#         # the -1 correct see https://github.com/dfm/emcee/issues/267#issuecomment-477556521
#         taus = 2 * cumsum(rho) - 1
#         window = auto_window(taus, c)
#         push!(out, taus[window])
#     end
#     converged = nsamples./out # probably converged if > 50
#     if warn && any(converged.<warnat)
#         Base.warn("Estimate of integrated autocorrelation likely not accurate!")
#     end
#     if any(isnan.(out)) || any(isnan.(converged))
#         # set both to -1 in this case (a hack)
#         out = out*false - 1
#         converged = converged*false -1
#     end
#     return out, converged
# end

# """
#     eff_samples(thetas::Array{<:Any,3}, c=5)

# Return scalars:
# - total mean number of effective samples
# - suggested thinning step (essentially equal to the mean τ)
# - mean τ-convergence estimate  (as ratio chain-length and τ,
#   should be greater than 50 or 100 or so)
# Return vectors for each θ component
# - total number of effective samples
# - τ auto correlation steps
# - per chain estimates whether τ-estimates are converged

# (the last two output are pass on from int_acorr)

# Example:

#     Neff, τ, conv, Neffs, τs, convs = eff_samples(theta)
# """
# function eff_samples(thetas::Array{<:Any,3}; c=5)
#     acorr, converged = int_acorr(thetas, c=c, warn=false)
#     ns = size(thetas,2)./acorr * size(thetas,3)
#     return (round.(Int, mean(ns)), round.(Int, size(thetas,2)*size(thetas,3)÷mean(ns)), mean(converged),
#             round.(Int, ns), acorr, converged)
# end

# """
#     samples_vs_tau(thetas::Array{<:Any,3})

# Returns samples vs τ as is plotted in
# https://emcee.readthedocs.io/en/latest/tutorials/autocorr/

# To plot do:

#     using Plots
#     N,taus = KissMCMC.samples_vs_tau(thetas);
#     plot(N,taus)
#     # additionally plot the N/50 line
#     plot!(N, N/50, ls=:dash, c=:black)
# """
# function samples_vs_tau(thetas::Array{<:Any,3})
#     N = round(Int, logspace(2,log10(size(thetas,2)), 10))
#     nthin = 1

#     taus = []
#     # converged = []
#     NN = []
#     for n in N
#         if length(1:nthin:n)>3
#             tau, conv = KissMCMC.int_acorr(thetas[:, 1:nthin:n, :], warn=false)
#             push!(taus, tau)
#             # push!(converged, conv)
#             push!(NN,length(1:nthin:n))
#         end
#     end
#     taus = hcat(taus...)'
#     return NN, taus
# end

# """
#     error_of_estimated_mean(thetas::Array{<:Any,3},
#                             Neff = eff_samples(thetas)[4])

# Calculates the error in the estimated mean of each theta.

# Return estimates of:
# - mean(theta)
# - error of mean
# - std(theta)

# Ref:
# 15.4.3 in https://mc-stan.org/docs/2_18/reference-manual/effective-sample-size-section.html
# """
# function error_of_estimated_mean(thetas::Array{<:Any,3},
#                                  Neff = eff_samples(thetas)[4])
#     mtheta = mean(mean(thetas,2),3)[:]
#     stheta = std(std(thetas,2),3)[:]
#     err_of_mean = stheta./sqrt(Neff)
#     return mtheta, err_of_mean, stheta
# end

# ## helper functions
# import Base.DFT
# function acor1d(x::AbstractVector, norm=true)
#     # TODO for speed:
#     #n = DFT.nextprod([2,3,4,5], length(x))

#     # Compute the FFT and then (from that) """
#     hasblob(pdf, theta) = length(pdf(theta))==2

# Test if the log-pdf returns a blob
# """
# hasblob(pdf, theta) = length(pdf(theta))==2
# the auto-correlation function
#     f = DFT.fft(x - mean(x))
#     acf = real.(DFT.ifft(f .* conj(f)))# TODO [1:length(x)]
#     acf ./= 4*length(x)

#     # Optionally normalize
#     if norm
#         acf ./= acf[1]
#     end

#     return acf[1:length(acf)÷2]
# end

# """
#     auto_window(taus, c)

# Gives the ±window over which the autocorr should be integrated
# to get a good estimate.  `c` defaults to 5.
# """
# function auto_window(taus, c)
#     for (i,t) in enumerate(taus)
#         i >= c*t && return i
#     end
#     return length(taus)-1
# end


text:
using Test
using StringParserPEG

include("../examples/calc1.jl")
@test transformed == 9

include("../examples/calc2.jl")
@test ast == 9

@test parse(Grammar("""start => '<<' & *(char) & '>>'; char  => r(.)r"""),"<<abc>>")[3] != nothing

@test parse(Grammar("""start => '<<' & *(char) & '>>'; char  => !('>>') & r(.)r"""),"<<abc>>")[3] == nothing

@test length(parse(Grammar("start => +(r([\\s\\S])r)"),"大家好!")[1].children) == 4

@test StringParserPEG.grammargrammar == Grammar(StringParserPEG.grammargrammar_string)

@test length(parse(Grammar("start => *(r(.*)r)"),"...")[1].children) == 1
@test length(parse(Grammar("start => *(*('.'))"),"...")[1].children) == 1


text:
# module TestAggregation

using Test, Dates, Temporal

@testset "Reducing / Collapsing" begin
    @testset "Periods" begin
        x = TS(rand(7))
        @test length(mondays(x)) == 1
        @test length(tuesdays(x)) == 1
        @test length(wednesdays(x)) == 1
        @test length(thursdays(x)) == 1
        @test length(fridays(x)) == 1
        @test length(saturdays(x)) == 1
        @test length(sundays(x)) == 1
        x = TS(rand(252*5))
        @test all(boy(boy(x).index)[2:end])
        @test all(eoy(eoy(x).index)[1:end-1])
        @test all(boq(boq(x).index)[2:end])
        @test all(eoq(eoq(x).index)[1:end-1])
        @test all(bom(bom(x).index)[2:end])
        @test all(eom(eom(x).index)[1:end-1])
        @test length(bow(x)) == length(mondays(x))
        @test length(eow(x)) == length(sundays(x))
    end
    @testset "Aggregation" begin
        X = TS(rand(252*5,4))
        x = X[:,1]
        m = collapse(x, eoy, fun=mean)
        prev_year = string(year(today())-1)
        @test round(m.values[end], digits=2) == round(mean(x[prev_year]), digits=2)
        m = collapse(x, (t)->eom(t,cal=true), fun=mean)
        prev_month = string(today()-Month(1))[1:7]
        @test mean(x[prev_month]) == m[prev_month].values[end]
        m = apply(X, 1, fun=mean)
        @test m.values[1] == mean(X[1].values)
        m = apply(X, 2, fun=mean)
        @test m[1] == mean(X.values[:,1])
    end
end

# end


text:
using ConstrainedDynamics

# Parameters
joint_axis = [1.0;0.0;0.0]

length1 = 1.0
width, depth = 0.1, 0.1
box = Box(width, depth, length1, length1)

p2 = [0.0;0.0;length1 / 2] # joint connection point

# Initial orientation
ϕ1 = π / 2
q1 = Quaternion(RotX(ϕ1))

# Links
origin = Origin{Float64}()
link1 = Body(box)

# Constraints
joint_between_origin_and_link1 = EqualityConstraint(Revolute(origin, link1, zeros(3), p2, joint_axis))

links = [link1]
constraints = [joint_between_origin_and_link1]
shapes = [box]

mech = Mechanism(origin, links, constraints, shapes = shapes)
setPosition!(mech,origin,link1,p2 = p2,Δq = q1)

simulate!(mech,save = true)
visualize!(mech)


text:
@symbol_func function kappa_P(cur_reactor::AbstractReactor)
  cur_kappa = 1.3

  cur_kappa *= cur_reactor.kappa_95

  cur_kappa
end


text:
immutable Shape
  vertices::AVec
end

"get an array of tuples of points on a circle with radius `r`"
function partialcircle(start_θ, end_θ, n = 20, r=1)
  @compat(Tuple{Float64,Float64})[(r*cos(u),r*sin(u)) for u in linspace(start_θ, end_θ, n)]
end

"interleave 2 vectors into each other (like a zipper's teeth)"
function weave(x,y; ordering = Vector[x,y])
  ret = eltype(x)[]
  done = false
  while !done
    for o in ordering
      try
          push!(ret, shift!(o))
      end
      # try
      #     push!(ret, shift!(y))
      # end
    end
    done = isempty(x) && isempty(y)
  end
  ret
end

"create a star by weaving together points from an outer and inner circle.  `n` is the number of arms"
function makestar(n; offset = -0.5, radius = 1.0)
    z1 = offset * π
    z2 = z1 + π / (n)
    outercircle = partialcircle(z1, z1 + 2π, n+1, radius)
    innercircle = partialcircle(z2, z2 + 2π, n+1, 0.4radius)
    Shape(weave(outercircle, innercircle)[1:end-2])
end

"create a shape by picking points around the unit circle.  `n` is the number of point/sides, `offset` is the starting angle"
function makeshape(n; offset = -0.5, radius = 1.0)
    z = offset * π
    Shape(partialcircle(z, z + 2π, n+1, radius)[1:end-1])
end

function makecross(; offset = -0.5, radius = 1.0)
    z2 = offset * π
    z1 = z2 - π/8
    outercircle = partialcircle(z1, z1 + 2π, 9, radius)
    innercircle = partialcircle(z2, z2 + 2π, 5, 0.5radius)
    Shape(weave(outercircle, innercircle,
                ordering=Vector[outercircle,innercircle,outercircle])[1:end-2])
end

const _shapes = @compat Dict(
    :ellipse    => makeshape(20),
    :rect       => makeshape(4, offset=-0.25),
    :diamond    => makeshape(4),
    :utriangle  => makeshape(3),
    :dtriangle  => makeshape(3, offset=0.5),
    :pentagon   => makeshape(5),
    :hexagon    => makeshape(6),
    :heptagon   => makeshape(7),
    :octagon    => makeshape(8),
    :cross      => makecross(offset=-0.25),
    :xcross     => makecross(),
    :vline      => Shape([(0,1),(0,-1)]),
    :hline      => Shape([(1,0),(-1,0)]),
  )

for n in [4,5,6,7,8]
  _shapes[symbol("star$n")] = makestar(n)
end

# -----------------------------------------------------------------------

immutable Font
  family::AbstractString
  pointsize::Int
  halign::Symbol
  valign::Symbol
  rotation::Float64
  color::Colorant
end

"Create a Font from a list of unordered features"
function font(args...)

  # defaults
  family = "Helvetica"
  pointsize = 14
  halign = :hcenter
  valign = :vcenter
  rotation = 0.0
  color = colorant"black"

  for arg in args
    T = typeof(arg)

    if arg == :center
      halign = :hcenter
      valign = :vcenter
    elseif arg in (:hcenter, :left, :right)
      halign = arg
    elseif arg in (:vcenter, :top, :bottom)
      valign = arg
    elseif T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      catch
        family = string(arg)
      end
    elseif typeof(arg) <: Integer
      pointsize = arg
    elseif typeof(arg) <: Real
      rotation = convert(Float64, arg)
    else
      warn("Unused font arg: $arg ($(typeof(arg)))")
    end
  end

  Font(family, pointsize, halign, valign, rotation, color)
end

"Wrap a string with font info"
immutable PlotText
  str::@compat(AbstractString)
  font::Font
end

function text(str, args...)
  PlotText(string(str), font(args...))
end

# -----------------------------------------------------------------------

immutable Stroke
  width
  color
  alpha
  style
end

function stroke(args...; alpha = nothing)
  # defaults
  # width = 1
  # color = colorant"black"
  # style = :solid
  width = nothing
  color = nothing
  style = nothing

  for arg in args
    T = typeof(arg)

    if arg in _allStyles
      style = arg
    elseif T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      end
    elseif typeof(arg) <: Real
      width = arg
    else
      warn("Unused stroke arg: $arg ($(typeof(arg)))")
    end
  end

  Stroke(width, color, alpha, style)
end

immutable Brush
  size  # fillrange, markersize, or any other sizey attribute
  color
  alpha
end

function brush(args...; alpha = nothing)
  # defaults
  # sz = 1
  # color = colorant"black"
  size = nothing
  color = nothing

  for arg in args
    T = typeof(arg)

    if T <: Colorant
      color = arg
    elseif T <: @compat Union{Symbol,AbstractString}
      try
        color = parse(Colorant, string(arg))
      end
    elseif typeof(arg) <: Real
      size = arg
    else
      warn("Unused brush arg: $arg ($(typeof(arg)))")
    end
  end

  Brush(size, color, alpha)
end

# -----------------------------------------------------------------------

"type which represents z-values for colors and sizes (and anything else that might come up)"
immutable ZValues
  values::Vector{Float64}
  zrange::Tuple{Float64,Float64}
end

function zvalues{T<:Real}(values::AVec{T}, zrange::Tuple{T,T} = (minimum(values), maximum(values)))
  ZValues(collect(float(values)), map(Float64, zrange))
end

# -----------------------------------------------------------------------

"represents a contour or surface mesh"
immutable Surface{M<:AMat}
  # x::AVec
  # y::AVec
  surf::M
end

Surface(f::Function, x, y) = Surface(Float64[f(xi,yi) for xi in x, yi in y])

Base.Array(surf::Surface) = surf.surf

# -----------------------------------------------------------------------

type OHLC{T<:Real}
  open::T
  high::T
  low::T
  close::T
end


text:
# A test case which covers most possible field types.
@eval Discord struct Foo
    a::String
    b::DateTime
    c::Snowflake
    d::Vector{String}
    e::Union{Int, Nothing}
    f::Union{Int, Missing}
    g::Union{Vector{String}, Nothing, Missing}
    h::Union{Foo, Missing}
    abcdefghij::Union{Int, Missing}  # 10 characters.
end
@eval Discord @boilerplate Foo :constructors :docs :lower :merge :mock

# A simple struct with merge.
@eval Discord struct Bar
    id::Int
    x::Int
    y::Union{Int, Missing}
end
@eval Discord @boilerplate Bar :merge

using Discord: Foo, Bar

@testset "Boilerplate" begin
    local f
    kwargs = Dict(
        :a => "foo",
        :b => "2018-10-08T05:20:22.782Z",
        :c => "1234567890",
        :d => ["a", "b", "c"],
        :e => 1,
        :f => 2,
        :g => ["a", "b", "c"],
    )
    kwargs[:h] = copy(kwargs)

    @testset "@constructors" begin
        f = Foo(; kwargs...)
        @test f.a == Foo(kwargs).a
        @test f.a == "foo"
        @test f.b == DateTime(2018, 10, 8, 5, 20, 22, 782)
        @test f.c == 1234567890
        @test f.d == ["a", "b", "c"]
        @test f.e == 1
        @test f.f == 2
        @test f.g == ["a", "b", "c"]
        @test f.h isa Foo && f.h.a == "foo"

        # Set e to nothing, f and g to missing.
        kwargs[:e] = nothing
        delete!(kwargs, :f)
        delete!(kwargs, :g)
        f = Foo(; kwargs...)
        @test f.e === nothing
        @test ismissing(f.f)
        @test ismissing(f.g)

        # Union{T, Nothing, Missing} works too.
        kwargs[:g] = nothing
        f = Foo(; kwargs...)
        @test f.g === nothing
    end

    @testset "@docs" begin
        docs = string(@doc Foo)
        # Variable names get padded to the longest one.
        @test occursin("$(rpad("a", 10)) :: String", docs)
        # Array{T,1} is replaced with Vector{T}.
        @test occursin("Vector{String}", docs)
        # Union{Missing, Nothing, T} is replaced with OptionalNullable{T]}.
        @test occursin("OptionalNullable{Vector{String}}", docs)
        # Union{Missing, T} is replaced with Optional{T}.
        @test occursin("Optional{Int}", docs)
        # Union{Nothing, T} is replaced with Nullable{T}.
        @test occursin("Nullable{Int}", docs)
    end

    @testset "@lower" begin
        kwargs[:b] = kwargs[:h][:b] = round(Int, datetime2unix(f.b))
        kwargs[:c] = kwargs[:h][:c] = snowflake(f.c)
        lowered = JSON.lower(f)
        # The result is always a Dict{Symbol, Any}.
        @test lowered isa Dict{Symbol, Any}
        @test lowered == Dict(
            :a => "foo",
            :b => kwargs[:b],
            :c => kwargs[:c],
            :d => ["a", "b", "c"],
            :e => nothing,
            :g => nothing,
            :h => kwargs[:h],
        )
    end

    @testset "@merge" begin
        # Anything except missing values should be taken from f2.
        f2 = Foo("bar", f.b, f.c, f.d, f.e, f.f, f.g, missing, missing)
        f3 = merge(f, f2)
        @test f3.a == f2.a
        @test f3.h == f.h
    end

    @testset "@mock" begin
        # Just make sure this doesn't error, it should cover most types.
        foreach(mock, subtypes(AbstractEvent))

        # We can also specify fields via keywords.
        m = mock(Message; content="foo")
        @test m.content == "foo"
    end
end


text:
using PicoQuant
using Test
using TestSetExtensions

include("test_utils.jl")

@testset "All the tests" begin
    @includetests ARGS
end


text:
"""
## formm!

This subroutine forms the derivatives of the invariants with respect to
stress in 2- or 3-d. See equation 6.25.

### Function
```julia
formm!(stress, m1, m2, m3)
```

### Arguments
```julia
* stress::Vector{Float64}    : Stress vector, see eq 6.25
* m1::Matrix{Float64}        : m1 matrix
* m2::Matrix{Float64}        : m2 matrix
* m3::Matrix{Float64}        : m3 matrix
```

### Return values
```julia
* m1::Matrix{Float64}        : m1 matrix
* m2::Matrix{Float64}        : m2 matrix
* m3::Matrix{Float64}        : m3 matrix
```
"""
function formm!(stress::Vector{Float64}, m1::Matrix{Float64},
  m2::Matrix{Float64}, m3::Matrix{Float64})
  #
  # This subroutine forms the derivatives of the invariants with respect to
  # stress in 2- or 3-d.
  #
  # Updates: m1, m2, m3
 nst = size(stress, 1)
 if nst == 4
   sx = stress[1]
   sy = stress[2]
   txy = stress[3]
   sz = stress[4]
   dx = (2.0*sx-sy-sz)/3.0
   dy = (2.0*sy-sz-sx)/3.0
   dz = (2.0*sz-sx-sy)/3.0
   sigm = (sx+sy+sz)/3.0
   m1[1,1:2] .= 1.0
   m1[2,1:2] .= 1.0
   m1[4,1:2] .= 1.0
   m1[1,4] = 1.0
   m1[4,4] = 1.0
   m1[2,4] = 1.0
   m1[:,:] = m1[:,:]/9.0/sigm
   m2[1,1] = 2.0/3.0
   m2[2,2] = 2.0/3.0
   m2[4,4] =  2.0/3.0
   m2[2,4] = -1.0/3.0
   m2[4,2] = -1.0/3.0
   m2[1,2] = -1.0/3.0
   m2[2,1] = -1.0/3.0
   m2[1,4] = -1.0/3.0
   m2[4,1] = -1.0/3.0
   m2[3,3] = 2.0
   m3[3,3] = -dz
   m3[1:2,3] .= txy/3.0
   m3[3,1:2] .= txy/3.0
   m3[3,4] = -2.0*txy/3.0
   m3[4,3] = -2.0*txy/3.0
   m3[1,1] = dx/3.0
   m3[2,4] = dx/3.0
   m3[4,2] = dx/3.0
   m3[2,2] = dy/3.0
   m3[1,4] = dy/3.0
   m3[4,1] = dy/3.0
   m3[4,4] = dz/3.0
   m3[1,2] = dz/3.0
   m3[2,1] = dz/3.0
 elseif nst == 6
   sx = stress[1]
   sy = stress[2]
   sz = stress[3]
   txy = stress[4]
   tyz = stress[5]
   tzx = stress[6]
   sigm = (sx+sy+sz)/3.0
   dx = sx-sigm
   dy = sy-sigm
   dz = sz-sigm
   m1[1:3,1:3] = 1.0/(3.0*sigm)
   for i in 1:3
     m2[i,i] = 2.0
     m2[i+3,i+3] = 6.0
   end
   m2[1,2] = -1.0
   m2[1,3] = -1.0
   m2[2,3] = -1.0
   m3[1,1] = dx
   m3[1,2] = dz
   m3[1,3] = dy
   m3[1,4] = txy
   m3[1,5] = -2.0*tyz
   m3[1,6] = tzx
   m3[2,2] = dy
   m3[2,3] = dx
   m3[2,4] = txy
   m3[2,5] = tyz
   m3[2,6] = -2.0*tzx
   m3[3,3] = dz
   m3[3,4] = -2.0*txy
   m3[3,5] = tyz
   m3[3,6] = tzx
   m3[4,4] = -3.0*dz
   m3[4,5] = 3.0*tzx
   m3[4,6] = 3.0*tyz
   m3[5,5] = -3.0*dx
   m3[5,6] = 3.0*txy
   m3[6,6] = -3.0*dy
   for i in 1:6
     for j in i+1:6
       m1[j,i] = m1[i,j]
       m2[j,i] = m2[i,j]
       m3[j,i] = m3[i,j]
     end
   end
   m1[:,:] = m1/3.0
   m2[:,:] = m2/3.0
   m3[:,:] = m3/3.0
 else
   println("nst size not recognised in formm.")
 end
 (m1, m2, m3)
end


text:
# Common univariate distributions.
distset = [
    Normal(1.23, 3.45),
    Exponential(2.12),
    Gamma(4.3),
    Gamma(2.2, 3.1,),
    Erlang(2, 3),
    Beta(2.123),
    LogNormal(4.5),
    DiscreteUniform(1, 10),
    Chisq(3),
    Uniform(1, 200),
    Binomial(101, 0.45),
    Categorical([0.2, 0.35, 0.15, 0.3])
]

for dist in distset
    println(dist)
    μ = mean(dist)
    σ = std(dist)
    # No convenience call.
    E = expectation(dist)
    @test E(x -> x) ≈ μ
    @test E(x -> x^2) - μ^2 ≈ σ^2
    @test E(x -> ((x - μ)/σ)^3) + 1. ≈ skewness(dist) + 1. # To avoid comparisons to 0.0 exactly.
    # Convenience call.
    @test expectation(x -> x, dist) ≈ μ
    @test expectation(x -> x^2, dist) - μ^2 ≈ σ^2
    @test expectation(x -> ((x - μ)/σ)^3, dist) + 1. ≈ skewness(dist) + 1.
    # Stress tests.
    # Many nodes.
    E2 = expectation(dist, n = 100)
    @test E2(x -> x) ≈ μ
    @test E2(x -> x^2) - μ^2 ≈ σ^2
    @test E2(x -> ((x - μ)/σ)^3) + 1. ≈ skewness(dist) + 1. # To avoid comparisons to 0.0 exactly.
end

# Linear operator behavior.
distset = [
    DiscreteUniform(1., 10.),
    Normal(1.45),
    Beta(2.1)
]

for dist in distset
    E = expectation(dist)
    h(x) = 2*x
    z = nodes(E)
    @test E * h.(z) ≈ E(x -> 2*x) # Right-multiplying.
    @test weights(2E) ≈ 2*weights(E) # Left-multiplying.
    @test 3E*z ≈ (3E) * z ≈ 3*(E * z) # Linearity.
end

# Error handling.
distset = [ # Noncompact dists
    LogNormal(Inf, Inf),
    Beta(Inf),
    Normal(0, Inf),
    Uniform(-Inf, Inf),
    Poisson(3)
]

for dist in distset
    @test_throws ArgumentError expectation(dist)
end

distset = [ # Compact dist
    Uniform(1,2),
    Arcsine(1,2),
    Beta(1,1),
    Truncated(LogNormal(1,1), 0., 10.),
    ]

for dist in distset
    E = expectation(dist)
    h(x) = 2*x
    z = nodes(E)
    @test E * h.(z) ≈ E(x -> 2*x) # Right-multiplying.
    @test weights(2E) ≈ 2*weights(E) # Left-multiplying.
    @test 3E*z ≈ (3E) * z ≈ 3*(E * z) # Linearity.

    nodeList = nodes(E);
    E = expectation(dist, nodeList)
    @test E * h.(z) ≈ E(x -> 2*x) # Right-multiplying.
    @test weights(2E) ≈ 2*weights(E) # Left-multiplying.
    @test 3E*z ≈ (3E) * z ≈ 3*(E * z) # Linearity.
end

# Other errors.
E = expectation(DiscreteUniform(1, 10))
@test_throws Exception E(x -> dot(x, ones(7))) # Non-applicable functions.
@test_throws MethodError (x -> 2*x).(nodes(E)) * E # Non-commutativity.

# Trapezoidal methods.
distset = [Beta()]

for dist in distset
    # Setup.
    x = support(dist)
    μ = mean(dist)
    σ = std(dist)
    # Regular grid.
    grid = range(minimum(x), stop = maximum(x), length = 100)
    E = expectation(dist, grid)
    @test E(x -> x) ≈ μ
    @test abs(E(x -> x^2) - μ^2 - σ^2) < 1e-4
    # Irregular grid.
    grid2 = unique([grid' range(minimum(x), stop = maximum(x), length = 77)'])
    E2 = expectation(dist, grid2)
    @test E2(x -> x) isa Number # no accuracy guarantees for the irregular grid
    @test abs(E2(x -> x^2) - μ^2 - σ^2) isa Number
    # Convenience method
    @test expectation(identity, dist, grid2) ≈ E2(x -> x)
end

# Quantile

distset = [
    Uniform(1., 2.)
]

for dist in distset
    E = expectation(dist, QuantileRange)
    h(x) = 2*x
    z = nodes(E)
    @test E * h.(z) ≈ E(x -> 2*x) # Right-multiplying.
    @test weights(2E) ≈ 2*weights(E) # Left-multiplying.
    @test 3E*z ≈ (3E) * z ≈ 3*(E * z) # Linearity.
end


text:
#trapla_path = "/home/scholzs/Documents/Research/Development/Projects/TrajectoryPlanning"
#push!(LOAD_PATH, trapla_path)

using Documenter, Trapla

pages = ["Home" => "index.md", "Gevrey" => "gevrey.md"]

# makedocs(sitename="TrajectoryPlanning.jl")

makedocs(
    sitename="Trapla.jl",
    format=Documenter.HTML(),
    modules=[Trapla],
    authors = "stephan.scholz@rwu.de",
    doctest = true,
    pages = pages
)

deploydocs(
    repo = "github.com/stephans3/Trapla.jl.git",
)


text:
# From Horowitz paper:
# choosing initial conditions uniformly
# from the simplex of species concentrations with total
# sum of concentrations fixed to 25
# (on average 1 mol/vol of each species).
function rand_concentrations(n_species::Integer, total_conc::Real)
    # see 'Uniform sampling from a simplex' from StackExchange
    # Generate a list of sorted numbers from 0 to 1
    lst = cat(dims=1, [0.0], sort(rand(n_species-1)), [1.0])
    # Take differences between consecutive elements
    return total_conc*(lst[2:end]-lst[1:end-1])
end


text:
# where is the package installed?
const _PATH_TO_SRC = dirname(pathof(@__MODULE__))
const _PATH_TO_DATABASE = joinpath(_PATH_TO_SRC, "database")

# IR dictionary keys -
const ir_master_reaction_table_key = "_master_reaction_table_key"
const ir_list_of_molecular_species_key = "_list_of_molecular_species_key"
const ir_list_of_reaction_tags_key = "_list_of_reaction_tags_key"
const ir_master_species_bounds_table_key = "_master_species_bounds_table_key"
const ir_sequence_section_table_key = "_sequence_section_table_key"

# using statement for external packages -
using DataFrames
using CSV
using Dates
using Logging
using WordTokenizers
using DelimitedFiles
using SQLite

# include package codes -
include("./base/Types.jl")

# Minerva lite -
include("./minerva/Types.jl")
include("./minerva/Scanner.jl")
include("./minerva/Parser.jl")
include("./minerva/Database.jl")
include("./minerva/grammar/type_assignment_grammar.jl")
include("./minerva/grammar/species_bound_grammar.jl")

include("./base/Extensions.jl")
include("./base/Factory.jl")
include("./base/Checks.jl")
include("./base/Files.jl")
include("./base/Parser.jl")
include("./base/Callbacks.jl")
include("./base/Sequence.jl")
include("./strategy/General.jl")

# language specific codes -
include("MakeJuliaModel.jl")

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
# the grid data structure contains 37 members
# (6) Surface structures
# sgeo, and the 5 ids
# (17) Volume structures
# vgeo, and the 16 ids
# (2) Volume to Surface index
# vmap⁻, vmap⁺
# (1) Topology
# topology
# (3) Operators
# grid.D (differentiation matrix), grid.Imat (integration matrix?), grid.ω (quadrature weights)
# information (1)
# activedofs (active degrees of freedom)
# Domain Boundary Information (1)
# elemtobndy (element to boundary)
# MPI stuff (6)
# interiorelems, exteriorelems, nabrtovmaprecv, nabrtovmapsend, vmaprecv, vmapsendv

# sgeo, the first index is for the ids
# the second index are the gauss-lobatto entries
# the third index is the face
# the last index is the element
# the ids are: (1, n1id), (2 n2id), (3 n3id) (4, sMid), (5, vMIid)
# grid.sgeo[1,:,:,:] are normals in the n1 direction
# norm(grid.vgeo[:, grid.MIid, :][grid.vmap⁻] - grid.sgeo[5,:,:,:]) is zero
# x[grid.vmap⁻[10]] is connected to x[grid.vmap⁺[10]]

# vgeo, first index is guass-lobatto points
# second index is ids
# third index is elements
# the are are 16 ids
# 1-9 are metric terms
# 10 is the mass matrix
# 11 is the inverse mass matrix
# 12 is MHid (horizontal mass matrix?)
# 13-15 are x1 x2 x3
# 16 is the vertical volume jacobian
# grid.vgeo[:, grid.ξ1x1id, :]

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
using FEM
using Test

@testset "FEM.jl" begin
    # Write your tests here.
end


text:
########################################
###### Vectorized "dot" operators ######
########################################

# For every binary operation like ^, there is a corresponding "dot" operation .^ that is
# automatically defined to perform ^ element-by-element on arrays. For example, [1,2,3] ^ 3
# is not defined, since there is no standard mathematical meaning to "cubing" a (non-square) array,
# but [1,2,3] .^ 3 is defined as computing the elementwise (or "vectorized") result [1^3, 2^3, 3^3].
# Similarly for unary operators like ! or √, there is a corresponding .√ that applies the operator elementwise.

# [1,2,3] .^ 3
println([1,2,3] .^ 3)  # [1, 8, 27]

max_stars_repo_name: mauricioriva/learn-julia

max_issues_repo_name: mauricioriva/learn-julia

max_forks_repo_name: mauricioriva/learn-julia


text:
module OwnTime

export owntime, totaltime, filecontains
export framecounts, frametotal, frames

using Printf
using Profile

const StackFrame = StackTraces.StackFrame

function countmap(iter)
    result = Dict{eltype(iter), Int64}()
    for i in iter
        if haskey(result, i)
            result[i] += 1
        else
            result[i] = 1
        end
    end
    result
end

mutable struct OwnTimeState
    lookups :: Dict{UInt64, Vector{StackFrame}}
end

const state = OwnTimeState(Dict{UInt64, Vector{StackFrame}}())

"""
    clear()

OwnTime has an internal cache for performance. Clear that cache.

See also: [`Profile.clear`](@ref)
"""
function clear()
    state.lookups = Dict{UInt64, Vector{StackFrame}}()
end

"""
    fetch()

Return a 2-tuple of the form `(instruction_pointers, is_profile_buffer_full)`.

This function is primarily for internal use.

See also: [`Profile.fetch`](@ref)
"""
function fetch()
    maxlen = Profile.maxlen_data()
    len = Profile.len_data()
    data = Vector{UInt}(undef, len)
    GC.@preserve data unsafe_copyto!(pointer(data), Profile.get_data_pointer(), len)
    return data, len == maxlen
end

"""
    backtraces(;warn_on_full_buffer=true)

Return an array of backtraces. A backtrace is an array of instruction pointers.

This function is primarily for internal use, try calling `OwnTime.stacktraces()` instead.

This function may give a warning if Julia's profiling buffer is full.
This warning can be disabled with the relative function parameter.

See also: [`backtrace`](@ref)
"""
function backtraces(;warn_on_full_buffer=true)
    profile_pointers, full_buffer = fetch()
    if warn_on_full_buffer && full_buffer
        @warn """The profile data buffer is full; profiling probably terminated
                 before your program finished. To profile for longer runs, call
                 `Profile.init()` with a larger buffer and/or larger delay."""
    end
    bts = Vector{UInt64}[]
    i = 1
    for j in 1:length(profile_pointers)
        # 0 is a sentinel value that indicates the start of a new backtrace.
        # See the source code for `tree!` in Julia's Profile package.
        if profile_pointers[j] == 0
            push!(bts, profile_pointers[i:j-1])
            i = j+1
        end
    end
    filter(!isempty, bts)
end

"""
    lookup(p::UInt64)

Similar to `StackTraces.lookup`, but internally caches lookup results for performance.

Cache can be cleared with `OwnTime.clear()`.

This function is for internal use.

See also: [`StackTraces.lookup`](@ref) [`OwnTime.clear`](@ref)
"""
function lookup(p::UInt64)
    if !haskey(state.lookups, p)
        state.lookups[p] = StackTraces.lookup(p)
    end
    state.lookups[p]
end

"""
    stacktraces([backtraces]; warn_on_full_buffer=true)

Return an array of `StackTrace`s. A `StackTrace` is an array of `StackFrame`s.

This function may take several minutes if you have a large profile buffer.

See also: [`stacktrace`](@ref), [`StackTraces.StackTrace`](@ref), [`StackTraces.StackFrame`](@ref)
"""
function stacktraces(;warn_on_full_buffer=true)
    bts = backtraces(warn_on_full_buffer=warn_on_full_buffer)
    stacktraces(bts)
end

function stacktraces(backtraces)
    map(backtraces) do backtrace
        filter(reduce(vcat, lookup.(backtrace))) do stackframe
            stackframe.from_c == false
        end
    end
end

struct FrameCounts
    counts :: Vector{Pair{StackFrame,Int64}}
    total :: Int64
end

framecounts(fcs::FrameCounts) = fcs.counts

frametotal(fcs::FrameCounts) = fcs.total

frames(fcs::FrameCounts) = map(fcs -> fcs.first, fcs.counts)

Base.getindex(fcs::FrameCounts, i) = framecounts(fcs)[i]
Base.iterate(fcs::FrameCounts) = iterate(framecounts(fcs))
Base.iterate(fcs::FrameCounts, s) = iterate(framecounts(fcs), s)
Base.length(fcs::FrameCounts) = length(framecounts(fcs))

function Base.show(io::IO, fcs::FrameCounts)
    for (i, (stackframe, count)) in enumerate(fcs)
        percent_of_time = round(count / frametotal(fcs) * 100)
        if percent_of_time >= 1
            @printf(io, "%4s %3d%% => %s\n", @sprintf("[%d]", i), percent_of_time, stackframe)
        end
    end
end

"""
    owntime([stacktraces]; stackframe_filter, warn_on_full_buffer=true)

Count the time spent on each `StackFrame` *excluding* its sub-calls.

If supplied, `stackframe_filter` should be a function that accepts a single `StackFrame`
and returns `true` if it should be included in the counts.

More advance filtering my by done by preprocessing the `StackTrace`s from
`OwnTime.stacktraces()` and then passing those `StackTrace`s to this function.

See also: [`OwnTime.stacktraces`](@ref) [`StackTraces.StackFrame`](@ref)
"""
function owntime(;stackframe_filter=stackframe -> true, warn_on_full_buffer=true)
    sts = stacktraces(warn_on_full_buffer=warn_on_full_buffer)
    owntime(sts; stackframe_filter=stackframe_filter)
end

function owntime(stacktraces; stackframe_filter=stackframe -> true)
    filtered_stacktraces = map(stacktraces) do stackframes
        filter(stackframe_filter, stackframes)
    end
    nonempty_stacktraces = filter(!isempty, filtered_stacktraces)
    framecounts = countmap(reduce(vcat, first.(nonempty_stacktraces), init=StackFrame[]))
    FrameCounts(sort(collect(framecounts), by=pair -> pair.second, rev=true), length(stacktraces))
end

"""
    totaltime([stacktraces]; stackframe_filter, warn_on_full_buffer=true)

Count the time spent on each `StackFrame` *including* its sub-calls.

If supplied, `stackframe_filter` should be a function that accepts a single `StackFrame`
and returns `true` if it should be included in the counts.

More advance filtering my by done by preprocessing the `StackTrace`s from
`OwnTime.stacktraces()` and then passing those `StackTrace`s to this function.

See also: [`OwnTime.stacktraces`](@ref) [`StackTraces.StackFrame`](@ref)
"""
function totaltime(;stackframe_filter=stackframe -> true, warn_on_full_buffer=true)
    sts = stacktraces(warn_on_full_buffer=warn_on_full_buffer)
    totaltime(sts; stackframe_filter=stackframe_filter)
end

function totaltime(stacktraces; stackframe_filter=stackframe -> true)
    filtered_stacktraces = map(stacktraces) do stackframes
        filter(stackframe_filter, stackframes)
    end
    framecounts = countmap(reduce(vcat, collect.(unique.(filtered_stacktraces)), init=StackFrame[]))
    FrameCounts(sort(collect(framecounts), by=pair -> pair.second, rev=true), length(stacktraces))
end

"""
    filecontains(needle)

A `StackFrame` filter that returns `true` if the given `StackFrame` has `needle` in its file path.

# Example
```julia-repl
julia> stackframe = stacktrace()[1]
top-level scope at REPL[6]:1

julia> stackframe.file
Symbol("REPL[6]")

julia> filecontains("REPL")(stackframe)
true
```

See also: [`StackTraces.StackFrame`](@ref)
"""
function filecontains(needle)
    function (stackframe)
        haystack = string(stackframe.file)
        occursin(needle, haystack)
    end
end

end # module

max_issues_repo_name: DevJac/OwnTime


text:
"""
    underapproximate(X::LazySet{N}, dirs::AbstractDirections;
                    [apply_convex_hull]::Bool=false) where {N}

Compute the underapproximation of a convex set by sampling support vectors.

### Input

- `X`                 -- set
- `dirs`              -- directions
- `apply_convex_hull` -- (optional, default: `false`) if `true`, post-process
                         the support vectors with a convex hull operation
### Output

The `VPolytope` obtained by taking the convex hull of the support vectors of `X`
along the directions determined by `dirs`.

### Notes

Since the support vectors are not always unique, this algorithm may return
a strict underapproximation even if the set can be exactly approximated using
the given template.
"""
function underapproximate(X::LazySet{N}, dirs::AbstractDirections;
                          apply_convex_hull::Bool=false) where {N}
    vinner = Vector{Vector{N}}(undef, length(dirs))
    underapproximate!(vinner, X, dirs; apply_convex_hull=apply_convex_hull)
end

# in-place version
function underapproximate!(vinner, X::LazySet{N}, dirs::AbstractDirections;
                           apply_convex_hull::Bool=false) where {N}
   @assert dim(X) == dim(dirs) "the dimension of the set, $(dim(X)), doesn't match " *
                               "the dimension of the template directions, $(dim(dirs))"

    @inbounds for (i, di) in enumerate(dirs)
        vinner[i] = σ(di, X)
    end
    if apply_convex_hull
        convex_hull!(vinner)
    end
    return VPolytope(vinner)
end

function underapproximate(X::LazySet, dirs::Type{<:AbstractDirections}; kwargs...)
    return underapproximate(X, dirs(dim(X)), kwargs...)
end


text:
module MaxPooling2DM
    using LoopVectorization, Polyester

    mutable struct MaxPooling2D
        save_layer::Any
        load_layer::Any
        activate::Any
        initialize::Any
        update::Any

        input_shape::Tuple
        output_shape::Tuple
        pool_size::Tuple{Int64, Int64}
        strides::Tuple{Int64, Int64}

        index::Array{Int64}
        output::Array{Float32}
        δ::Array{Float32}

        function MaxPooling2D(;input_shape::Tuple, pool_size::Tuple{Int64, Int64}, strides::Tuple{Int64,Int64}=(pool_size[1],pool_size[2]))
            pool_num_per_col = (input_shape[1]-pool_size[1])÷strides[1]+1
            pool_num_per_row = (input_shape[2]-pool_size[2])÷strides[2]+1
            output_shape = (pool_num_per_col, pool_num_per_row, input_shape[3])
            new(save_MaxPooling2D, load_MaxPooling2D, activate_MaxPooling2D, init_MaxPooling2D, update_MaxPooling2D, input_shape, output_shape, pool_size, strides)
        end
    end

    function init_MaxPooling2D(layer::MaxPooling2D, mini_batch::Int64)
        layer.index = zeros(Int64, 2, layer.output_shape..., mini_batch)
        layer.output = zeros(Float32, layer.output_shape..., mini_batch)
        layer.δ = zeros(Float32, layer.input_shape..., mini_batch)
    end

    function activate_MaxPooling2D(layer::MaxPooling2D, input::Array{Float32})
        x, y = layer.strides
        # Waiting for the update of LoopVectorization
        #=@avxt for i in axes(layer.output, 1), j in axes(layer.output, 2), c in axes(layer.output, 3), b in axes(layer.output, 4)
            s = -Inf32
            for p₁ in 1:layer.pool_size[1], p₂ in 1:layer.pool_size[2]
                s = ifelse(input[(i-1)*x+p₁, (j-1)*y+p₂, c, b]>=layer.output[i, j, c, b], input[(i-1)*x+p₁, (j-1)*y+p₂, c, b], s)
                layer.index[1, i, j, c, b] = ifelse(input[(i-1)*x+p₁, (j-1)*y+p₂, c, b]>=s, (i-1)*x+p₁, layer.index[1, i, j, c, b])
                layer.index[2, i, j, c, b] = ifelse(input[(i-1)*x+p₁, (j-1)*y+p₂, c, b]>=s, (j-1)*y+p₂, layer.index[2, i, j, c, b])
            end
            layer.output[i, j, c, b] = s
        end=#

        @batch for i in axes(layer.output, 1), j in axes(layer.output, 2), c in axes(layer.output, 3), b in axes(layer.output, 4)
            temp = findmax(view(input, (i-1)*x+1:(i-1)*x+layer.pool_size[1], (j-1)*y+1:(j-1)*y+layer.pool_size[2], c, b))
            layer.output[i, j, c, b] = temp[1]
            layer.index[1, i, j, c, b] = (i-1)*x+temp[2][1]
            layer.index[2, i, j, c, b] = (j-1)*y+temp[2][2]
        end
    end

    function update_MaxPooling2D(layer::MaxPooling2D, optimizer::String, input::Array{Float32}, δₗ₊₁::Array{Float32}, α::Float64, parameters::Tuple, direction::Int64=0)
        @avxt layer.δ .= 0.0f0
        @batch for i in axes(δₗ₊₁, 1), j in axes(δₗ₊₁, 2), c in axes(δₗ₊₁, 3), b in axes(δₗ₊₁, 4)
            layer.δ[layer.index[1, i, j, c, b], layer.index[2, i, j, c, b], c, b] += δₗ₊₁[i, j, c, b]
        end
    end

    function save_MaxPooling2D(layer::MaxPooling2D, file::Any, id::String)
        write(file, id, "MaxPooling2D")
        write(file, id*"input_shape", collect(layer.input_shape))
        write(file, id*"pool_size", collect(layer.pool_size))
        write(file, id*"strides", collect(layer.strides))
    end

    function load_MaxPooling2D(layer::MaxPooling2D, file::Any, id::String)
    end

    function get_args(file::Any, id::String)
        input_shape = tuple(read(file, id*"input_shape")...)
        pool_size = tuple(read(file, id*"pool_size")...)
        strides = tuple(read(file, id*"strides")...)
        return (input_shape=input_shape, pool_size=pool_size, strides=strides)
    end
end


text:
### Exports ###

export NConfig, NState, Neurons
export show, size, nupdate!

### Types ###

mutable struct NConfig{S}
  size::S
  a::Float32
  b::Float32
  c::Float32
  d::Float32
  thresh::Float32
  θRate::Float32
  traceRate::Float32
  boostRate::Float32
end
mutable struct NState{RT<:AbstractArray{Float32}, BT<:AbstractArray{Bool}}
  V::RT
  U::RT
  I::RT
  F::BT
  T::RT
  B::RT
  θ::RT
end
mutable struct Neurons{C<:NConfig, S<:NState} <: AbstractNeurons
  conf::C
  state::S
end

NCPUCont{C,S} = CPUContainer{Neurons{C,S}}

### Methods ###

function reinit!(neurons::Neurons)
  fill!(neurons.state.V, neurons.conf.c)
  fill!(neurons.state.U, 0f0)
  fill!(neurons.state.I, 0f0)
  fill!(neurons.state.F, false)
  fill!(neurons.state.T, 0f0)
  fill!(neurons.state.B, 0f0)
  fill!(neurons.state.θ, 0f0)
end
function Base.show(io::IO, neurons::Neurons)
  print(io, "Neurons ($(neurons.conf.size))")
end

Base.size(neurons::Neurons) = neurons.conf.size

@inline function _nupdate!(i, a, b, c, d, thresh, θRate, traceRate, boostRate, V, U, I, F, T, B, θ)
  # Update V and U
  V[i] += 0.5f0 * ((0.04f0 * (V[i] * V[i])) + (5f0 * V[i]) + 140f0 - U[i] + I[i])
  V[i] += 0.5f0 * ((0.04f0 * (V[i] * V[i])) + (5f0 * V[i]) + 140f0 - U[i] + I[i])
  U[i] += a * ((b * V[i]) - U[i])
  #V[i] += -4f0 * (B[i] - 0.8f0)

  # Spike and reset
  F[i] = (V[i] >= thresh + θ[i])
  V[i] += (F[i] * (c - V[i]))
  U[i] += F[i] * d

  # Update θ firing threshold
  θ[i] += F[i] + (θRate * -θ[i] * !F[i])

  # Update traces
  T[i] += F[i] + (traceRate * -T[i] * !F[i])
  #local traceT = fasttanh(T[i])
  #T[i] += traceRate * (((.75f0 - traceT) * F[i]) + ((0.00f0 - traceT) * !F[i]))

  # Update boosts
  # B[i] += boostRate * (T[i] - B[i])

  # Clear inputs
  I[i] = 0f0
end

function nupdate!(neurons::Neurons)
  conf = neurons.conf
  state = neurons.state
  @inbounds @simd for i in eachindex(state.V)
    _nupdate!(i, conf.a, conf.b, conf.c, conf.d, conf.thresh, conf.θRate, conf.traceRate, conf.boostRate, state.V, state.U, state.I, state.F, state.T, state.B, state.θ)
  end
end
nupdate!(ncont::NCPUCont{C,S} where {C,S}) = nupdate!(transient(ncont))
cycle_neurons!(neurons) = nupdate!(neurons)

### Includes ###

include("generic/core.jl")
include("conv/core.jl")


text:
using AtomsBase
using Unitful
using Test

@testset "atomic systems" begin
    @testset "Atom construction" begin
        at = Atom(:Si, zeros(3) * u"m", extradata=42)

        @test n_dimensions(at) == 3
        @test position(at) == zeros(3) * u"m"
        @test velocity(at) == zeros(3) * u"bohr/s"
        @test at.atomic_symbol == :Si
        @test at.atomic_number == 14
        @test hasproperty(at, :atomic_mass)
        @test hasproperty(at, :atomic_symbol)
        @test hasproperty(at, :atomic_number)
        @test hasproperty(at, :extradata)
        @test at.extradata == 42

        @test propertynames(at) == (:position, :velocity, :atomic_symbol,
                                    :atomic_number, :atomic_mass, :extradata)

        # Test update constructor
        newatom = Atom(at; extradata=43, atomic_number=15)
        @test propertynames(at) == propertynames(newatom)
        @test newatom.extradata == 43
        @test newatom.atomic_number == 15

        newatom = Atom(; atom=at, extradata=43, atomic_number=15)
        @test propertynames(at) == propertynames(newatom)
        @test newatom.extradata == 43
        @test newatom.atomic_number == 15
    end

    @testset "flexible atomic systems" begin
        box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"
        bcs = [Periodic(), DirichletZero(), DirichletZero()]
        atoms = [:Si => [0.0, 1.0, 1.5]u"Å",
                 :C  => [0.0, 0.8, 1.7]u"Å",
                 Atom(:H, zeros(3) * u"Å", ones(3) * u"bohr/s")]
        system = atomic_system(atoms, box, bcs)
        @test length(system) == 3
        @test atomic_symbol(system) == [:Si, :C, :H]
        @test boundary_conditions(system) == [Periodic(), DirichletZero(), DirichletZero()]
        @test position(system) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7], [0.0, 0.0, 0.0]]u"Å"
        @test velocity(system) == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [1.0, 1.0, 1.0]]u"bohr/s"

        # Test update constructor
        newatoms  = [system[1], system[2]]
        newsystem = FlexibleSystem(system; atoms=newatoms,
                                   boundary_conditions=[Periodic(), Periodic(), Periodic()])
        @test length(newsystem) == 2
        @test atomic_symbol(newsystem) == [:Si, :C]
        @test boundary_conditions(newsystem) == [Periodic(), Periodic(), Periodic()]
        @test position(newsystem) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7]]u"Å"
    end

    @testset "isolated_system" begin
        system = isolated_system([:Si => [0.0, 1.0, 1.5]u"Å",
                                  :C  => [0.0, 0.8, 1.7]u"Å",
                                  Atom(:H, zeros(3) * u"Å")])
        @test length(system) == 3
        @test atomic_symbol(system) == [:Si, :C, :H]
        @test boundary_conditions(system) == [DirichletZero(), DirichletZero(), DirichletZero()]
        @test position(system) == [[0.0, 1.0, 1.5], [0.0, 0.8, 1.7], [0.0, 0.0, 0.0]]u"Å"
        @test velocity(system) == [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 0.0]]u"bohr/s"
        @test bounding_box(system) == infinite_box(3)
    end

    @testset "periodic_system" begin
        box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"
        atoms = [:Si => [0.0, -0.125, 0.0],
                 :C  => [0.125, 0.0, 0.0],
                 Atom(:H, zeros(3) * u"Å")]
        system = periodic_system(atoms, box; fractional=true)

        @test length(system) == 3
        @test atomic_symbol(system) == [:Si, :C, :H]
        @test boundary_conditions(system) == [Periodic(), Periodic(), Periodic()]
        @test position(system) == [[0.0, -0.625, 0.0], [1.25, 0.0, 0.0], [0.0, 0.0, 0.0]]u"Å"
    end

    @testset "no stackoverflow" begin
        box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"
        atoms = Any[:Si => [0.0, -0.125, 0.0],
                    :C  => [0.125, 0.0, 0.0],
                    Atom(:H, zeros(3) * u"Å")]
        system = periodic_system(atoms, box; fractional=true)
        @test length(system) == 3

        system = isolated_system(Any[:Si => [0.0, 1.0, 1.5]u"Å",
                                     :C  => [0.0, 0.8, 1.7]u"Å",
                                     Atom(:H, zeros(3) * u"Å")])
        @test length(system) == 3
    end
end


text:
using Test, Statistics
using Knet.Ops20: batchnorm, bnmoments, bnparams
using Knet.KnetArrays: KnetArray
using AutoGrad: gradcheck
using CUDA: CUDA, functional

@testset "batchnorm" begin

    # Random.seed!(42)

    # utils
    std2(x) = let x_mu = x .- mean(x)
        mean(x_mu .* x_mu)
    end

    # gradcheck functions
    bn3(a) = batchnorm(a[1], nothing, a[2]; training=true)
    bn1(a) = batchnorm(a; training=true)
    bn3ts(a) = batchnorm(a[1], bnmoments(), a[2]; training=false)
    bn1ts(a) = batchnorm(a, bnmoments(); training=false)

    types = (Float64,) #TODO: [Float32, Float64]
    sizes = Dict([2=>(5,10), 4=>(3,4,5,3), 5=>(4,3,4,5,2)])
    dims = [2, 4, 5]
    gpu_av = CUDA.functional()

    for d in dims
        for et in types
            sz = sizes[d]
            tn(str) = string(str, d, "d")
            # CPU buffers
            ax = 5rand(et, sz...);
            C = sz[end-1]
            aw = bnparams(et, C)
            # GPU buffers
            kax = gpu_av ? KnetArray{et}(5rand(sz...)) : nothing
            kaw = gpu_av ? KnetArray{et}(bnparams(C)) : nothing

            # Array types
            at = typeof(ax)
            kat = typeof(kax)

            @testset "{$et, $d}" begin
                @testset "cpu-stat" begin
                    # The unit gaussian output
                    @test abs(mean(batchnorm(ax; training=true))) < 1e-3
                    @test isapprox(std2(batchnorm(ax; training=true)), 1.0; rtol=1e-2)
                end

                @testset "cpu-grads" begin
                    @test gradcheck(bn1, ax)
                    @test gradcheck(bn3, (ax, aw))
                end

                if gpu_av
                    @testset "gpu-stats" begin
                        # isapprox 0 didn't work
                        @test abs(mean(batchnorm(kax; training=true))) < 1e-3
                        @test isapprox(std2(batchnorm(kax; training=true)), 1.0; rtol=1e-2)
                    end

                    @testset "gpu-grads" begin
                        @test gradcheck(bn1, kax)
                        @test gradcheck(bn3, (kax, kaw))
                    end

                    @testset "dev-consistency" begin
                        mc = bnmoments() #cpu
                        mg = bnmoments() #gpu
                        y1 = batchnorm(ax, mc)
                        y2 = batchnorm(kat(ax), mg) #use the same array
                        @test isapprox(mc.mean, at(mg.mean))
                        @test isapprox(mc.var, at(mg.var))
                        @test isapprox(y1, at(y2))
                    end
                end

                # Moments should be updates
                @testset "test-moments" begin
                    mc = bnmoments()
                    batchnorm(ax, mc)
                    m, v = copy(mc.mean), copy(mc.var)
                    batchnorm(.1ax.+3, mc)
                    @test isapprox(m, mc.mean; rtol=1e-20)
                    @test isapprox(v, mc.var; rtol=1e-20)
                    if gpu_av
                        mg = bnmoments()
                        batchnorm(kax, mg)
                        m, v = copy(mg.mean), copy(mg.var)
                        batchnorm(.1kax.+3, mg)
                        @test isapprox(m, mg.mean; rtol=1e-20)
                        @test isapprox(v, mg.var; rtol=1e-20)
                    end
                end

                @testset "training-moments" begin
                    mc = bnmoments()
                    batchnorm(ax, mc; training=true)
                    m, v = copy(mc.mean), copy(mc.var)
                    batchnorm(.1ax.+3, mc; training=true)
                    @test ~isapprox(m, mc.mean; rtol=1e-20)
                    @test ~isapprox(v, mc.var; rtol=1e-20)
                    if gpu_av
                        mg = bnmoments()
                        batchnorm(kax, mg; training=true)
                        m, v = copy(mg.mean), copy(mg.var)
                        batchnorm(.1kax.+3, mg; training=true)
                        @test ~isapprox(m, mg.mean; rtol=1e-20)
                        @test ~isapprox(v, mg.var; rtol=1e-20)
                    end
                end
                # TODO: Remove this if after 2d mode test backward supported
                if d > 2
                    @testset "cpu-grads-testing" begin
                        m1 = bnmoments()
                        @test gradcheck(bn1ts, ax)
                        @test gradcheck(bn3ts, (ax, aw))
                    end

                    if gpu_av
                        @testset "gpu-grads-testing" begin
                            @test gradcheck(bn1ts, kax)
                            @test gradcheck(bn3ts, (kax, kaw))
                        end
                    end
                end
            end #end of {dim, type} testset
        end
    end
end #end of testset

# suppress the return
nothing


text:
const libceed_envvar = "JULIA_LIBCEED_LIB"

if isfile("deps.jl")
   rm("deps.jl")
end

if haskey(ENV, libceed_envvar)
   ceedpath = ENV[libceed_envvar]
   if !isfile(ceedpath)
      error("No library file found at $libceed_envvar, given by $ceedpath")
   end
   @info "Using libCEED library specified by $libceed_envvar at $ceedpath"
   open("deps.jl", write=true) do f
      println(f, "const libceed = \"$(escape_string(ceedpath))\"")
   end
else
   @info "Using prebuilt libCEED binaries provided by libCEED_jll"
   open("deps.jl", write=true) do f
      println(f, "using libCEED_jll")
   end
end

max_stars_repo_name: DiffeoInvariant/libCEED

max_issues_repo_name: DiffeoInvariant/libCEED

max_forks_repo_name: DiffeoInvariant/libCEED


text:
# # Example: tagged paired-end reads
#
# This is an example generated from this source
# file: [`tg-example.jl`](@__REPO_ROOT_URL__/examples/sequencing/tg-example.jl)
# You are seeing the
#md # online documentation version. The corresponding notebook can be found
#md # here: [`tg-example.ipynb`](@__NBVIEWER_ROOT_URL__/man/sequencing/examples/tg-example.ipynb)
#nb # jupyter notebook version. The corresponding online documentation page can
#nb # be found here: [`tg-example.html`](https://bioinfologics.github.io/Pseudoseq.jl/latest/examples/sequencing/tg-example),
# and the script can be found here: [`tg-example.jl`](./tg-example.jl)

# Let's see how we might simulate something like an 10x sequencing experiment.
#
# For this simulation script we will:
#
# 1. Create a pool of 5000 copies of a reference genome.
# 2. Fragment the DNA molecules in the pool, to an average length of 40,000bp.
# 3. Tag the long molecules in the pool randomly with a set of 1,000,000 tags.
# 4. Fragment the molecules in the pool to an average length of 700bp.
# 5. Subsample the molecules in the pool to achieve approximatly 50x coverage.
# 6. Create a set of 250bp paired-end reads.
# 7. Apply errors to the paired-end reads at a rate of 0.001 (.1%).
# 8. Generate an output FASTQ file.

using Pseudoseq.Sequencing

# Let's start with a pool of 5000 copies of a genome contained in a FASTA file:

m = Molecules("ecoli-ref.fasta")

# We want to first have 5000 full copies of the genome, so we will make an amplifier.

amp = amplify(5000)

# Next we want a fragmenter to make a pool of shorter DNA molecules, with an average
# length of 40,000bp.

firstfrag = fragment(40000bp)

# Ok, now we want to assign random tags to these large fragments.
# Once you tag a molecule any other molecules that are derived from that tagged
# molecule will inherit the same tag.
# Our tagger will have a pool of 1,000,000 distinct tags.
# Which fragment gets a certain tag is random.
# The size of the tag pool, and the number of fragments in your universe will
# determine how likely it is that any two fragments get the same tag.

tagger = tag(1000000)

# We'll want to fragment the molecules again to even shorter lengths.

secondfrag = fragment(700bp)

# We will want a filter to make sure no extreme sized molecules make their way through.

size_filter = select(x -> 900 >= length(x) >= 450)

# Next we create a subsampler which will randomly sample molecules, to give us
# a desired expected coverage. We will want our reads to be 250bp in length.

ssmpl = subsample(50X, 2 * 250bp)

# We then want a read-maker that will give us paired end reads. We want our
# reads to be 250bp in length.

readmaker = makereads(2 * 250bp)

# Once we have reads, we will mark positions in the reads that are incorrectly
# detected by the sequencer: errors.
# We will construct a `FixedProbSubstitutions` function with a per base error
# probability of 0.001 and pass it to the `make_substitutions` method.
# This will make errors fall totally randomly over each read.

errmaker = make_substitutions(FixedProbSubstitutions(0.001))

# Now we can push our molecules through a pipeline of these processors, and out
# to a FASTQ file:

m |> amp |> firstfrag |> tagger |> secondfrag |> size_filter |> ssmpl |> readmaker |> errmaker |> generate("tagged_reads.fastq")

# You can also compose the processors together into one whole function.
# Typing \circ in the julia repl and then hitting tab gives you the circular
# composition symbol. Note how pipelining above progresses from left to right,
# but composition is right to left in order.

my_protocol = errmaker ∘ readmaker ∘ ssmpl ∘ size_filter ∘ secondfrag ∘ tagger ∘ firstfrag ∘ amp

m |> my_protocol |> generate("tagged_reads.fastq")


text:
"""
Implementation of Grid Partial Rejection Sampling of [MoKr20](@cite).
"""

"""
    AbstractCellGridPRS

Abstract type describing a cell used Grid Partial Rejection Sampling of [MoKr20](@cite).
"""
abstract type AbstractCellGridPRS end

"""
    GraphCellGridPRS{T} <: AbstractCellGridPRS

Mutable struct describing a graph node with fields

- `window` of type [`PRS.GraphNode`](@ref)`{T}`
- `value` of type `T`
"""
mutable struct GraphCellGridPRS{T} <: AbstractCellGridPRS
    window::GraphNode{T}
    value::T
end

"""
    SpatialCellGridPRS{T<:Vector{Float64}} <: AbstractCellGridPRS

Mutable struct describing a spatial cell with fields

- `window` of type `Union{`[`PRS.RectangleWindow`](@ref)`,`[`PRS.RectangleWindow`](@ref)`}`
- `value` of type `Vector{T}`
"""
mutable struct SpatialCellGridPRS{T<:Vector{Float64}} <: AbstractCellGridPRS
    window::Union{RectangleWindow,SquareWindow}
    value::Vector{T}
end

"""
    window(cell::AbstractCellGridPRS) = cell.window
"""
window(cell::AbstractCellGridPRS) = cell.window

"""
    dimension(cell::AbstractCellGridPRS) = dimension(window(cell))
"""
dimension(cell::AbstractCellGridPRS) = dimension(window(cell))

"""
    Base.iterate(cell::AbstractCellGridPRS, state=1) = iterate(cell.value, state)
"""
Base.iterate(cell::AbstractCellGridPRS, state=1) = iterate(cell.value, state)

"""
    Base.isempty(cell::SpatialCellGridPRS) = isempty(cell.value)
"""
Base.isempty(cell::SpatialCellGridPRS) = isempty(cell.value)

"""
    generate_sample!(
        rng::Random.AbstractRNG,
        cell::AbstractCellGridPRS,
        pp::AbstractPointProcess
    )

Generate an exact sample of `pp` in `cell.window` and save it in `cell.value`
"""
function generate_sample!(
    rng::Random.AbstractRNG,
    cell::AbstractCellGridPRS,
    pp::AbstractPointProcess
)
    cell.value = generate_sample(rng, pp, cell.window)
end

"""
    generate_sample!(
        rng::Random.AbstractRNG,
        cells::Vector{T},
        indices,
        pp::AbstractPointProcess
    ) where {T<:AbstractCellGridPRS}

Apply [`PRS.generate_sample!`](@ref) to each cell of `cells` indexed by `indices`.
"""
function generate_sample!(
    rng::Random.AbstractRNG,
    cells::Vector{T},
    indices,
    pp::AbstractPointProcess
) where {T<:AbstractCellGridPRS}
    for i in indices
        generate_sample!(rng, cells[i], pp)
    end
end

"""
    generate_sample_grid_prs(
        [rng::Random.AbstractRNG,]
        pp::AbstractPointProcess{T}
    )::Vector{T} where {T}

Generate an exact sample from `pp` using grid Partial Rejection Sampling (grid PRS) of [MoKr20](@cite).
"""
function generate_sample_grid_prs(
    rng::Random.AbstractRNG,
    pp::AbstractPointProcess{T}
)::Vector{T} where {T}
    g = weighted_interaction_graph(rng, pp)
    cells = initialize_cells(pp)
    resample_indices = Set(1:length(cells))

    while !isempty(resample_indices)
        generate_sample!(rng, cells, resample_indices, pp)
        resample_indices = find_cells_to_resample_indices!(rng, g, cells, pp)
    end
    return vcat(getfield.(cells, :value)...)
end

function generate_sample_grid_prs(pp::AbstractPointProcess, args...)
    return generate_sample_grid_prs(Random.default_rng(), pp, args...)
end

@doc raw"""
    gibbs_interaction(
        pp::AbstractPointProcess{T},
        cell1::AbstractCellGridPRS{T},
        cell2::AbstractCellGridPRS{T}
    )::Real where {T}

Compute the pairwise Gibbs interaction `pp` between ``C_1``=`cell1` and ``C_2``=`cell2`.

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function gibbs_interaction end

@doc raw"""
    find_bad_cells_indices!(
        [rng::Random.AbstractRNG,]
        g::SWG.SimpleWeightedGraph{T,U},
        cells::Vector{V},
        pp::AbstractPointProcess
    )::Set{T} where {T,U,V<:AbstractCellGridPRS}

Identify bad events and return the corresponding cells' index.
An event ``\{i,j\}`` is said to be \"bad\"

```math
    \left\{U_{ij} > \exp \left[ -\sum_{x \in C_i} \sum_{y \in C_j} V(x,y) \right] \right\}
```

where ``U_{ij}`` is the weight of the edge ``\{i,j\}`` in the interaction graph ``g`` created by [`PRS.weighted_interaction_graph`](@ref) and ``V`` the Gibbs potential discribing the pairwise Gibbs interaction of `pp`

**Note** when a bad event occurs, the corresponding ``U_{ij}`` is resampled hence the "!"

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function find_bad_cells_indices!(
    rng::Random.AbstractRNG,
    g::SWG.SimpleWeightedGraph{T,U},
    cells::Vector{V},
    pp::AbstractPointProcess
)::Set{T} where {T,U,V<:AbstractCellGridPRS}
    bad = Set{T}()
    for e in LG.edges(g)
        i, j = Tuple(e)
        if e.weight > gibbs_interaction(pp, cells[i], cells[j])
            union!(bad, i, j)
            # resample U_ij associated to the bad event
            @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng, U)
        end
    end
    return bad
end

@doc raw"""
    find_cells_to_resample_indices!(
        rng::Random.AbstractRNG,
        g::SWG.SimpleWeightedGraph{T,U},
        cells::Vector{V},
        pp::AbstractPointProcess
    )::Set{T} where {T,U,V<:AbstractCellGridPRS}

Identify the set of events to be resampled as constructed by Algorithm 5 in [GuJeLi19](@cite) as part of the Partial Rejection Sampling (PRS) method.
Return the indices of the variables (here cells) involved in the corresponding events.

This function is used as a subroutine of the grid PRS methodology of [MoKr20](@cite), see [`PRS.generate_sample_grid_prs`](@ref).

**Note** if the event associated to the edge ``\{i,j\}`` of `g` is selected to be resampled, the uniform random variable encoded as the weight of the correspond edge is resampled (hence the "!")

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function find_cells_to_resample_indices!(
    rng::Random.AbstractRNG,
    g::SWG.SimpleWeightedGraph{T,U},
    cells::Vector{V},
    pp::AbstractPointProcess
)::Set{T} where {T,U,V<:AbstractCellGridPRS}
    R = find_bad_cells_indices!(rng, g, cells, pp)
    ∂R, ∂R_tmp = copy(R), empty(R)
    while !isempty(∂R)
        for i in ∂R
            isempty(cells[i]) && continue
            for j in LG.neighbors(g, i)
                if j in R
                    if is_inner_interaction_possible(pp, cells[i], cells[j])
                        @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng, U)
                    end
                elseif is_outer_interaction_possible(pp, cells[i], cells[j])
                    union!(R, j)
                    union!(∂R_tmp, j)
                    @inbounds g.weights[i, j] = g.weights[j, i] = rand(rng, U)
                end
            end
        end
        ∂R, ∂R_tmp = ∂R_tmp, empty(∂R)
    end
    return R
end

## Spatial point processes

"""
    weighted_interaction_graph(
        rng::Random.AbstractRNG,
        pp::AbstractSpatialPointProcess;
    )::SWG.SimpleWeightedGraph

Construct the weighted interaction graph ([King graph](https://en.wikipedia.org/wiki/King%27s_graph)) used in [`PRS.generate_sample_grid_prs`](@ref), to generate exact samples from [`PRS.AbstractSpatialPointProcess`](@ref)

The `pp.window` is divided into cells of length the interaction range `pp.r`.
Each cell represents a vertex of the interaction (king) graph and each edge carries a uniform random varialble.

**See also**

- Figure 4 [MoKr20](@cite)
"""
function weighted_interaction_graph(
    rng::Random.AbstractRNG,
    spp::AbstractSpatialPointProcess
)::SWG.SimpleWeightedGraph
    window_ = window(spp)
    if allequal(window_.w)
        nb_cells_x = Int(cld(window_.w[1], spp.r))
        return uniform_weighted_graph(rng, king_graph(nb_cells_x))
    else
        nb_cells_x, nb_cells_y = Int(cld.(window_.w, spp.r))
        return uniform_weighted_graph(rng, king_graph(nb_cells_x, nb_cells_y))
    end
end

@doc raw"""
    initialize_cells(
        spp::AbstractSpatialPointProcess{T},
    )::Vector{SpatialCellGridPRS{T}} where {T}

The `spp.window` ([`PRS.RectangleWindow`](@ref) or [`PRS.SquareWindow`](@ref)) is divided into [`PRS.SpatialCellGridPRS`](@ref) of length the interaction range `pp.r` following the construction of [MoKr20](@cite) in their grid Partial Rejection Sampling (grid PRS) methodology.

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function initialize_cells(
    spp::AbstractSpatialPointProcess{T},
)::Vector{SpatialCellGridPRS{T}} where {T}
    dimension(spp) != 2 && throw(DomainError(spp, "must be 2D point process for now"))
    window_ = window(spp)
    nb_cells_x = nb_cells_y = 1
    if allequal(window_.w)  # SquareWindow
        nb_cells_x = nb_cells_y = Int(cld(window_.w[1], spp.r))
    else
        nb_cells_x, nb_cells_y = Int(cld.(Int, window_.w, spp.r))
    end

    cells = Vector{SpatialCellGridPRS{T}}(undef, nb_cells_x * nb_cells_y)
    for i in eachindex(cells)
        c_ij = divrem(i-1, nb_cells_x)  # coordinate (i, j) of the cell in the grid
        c_xy = @. window_.c + spp.r * c_ij
        win_i = rectangle_square_window(c_xy, @. min(spp.r, window_.w - c_xy))
        cells[i] = SpatialCellGridPRS(win_i, T[])
    end
    return cells
end

@doc raw"""
    is_inner_interaction_possible(
        spp::AbstractSpatialPointProcess,
        cell_i::SpatialCellGridPRS,
        cell_j::SpatialCellGridPRS
    ) = false

Assume `cell_i` and `cell_j` are neighboring cells in the weighted interaction graph constructed by [`PRS.weighted_interaction_graph`](@ref) from `spp` and already identified in the set of variables to be resampled in [`PRS.generate_sample_grid_prs`](@ref)
Since the configuration of points in the corresponding cells and the uniform random variable associated to the event `\{i,j\}` are considered fixed, there is no degree of freedom to make the interaction between `cell_i` and `cell_j` possible.

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function is_inner_interaction_possible(
        spp::AbstractSpatialPointProcess,
        cell_i::SpatialCellGridPRS,
        cell_j::SpatialCellGridPRS
)::Bool
    return false
end

@doc raw"""
    is_outer_interaction_possible(
        spp::AbstractSpatialPointProcess,
        cell_i::SpatialCellGridPRS,
        cell_j::SpatialCellGridPRS
    )::Bool

Assume `cell_i` and `cell_j` are neighboring cells in the weighted interaction graph constructed by [`PRS.weighted_interaction_graph`](@ref) from `spp`.
Given the configuration of points in `cell_i`, check whether a realization of `spp` in `cell_j` can induce a bad event.

This is a subroutine of [`PRS.generate_sample_grid_prs`](@ref).
"""
function is_outer_interaction_possible(
        spp::AbstractSpatialPointProcess,
        cell_i::SpatialCellGridPRS,
        cell_j::SpatialCellGridPRS
)::Bool
    isempty(cell_i) && return false

    win_i, win_j = window(cell_i), window(cell_j)
    # If i-j is a vertical or horizontal neighborhood
    any(win_i.c .== win_j.c) && return true

    # If i-j is a diagonal neighborhood
    c_ij = any(win_i.c .< win_j.c) ? win_j.c : win_i.c
    return any(Distances.euclidean(c_ij, x) < spp.r for x in cell_i)
end


text:
using Luxor
using Test

function checkpoints(points, hull)
    flag  = true
    for pt in points
        if ! isinside(pt, offsetpoly(hull, 0.1), allowonedge=true)
            flag = false
        end
    end
    return flag
end

function testpolyhull1(fname)
    Drawing(800, 1200, fname)
    origin()
    setopacity(0.5)
    tiles = Tiler(600, 600, 4, 4)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            if n == 1
                pts = box(O, 100, 100, vertices=true)
            elseif n == 2
                pts = ngon(O, 60, 3, vertices=true)
            elseif n == 3
                pts = randompointarray(BoundingBox(box(O, 100, 100, vertices=true)), 25)
            elseif n == 4
                pts = star(O, 60, 5, 0.2, vertices=true)
            else
                pts = ngon(O, 60, n, vertices=true)
            end
            hull = polyhull(pts)
            @test checkpoints(pts, hull)
            sethue("black")
            setline(10)
            poly(pts, :fill, close=true)
            sethue("red")
            setline(2)
            poly(hull, :stroke, close=true)
        end
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

function testpolyhull2(fname)
    Drawing(800, 1200, fname)
    origin()
    setopacity(0.5)
    points = randompointarray(BoundingBox(), 30)
    while length(points) > 2
        hull = polyhull(points)
        randomhue()
        prettypoly(hull, close=true, :fill)
        # discard this hull and repeat
        points = setdiff(points, hull)
    end
    @test finish() == true
    println("...finished test: output in $(fname)")
end

testpolyhull1("polyhull-test-1.png")
testpolyhull2("polyhull-test-2.png")


text:
# First load the package
using PreprocessData

# List all registered datasets
listdatasets()

# Looking for a classification dataset, list them.
listdatasets(Classification)

# Show information about Abalone dataset
info(PreprocessData.Abalone())

# Export Abalone dataset for easy usage
using PreprocessData: Abalone

# load dataset to see how it looks like
load(Abalone())

# show column names, if there are any
getheader(Abalone())

# split dataset
train, test = split_traintest(Abalone(), header=true)

# normalization
# compute mean and standard deviation
mean,std=meanstd(train)
normalize!(Std,train)
normalize!(Std,test,mean,std)

# show all possible classes of Abalone dataset
classes(Abalone())

# binarize labels - numbers 1 to 9 are positive labels, others are negative
bin_train=binarize(train[:,end],collect(1:9))
bin_test=binarize(test[:,end],collect(1:9))

# delete files containing Abalone dataset
remove(Abalone())

# check whether Abalone dataset is downloaded
isdownloaded(Abalone())

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
using NoteSequences.PerformanceRepr

const pe = PerformanceEvent

@testset "PerformanceEvent encoding/decoding" begin
    perfencoder = PerformanceOneHotEncoding(num_velocitybins=16)
    pairs = [
        (pe(NOTE_ON, 60), 61),
        (pe(NOTE_ON, 0), 1),
        (pe(NOTE_ON, 22), 23),
        (pe(NOTE_ON, 127), 128),
        (pe(NOTE_OFF, 72), 201),
        (pe(NOTE_OFF, 0), 129),
        (pe(NOTE_OFF, 22), 151),
        (pe(NOTE_OFF, 127), 256),
        (pe(TIME_SHIFT, 10), 266),
        (pe(TIME_SHIFT, 1), 257),
        (pe(TIME_SHIFT, 72), 328),
        (pe(TIME_SHIFT, 100), 356),
        (pe(VELOCITY, 5), 361),
        (pe(VELOCITY, 1), 357),
        (pe(VELOCITY, 16), 372)
    ]

    for (event, index) in pairs
        @test index == encode_event(event, perfencoder)
        @test event == decode_event(index, perfencoder)
    end
end

@testset "Performance steps and setlength!" begin
    @testset "Add length" begin
        perf = Performance(100)
        setlength!(perf, 42)
        @test perf.numsteps == 42
        @test perf.events == [pe(TIME_SHIFT, 42)]

        setlength!(perf, 142)
        @test perf.numsteps == 142
        @test perf.events == [pe(TIME_SHIFT, 100), pe(TIME_SHIFT, 42)]

        setlength!(perf, 200)
        @test perf.numsteps == 200
        @test perf.events == [pe(TIME_SHIFT, 100), pe(TIME_SHIFT, 100)]
    end

    @testset "Remove length" begin
        perf = Performance(100)
        events = [
            pe(NOTE_ON, 60),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 60),
            pe(NOTE_ON, 64),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 64),
            pe(NOTE_ON, 67),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 67)
        ]
        append!(perf, events)
        setlength!(perf, 200)
        result_events = [
            pe(NOTE_ON, 60),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 60),
            pe(NOTE_ON, 64),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 64),
            pe(NOTE_ON, 67)
        ]

        @test perf.events == result_events

        setlength!(perf, 50)
        result_events = [
            pe(NOTE_ON, 60),
            pe(TIME_SHIFT, 50)
        ]

        @test perf.events == result_events
    end

    @testset "numsteps" begin
        perf = Performance(100)
        events = [
            pe(VELOCITY, 32),
            pe(NOTE_ON, 60),
            pe(TIME_SHIFT, 100),
            pe(NOTE_OFF, 60)
        ]
        append!(perf, events)

        @test perf.numsteps == 100
    end
end