text:
"""
    struct ParamIntervalInput
        theta_init::Vector{Float64} # initial parameters vector
        theta_num::Int # number of the parameter for identifiability analysis
        scan_func::Function # scan function
        loss_func::Function # loss function
        loss_crit::Float64 # loss function maximum value, "identifiability level"
        scale::Vector{Symbol}
        theta_bounds::Vector{Tuple{Float64, Float64}} # search bounds for id parameter
        scan_bounds::Tuple{Float64,Float64}
        scan_tol::Float64 # fitting tolerance for local optimizer (default - `1e-3`)
        loss_tol::Float64 # constraints tolerance
        local_alg::Symbol # local fitting algorithm (default - `:LN_NELDERMEAD`)
        fitter_options::Any
    end

Structure storing input data for parameter identification
"""
struct ParamIntervalInput
    theta_init::Vector{Float64} # initial parameters vector
    theta_num::Int # number of the parameter for analysis
    scan_func::Function
    loss_func::Function # loss function
    loss_crit::Float64 # loss function maximum value, "identifiability level"
    scale::Vector{Symbol}
    theta_bounds::Vector{Tuple{Float64, Float64}} # search bounds for id parameter
    scan_bounds::Tuple{Float64,Float64}
    scan_tol::Float64 # fitting tolerance for local optimizer (default - 1e-3)
    loss_tol::Float64 # constraints tolerance
    local_alg::Symbol # local fitting algorithm (default - :LN_NELDERMEAD)
    fitter_options::Any
end

"""
    struct ParamInterval
        input::ParamIntervalInput
        loss_init::Float64
        method::Symbol
        result::Tuple{EndPoint, EndPoint}
    end

Structure storing result of parameter identification
"""
struct ParamInterval
    input::ParamIntervalInput
    loss_init::Float64
    method::Symbol
    result::Tuple{EndPoint, EndPoint}
end

"""
    get_interval(
        theta_init::Vector{Float64},
        theta_num::Int,
        loss_func::Function,
        method::Symbol;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bounds::Tuple{Float64,Float64} = unscaling.(
            (-9.0, 9.0),
            scale[theta_num]
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 1e-3,
        local_alg::Symbol = :LN_NELDERMEAD,
        kwargs...
        )
Computes confidence interval for single component `theta_num` of parameter vector.

## Return
`ParamInterval` structure storing all input data and estimated confidence interval.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `theta_num`: index of vector component for identification: `theta_init(theta_num)`.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` for profile likelihood-based (PL) identification. Usually we use log-likelihood for PL analysis: ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently the following methods are implemented: `:CICO_ONE_PASS`, `:LIN_EXTRAPOL`, `:QUADR_EXTRAPOL`.

## Keyword arguments
- `loss_crit`: critical level of loss function. Confidence interval's endpoint value is the intersection point of profile likelihood and `loss_crit` level.
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct, :log, :logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `theta_bounds`: vector of tuple `(lower_bound, upper_bound)` for each parameter. Bounds define the ranges for possible parameter values. Default bounds are `(-Inf,Inf)`.
- `scan_bounds`: scan bounds tuple for `theta_num` parameter. Should be within the `theta_bounds` for `theta_num` parameter. Default is `(1e-9,1e9)`.
- `scan_tol`: Absolute tolerance for `theta_num` parameter used as termination criterion.
- `loss_tol`: Absolute tolerance controlling `loss_func` closenes to `loss_crit` (termination criterion). Currently doesn't work for `:CICO_ONE_PASS` method because of limitation in `LN_AUGLAG` interface.
- `local_alg`: algorithm of optimization. Derivative-free and gradient-based algorithms form NLopt package.
- `kwargs...`: the additional `method` specific keyword arguments.

"""
function get_interval(
    theta_init::Vector{Float64},
    theta_num::Int,
    loss_func::Function,
    method::Symbol;

    loss_crit::Float64 = 0.0,
    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_bounds::Tuple{Float64,Float64} = unscaling.(
        (-9.0, 9.0),
        scale[theta_num]
        ),
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3,
    local_alg::Symbol = :LN_NELDERMEAD,
    kwargs... # other options for get_right_endpoint
    )
    # both endpoints
    endpoints = [get_endpoint(
        theta_init,
        theta_num,
        loss_func,
        method,
        [:left,:right][i]; # method
        loss_crit = loss_crit,
        scale = scale,
        theta_bounds = theta_bounds,
        scan_bound = scan_bounds[i],
        scan_tol = scan_tol,
        loss_tol = loss_tol,
        local_alg = local_alg,
        kwargs... # options for local fitter
        ) for i in 1:2]

    input = ParamIntervalInput(
        theta_init,
        theta_num,
        (x)->x[theta_num],
        loss_func,
        loss_crit,
        scale,
        theta_bounds,
        scan_bounds,
        scan_tol,
        loss_tol,
        local_alg,
        kwargs
        )

    ParamInterval(
        input,
        loss_func(theta_init),
        method,
        Tuple(endpoints)
        )
end

"""
    get_interval(
        theta_init::Vector{Float64},
        scan_func::Function,
        loss_func::Function,
        method::Symbol;

        loss_crit::Float64 = 0.0,
        scale::Vector{Symbol} = fill(:direct, length(theta_init)),
        theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
            fill((-Inf, Inf), length(theta_init)),
            scale
            ),
        scan_bounds::Tuple{Float64,Float64} = unscaling.(
            (-9.0, 9.0),
            :direct
            ),
        scan_tol::Float64 = 1e-3,
        loss_tol::Float64 = 1e-3,
        local_alg::Symbol = :LN_NELDERMEAD,
        kwargs...
        )
Computes confidence interval for function of parameters `scan_func`.

## Return
`ParamInterval` structure storing all input data and estimated confidence interval.

## Arguments
- `theta_init`: starting values of parameter vector ``\\theta``. The starting values should not necessary be the optimum values of `loss_func` but `loss_func(theta_init)` should be lower than `loss_crit`.
- `scan_func`: scan function of parameters.
- `loss_func`: loss function ``\\Lambda\\left(\\theta\\right)`` the profile of which is analyzed. Usually we use log-likelihood for profile analysis in form ``\\Lambda( \\theta ) = - 2 ln\\left( L(\\theta) \\right)``.
- `method`: computational method to estimate confidence interval's endpoint. Currently supports only `:CICO_ONE_PASS` method.

## Keyword arguments
- `loss_crit`: critical level of loss function. Confidence interval's endpoint value is the intersection point of profile likelihood and `loss_crit` level.
- `scale`: vector of scale transformations for each parameters' component. Possible values: `:direct, :log, :logit`. This option can speed up the optimization, especially for wide `theta_bounds`. The default value is `:direct` (no transformation) for all parameters.
- `theta_bounds`: vector of tuple `(lower_bound, upper_bound)` for each parameter. Bounds define the ranges for possible parameter values. Default bounds are `(-Inf,Inf)`.
- `scan_bounds`: scan bounds tuple for `scan_func` values. Default is `(1e-9, 1e9)` .
- `scan_tol`: Absolute tolerance for `theta_num` parameter used as termination criterion.
- `loss_tol`: Absolute tolerance controlling `loss_func` closenes to `loss_crit` (termination criterion). Currently doesn't work for `:CICO_ONE_PASS` method because of limitation in `LN_AUGLAG` interface.
- `local_alg`: algorithm of optimization. Derivative-free and gradient-based algorithms form NLopt package.
- `kwargs...`: the additional `method` specific keyword arguments.
"""
function get_interval(
    theta_init::Vector{Float64},
    scan_func::Function,
    loss_func::Function,
    method::Symbol;

    loss_crit::Float64 = 0.0,
    scale::Vector{Symbol} = fill(:direct, length(theta_init)),
    theta_bounds::Vector{Tuple{Float64,Float64}} = unscaling.(
        fill((-Inf, Inf), length(theta_init)),
        scale
        ),
    scan_bounds::Tuple{Float64,Float64} = (-1e9, 1e9), # log scan bound is not implemented
    scan_tol::Float64 = 1e-3,
    loss_tol::Float64 = 1e-3,
    local_alg::Symbol = :LN_NELDERMEAD,
    kwargs... # other options for get_right_endpoint
    )
    # both endpoints
    endpoints = [get_endpoint(
        theta_init,
        scan_func,
        loss_func,
        method,
        [:left,:right][i]; # method
        loss_crit = loss_crit,
        scale = scale,
        theta_bounds = theta_bounds,
        scan_bound = scan_bounds[i],
        scan_tol = scan_tol,
        loss_tol = loss_tol,
        local_alg = local_alg,
        kwargs... # options for local fitter
        ) for i in 1:2]

    input = ParamIntervalInput(
        theta_init,
        0,
        scan_func,
        loss_func,
        loss_crit,
        scale,
        theta_bounds,
        scan_bounds,
        scan_tol,
        loss_tol,
        local_alg,
        kwargs
        )

    ParamInterval(
        input,
        loss_func(theta_init),
        method,
        Tuple(endpoints)
        )
end


text:
"Context manager to reduce memory usage"
mutable struct Context
    predicates::Dict{Int,Dict{String,Predicate}}
    constants::Dict{Type,Dict{String,Constant}}
    variables::Dict{Type,Dict{String,Variable}}
    functors::Dict{Int,Dict{String,Functor}}
    types::Dict{String,Type}
    propositions::Dict{String,Proposition}

    function Context()
        return new(Dict(), Dict(), Dict(), Dict(), Dict(), Dict())
    end
end

const global_context = Context()

"constructs a new type within context"
function c_type!(name::String, context::Context)
    if !haskey(context.types, name)
        context.types[name] = HerbLanguage.Type(name)
    end

    context.types[name]
end

c_type(name::String) = c_type!(name, global_context)

"Creates constant within the context"
function c_const!(name::String, type::Type, context::Context)
    if !haskey(context.constants, type)
        context.constants[type] = Dict{String, Constant}()
    end

    if !haskey(context.constants[type], name)
        context.constants[type][name] = Constant(name, type)
    end

    context.constants[type][name]
end

c_const(name::String, type::Type) = c_const!(name, type, global_context)
c_const!(name::String, context::Context) = c_const!(name, c_type!("thing", context), context)
c_const(name::String) = c_const!(name, c_type(name), global_context)

"Creates a variable within a context"
function c_var!(name::String, type::Type, context::Context)
    if !haskey(context.variables, type)
        context.variables[type] = Dict{String,Variable}()
    end

    if !haskey(context.variables[type], name)
        context.variables[type][name] = Variable(name, type)
    end

    context.variables[type][name]
end

c_var(name::String, type::String) = c_var!(name, type, global_context)
c_var!(name::String, context::Context) = c_var!(name, c_type!("thing", context), context)
c_var(name::String) = c_var!(name, c_type("thing"), global_context)

"Constructs a predicate within context"
function c_pred!(name::String, arity::Int, arg_types::Vector{Type}, context::Context)
    if !haskey(context.predicates, arity)
        context.predicates[arity] = Dict{String,Predicate}()
    end

    if !haskey(context.predicates[arity], name)
        context.predicates[arity][name] = Predicate(name, arity, arg_types)
    end

    context.predicates[arity][name]
end

c_pred!(name::String, arity::Int, context::Context) = c_pred!(name, arity, fill(c_type!("thing", context), (arity,)), context)
c_pred!(name::String, arg_types::Vector{Type}, context::Context) = c_pred!(name, length(arg_types), arg_types, context)
c_pred(name::String, arity::Int) = c_pred!(name, arity, fill(c_type("thing"), (arity,)), global_context)
c_pred(name::String, arg_types::Vector{Type}) = c_pred!(name, length(arg_types), arg_types, global_context)

"Construct functor"
function c_functor!(name::String, arity::Int, context::Context)
    if !haskey(context.functors, arity)
        context.functors[arity] = Dict{String, Functor}()
    end

    if !haskey(context.functors[arity], name)
        context.functors[arity][name] = Functor(name, arity)
    end

    context.functors[arity][name]
end

c_functor(name::String, arity::Int) = c_functor!(name, arity, global_context)

"Construct proposition"
function c_prop!(name::String, context::Context)
    if !haskey(context.propositions, name)
        context.propositions[name] = Proposition(name)
    end

    context.propositions[name]
end

c_prop(name) = c_prop!(name, global_context)

ext: jl


text:
"""
For a given grid size, this function runs the greedy_multi_path routing algorithm
on randomly placed end-user pairs.
"""

using QuNet
using LightGraphs, SimpleWeightedGraphs, GraphPlot, MetaGraphs, Plots, Colors, Statistics
using LaTeXStrings
using JLD
using Parameters

function asymptotic_costs(gridsize::Int64)
    N = 1000
    G = GridNetwork(gridsize, gridsize)
    T = QuNet.TemporalGraph(G, 5, memory_costs = unit_costvector())
    p, dum1, dum2, dum3 = net_performance(T, N, 1, max_paths=4)
    return p["loss"], p["Z"]
end

datafile = "data/timedepth"

# Params
# 15
# 50
# 10
# 1000 trials
num_trials = 100::Int64
max_depth = 15::Int64
num_pairs = 50::Int64
grid_size = 10::Int64
async_cost = 1000

generate_new_data = true
if generate_new_data == true

    perf_data = []
    perf_err = []
    path_data = []
    path_err = []

    for i in 1:max_depth
        println("Collecting for time depth $i")
        G = GridNetwork(grid_size, grid_size)
        # Create a Temporal Graph from G with timedepth i
        T = QuNet.TemporalGraph(G, i, memory_costs = zero_costvector())
        # Get data
        p, p_e, pat, pat_e = net_performance(T, num_trials, num_pairs, max_paths=4,
        async_cost=async_cost)
        push!(perf_data, p)
        push!(perf_err, p_e)
        push!(path_data, pat)
        push!(path_err, pat_e)
    end

    # Collect data for horizontal lines:
    println("Collecting data for asymptote")
    as = asymptotic_costs(grid_size)
    e_as = ones(length(1:max_depth)) * as[1]
    f_as = ones(length(1:max_depth)) * as[2]

    # Save data
    d = Dict{Symbol, Any}()
    @pack! d = num_trials, max_depth, num_pairs, grid_size, perf_data, perf_err, path_data, path_err, e_as, f_as
    save("$datafile.jld", "data", d)

else
    # Load data
    if !isfile("$datafile.jld")
        error("$datafile.jld not found in working directory")
    end
    d = load("$datafile.jld")["data"]
    @unpack num_trials, max_depth, num_pairs, grid_size, perf_data, perf_err, path_data, path_err, e_as, f_as = d
end

asymptotic_costs(10)

# Get values for x axis
x = collect(1:max_depth)

# Extract from performance data
loss = collect(map(x->x["loss"], perf_data))
z = collect(map(x->x["Z"], perf_data))
loss_err = collect(map(x->x["loss"], perf_err))
z_err = collect(map(x->x["Z"], perf_err))

# Extract from path data
P0 = [path_data[i][1]/num_pairs for i in 1:max_depth]
P1 = [path_data[i][2]/num_pairs for i in 1:max_depth]
P2 = [path_data[i][3]/num_pairs for i in 1:max_depth]
P3 = [path_data[i][4]/num_pairs for i in 1:max_depth]
P4 = [path_data[i][5]/num_pairs for i in 1:max_depth]

P0e = [path_err[i][1]/num_pairs for i in 1:max_depth]
P1e = [path_err[i][2]/num_pairs for i in 1:max_depth]
P2e = [path_err[i][3]/num_pairs for i in 1:max_depth]
P3e = [path_err[i][4]/num_pairs for i in 1:max_depth]
P4e = [path_err[i][5]/num_pairs for i in 1:max_depth]

# Plot
# after seriestype: marker = (5)
plot(x, loss, ylims=(0,1), xlims=(0,15), seriestype = :scatter, yerror = loss_err, label=L"$\eta$",
legend=:right, xguidefontsize=14, tickfontsize=12, legendfontsize=10, fontfamily="computer modern",
markersize=5, color=:DodgerBlue, markershape =:utriangle)
plot!(x, z, seriestype = :scatter, yerror = z_err, label=L"$F$", markersize=5, color=:Crimson)

# Plot asymptote
plot!(x, e_as, linestyle=:dot, color=:blue, linewidth=2, label=L"$\textrm{Asymptotic } \eta$")
plot!(x, f_as, linestyle=:dot, color=:red, linewidth=2, label=L"$\textrm{Asymptotic } F$")
xaxis!(L"$\textrm{Time Depth of Tempral Meta-Graph}$")

savefig("plots/cost_temporal.png")
savefig("plots/cost_temporal.pdf")

plot(x, P0, ylims=(0,1), xlims=(0,15), seriestype=:scatter, yerr = P0e, label=L"$P_0$", legend= :right,
xguidefontsize=14, tickfontsize=12, legendfontsize=10, fontfamily="computer modern",
markersize=5, color_palette = palette(:plasma, 10))
plot!(x, P0, linewidth=1, label=false)
plot!(x, P1, seriestype=:scatter, yerr = P1e, label=L"$P_1$", markersize=5)
plot!(x, P1, linewidth=1, label=false)
plot!(x, P2, seriestype=:scatter, yerr = P2e, label=L"$P_2$", markersize=5)
plot!(x, P2, linewidth=1, label=false)
plot!(x, P3, seriestype=:scatter, yerr = P3e, label=L"$P_3$", markersize=5)
plot!(x, P3, linewidth=1, label=false)
plot!(x, P4, seriestype=:scatter, yerr = P4e, label=L"$P_4$", markersize=5)
plot!(x, P4, linewidth=1, label=false)
xaxis!(L"$\textrm{Time Depth of Temporal Meta-Graph}$")

savefig("plots/path_temporal.png")
savefig("plots/path_temporal.pdf")

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
import ..SamplingConstraints:
    SamplingConstraint,
    SequentialSamplingConstraint,
    OrderedSamplingAlgorithm,
    StartToEnd,
    StrictlyIncreasing,
    sequence_exists,
    TruncateMinimum,
    TruncateRange,
    TruncateQuantiles,
    constrain

export resample

import IntervalArithmetic:
    interval

function _draw!(s, x, c::StrictlyIncreasing{<:StartToEnd}, mins, maxs)
    L = length(x)

    # TODO: add slight margin?
    for i = 1:L
        if i == 1
            hi = minimum(maxs[2:end])
            truncated_distribution = truncate(x[i], TruncateMaximum(hi))
            s[i] = resample(truncated_distribution)
        end

        if 1 < i < L
            lo = max(s[i - 1], mins[i])
            hi = min(maxs[i], minimum(maxs[i+1:end]))

            lo <= hi || error("Truncation range invalid for point $i. Got lo < hi ($lo < $hi), which should be impossible.")

            truncated_distribution = truncate(x[i], TruncateRange(lo, hi))
            s[i] = resample(truncated_distribution)
        end

        if i == L
            lo = max(s[i - 1], mins[i])
            truncated_distribution = truncate(x[i], TruncateMinimum(lo))
            s[i] = resample(truncated_distribution)

        end
    end

    return s
end

"""
    _draw(x, c::StrictlyIncreasing{StartToEnd}, mins, maxs)

Sample `x` in a strictly increasing manner, given pre-computed minimum and maximum
values for each distribution in `c`.

Implicitly assumes a strictly increasing sequence exists, but does not check that condition.
"""
function _draw(x, c::StrictlyIncreasing{<:StartToEnd}, mins, maxs)
    L = length(x)
    samples = zeros(Float64, L) # a vector to hold the element-wise samples

    _draw!(samples, x, c, mins, maxs)
end


text:
# ---
# title: 769. Max Chunks To Make Sorted
# id: problem769
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/max-chunks-to-make-sorted/description/>
# hidden: true
# ---
#
# Given an array `arr` that is a permutation of `[0, 1, ..., arr.length - 1]`,
# we split the array into some number of "chunks" (partitions), and individually
# sort each chunk.  After concatenating them, the result equals the sorted
# array.
#
# What is the most number of chunks we could have made?
#
# **Example 1:**
#
#
#
#     Input: arr = [4,3,2,1,0]
#     Output: 1
#     Explanation:
#     Splitting into two or more chunks will not return the required result.
#     For example, splitting into [4, 3], [2, 1, 0] will result in [3, 4, 0, 1, 2], which isn't sorted.
#
#
# **Example 2:**
#
#
#
#     Input: arr = [1,0,2,3,4]
#     Output: 4
#     Explanation:
#     We can split into two chunks, such as [1, 0], [2, 3, 4].
#     However, splitting into [1, 0], [2], [3], [4] is the highest number of chunks possible.
#
#
# **Note:**
#
#   * `arr` will have length in range `[1, 10]`.
#   * `arr[i]` will be a permutation of `[0, 1, ..., arr.length - 1]`.
#
#
#
#
## @lc code=start
using LeetCode

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


text:
# Descriptor of a cartesian grid

"""
    struct CartesianDescriptor{D,T,F<:Function}
      origin::Point{D,T}
      sizes::NTuple{D,T}
      partition::NTuple{D,Int}
      map::F
    end

Struct that stores the data defining a Cartesian grid.
"""
struct CartesianDescriptor{D,T,F<:Function} <: GridapType
  origin::Point{D,T}
  sizes::NTuple{D,T}
  partition::NTuple{D,Int}
  map::F
  isperiodic::NTuple{D,Bool}
  @doc """
      CartesianDescriptor(
        origin::Point{D},
        sizes::NTuple{D},
        partition;
        map::Function=identity,
        isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D

  `partition` is a 1D indexable collection of arbitrary type.
  """
  function CartesianDescriptor(
    origin::Point{D},
    sizes::NTuple{D},
    partition;
    map::Function=identity,
    isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D

    for i in 1:D
      if isperiodic[i]
        @assert (partition[i] > 2) "A minimum of 3 elements is required in any
          periodic direction"
      end
    end

    T = eltype(sizes)
    F = typeof(map)
    new{D,T,F}(origin,sizes,Tuple(partition),map,isperiodic)
  end
end

"""
    CartesianDescriptor(
      domain,
      partition;
      map::Function=identity,
      isperiodic::NTuple{D,Bool}=tfill(false,Val{D}))

`domain` and `partition` are 1D indexable collections of arbitrary type.
"""
function CartesianDescriptor(
  domain,
  partition;
  map::Function=identity,
  isperiodic::NTuple=tfill(false,Val{length(partition)}()))

  D = length(partition)
  limits = [(domain[2*d-1],domain[2*d]) for d in 1:D]
  sizes = Tuple([(limits[d][2]-limits[d][1])/partition[d] for d in 1:D])
  origin = Point([ limits[d][1] for d in 1:D]...)
  CartesianDescriptor(origin,sizes,partition;map=map,isperiodic=isperiodic)
end

"""
    CartesianDescriptor(
      pmin::Point{D},
      pmax::Point{D},
      partition;
      map::Function=identity,
      isperiodic::NTuple{D,Bool}=tfill(false,Val{D})) where D

`partition` is a 1D indexable collection of arbitrary type.
"""
function CartesianDescriptor(
  pmin::Point{D},
  pmax::Point{D},
  partition;
  map::Function=identity,
  isperiodic::NTuple{D,Bool}=tfill(false,Val{D}())) where D

  T = eltype(pmin)
  domain = zeros(T,2*D)
  for d in 1:D
    domain[2*(d-1)+1] = pmin[d]
    domain[2*(d-1)+2] = pmax[d]
  end
  CartesianDescriptor(domain,partition;map=map,isperiodic=isperiodic)
end

# Deprecated signatures for backwards compatibility

@deprecate CartesianDescriptor(origin::Point{D}, sizes::NTuple{D}, partition, map::Function) where D CartesianDescriptor(origin, sizes, partition; map=map)

@deprecate CartesianDescriptor(domain,partition,map::Function) CartesianDescriptor(domain,partition;map=map)

@deprecate CartesianDescriptor(pmin::Point{D}, pmax::Point{D}, partition, map::Function) where D CartesianDescriptor(pmin, pmax, partition; map=map)

# Coordinates

struct CartesianCoordinates{D,T,F} <: AbstractArray{Point{D,T},D}
  data::CartesianDescriptor{D,T,F}
end

Base.size(a::CartesianCoordinates) = a.data.partition .+ 1

Base.IndexStyle(::Type{<:CartesianCoordinates}) = IndexCartesian()

function Base.getindex(a::CartesianCoordinates{D,T}, I::Vararg{Integer,D}) where {D,T}
  p = zero(Mutable(Point{D,T}))
  x0 = a.data.origin
  dx = a.data.sizes
  @inbounds for d in 1:D
    p[d] =  x0[d] + (I[d]-1)*dx[d]
  end
  a.data.map(Point(p))
end

# Cell nodes

struct CartesianCellNodes{D} <: AbstractArray{Vector{Int32},D}
  partition::NTuple{D,Int}
  function CartesianCellNodes(partition)
    D = length(partition)
    new{D}(partition)
  end
end

Base.size(a::CartesianCellNodes) = a.partition

Base.IndexStyle(::Type{<:CartesianCellNodes}) = IndexCartesian()

function array_cache(a::CartesianCellNodes{D}) where D
  zeros(Int32,2^D)
end

#function getindex!(cache,a::CartesianCellNodes,i::Integer)
#  cis = CartesianIndices(size(a))
#  ci = cis[i]
#  getindex!(cache,a,Tuple(ci)...)
#end
#
#function getindex!(v,a::CartesianCellNodes{1},i::Integer)
#  v[1] = i
#  v[2] = i+1
#  v
#end

function getindex!(v,a::CartesianCellNodes{D},i::Vararg{Integer,D}) where D
  nodes = LinearIndices(size(a).+1)
  lnodes = CartesianIndices(tfill(Int32(2),Val{D}()))
  j = i .- 1
  for (i,lnode) in enumerate(lnodes)
    k = Tuple(lnode) .+ j
    v[i] = nodes[ k... ]
  end
  v
end

function Base.getindex(a::CartesianCellNodes{D},i::Vararg{Integer,D}) where D
  cache = array_cache(a)
  getindex!(cache,a,i...)
end

# Cartesian Grid

"""
    struct CartesianGrid{D,T,F} <: Grid{D,D}
      # private fields
    end
"""
struct CartesianGrid{D,T,F} <: Grid{D,D}
  node_coords::CartesianCoordinates{D,T,F}
  cell_node_ids::CartesianCellNodes{D}
  cell_type::Fill{Int8,1,Tuple{Base.OneTo{Int}}}
  @doc """
      CartesianGrid(desc::CartesianDescriptor)
  """
  function CartesianGrid(desc::CartesianDescriptor{D,T,F}) where {D,T,F}
    node_coords = CartesianCoordinates(desc)
    cell_node_ids = CartesianCellNodes(desc.partition)
    cell_type = Fill(Int8(1),length(cell_node_ids))
    new{D,T,F}(node_coords,cell_node_ids,cell_type)
  end
end

OrientationStyle(::Type{<:CartesianGrid}) = Oriented()

"""
    get_cartesian_descriptor(grid::CartesianGrid)

Get the descriptor of the Cartesian grid
"""
function get_cartesian_descriptor(a::CartesianGrid)
  a.node_coords.data
end

get_node_coordinates(g::CartesianGrid) = g.node_coords

get_cell_type(g::CartesianGrid) = g.cell_type

get_cell_node_ids(g::CartesianGrid) = g.cell_node_ids

function get_cell_map(grid::CartesianGrid{D,T,typeof(identity)} where {D,T})
  CartesianMap(grid.node_coords.data)
end

function get_reffes(g::CartesianGrid{D}) where D
  p = Polytope(tfill(HEX_AXIS,Val{D}()))
  order = 1
  reffe = LagrangianRefFE(Float64,p,order)
  [reffe,]
end

get_reffes(g::CartesianGrid{1}) = [SEG2,]

get_reffes(g::CartesianGrid{2}) = [QUAD4,]

get_reffes(g::CartesianGrid{3}) = [HEX8,]

"""
    CartesianGrid(args...;kwargs...)

Same args needed to construct a `CartesianDescriptor`
"""
function CartesianGrid(args...;kwargs...)
  desc = CartesianDescriptor(args...;kwargs...)
  CartesianGrid(desc)
end

# Cell map

struct CartesianMap{D,T,L} <: AbstractArray{AffineMap{D,D,T,L},D}
  data::CartesianDescriptor{D,T,typeof(identity)}
  function CartesianMap(des::CartesianDescriptor{D,T}) where {D,T}
    L = D*D
    new{D,T,L}(des)
  end
end

Base.size(a::CartesianMap) = a.data.partition

Base.IndexStyle(::Type{<:CartesianMap}) = IndexCartesian()

function Base.getindex(a::CartesianMap{D,T},I::Vararg{Integer,D}) where {D,T}
  p = zero(Mutable(Point{D,T}))
  x0 = a.data.origin
  dx = a.data.sizes
  @inbounds for d in 1:D
    p[d] =  x0[d] + (I[d]-1)*dx[d]
  end
  origin =  Point(p)
  grad = diagonal_tensor(VectorValue(dx))
  AffineMap(grad,origin)
end

function lazy_map(::typeof(∇),a::CartesianMap)
  dx = a.data.sizes
  grad = diagonal_tensor(VectorValue(dx))
  j = ConstantField(grad)
  Fill(j,size(a))
end

function lazy_map(::typeof(∇),a::LazyArray{<:Fill{<:Reindex{<:CartesianMap}}})
  i_to_map = a.maps.value.values
  j_to_i = a.args[1]
  i_to_grad = lazy_map(∇,i_to_map)
  lazy_map(Reindex(i_to_grad),j_to_i)
end


text:
using BHAtp, Statistics

ProjDir = @__DIR__
!isdir(joinpath(ProjDir, "plots")) && mkdir(joinpath(ProjDir, "plots"))

# Actual data

incs = 45:5:55
wobs = 20:20:60
tpf = DataFrame()
for i in wobs
  for j in incs
    append!(tpf, DataFrame(incl=j, wob=i))
  end
end

tpf[!, :tp] = [1.55, 1.70, 1.85, 1.65, 1.83, 2.00, 1.75, 1.95, 2.15]
tpf[!, :ft] = [2200, 2400, 2600, 2400, 2600, 2800, 2650, 2750, 2950]

actual = DataFrame(
  incls = [49.75, 50.75, 51.50, 52.80, 54, 55],
  atp = [1.05, 1.07, 1.20, 1.15, 1.60, 1.22]
)

name = "figure_07_2"

figure_07_2 = @pgf GroupPlot({
    group_style = {
      group_size = "1 by 1",
      vertical_sep = "0pt",
      xticklabels_at = "edge bottom"
    }
  },
  Axis({
      height = "6cm", width = " 8cm",
      xmin = 42, xmax = 56,
      ymin = 0, ymax = 2.5,
      grid = "major",
      xlabel = "Inclination [°]",
      ylabel = "TP [°/100ft]",
      title = "Figure 7:Actual vs. Isotropic TP",
      legend_pos = "south west",
    },
    Plot({ color="blue", mark = "." },
      Table(tpf.incl[1:3], tpf.tp[1:3])),
    LegendEntry("TP wob=20"),
    Plot({ color="gray", mark = "." },
      Table(tpf.incl[4:6], tpf.tp[4:6])),
    LegendEntry("TP wob=40"),
    Plot({ color="green", mark = "." },
      Table(tpf.incl[7:9], tpf.tp[7:9])),
    LegendEntry("TP wob=60"),
    Plot({ thick, color="black", mark = "+", },
      Table(actual.incls, actual.atp)),
    LegendEntry("TP, actual")))

savefig(joinpath(ProjDir, "plots", name), figure_07_2, ProjDir)


text:
using AbstractGPs
using ConjugateComputationVI
using Distributions
using Optim
using ParameterHandling
using Plots
using Random
using StatsBase
using StatsFuns
using TemporalGPs
using Zygote

using ConjugateComputationVI: GaussHermiteQuadrature, UnivariateFactorisedLikelihood

# Adjoint for the Poisson logpdf.
# Someone should implement this in the correct place.
# log(λ^x exp(-λ) / x!) =
# x log(λ) - λ - log(x!)
# dλ = x / λ - 1
Zygote.@adjoint function StatsFuns.poislogpdf(λ::Float64, x::Union{Float64, Int})
    function poislogpdf_pullback(Δ::Real)
        return Δ * (x / λ - 1), nothing
    end
    return StatsFuns.poislogpdf(λ, x), poislogpdf_pullback
end

function Zygote._pullback(ctx::Zygote.AContext, ::typeof(logpdf), d::Poisson, x)
    out, pb = Zygote._pullback(ctx, poislogpdf, d.λ, x)
    function logpdf_pullback(Δ::Real)
        _, Δλ, Δx = pb(Δ)
        return nothing, (λ=Δλ, ), Δx
    end
    return out, logpdf_pullback
end

# Specify a model.
θ_init = (scale=positive(1.0), stretch=positive(1.0));
θ_init_flat, unflatten = ParameterHandling.flatten(θ_init);

function build_latent_gp(θ::AbstractVector{<:Real})
    return build_latent_gp(ParameterHandling.value(unflatten(θ)))
end
function build_latent_gp(θ::NamedTuple)
    gp_naive = GP(θ.scale * AbstractGPs.transform(Matern52Kernel(), θ.stretch))
    gp = to_sde(gp_naive, SArrayStorage(Float64))
    lik = UnivariateFactorisedLikelihood(f -> Poisson{Float64}(exp(f)))
    return LatentGP(gp, lik, 1e-9)
end

# Generate synthetic data.
x = RegularSpacing(0.0, 0.001, 10_000);
y = rand(build_latent_gp(θ_init_flat)(x)).y;

f_approx_post, results_summary = ConjugateComputationVI.optimize_elbo(
    build_latent_gp,
    GaussHermiteQuadrature(10),
    x,
    y,
    θ_init_flat,
    BFGS(
        alphaguess = Optim.LineSearches.InitialStatic(scaled=true),
        linesearch = Optim.LineSearches.BackTracking(),
    ),
    Optim.Options(
        show_trace = true,
        iterations=25,
        f_calls_limit=50,
    ),
);

function approx_post_marginal_samples(x::AbstractVector, N::Int)
    return map(marginals(f_approx_post(x))) do latent_marginal
        f = rand(latent_marginal, N)
        return exp.(f)
    end
end

function approx_post_95_CI(x::AbstractVector, N::Int)
    samples = approx_post_marginal_samples(x, N)
    return quantile.(samples, Ref((0.025, 0.5, 0.975)))
end

# Plot the predictions.
qs = approx_post_95_CI(x, 100);
p1 = plot(
    x, getindex.(qs, 1);
    linealpha=0,
    fillrange=getindex.(qs, 3),
    label="95% CI",
    fillalpha=0.3,
    color=:blue,
);
scatter!(p1, x, y; markersize=1, label="Observations");

p2 = plot(f_approx_post(x, 1e-6); ribbon_scale=3, color=:blue, label="approx post latent");

plot(p1, p2; layout=(2, 1))


text:
# @enum Side NoneSide UpperSide LowerSide
type NoneSide <: BC_Side end
type UpperSide <: BC_Side end
type LowerSide <: BC_Side end

type DefaultBoundaryCondition <: BoundaryConditionType end
type NeumannBCType <: BoundaryConditionType end
type DirichletBCType <: BoundaryConditionType end

type BoundaryCondition{B <: BoundaryConditionType, S <: BC_Side}
  condType::B
  value::Float64
  side::S
end

typealias NeumannBC{S} BoundaryCondition{NeumannBCType, S}
build_NeumannBC(value::Float64, side::BC_Side) = BoundaryCondition(NeumannBCType(), value, side)

set_time!(bc::NeumannBC, ::Float64) = bc # returns self

apply_before_applying!(bc::NeumannBC, L::TridiagonalOperator) = apply_before_applying!(bc, bc.side, L)

function apply_before_applying!(bc::NeumannBC, ::LowerSide, L::TridiagonalOperator)
  set_first_row!(L, -1.0, 1.0)
  return L
end

function apply_before_applying!(bc::NeumannBC, ::UpperSide, L::TridiagonalOperator)
  set_last_row!(L, -1.0, 1.0)
  return L
end

apply_after_applying!(bc::NeumannBC, u::Vector{Float64}) = apply_after_applying!(bc, bc.side, u)

function apply_after_applying!(bc::NeumannBC, ::LowerSide, u::Vector{Float64})
  u[1] = u[2] - bc.value
  return u
end

function apply_after_applying!(bc::NeumannBC, ::UpperSide, u::Vector{Float64})
  u[end] = u[end-1] + bc.value
  return u
end

apply_before_solving!(bc::NeumannBC, L::TridiagonalOperator, rhs::Vector{Float64}) = apply_before_solving!(bc, bc.side, L, rhs)

function apply_before_solving!(bc::NeumannBC, ::LowerSide, L::TridiagonalOperator, rhs::Vector{Float64})
  set_first_row!(L, -1.0, 1.0)
  rhs[1] = bc.value
  return L, rhs
end

function apply_before_solving!(bc::NeumannBC, ::UpperSide, L::TridiagonalOperator, rhs::Vector{Float64})
  set_last_row!(L, -1.0, 1.0)
  rhs[end] = bc.value
  return L, rhs
end

apply_after_solving!(bc::NeumannBC, ::Vector{Float64}) = bc # returns self


text:
struct CompiledForward{T<:AbstractTape} <: AbstractTape
    tape::T
    forward_exec::Vector{FunctionWrapper{Nothing, Tuple{}}}
end

Base.show(io::IO, t::CompiledForward) = print(io, typeof(t).name, "($(t.tape.func))")

function CompiledForward(t::T) where T<:AbstractTape
    CompiledForward{T}(t, [FunctionWrapper{Nothing, Tuple{}}(ForwardExecutor(instruction)) for instruction in t.tape])
end

Base.length(ct::CompiledForward) = length(ct.tape)

@inline func_hook(ct::CompiledForward) = func_hook(ct.tape)

@inline input_hook(ct::CompiledForward) = input_hook(ct.tape)

@inline output_hook(ct::CompiledForward) = output_hook(ct.tape)

function forward_pass!(compiled_tape::CompiledForward)
    for wrapper in compiled_tape.forward_exec
        wrapper()
    end
    nothing
end

function compile(f, args...)
    tp = GradientTape(f, (args))
    cp = CompiledForward(tp)
end

value(x::TrackedArray) = x.value
value(x::Tuple) = value.(x)
value(x) = x

function (tp::CompiledForward)(in...)
  ReverseDiff.seeded_forward_pass!(tp, in)
  tp.tape.output |> value
end


text:
import JuliaInterpreter
import JuliaInterpreter: enter_call, step_expr!, @lookup, Frame, SSAValue, SlotNumber

getexpr(frame::Frame, pc::Int) = frame.framecode.src.code[pc]
current_expr(frame::Frame) = getexpr(frame, frame.pc)

"""
Split JuliaInterpreter call expression into a tuple of 2 elements:

 * callable (e.g. function or callable struct) & its arguments
 * vars on the tape corresponding to these function & args

If arguments include free parameters (not SlotNumber or SSAValue), these are recorded
to the tape as constants
"""
function split_int_call!(tape::Tape, frame::Frame, frame_vars::Dict, ex)
    arr = Meta.isexpr(ex, :(=)) ? ex.args[2].args : ex.args
    # for whatever reason JuliaInterpreter wraps some nodes in the original code into QuoteNode
    arr = [isa(x, QuoteNode) ? x.value : x for x in arr]
    cfargs = [x isa Symbol ? x : @lookup(frame, x) for x in arr]
    cvars = Vector{Int}(undef, length(cfargs))
    for (i, x) in enumerate(arr)
        # if isa(x, JuliaInterpreter.SlotNumber) || isa(x, JuliaInterpreter.SSAValue)
        if haskey(frame_vars, x)
            cvars[i] = frame_vars[x]
        else
            val = x isa Symbol ? x : @lookup(frame, x)
            id = record!(tape, Constant, val)
            cvars[i] = id
            if val != x
                # if constant appeared to be a SlotNumber or SSAValue
                # store its mapping into frame_vars
                frame_vars[x] = id
            end
        end
    end
    return cfargs, cvars
end

"""
Given a Frame and current expression, extract LHS location (SlotNumber or SSAValue)
"""
get_location(frame::Frame, ex) = Meta.isexpr(ex, :(=)) ? ex.args[1] : JuliaInterpreter.SSAValue(frame.pc)

is_int_call_expr(ex) = Meta.isexpr(ex, :call) || (Meta.isexpr(ex, :(=)) && Meta.isexpr(ex.args[2], :call))
is_int_assign_expr(ex) = Meta.isexpr(ex, :(=)) && (isa(ex.args[2], SlotNumber) || isa(ex.args[2], SSAValue))

is_interesting_expr(ex) = is_int_call_expr(ex) || is_int_assign_expr(ex) || Meta.isexpr(ex, :return)

function itrace!(tape::Tape, fargs, vars; primitives)
    frame = enter_call(fargs...)
    frame_vars = Dict{Any, Int}(JuliaInterpreter.SlotNumber(i) => v for (i, v) in enumerate(vars))
    # f might be a callable struct, so we need to record it and add to frame_vars
    # f_id = record!(tape, Constant, f)  # should not be constant! should caller record it?
    # frame_vars[SlotNumber(1)] = f
    is_interesting_expr(current_expr(frame)) || step_expr!(frame)  # skip non-call expressions
    ex = current_expr(frame)
    while !Meta.isexpr(ex, :return)
        if is_int_assign_expr(ex)
            lhs, rhs = ex.args
            frame_vars[lhs] = frame_vars[rhs]
            step_expr!(frame)
        elseif is_int_call_expr(ex)
            # read as "current function & arguments", "current variables"
            cfargs, cvars = split_int_call!(tape, frame, frame_vars, ex)
            cf = cfargs[1]
            loc = get_location(frame, ex)
            # there are several special cases such as NamedTuples and constructors
            # we replace these with calls to special helper functions
            if cf isa UnionAll && cf <: NamedTuple
                # replace cf with namedtuple function, adjust arguments
                names = collect(cf.body.parameters)[1]
                cf = namedtuple
                insert!(cfargs, 2, names)
                names_var_id = record!(tape, Constant, names)
                cvars = [names_var_id; cvars]
            elseif cf isa DataType
                # constructor, replace with a call to __new__ which we know how to differentiate
                T = cf
                cf = __new__
                insert!(cfargs, 2, T)
                T_var_id = record!(tape, Constant, T)
                cvars = [T_var_id; cvars]
            elseif cf == Base.tuple
                cf = __tuple__
            elseif cf == Base.getfield || (cf == Base.getindex && isa(tape[cvars[2]].val, NamedTuple))
                # similar to constuctors, there's a special case for __getfield__ in backprop
                cf = __getfield__
            end
            # if current function is a primitive of a built-in, write it to the tape
            # otherwise recurse into the current function
            if cf in primitives || isa(cf, Core.Builtin) || isa(cf, Core.IntrinsicFunction)
                step_expr!(frame)
                retval = @lookup(frame, loc)
                ret_id = record!(tape, Call, retval, cf, cvars[2:end])
                frame_vars[loc] = ret_id  # for slots it may overwrite old mapping
            else
                try
                    retval, ret_id = itrace!(tape, cfargs, cvars; primitives=primitives)
               catch
                    println("Failed to trace through function $cf")
                    rethrow()
                end
                frame_vars[loc] = ret_id  # for slots it may overwrite old mapping
                step_expr!(frame)  # can we avoid this double execution?
            end
        else
            step_expr!(frame)
        end
        ex = current_expr(frame)
    end
    retval = @lookup(frame, ex.args[1])
    ret_id = frame_vars[ex.args[1]]
    return retval, ret_id  # return var ID of a result variable
end

"""
Trace function f with arguments args using JuliaInterpreter
"""
function itrace(f, args...; primitives=PRIMITIVES, optimize=true)
    tape = Tape(guess_device(args))
    # record arguments as input variables
    fargs = Vector(undef, length(args) + 1)
    vars = Vector{Int}(undef, length(args) + 1)
    for (i, arg) in enumerate([f, args...])
        id = record!(tape, Input, arg)
        fargs[i] = arg
        vars[i] = i
    end
    val, resultid = itrace!(tape, fargs, vars; primitives=primitives)
    tape.resultid = resultid
    if optimize
        tape = simplify(tape)
    end
    return val, tape
end


text:
using DiffEqEnvironments
using DifferentialEquations
using ReinforcementLearning
using Test

"""
This test defines an autonomous system and checks
whether the integration and terminal condition work
when stepping through the environment.
"""

# Define ODEProblem of scalar autonomous system
ode(s,a,t) = 1.01 * s
s0 = 1 / 2
tspan = (1f0, 2f0)
problem = ODEProblem(ode, s0, tspan)

# Solve ODE Problem for ground truth
sol = solve(problem, Tsit5())

# Define DiffEqEnv
r = ASReward((a, s) -> sum(abs2, s))
n_actions = 1
dt = 0.1
env = DiffEqEnv(problem, r, n_actions, dt; T=Float32, solver=Tsit5())

#= ==== Testing starts here ==== =#
# Test initialization
@test env.t == problem.tspan[1]
@test env.state[1] == s0

# Step through env until done
while !env.done
    env(0f0) # action doesn't matter for autonomous system
end

@test env.steps == (tspan[2] - tspan[1]) / dt
@test env.t ≈ problem.tspan[2]
@test env.state[1] ≈ sol.u[end]

#= ==== Reset and try again ==== =#
# Test initialization
RLBase.reset!(env)
@test env.t == problem.tspan[1]
@test env.state[1] == s0

# Step through env until done
while !env.done
    env(0f0) # action doesn't matter for autonomous system
end

@test env.steps == (tspan[2] - tspan[1]) / dt
@test env.t ≈ problem.tspan[2]
@test env.state[1] ≈ sol.u[end]

#= = Remake env using different types T = =#
env64 = DiffEqEnv(problem, r, n_actions, dt, T=Float64)
@test env64.state[1] == s0

env16 = DiffEqEnv(problem, r, n_actions, dt, T=Float16)
@test env16.state[1] == s0

#= = Run random policy using ReinforcementLearning.jl = =#
RLBase.reset!(env)
hook = TotalRewardPerEpisode()
run(
    Agent(
        ;policy=RandomPolicy(env),
        trajectory=VectCompactSARTSATrajectory(
            state_type=Any,
            action_type=Any,
            reward_type=Real,
            terminal_type=Bool,
        ),
    ),
    env,
    StopAfterEpisode(10),
    hook
)


text:
module TestCSetDataStructures
using Test

using Catlab: @present, generator
using Catlab.Theories: compose, id
using Catlab.CSetDataStructures

# Discrete dynamical systems
############################

@present TheoryDDS(FreeSchema) begin
  X::Ob
  Φ::Hom(X,X)
end

const AbstractDDS = AbstractCSetType(TheoryDDS)
const DDS = CSetType(TheoryDDS, index=[:Φ])
@test DDS <: AbstractDDS
@test DDS <: CSet

dds = DDS()
@test keys(tables(dds)) == (:X,)
@test keys(dds.indices) == (:Φ,)
@test nparts(dds, :X) == 0
@test add_part!(dds, :X) == 1
@test nparts(dds, :X) == 1
@test incident(dds, 1, :Φ) == []

set_subpart!(dds, 1, :Φ, 1)
@test subpart(dds, 1, :Φ) == 1
@test incident(dds, 1, :Φ) == [1]

@test add_part!(dds, :X, Φ=1) == 2
@test add_part!(dds, :X, Φ=1) == 3
@test subpart(dds, :Φ) == [1,1,1]
@test subpart(dds, [2,3], :Φ) == [1,1]
@test incident(dds, 1, :Φ) == [1,2,3]

@test has_part(dds, :X)
@test !has_part(dds, :nonpart)
@test has_part(dds, :X, 3)
@test !has_part(dds, :X, 4)
@test has_part(dds, :X, 1:5) == [true, true, true, false, false]

@test has_subpart(dds, :Φ)
@test !has_subpart(dds, :nonsubpart)
@test_throws ArgumentError subpart(dds, 1, :nonsubpart)
@test_throws ArgumentError incident(dds, 1, :nonsuppart)
@test_throws ArgumentError set_subpart!(dds, 1, :nonsubpart, 1)

# Deletion.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[2,3,3])
rem_part!(dds, :X, 2)
@test nparts(dds, :X) == 2
@test subpart(dds, :Φ) == [0,2]
@test incident(dds, 1, :Φ) == []
@test incident(dds, 2, :Φ) == [2]
rem_part!(dds, :X, 2)
@test nparts(dds, :X) == 1
@test subpart(dds, :Φ) == [0]
rem_part!(dds, :X, 1)
@test nparts(dds, :X) == 0

dds = DDS()
add_parts!(dds, :X, 4, Φ=[2,3,3,4])
@test_throws ErrorException rem_parts!(dds, :X, [4,1])
rem_parts!(dds, :X, [1,4])
@test subpart(dds, :Φ) == [1,1]
@test incident(dds, 1, :Φ) == [1,2]
@test incident(dds, 2, :Φ) == []

# Pretty printing.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[2,3,3])
s = sprint(show, dds)
@test startswith(s, "CSet")
@test occursin("X = 1:3", s)
@test occursin("Φ : X → X = ", s)

s = sprint(show, MIME"text/plain"(), dds)
@test startswith(s, "CSet")
@test occursin("X = 1:3", s)
@test occursin("│ X │", s)

s = sprint(show, MIME"text/html"(), dds)
@test startswith(s, "
")
@test occursin("
", s)
@test endswith(rstrip(s), "
")

# Special case of pretty print: empty table.
empty_dds = DDS()
@test !isempty(sprint(show, empty_dds))
@test !isempty(sprint(show, MIME"text/plain"(), empty_dds))
@test !isempty(sprint(show, MIME"text/html"(), empty_dds))

# Error handling when adding parts.
dds = DDS()
add_parts!(dds, :X, 3, Φ=[1,1,1])
@test_throws AssertionError add_part!(dds, :X, Φ=5)
@test nparts(dds, :X) == 3
@test subpart(dds, :Φ) == [1,1,1]
@test_throws AssertionError add_parts!(dds, :X, 2, Φ=[3,6])
@test nparts(dds, :X) == 3
@test incident(dds, 3, :Φ) == []

# Incidence without indexing.
UnindexedDDS = CSetType(TheoryDDS)
dds = UnindexedDDS()
add_parts!(dds, :X, 4, Φ=[3,3,4,4])
@test isempty(keys(dds.indices))
@test incident(dds, 3, :Φ) == [1,2]
@test incident(dds, 4, :Φ) == [3,4]

# Dendrograms
#############

# Strictly speaking, this data structure is not a single dendrogram but a forest
# of dendrograms (whose roots are the elements fixed under the `parent` map) and
# in order to be valid dendrograms, there must be no nontrivial cycles and the
# `height` map must satisfy `compose(parent, height) ≥ height` pointwise.

@present TheoryDendrogram(FreeSchema) begin
  X::Ob
  R::Data
  parent::Hom(X,X)
  height::Attr(X,R)
end

const AbstractDendrogram = AbstractACSetType(TheoryDendrogram)
const Dendrogram = ACSetType(TheoryDendrogram, index=[:parent])
@test Dendrogram <: AbstractDendrogram
@test Dendrogram <: ACSet
@test Dendrogram{Real} <: AbstractDendrogram{Real}
@test_throws ErrorException CSetType(TheoryDendrogram)

d = Dendrogram{Int}()
add_parts!(d, :X, 3, height=0)
add_parts!(d, :X, 2, height=[10,20])
set_subpart!(d, 1:3, :parent, 4)
set_subpart!(d, [4,5], :parent, 5)

@test nparts(d, :X) == 5
@test subpart(d, [1,2,3], :parent) == [4,4,4]
@test subpart(d, 4, :parent) == 5
@test subpart(d, :, :parent) == [4,4,4,5,5]
@test incident(d, 4, :parent) == [1,2,3]
@test has_subpart(d, :height)
@test subpart(d, [1,2,3], :height) == [0,0,0]
@test subpart(d, 4, :height) == 10
@test subpart(d, :, :height) == [0,0,0,10,20]

# Chained accessors.
@test subpart(d, 3, [:parent, :parent]) == 5
@test subpart(d, 3, [:parent, :height]) == 10
@test incident(d, 5, [:parent, :parent]) == [1,2,3,4,5]
@test incident(d, 10, [:parent, :height]) == [1,2,3]

X, parent, height = TheoryDendrogram[[:X, :parent, :height]]
@test subpart(d, 3, parent) == 4
@test subpart(d, 3, compose(parent, height)) == 10
@test subpart(d, 3, id(X)) == 3
@test incident(d, 10, compose(parent, height)) == [1,2,3]

# Indexing syntax.
@test d[3, :parent] == 4
@test d[3, [:parent, :height]] == 10
@test d[:, :parent] == [4,4,4,5,5]

# Copying parts.
d2 = Dendrogram{Int}()
copy_parts!(d2, d, X=[4,5])
@test nparts(d2, :X) == 2
@test subpart(d2, [1,2], :parent) == [2,2]
@test subpart(d2, [1,2], :height) == [10,20]

du = disjoint_union(d, d2)
@test nparts(du, :X) == 7
@test subpart(du, :parent) == [4,4,4,5,5,7,7]
@test subpart(du, :height) == [0,0,0,10,20,10,20]

# Pretty printing of data attributes.
s = sprint(show, d)
@test startswith(s, "ACSet")
@test occursin("R = Int64", s)
@test occursin("height : X → R = ", s)

s = sprint(show, MIME"text/plain"(), d)
@test startswith(s, "ACSet")

# Allow type inheritance for data attributes.
d_abs = Dendrogram{Number}()
add_parts!(d_abs, :X, 2, height=[10.0, 4])
@test subpart(d_abs, :height) == [10.0, 4]

# Dendrograms with leaves
#------------------------

@present TheoryLDendrogram <: TheoryDendrogram begin
  L::Ob
  leafparent::Hom(L,X)
end

const LDendrogram = ACSetType(TheoryLDendrogram, index=[:parent, :leafparent])

# Copying between C-sets and C′-sets with C != C′.
ld = LDendrogram{Int}()
copy_parts!(ld, d)
@test nparts(ld, :L) == 0
@test subpart(ld, :parent) == subpart(d, :parent)

add_parts!(ld, :L, 3, leafparent=[2,3,4])
@test subpart(ld, :leafparent) == [2,3,4]
d′ = Dendrogram{Int}()
copy_parts!(d′, ld)
@test d′ == d

# Labeled sets
##############

# The simplest example of a C-set with a data attribute, to test data indexing.

@present TheoryLabeledSet(FreeSchema) begin
  X::Ob
  Label::Data
  label::Attr(X,Label)
end

# Labeled sets with index
#------------------------

const IndexedLabeledSet = ACSetType(TheoryLabeledSet, index=[:label])

lset = IndexedLabeledSet{Symbol}()
@test keys(lset.indices) == (:label,)
add_parts!(lset, :X, 2, label=[:foo, :bar])
@test subpart(lset, :, :label) == [:foo, :bar]
@test incident(lset, :foo, :label) == [1]
@test isempty(incident(lset, :nonkey, :label))

add_part!(lset, :X, label=:foo)
@test incident(lset, :foo, :label) == [1,3]
set_subpart!(lset, 1, :label, :baz)
@test subpart(lset, 1, :label) == :baz
@test incident(lset, [:foo,:baz], :label) == [[3],[1]]
set_subpart!(lset, 3, :label, :biz)
@test incident(lset, :foo, :label) == []

# Deletion with indexed data attribute.
lset = IndexedLabeledSet{Symbol}()
add_parts!(lset, :X, 3, label=[:foo, :foo, :bar])
rem_part!(lset, :X, 1)
@test subpart(lset, :label) == [:bar, :foo]
@test incident(lset, [:foo, :bar], :label) == [[2], [1]]

# Deletion with unitialized data attribute.
lset = IndexedLabeledSet{Symbol}()
add_part!(lset, :X)
rem_part!(lset, :X, 1)
@test nparts(lset, :X) == 0

# Pretty-printing with unitialized data attribute.
lset = IndexedLabeledSet{Symbol}()
add_part!(lset, :X)
@test occursin("#undef", sprint(show, lset))
@test occursin("#undef", sprint(show, MIME"text/plain"(), lset))
@test occursin("#undef", sprint(show, MIME"text/html"(), lset))

# Labeled sets with unique index
#-------------------------------

const UniqueIndexedLabeledSet = ACSetType(TheoryLabeledSet,
                                          unique_index=[:label])

lset = UniqueIndexedLabeledSet{Symbol}()
@test keys(lset.indices) == (:label,)
add_parts!(lset, :X, 2, label=[:foo, :bar])
@test subpart(lset, :, :label) == [:foo, :bar]
@test incident(lset, :foo, :label) == 1
@test incident(lset, [:foo,:bar], :label) == [1,2]
@test incident(lset, :nonkey, :label) == 0

set_subpart!(lset, 1, :label, :baz)
@test subpart(lset, 1, :label) == :baz
@test incident(lset, :baz, :label) == 1
@test incident(lset, :foo, :label) == 0

@test_throws ErrorException set_subpart!(lset, 1, :label, :bar)

end


text:
"""
    L5_code_length

The number of bits in the L5 signal for a 1ms segment = `10230`
"""
const L5_code_length = 10230  # chips

"""
    nh20_code_length

The number of bits in the Neuman sequence in a 1s segment = `20`
"""
const nh20_code_length = 20  # chips

"""
    nh10_code_length

The number of bits in the Neuman sequence in a 1s segment = `10`
"""
const nh10_code_length = 10  # chips

"""
    L5_chipping_rate

The chipping rate of the L5 signal = `10.23e6Hz`
"""
const L5_chipping_rate = 10.23e6  # Hz

"""
    nh_chipping_rate

The chipping rate of the Neuman sequence = `1000Hz`
"""
const nh_chipping_rate = 1000.  # Hz

"""
   L5_db_chipping_rate

The chipping rate of the databits for L5I = `100Hz`
"""
const L5_db_chipping_rate = 100.  # Hz

"""
    XA_XB_code_length

The length of the XA and XB codes used to construct the L5 signal = `8190`
"""
const XA_XB_code_length = 8190  # chips

"""
    XA_taps

The taps for the XA shift register.

Taps used = `[9, 10, 12, 13]`
"""
const XA_taps = [9, 10, 12, 13]

"""
    XB_taps

The taps for the XB shift register.

Taps used = `[1, 3, 4, 6, 7, 8, 12, 13]`
"""
const XB_taps = [1, 3, 4, 6, 7, 8, 12, 13]

"""
    XB_Q_intial_conditions

The initial conditions of the XB shift register per PRN for Q5 only. Used with
the dataless signal.
"""
const XB_Q_intial_conditions = Dict( 1 => [1, 0, 0, 1, 0, 1, 1, 0, 0, 1, 1, 0, 0],
                                     2 => [0, 1, 0, 0, 0, 1, 1, 1, 1, 0, 1, 1, 0],
                                     3 => [1, 1, 1, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1],
                                     4 => [0, 0, 1, 1, 1, 0, 1, 1, 0, 1, 0, 1, 0],
                                     5 => [0, 0, 1, 1, 1, 1, 0, 1, 1, 0, 0, 1, 0],
                                     6 => [0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0, 1],
                                     7 => [1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1],
                                     8 => [0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 0],
                                     9 => [1, 0, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 1],
                                    10 => [0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 1, 1, 0],
                                    11 => [0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 1],
                                    12 => [0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0, 1],
                                    13 => [0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1],
                                    14 => [1, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
                                    15 => [1, 0, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 1],
                                    16 => [1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1],
                                    17 => [1, 1, 1, 0, 0, 1, 1, 0, 0, 1, 0, 0, 0],
                                    18 => [1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0],
                                    19 => [0, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0, 1, 1],
                                    20 => [1, 1, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 1],
                                    21 => [0, 1, 1, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0],
                                    22 => [0, 0, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 0],
                                    23 => [1, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1],
                                    24 => [0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 1],
                                    25 => [0, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1],
                                    26 => [0, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 0, 0],
                                    27 => [1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0],
                                    28 => [1, 1, 1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0],
                                    29 => [0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 0],
                                    30 => [1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 1],
                                    31 => [0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 1, 0, 1],
                                    32 => [1, 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 1, 0])

"""
    XB_I_intial_conditions

The initial conditions of the XB shift register per PRN for I5 only. Used with
the data signal.
"""
const XB_I_intial_conditions = Dict( 1 => [0, 1, 0, 1, 0, 1, 1, 1, 0, 0, 1, 0, 0],
                                     2 => [1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1],
                                     3 => [0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
                                     4 => [1, 0, 1, 1, 0, 0, 0, 1, 0, 0, 1, 1, 0],
                                     5 => [1, 1, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, 1],
                                     6 => [0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 1, 0],
                                     7 => [1, 0, 1, 0, 0, 1, 0, 0, 1, 1, 1, 1, 1],
                                     8 => [1, 0, 1, 1, 1, 1, 0, 1, 0, 0, 1, 0, 0],
                                     9 => [1, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1],
                                    10 => [0, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 0],
                                    11 => [0, 0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 1, 0],
                                    12 => [1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0, 0, 1],
                                    13 => [0, 0, 0, 1, 1, 1, 0, 0, 1, 1, 1, 0, 0],
                                    14 => [0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 1],
                                    15 => [0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0, 1, 0],
                                    16 => [0, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1],
                                    17 => [0, 1, 0, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
                                    18 => [1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 1, 1, 0],
                                    19 => [1, 1, 0, 0, 1, 0, 0, 0, 1, 1, 1, 1, 1],
                                    20 => [0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 1, 0, 1],
                                    21 => [0, 0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
                                    22 => [1, 1, 1, 0, 1, 1, 1, 1, 0, 1, 1, 1, 1],
                                    23 => [1, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 0],
                                    24 => [1, 1, 0, 0, 0, 1, 0, 1, 1, 0, 1, 0, 0],
                                    25 => [1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1],
                                    26 => [1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 0],
                                    27 => [0, 1, 0, 1, 0, 1, 1, 0, 1, 1, 1, 1, 0],
                                    28 => [0, 1, 1, 1, 1, 0, 1, 0, 1, 0, 1, 1, 0],
                                    29 => [0, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1],
                                    30 => [1, 0, 0, 0, 0, 1, 0, 1, 1, 0, 1, 1, 1],
                                    31 => [0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 1, 1, 0],
                                    32 => [0, 0, 0, 0, 0, 1, 0, 1, 1, 1, 0, 0, 1])

"""
    nh20

A 20-bit Neuman-Hofman code XORed to the L5Q code. `nh20` is 20ms in length,
where every bit is XORed to each 1ms L5Q code sequence.
"""
const nh20 = [0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 0, 0, 1, 1, 1, 0]

"""
    nh10

A 10-bit Neuman-Hofman code XORed to the L5I code. `nh10` is 10ms in length,
where every bit is XORed to each 1ms L5I code sequence.
"""
const nh10 = [0, 0, 0, 0, 1, 1, 0, 1, 0, 1]

"""
    xaregister()

Generates the XA register 8091 bit sequence. Restarts after the first 8090 bits,
where it repeats a second time until the 2040th state where it is reset again
for the start of the next 1ms.

Arguments: `None`

Returns:

- `xa_code::Vector`: vector containg the XA register result for 1ms period of
                     time
"""
function xaregister()
    xa_register = ones(Int, 13)
    xa_code = zeros(Int, L5_code_length)
    @inbounds for i in 1:XA_XB_code_length
        xa_code[i] = xa_register[end]
        xa_register[end] = xor(xa_register[XA_taps]...)
        xa_register = circshift(xa_register, 1)
    end
    xa_code[XA_XB_code_length+1:end] = xa_code[1:2040]
    return xa_code
end

"""
    xbregister(prn::Int, intial_conditions)

Generates the XB code sequence for 1ms for the specified PRN. Register starts
with initial conditions defined by `intial_conditions[prn]`.

Required Arguments:

- `prn::Int`: the pseudorandom noise (PRN) number, which is also related to the
              satellite vehicle
- `initial_conditions::Dict`: dictionary containing the initial states of the
                              register for each PRN

Returns:

- `xb_code::Vector`: vector containing the XB code sequence for a given PRN
                     and initial conditions
"""
function xbregister(prn::Int, intial_conditions)
    xb_register = deepcopy(intial_conditions[prn])
    xb_code = zeros(eltype(intial_conditions[prn]),
    	            L5_code_length)
    @inbounds for i in 1:L5_code_length
        xb_code[i] = xb_register[end]
        xb_register[end] = xor(xb_register[XB_taps]...)
        xb_register = circshift(xb_register, 1)
    end
    return xb_code
end

"""
    xa_code

Constant output of the XA register for the L5Q signal.
"""
const xa_code = xaregister()

"""
    gen_xb_codes()

Generates the output of the XB registers for all PRNs in `XB_intial_conditions.`

Arguments: `None`

Returns:

- `xb_codes::Dict`: dictionary containing the output of the XB register for all
                    PRN numbers
"""
    function gen_xb_codes()
	xb_codes = typeof(XB_Q_intial_conditions)()
	for prn in keys(XB_Q_intial_conditions)
		xb_codes[prn] = xbregister(prn, XB_Q_intial_conditions)
	end
	return xb_codes
end

# """
#     xb_codes

# Dictionary containing the output of the XB register for
# all PRNs.
# """
# const xb_codes = gen_xb_codes(XB_intial_conditions)

"""
    gen_l5q_codes()

Generates the L5Q codes for each PRN listed in `XB_Q_intial_conditions.`

Arguments: `None`

Returns:

- `l5q_codes::Dict`: dictionary containing L5Q codes for all PRN numbers
"""
function gen_l5q_codes()
	l5q_codes = typeof(XB_Q_intial_conditions)()
	for prn in keys(XB_Q_intial_conditions)
		l5q_codes[prn] = xor.(xa_code, xbregister(prn, XB_Q_intial_conditions))
	end
	return l5q_codes
end

"""
    gen_l5i_codes()

Generates the L5I codes for each PRN listed in `XB_I_intial_conditions.`

Arguments: `None`

Returns:

- `l5i_codes::Dict`: dictionary containing L5I codes for all PRN numbers
"""
function gen_l5i_codes()
    l5i_codes = typeof(XB_I_intial_conditions)()
    for prn in keys(XB_I_intial_conditions)
        l5i_codes[prn] = xor.(xa_code, xbregister(prn, XB_I_intial_conditions))
    end
    return l5i_codes
end

"""
    l5q_codes

Dictionary containing the L5Q codes for all PRNs.
"""
const l5q_codes = gen_l5q_codes()

"""
    l5i_codes

Dictionary containing the L5I codes for all PRNs.
"""
const l5i_codes = gen_l5i_codes()

"""
	define_l5_code_type(t_length=missing; channel="both",
						prns=collect(keys(l5i_codes)), sig_freq=L5_freq)

Generates a `SignalType` struct which defines the L5 codes that can be used
to define and generate L5 signals. The databits generated are random.

Optional Arguments:

- `t_length:Float64`: the length of the planned signal in seconds
    * this argument is used only to generate a `Vector` of databits so they do
      not repeat throughout the duration of the signal, once it is generated
    * `(default = missing)`
	* if given, databits are included in signal type
- `channel`: channel to include (either `"I"`, `"Q"`, or `"both"`)
	* `(default = "both")`
- `prns`: codes to define for signal type `(default = keys(l5i_codes))`
	* default is to define all PRN codes and assume that databits for each
	  are different
- `sig_freq`: carrier frequency in Hz `(default = L5_freq)`
- `databits`: optional databits `(default = missing)`, where random bits are used
- `B`: receiver bandwidth `(default = missing)`

Returns:

- `siganl_type::SignalType`: can be used in `definesignal` function to create
                             a signal, which can be used with `generatesignal!`
                             to generate the signal
"""
function define_l5_code_type(t_length=missing; channel="both",
	                         prns=collect(keys(l5i_codes)), sig_freq=L5_freq,
							 databits=missing, B=missing)
	if (channel == "I") || (channel == "both")
        if ~ismissing(t_length)
            # Generate random databits
			if ismissing(databits)
				databits = random_databits(L5_db_chipping_rate, t_length;
				                           prns=prns)
	        end
            # Include databits in I channel codes
            I_codes = definecodetype([copy_dictionary(l5i_codes, prns), nh10],
			                         [L5_chipping_rate, nh_chipping_rate];
                                     databits=[databits, L5_db_chipping_rate])
	    else
	        # Define I channel codes without databits
	        I_codes = definecodetype([copy_dictionary(l5i_codes, prns), nh10],
			                         [L5_chipping_rate, nh_chipping_rate])
	    end
		if channel == "I"
			signal_type = definesignaltype(I_codes, sig_freq, "I"; name="L5",
                                           B=B)
		else  # means that `channel` == "both"
			Q_codes = definecodetype([copy_dictionary(l5q_codes, prns), nh20],
									 [L5_chipping_rate, nh_chipping_rate])
			signal_type = definesignaltype(I_codes, Q_codes, sig_freq; name="L5",
                                           B=B)
		end
	elseif channel == "Q"
		Q_codes = definecodetype([copy_dictionary(l5q_codes, prns), nh20],
								 [L5_chipping_rate, nh_chipping_rate])
		signal_type = definesignaltype(Q_codes, sig_freq, "Q"; name="L5", B=B)
	else
		error("Invalid channel specifed.")
	end
    return signal_type
end


text:
include("../../matplotlib_settings.jl")

set_config(:tab_compensated_summation, false)

using ChargedParticleDynamics.GuidingCenter4d: plot_integral_error, plot_poincare_surface

const coll_title = "Guiding Center 4D Symmetric 2nd Poincaré Integral Invariant"

const XMIN = -0.25
const XMAX = +0.25
const YMIN = -0.25
const YMAX = +0.25

include("../guiding_center_4d_poincare_invariant_2nd.jl")


text:
using Stage

# test the macros
@stage function add(a, b)
  sleep(10)
  a + b
end
@stage function mult(a, b)
  sleep(2)
  a * b
end
@stage function sub(a, b)
  sleep(2)
  a - b
end

@debug "initiating test"
res2 = sub("sub2", sub("sub", mult("mult", 4, 3), 2), 2)
res_future = add("add", mult("mult2", 1, 2), 3)
@debug "mult result: $(fetch(res2))"
@debug "madd result: $(fetch(res_future))"

res3 = sub("s1", mult("m1", 4, 3), mult("m2", 5, 6))
@expect fetch(res3) == -18
@debug "res3 = $(fetch(res3)) [should be -18 and should have taken approx. 4 seconds]"

res4 = sub("sub2", mult("m1", 4, 3), mult("m2", 5, 6))
@expect fetch(res4) == 8
@debug "res4 = $(fetch(res4)) [should be 8 and should have taken approx. 0 seconds]"
@debug "completed test"


text:
#####
##### helper functions
#####

"function wavespeed_1D(rho,rhou,E)
    one-dimensional wavespeed (for DG penalization terms)"
function wavespeed_1D(rho,rhou,E)
    p = pfun_nd(rho,(rhou,),E)
    cvel = @. sqrt(γ*p/rho)
    return @. abs(rhou/rho) + cvel
end

#####
##### one-dimensional fluxes
#####

"function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                        rhoR,uR,betaR,rhologR,betalogR)"
function euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                        rhoR,uR,betaR,rhologR,betalogR)

    rholog = logmean.(rhoL,rhoR,rhologL,rhologR)
    betalog = logmean.(betaL,betaR,betalogL,betalogR)

    # arithmetic avgs
    rhoavg = (@. .5*(rhoL+rhoR))
    uavg   = (@. .5*(uL+uR))

    unorm = (@. uL*uR)
    pa    = (@. rhoavg/(betaL+betaR))
    f4aux = (@. rholog/(2*(γ-1)*betalog) + pa + .5*rholog*unorm)

    FxS1 = (@. rholog*uavg)
    FxS2 = (@. FxS1*uavg + pa)
    FxS3 = (@. f4aux*uavg)

    return (FxS1,FxS2,FxS3)
end

#####
##### shim functions for 1D
#####

# dispatch to n-dimensional constitutive routines, with optional entropy scaling
function primitive_to_conservative(rho,u,p)
   rho,rhou,E = primitive_to_conservative_nd(rho,tuple(u),p)
   return rho,first(rhou),E
end
function v_ufun(rho,rhou,E)
    v1,vU,vE = v_ufun_nd(rho,tuple(rhou),E)
    return scale_entropy_output(v1,first(vU),vE)
end
function u_vfun(v1,vU,vE)
    v1,vU,vE = scale_entropy_input(v1,vU,vE)
    rho,rhoU,E = u_vfun_nd(v1,tuple(vU),vE)
    return rho,first(rhoU),E
end
function conservative_to_primitive_beta(rho,rhou,E)
    rho,U,beta = conservative_to_primitive_beta_nd(rho,tuple(rhou),E)
    return rho,first(U),beta
end

Sfun(rho,rhou,E) = Sfun_nd(rho,tuple(rhou),E)
pfun(rho,rhou,E) = pfun_nd(rho,tuple(rhou),E)
betafun(rho,rhou,E) = betafun_nd(rho,tuple(rhou),E)

"function euler_fluxes_1D(rhoL,uL,betaL,rhoR,uR,betaR)"
function euler_fluxes(rhoL,uL,betaL,rhoR,uR,betaR)
    rhologL,betalogL,rhologR,betalogR = map(x->log.(x),(rhoL,betaL,rhoR,betaR))
    return euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,
                           rhoR,uR,betaR,rhologR,betalogR)
end
euler_fluxes(rhoL,uL,betaL,rhologL,betalogL,rhoR,uR,betaR,rhologR,betalogR) =
    euler_fluxes_1D(rhoL,uL,betaL,rhologL,betalogL,rhoR,uR,betaR,rhologR,betalogR)

    euler_fluxes(varL,varR) = euler_fluxes(varL...,varR...)
    euler_fluxes(varL,varR,logL,logR) = euler_fluxes(varL...,logL...,varR...,logR...)

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
"""
```
gr(r::Vector{Float64}, count::Vector{Float64}, density::Float64, binstep::Float64)
```

Computes the radial distribution function from the count data and
the density.

This is exactly a conventional g(r)
if a single atom was chosen as the solute and solvent selections.

Returns both the g(r) and the kb(r)

"""
function gr(r::Vector{Float64}, count::Vector{Float64}, density::Float64, binstep::Float64)
    nbins = length(r)
    gr = zeros(nbins)
    kb = zeros(nbins)
    for i = 1:nbins
        gr[i] = (count[i] / sphericalshellvolume(i, binstep)) / density
        if i == 1
            kb[i] = 4π * (gr[i] - 1) * r[i]^2 * binstep
        else
            kb[i] = kb[i-1] + 4π * (gr[i] - 1) * r[i]^2 * binstep
        end
    end
    @. kb = units.Angs3tocm3permol * kb
    return gr, kb
end

"""

```
gr(R::Result) = gr(R.d,R.rdf_count,R.density.solvent_bulk,R.options.binstep)
```

If a Result structure is provided without further details, use the rdf count
and the bulk solvent density.
"""
gr(R::Result) = gr(R.d, R.rdf_count, R.density.solvent_bulk, R.options.binstep)


text:
struct DiscreteDistribution
    p::Vector{Float64}
    P::Vector{Float64} # p(k<=i)
end

function DiscreteDistribution(w)
    @assert all(w .>= 0) "All weights should be non-negative"
    p = w / sum(w)
    P = cumsum(p)
    return DiscreteDistribution(p, P)
end

rand(d::DiscreteDistribution) = searchsortedfirst(d.P, rand())

"""
    randint([rng,] probs::Vector{Float64})

Generate a random integer `i` between `1` and `length(probs)`.
`i` is sampled with probability probs[i].
The vector `probs` is assumed to sum to one.
"""
function randint(rng::AbstractRNG, probs::Vector{Float64})
    r = rand(rng)
    c = 0.
    i = 0
    for p in probs
        i += 1
        c += p
        r <= c && return i
    end
    return 0 #not supposed to arrive here
end

randint(probs::Vector{Float64}) = randint(GLOBAL_RNG, probs)

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
# -----------------------------------------------------------------------------
##### CYLINDER/SHELL #####-----------------------------------------------------

function shellregioncylinder(dataobject::HydroDataType;
                            radius::Array{<:Real,1}=[0.,0.],
                            height::Real=0.,
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            direction::Symbol=:z,
                            cell::Bool=true,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    radius_in  = radius[1]
    radius_out = radius[2]
    if radius_in == 0. || radius_out == 0. || height == 0. || in(0., center)
        error("[Mera]: given radius, height or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift, height_shift = prep_cylindrical_shellranges(dataobject.info, center, radius_in, radius_out, height, range_unit, verbose)

    if inverse == false
        if isamr
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                >= (cell_shift(p.level, radius_in_shift,cell))  &&

                                  get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                <= (cell_shift(p.level, radius_out_shift,cell))  &&

                                get_height_cylinder(p.cz, p.level, cz_shift)
                                <= (cell_shift(p.level, height_shift,cell)),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                >= (cell_shift(lmax, radius_in_shift,cell))  &&

                                  get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                <= (cell_shift(lmax, radius_out_shift,cell))  &&

                                get_height_cylinder(p.cz, lmax, cz_shift)
                                <= (cell_shift(lmax, height_shift,cell)),
                                dataobject.data)
        end

    elseif inverse == true
        if isamr
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                < (cell_shift(p.level, radius_in_shift,cell))  ||

                                  get_radius_cylinder(p.cx, p.cy, p.level, cx_shift, cy_shift)
                                > (cell_shift(p.level, radius_out_shift,cell))  ||

                                get_height_cylinder(p.cz, p.level, cz_shift)
                                > (cell_shift(p.level, height_shift,cell)),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                < (cell_shift(lmax, radius_in_shift,cell))  ||

                                  get_radius_cylinder(p.cx, p.cy, lmax, cx_shift, cy_shift)
                                > (cell_shift(lmax, radius_out_shift,cell))  ||

                                get_height_cylinder(p.cz, lmax, cz_shift)
                                > (cell_shift(lmax, height_shift,cell)),
                                dataobject.data)
        end

        ranges = dataobject.ranges

    end

        printtablememory(sub_data, verbose)

        hydrodata = HydroDataType()
        hydrodata.data = sub_data
        hydrodata.info = dataobject.info
        hydrodata.lmin = dataobject.lmin
        hydrodata.lmax = dataobject.lmax
        hydrodata.boxlen = dataobject.boxlen
        hydrodata.ranges = ranges
        hydrodata.selected_hydrovars = dataobject.selected_hydrovars
        hydrodata.used_descriptors = dataobject.used_descriptors
        hydrodata.smallr = dataobject.smallr
        hydrodata.smallc = dataobject.smallc
        hydrodata.scale = dataobject.scale
        return hydrodata

end

# -----------------------------------------------------------------------------
##### SPHERE/SHELL #####-------------------------------------------------------

function shellregionsphere(dataobject::HydroDataType;
                            radius::Array{<:Real,1}=[0.,0.],
                            center::Array{<:Any,1}=[0.,0.,0.],
                            range_unit::Symbol=:standard,
                            cell::Bool=true,
                            inverse::Bool=false,
                            verbose::Bool=verbose_mode)

    printtime("", verbose)

    radius_in  = radius[1]
    radius_out = radius[2]
    if radius_in == 0. || radius_out == 0. || in(0., center)
        error("[Mera]: given inner and outer radius or center should be != 0.")
    end

    boxlen = dataobject.boxlen
    scale = dataobject.scale
    lmax = dataobject.lmax
    isamr = checkuniformgrid(dataobject, lmax)

    # convert given ranges and print overview on screen
    ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift = prep_spherical_shellranges(dataobject.info, center, radius_in, radius_out, range_unit, verbose)

    if inverse == false
        if isamr
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)
                                >= cell_shift(p.level, radius_in_shift, cell) &&

                                get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)
                                <= cell_shift(p.level, radius_out_shift,cell),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)
                                >= cell_shift(lmax, radius_in_shift, cell) &&

                                get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)
                                <= cell_shift(lmax, radius_out_shift,cell),
                                dataobject.data)
        end

    elseif inverse == true
        if isamr
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)
                                < cell_shift(p.level, radius_in_shift, cell) ||

                                get_radius_sphere(p.cx, p.cy, p.cz, p.level, cx_shift, cy_shift, cz_shift)
                                > cell_shift(p.level, radius_out_shift, cell),
                                dataobject.data)
        else # for uniform grid
            sub_data = filter(p-> get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)
                                < cell_shift(lmax, radius_in_shift, cell) ||

                                get_radius_sphere(p.cx, p.cy, p.cz, lmax, cx_shift, cy_shift, cz_shift)
                                > cell_shift(lmax, radius_out_shift, cell),
                                dataobject.data)
        end
        ranges = dataobject.ranges

    end

    printtablememory(sub_data, verbose)

    hydrodata = HydroDataType()
    hydrodata.data = sub_data
    hydrodata.info = dataobject.info
    hydrodata.lmin = dataobject.lmin
    hydrodata.lmax = dataobject.lmax
    hydrodata.boxlen = dataobject.boxlen
    hydrodata.ranges = ranges
    hydrodata.selected_hydrovars = dataobject.selected_hydrovars
    hydrodata.used_descriptors = dataobject.used_descriptors
    hydrodata.smallr = dataobject.smallr
    hydrodata.smallc = dataobject.smallc
    hydrodata.scale = dataobject.scale
    return hydrodata

end


text:
include("main.jl")
include("bufferinfo.jl")
include("utility.jl")


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ f7e8e8ce-7dfc-11eb-0677-cb442926e519
using Distributions

# ╔═╡ e9b0055a-7dfc-11eb-2270-e30e833cee00
using StatsPlots

# ╔═╡ ef30e082-7dfd-11eb-3116-b590d75e8735
begin
	k = @bind K html"<input type='range' min='1.0' max='20.0' step='1.0' value='10'>"
	P = @bind p html"<input type='range' min='0.0' max='1.0' step='0.01' value='0.5'>"
	md"""**distribution parameters:**
	$$K$$: $k
	$$p$$: $P
	"""
end

# ╔═╡ 7f8de914-7dfd-11eb-100b-a33c8cfbd43d
x = [0:1:K;]

# ╔═╡ b887bd70-7e03-11eb-2ed9-4127c977a6d7
md"""
 $$K$$ = $K

 $$p$$ = $p
"""

# ╔═╡ c7a6a7aa-7dfe-11eb-1110-f313180ab599
begin
	m = K * p
	v = m * (1 - p)
	md"""
	Mean: $$K * p$$ = $m

	Variance: $$K * p * (1-p)$$ = $v
	"""
end

# ╔═╡ 821bd6d6-7dfe-11eb-1e6a-c5e52652461b
B = Binomial(K, p)

# ╔═╡ 7d111506-7dff-11eb-1b4d-f5b37e333a55
plot(x, pdf.(B, x),
	t = "scatter", label = "Binomial($K, $p)",
	ylims = (0, maximum(pdf.(B, x)) + 0.1))

# ╔═╡ Cell order:
# ╠═7f8de914-7dfd-11eb-100b-a33c8cfbd43d
# ╟─ef30e082-7dfd-11eb-3116-b590d75e8735
# ╟─b887bd70-7e03-11eb-2ed9-4127c977a6d7
# ╟─c7a6a7aa-7dfe-11eb-1110-f313180ab599
# ╠═821bd6d6-7dfe-11eb-1e6a-c5e52652461b
# ╟─7d111506-7dff-11eb-1b4d-f5b37e333a55
# ╟─f7e8e8ce-7dfc-11eb-0677-cb442926e519
# ╟─e9b0055a-7dfc-11eb-2270-e30e833cee00


text:
"""
    centre_of_gravity(img)

Given a matrix (i.e. an image) return the centre of gravity.
This is an index position (not rounded).
"""
function centre_of_gravity(img)
    tot = sum(img)
    cog_x = zero(eltype(img))
    for i in axes(img,1)
        cog_xi = zero(eltype(img))
        for j in axes(img,2)
            cog_xi += img[i,j]
        end
        cog_x += cog_xi * i
    end
    cog_x/=tot
    cog_y = zero(eltype(img))
    for i in axes(img,2)
        cog_yi = zero(eltype(img))
        for j in axes(img,1)
            cog_yi += img[j,i]
        end
        cog_y += cog_yi * i
    end
    cog_y/=tot
    return (cog_x, cog_y)
end
export centre_of_gravity

"""
    profile2image(profile, image)

Given an image and a vector profile (could be contrast, any
1D function), return a new image with the axes of `image`
and the values of `profile` interpolated radially from the
centre.
"""
function profile2image(profile, image)
    img_out = similar(image)

    I = origin(image)
    if I[1] == I[2] == 1
        image = centered(image)
    end
    dx = axes(image,1)
    dy = axes(image,2)
    dr = sqrt.(
        dx.^2 .+ (dy').^2
    )

    img_out .= profile.(dr)
    return img_out
end
export profile2image

"""
    lookup_coord(img, [125.24, -129.1], 12.43)

Given an image, a coordinate vector in milliarcseconds, and the platescale
in milliarcseconds, return the pixel at that coordinate using bilinear interpolation.
"""
## Function for looking up a single pixel in an image, with bi-linear interpolation
function lookup_coord(image::AbstractArray, x,y, platescale)
    return lookup_coord(image, (x,y), platescale)
end

function lookup_coord(image::AbstractArray, vec, platescale)
    ix = vec[1] / platescale
    iy = vec[2] / platescale
    axx = axes(image, 1)
    axy = axes(image, 2)

    # floor & ceil can actually fail with an InexactError if ix or iy fall outside
    # the range representable by Int64. In these absurd cases, we can safely
    # return 0 as the position is necessarily outside the range of indices possible
    # to store in an array on this system.
    if !(typemin(Int) < ix < typemax(Int)) || !(typemin(Int) < iy < typemax(Int))
        return 0
    end

    ix1 = floor(Int, ix)
    ix2 = ceil(Int, ix)
    iy1 = floor(Int, iy)
    iy2 = ceil(Int, iy)
    ix2 += (ix1 == ix2)
    iy2 += (iy1 == iy2)

    # Bi-linear interpolation, if in-bounds
    tot = zero(eltype(image))
    if (ix1 ∈ axx) & (ix2 ∈ axx) & (iy1 ∈ axy) & (iy2 ∈ axy)
        top_ave =
            image[ix1, iy2] * (ix2 - ix) +
            image[ix2, iy2] * (ix - ix1)
        bot_ave =
            image[ix1, iy1] * (ix2 - ix) +
            image[ix2, iy1] * (ix - ix1)
        tot += top_ave * (iy - iy1) + bot_ave * (iy2 - iy)
    end
    return tot
end
export lookup_coord

using Statistics

"""
    imgsep(img)

Given an array, return the distance of each pixel
to the centre.
"""
function imgsep(
    img,
    centrex=mean(axes(img,1)),
    centrey=mean(axes(img,2)),
)
    xs = axes(img,1) .- centrex
    ys = axes(img,2) .- centrey
    r = sqrt.(xs.^2 .+ ys'.^2)
    return r
end
export imgsep

"""
    fakeimage((100,100), [1.2])

Given desired image dimensions and gaussian blur level, create
a random image with a noise distribution a bit like a real direct
image.
"""
function fakeimage(dims,gaus_size=1.2)
    image = zeros(dims)
    cont, views = contrast!_prep(image,1:2:ceil(Int,maximum(size(image))*√2))
    for (i,v) in pairs(views)
        v .=1/ i^2*randn(length(v))
    end
    blurred = imfilter(image, Kernel.gaussian(gaus_size),NA())
    return blurred
end
export fakeimage

max_stars_repo_name: sefffal/SpaceImages


text:
abstract type OneDimQuadrature end

struct GaussLegendreQuad{T1,T2,T3} <: OneDimQuadrature
    n::T1
    nodes::T2
    weights::T3

    function GaussLegendreQuad(n; wType = Float64)
        @assert n > 0
        quad = gausslegendre(n)
        nodes = wType.(quad[1])
        weights = wType.(quad[2])
        return new{typeof(n),typeof(nodes),typeof(weights)}(n, nodes, weights)
    end
end

struct GaussHermiteQuad{T1,T2,T3} <: OneDimQuadrature
    n::T1
    nodes::T2
    weights::T3

    function GaussHermiteQuad(n; wType = Float64)
        @assert n > 0
        quad = gausshermite(n)
        nodes = wType.(quad[1])
        weights = wType.(quad[2])
        return new{typeof(n),typeof(nodes),typeof(weights)}(n, nodes, weights)
    end
end

max_stars_repo_name: d-netto/PolyChaosODE

max_issues_repo_name: d-netto/PolyChaosODE

max_forks_repo_name: d-netto/PolyChaosODE


text:
#=
Copyright 2022 INSIGNEO Institute for in silico Medicine

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

   http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.
=#

using openBF

parsed_args = openBF.parseCommandline()

input_filename = parsed_args["input_filename"]
verbose = parsed_args["verbose"]
out_files = parsed_args["out_files"]
conv_ceil = parsed_args["conv_ceil"]

openBF.runSimulation(input_filename, verbose=verbose, out_files=out_files, conv_ceil=conv_ceil)

rm("main.jl")


text:
module WASMVector

import ExportWebAssembly
export WVector, S_str, asvec
using LLVM
using LLVM.Interop

struct Ptr32{T}
    address::Int32
end

@inline malloc(len) = ExportWebAssembly.@extern(:malloc, Int32, Tuple{Int32}, Int32(len))

struct WVector{T} <: AbstractVector{T}
    ptr::Ptr32{T}
    len::Int32
    maxlen::Int32
end
WVector{T}(ptr::Ptr32{T}, len::Integer) where {T} = WVector{T}(ptr, Int32(len), Int32(len))
WVector{T}(addr::Int32, len::Integer) where {T} = WVector{T}(Ptr32{T}(addr), len)
WVector{T}(addr::Int32, len::Integer, maxlen::Integer) where {T} = WVector{T}(Ptr32{T}(addr), len, maxlen)

@inline function WVector{T}(len::Int32) where {T}
    #addr = malloc(27)
    #WVector{T}(addr, Int32(27))
    multiplier = Int32(4)
    maxlen = multiplier * len
    nbytes = Int32(sizeof(T))
    addr = malloc(nbytes * multiplier)
    WVector{T}(addr, len, len * multiplier)
end

#function WVector(args...)
#    n = length(args)
#    T = promote_type(typeof.(args)...)
#    res = WVector{T}(n)
#    for i in 1:n
#        res[i] = args[i]
#    end
#    res
#end

# function WVector(x::Vector{T}) where {T}
#     res = WVector{T}(length(x))
#     res .= x
#     res
# end

@inline function Base.getindex(x::WVector{T}, i::Integer) where {T}
    # @assert 0 < i <= x.len
    myunsafe_load(x.ptr, Int32(i))
end

@generated function myunsafe_load(p::Ptr32{T}, i::Int32) where {T}
    lt = ExportWebAssembly.llvmtype(T)
    runstr = """
        %pa = inttoptr i32 %0 to $lt*
        %pe = getelementptr $lt, $lt* %pa, i32 %1
        %ret = load $lt, $lt* %pe
        ret $lt %ret
        """
    quote
        Base.llvmcall($runstr, $T, Tuple{Int32, Int32}, p.address, i - Int32(1))
    end
end
# function unsafe_load(p::Ptr32{Float64}, i::Int32)
#     runstr = """
#         %pa = inttoptr i32 %0 to double*
#         %pe = getelementptr double, double* %pa, i32 %1
#         %ret = load double, double* %pe
#         ret double %ret
#         """
#     Base.llvmcall(runstr, Float64, Tuple{Int32, Int32}, p.address, i - Int32(1))
# end

@inline function Base.setindex!(x::WVector{T}, val, i::Integer) where {T}
    # @assert 0 < i <= x.len
    # unsafe_store!(x.ptr, convert(T, val), Int32(i))
    unsafe_store!(x.ptr, val, Int32(i))
    x
end

@generated function Base.unsafe_store!(p::Ptr32{T}, x::T, i::Int32) where {T}
    lt = ExportWebAssembly.llvmtype(T)
    runstr = """
        %pa = inttoptr i32 %0 to $lt*
        %pe = getelementptr $lt, $lt* %pa, i32 %2
        store $lt %1, $lt* %pe
        ret void
        """
    quote
        Base.llvmcall($runstr, Nothing, Tuple{Int32, $T, Int32}, p.address, x, i - Int32(1))
    end
end

Base.size(x::WVector) = (x.len,)

Base.IndexStyle(::Type{T}) where T <: WVector = IndexLinear()

Base.similar(A::WVector{T}, ::Type{T}, dims::Dims) where {T} =
    WVector{T}(dims[1])

# Array with AssemblyScript layout for use with wasm-ffi
function asvec(T, addr)
    x = WVector{Int32}(addr, Int32(4))
    WVector{T}(addr + Int32(16), x[3] ÷ Int32(sizeof(T)))
end

macro S_str(s)
    T_i32 = LLVM.IntType(32, JuliaContext())
    llvmf, _ = create_function(T_i32)
    len = length(s) + 1
    # generate IR
    Builder(JuliaContext()) do builder
        entry = BasicBlock(llvmf, "entry", JuliaContext())
        position!(builder, entry)
        s = globalstring!(builder, s)
        i = ptrtoint!(builder, s, T_i32)
        ret!(builder, i)
    end
    e = call_function(llvmf, Int32, Tuple{})
    e
    # quote
    #     WVector{Int8}($e, $len)
    # end
end

end # module

# module VectorStrings

# using ..WASMVector

# struct VectorString <: AbstractString
#     data::WVector{Int8}
# end

# VectorString(s::String) = VectorString(WVector(Vector{Int8}(s)))

# Base.ncodeunits(s::VectorString) = length(s.data)
# Base.codeunit(s::VectorString) = Int8

# Base.codeunit(s::VectorString, i::Integer) = s.data[i]

# Base.isvalid(s::VectorString, i::Int) = checkbounds(Bool, s, i) && thisind(s, i) == i

# Base.@propagate_inbounds function Base.next(s::VectorString, i::Int)
#     b = Base.codeunit(s, i)
#     u = Int32(b) << 24
#     Base.between(b, 0x80, 0xf7) || return reinterpret(Char, u), i+1
#     return Base.next_continued(s, i, u)
# end
# function Base.next_continued(s::VectorString, i::Int, u::Int32)
#     u < 0xc0000000 && (i += 1; @goto ret)
#     n = ncodeunits(s)
#     # first continuation byte
#     (i += 1) > n && @goto ret
#     @inbounds b = codeunit(s, i)
#     b & 0xc0 == 0x80 || @goto ret
#     u |= Int32(b) << 16
#     # second continuation byte
#     ((i += 1) > n) | (u < 0xe0000000) && @goto ret
#     @inbounds b = codeunit(s, i)
#     b & 0xc0 == 0x80 || @goto ret
#     u |= Int32(b) << 8
#     # third continuation byte
#     ((i += 1) > n) | (u < 0xf0000000) && @goto ret
#     @inbounds b = codeunit(s, i)
#     b & 0xc0 == 0x80 || @goto ret
#     u |= Int32(b); i += 1
# @label ret
#     return reinterpret(Char, u), i
# end
# Base.thisind(s::VectorString, i::Int) = Base._thisind_str(s, i)

# Base.string(ss::VectorString...) = VectorString(vcat((x.data for x in ss)...))
# # macro S_str(s, flags...) VectorString(s) end
# # macro S_mstr(s, flags...) VectorString(s) end

# export VectorString

# end


text:
"""

`union(dim, arrs...)`

Take uinon of arrays of type `LabeledArray`/`DictArray`/`AbstractArrayrapper`. First start with the first array and remove any duplicate slices in there. Then concat the second array along the specified direction. However, skip any duplicate slices. Proceed similarly for the rest of the arrays.

##### Arguments

* `dim` : direction along which to concat.
* `arrs...` : arrays to unite.

##### Examples

```julia
julia> union(1, nalift([1 2 3;4 5 6]), nalift([1 2 3;5 5 6]))
3x3 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)
 Nullable(4)  Nullable(5)  Nullable(6)
 Nullable(5)  Nullable(5)  Nullable(6)

julia> union(2, nalift([1 2 3;4 5 6]), nalift([1 2 3;5 5 6]))
2x4 DataCubes.AbstractArrayWrapper{Nullable{Int64},2,Array{Nullable{Int64},2}}:
 Nullable(1)  Nullable(2)  Nullable(3)  Nullable(1)
 Nullable(4)  Nullable(5)  Nullable(6)  Nullable(5)

julia> union(1, darr(a=[:x,:y,:z]), darr(a=[:x,:x,:u]), darr(a=[:v,:u,:y]))
5 DictArray

a
--
x
y
z
u
v

julia> union(1, larr(a=[1 2 3;4 5 6], axis1=[:X,:Y]), larr(a=[1 2 3;4 3 2], axis1=[:Z,:Y]))
4 x 3 LabeledArray

  |1 |2 |3
--+--+--+--
  |a |a |a
--+--+--+--
X |1 |2 |3
Y |4 |5 |6
Z |1 |2 |3
Y |4 |3 |2
```

"""
Base.union(dim::Int, arr0::LabeledArray, arr_rest::LabeledArray...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  elems_sofar = Set()
  isdefaxis = isa(permuted_arrs[1].axes[ndimsarr], DefaultAxis)
  result_arrs = if isdefaxis
    map(permuted_arrs) do arr
      arrlastsize = size(arr, ndimsarr)
      coords_toshow = falses(arrlastsize)
      for i in 1:arrlastsize
        full_coords = (remainingcoords..., i)
        oneslice = getindexvalue(arr.data, full_coords...)
        if !(oneslice in elems_sofar)
          push!(elems_sofar, oneslice)
          coords_toshow[i] = true
        end
      end
      arr[remainingcoords..., coords_toshow]
    end
  else
    map(permuted_arrs) do arr
      arrlastsize = size(arr, ndimsarr)
      coords_toshow = falses(arrlastsize)
      for i in 1:arrlastsize
        full_coords = (remainingcoords..., i)
        oneslice = getindexvalue(arr.data, full_coords...)
        labels = arr.axes[ndimsarr][i]
        if !((oneslice,labels) in elems_sofar)
          push!(elems_sofar, (oneslice,labels))
          coords_toshow[i] = true
        end
      end
      arr[remainingcoords..., coords_toshow]
    end
  end
  ipermutedims_if_necessary(cat(ndimsarr, result_arrs...), newdirs)
end

Base.union(dim::Int, arr0::DictArray, arr_rest::DictArray...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  elems_sofar = Set()
  result_arrs = map(permuted_arrs) do arr
    arrlastsize = size(arr, ndimsarr)
    coords_toshow = falses(arrlastsize)
    for i in 1:arrlastsize
      full_coords = (remainingcoords..., i)
      oneslice = getindexvalue(arr, full_coords...)
      if !(oneslice in elems_sofar)
        push!(elems_sofar, oneslice)
        coords_toshow[i] = true
      end
    end
    arr[remainingcoords..., coords_toshow]
  end
  ipermutedims_if_necessary(cat(ndimsarr, result_arrs...), newdirs)
end

Base.union{T,U}(dim::Int,
                arr0::AbstractArrayWrapper{Nullable{T}},
                arr_rest::AbstractArrayWrapper{Nullable{U}}...) = begin
  arrs = (arr0, arr_rest...)
  # first rearrange the axes so that the calculatioin can be more convenient.
  # this might be later optimized.
  # all arrays will be of the same dimensions.
  ndimsarr = ndims(arrs[1])
  remainingdirs = filter(x->x != dim, 1:ndimsarr)
  remainingcoords = ntuple(d->Colon(), length(remainingdirs))
  newdirs = (remainingdirs...,dim)
  permuted_arrs = map(arr->permutedims_if_necessary(arr, newdirs), arrs)
  elems_sofar = Set()
  result_arrs = map(permuted_arrs) do arr
    arrlastsize = size(arr, ndimsarr)
    coords_toshow = falses(arrlastsize)
    for i in 1:arrlastsize
      full_coords = (remainingcoords..., i)
      oneslice = getindex(arr, full_coords...)
      if !(oneslice in elems_sofar)
        push!(elems_sofar, oneslice)
        coords_toshow[i] = true
      end
    end
    arr[remainingcoords..., coords_toshow]
  end
  ipermutedims_if_necessary(cat(ndimsarr, result_arrs...), newdirs)
end


text:
module OutlierDetectionPython
    using OutlierDetectionInterface
    const OD = OutlierDetectionInterface

    include("utils.jl")
    include("models.jl")

    UUID = "2449c660-d36c-460e-a68b-92ab3c865b3e"
    MODELS = [:ABODDetector,
              :CBLOFDetector,
              :COFDetector,
              :COPODDetector,
              :HBOSDetector,
              :IForestDetector,
              :KNNDetector,
              :LMDDDetector,
              :LODADetector,
              :LOFDetector,
              :LOCIDetector,
              :MCDDetector,
              :OCSVMDetector,
              :PCADetector,
              :RODDetector,
              :SODDetector,
              :SOSDetector]

    for model in MODELS
        @eval begin
            OD.@default_frontend $model
            OD.metadata_pkg($model, package_name=string(@__MODULE__), package_uuid=$UUID,
                            package_url="https://github.com/OutlierDetectionJL/$(@__MODULE__).jl",
                            is_pure_julia=false, package_license="MIT", is_wrapper=false)
            OD.load_path(::Type{$model}) = string($model)
            export $model
        end
    end
end