text:
## SPLITTING DATA SETS

# Helper function for partitioning in the non-stratified case
function _partition(rows, fractions, ::Nothing)
    # container for the row selections (head:tail)
    n_splits = length(fractions) + 1
    heads    = zeros(Int, n_splits)
    tails    = zeros(Int, n_splits)
    n_rows   = length(rows)
    head     = 1
    for (i, p) in enumerate(fractions)
        n = round(Int, p * n_rows)
        iszero(n) && (@warn "A split has only one element."; n = 1)
        # update tail
        tail = head + n - 1
        # store
        heads[i] = head
        tails[i] = tail
        # update head
        head = tail + 1
    end
    if head > n_rows
        @warn "Last vector in the split has only one element."
        head = n_rows
    end
    heads[end] = head
    tails[end] = n_rows
    return tuple((rows[h:t] for (h, t) in zip(heads, tails))...)
end

_make_numerical(v::AbstractVector) =
    throw(ArgumentError("`stratify` must have `Count`, `Continuous` "*
                        "or `Finite` element scitpye. Consider "*
                        "`coerce(stratify, Finite)`. "))
_make_numerical(v::AbstractVector{<:Union{Missing,Real}}) = v
_make_numerical(v::AbstractVector{<:Union{Missing,CategoricalValue}}) =
                                int.(v)

# Helper function for partitioning in the stratified case
function _partition(rows, fractions, raw_stratify::AbstractVector)
    stratify = _make_numerical(raw_stratify)
    length(stratify) == length(rows) ||
        throw(ArgumentError("The stratification vector must "*
                            "have as many entries as " *
                            "the rows to partition."))
    uv    = unique(stratify)
    # construct table (n_classes * idx_of_that_class)
    # NOTE use of '===' is important to handle missing.
    idxs  = [[i for i in rows if stratify[rows[i]] === v] for v in uv]

    # number of occurences of each class and proportions
    nidxs = length.(idxs)
    props = length.(idxs) ./ length(rows)

    n_splits = length(fractions) + 1
    n_rows   = length(rows)

    ns_props = round.(Int, n_rows * fractions * props')
    ns_props = vcat(ns_props, nidxs' .- sum(ns_props, dims=1))

    # warn if anything is >= 1
    if !all(e -> e > 1, ns_props)
        @warn "Some splits have a single or no representative of some class."
    end
    # container for the rows
    split_rows = []
    heads      = ones(Int, length(uv))
    for r in 1:size(ns_props, 1)
        tails   = heads .+ ns_props[r, :] .- 1
        # take chunks of the indices corresponding to the current fraction
        indices = vcat((idxs[i][heads[i]:tails[i]] for i in eachindex(uv))...)
        # rearrange by order of appearance
        indices = sort(indices)
        push!(split_rows, rows[indices])
        heads .= tails .+ 1
    end
    if !all(sl -> sl > 1, length.(split_rows))
        @warn "Some splits have a single or no representative of some class."
    end
    return tuple(split_rows...)
end

"""
    partition(rows::AbstractVector{Int}, fractions...;
              shuffle=nothing, rng=Random.GLOBAL_RNG, stratify=nothing)

Splits the vector `rows` into a tuple of vectors whose lengths are
given by the corresponding `fractions` of `length(rows)` where valid
fractions are in (0,1) and sum up to less than 1. The last
fraction is not provided, as it is inferred from the preceding
ones. So, for example,

    julia> partition(1:1000, 0.8)
    ([1,...,800], [801,...,1000])

    julia> partition(1:1000, 0.2, 0.7)
    ([1,...,200], [201,...,900], [901,...,1000])

## Keywords

* `shuffle=nothing`:        if set to  `true`, shuffles the rows before taking fractions.
* `rng=Random.GLOBAL_RNG`:  specifies the random number generator to be used, can be an integer
                            seed. If specified, and `shuffle === nothing` is interpreted as true.
* `stratify=nothing`:       if a vector is specified, the partition will match the stratification
                            of the given vector. In that case, `shuffle` cannot be `false`.
"""
function partition(rows::AbstractVector{Int}, fractions::Real...;
                   shuffle::Union{Nothing,Bool}=nothing, rng=Random.GLOBAL_RNG,
                   stratify::Union{Nothing,AbstractVector}=nothing)
    # if rows is a unitrange, collect it
    rows = collect(rows)
    # check the fractions
    if !all(e -> 0 < e < 1, fractions) || sum(fractions) >= 1
        throw(DomainError(fractions, "Fractions must be in (0, 1) with sum < 1."))
    end
    # check the rng & adjust shuffling
    if rng isa Integer
        rng = MersenneTwister(rng)
    end
    if rng != Random.GLOBAL_RNG && shuffle === nothing
        shuffle = true
    end
    shuffle !== nothing && shuffle && shuffle!(rng, rows)
    return _partition(rows, collect(fractions), stratify)
end

"""

    t1, t2, ...., tk = unnpack(table, f1, f2, ... fk;
                               wrap_singles=false,
                               shuffle=false,
                               rng::Union{AbstractRNG,Int}=nothing)
)

Split any Tables.jl compatible `table` into smaller tables (or
vectors) `t1, t2, ..., tk` by making selections *without replacement*
from the column names defined by the filters `f1`, `f2`, ...,
`fk`. A *filter* is any object `f` such that `f(name)` is `true`
or `false` for each column `name::Symbol` of `table`.

Whenever a returned table contains a single column, it is converted to
a vector unless `wrap_singles=true`.

Scientific type conversions can be optionally specified (note
semicolon):

    unpack(table, t...; wrap_singles=false, col1=>scitype1, col2=>scitype2, ... )

If `shuffle=true` then the rows of `table` are first shuffled, using
the global RNG, unless `rng` is specified; if `rng` is an integer, it
specifies the seed of an automatically generated Mersenne twister. If
`rng` is specified then `shuffle=true` is implicit.

### Example

```
julia> table = DataFrame(x=[1,2], y=['a', 'b'], z=[10.0, 20.0], w=[:A, :B])
julia> Z, XY = unpack(table, ==(:z), !=(:w);
               :x=>Continuous, :y=>Multiclass)
julia> XY
2×2 DataFrame
│ Row │ x       │ y            │
│     │ Float64 │ Categorical… │
├─────┼─────────┼──────────────┤
│ 1   │ 1.0     │ 'a'          │
│ 2   │ 2.0     │ 'b'          │

julia> Z
2-element Array{Float64,1}:
 10.0
 20.0
```
"""
function unpack(X, tests...;
                wrap_singles=false,
                shuffle=nothing,
                rng=nothing, pairs...)

    shuffle, rng = shuffle_and_rng(shuffle, rng)

    shuffle && (X = selectrows(X, Random.shuffle(rng, 1:nrows(X))))

    if isempty(pairs)
        Xfixed = X
    else
        Xfixed = coerce(X, pairs...)
    end

    unpacked = Any[]
    names_left = schema(Xfixed).names |> collect
    history = ""
    counter = 1
    for c in tests
        names = filter(c, names_left)
        filter!(!in(names), names_left)
        history *= "selection $counter: $names\n remaining: $names_left\n"
        isempty(names) &&
            error("Empty column selection encountered at selection $counter"*
                  "\n$history")
        length(names) == 1 && !wrap_singles && (names = names[1])
        push!(unpacked, selectcols(Xfixed, names))
        counter += 1
    end
    return Tuple(unpacked)
end

## RESTRICTING TO A FOLD

struct FoldRestrictor{i,N}
    f::NTuple{N,Vector{Int}}
end
(r::FoldRestrictor{i})(X) where i = selectrows(X, (r.f)[i])

"""
    restrict(X, folds, i)

The restriction of `X`, a vector, matrix or table, to the `i`th fold
of `folds`, where `folds` is a tuple of vectors of row indices.

The method is curried, so that `restrict(folds, i)` is the operator
on data defined by `restrict(folds, i)(X) = restrict(X, folds, i)`.

### Example

    folds = ([1, 2], [3, 4, 5],  [6,])
    restrict([:x1, :x2, :x3, :x4, :x5, :x6], folds, 2) # [:x3, :x4, :x5]

See also [`corestrict`](@ref)

"""
restrict(f::NTuple{N}, i) where N = FoldRestrictor{i,N}(f)
restrict(X, f, i) = restrict(f, i)(X)

## RESTRICTING TO A FOLD COMPLEMENT

"""
    complement(folds, i)

The complement of the `i`th fold of `folds` in the concatenation of
all elements of `folds`. Here `folds` is a vector or tuple of integer
vectors, typically representing row indices or a vector, matrix or
table.

    complement(([1,2], [3,], [4, 5]), 2) # [1 ,2, 4, 5]

"""
complement(f, i) = reduce(vcat, collect(f)[Not(i)])

struct FoldComplementRestrictor{i,N}
    f::NTuple{N,Vector{Int}}
end
(r::FoldComplementRestrictor{i})(X) where i =
    selectrows(X, complement(r.f, i))

"""
    corestrict(X, folds, i)

The restriction of `X`, a vector, matrix or table, to the *complement*
of the `i`th fold of `folds`, where `folds` is a tuple of vectors of
row indices.

The method is curried, so that `corestrict(folds, i)` is the operator
on data defined by `corestrict(folds, i)(X) = corestrict(X, folds, i)`.

### Example

    folds = ([1, 2], [3, 4, 5],  [6,])
    corestrict([:x1, :x2, :x3, :x4, :x5, :x6], folds, 2) # [:x1, :x2, :x6]

"""
corestrict(f::NTuple{N}, i) where N = FoldComplementRestrictor{i,N}(f)
corestrict(X, f, i) = corestrict(f, i)(X)

## to be replaced (not used anywhere):
## ACCESSORS FOR JULIA NDSPARSE ARRAYS (N=2)

# nrows(::Val{:sparse}, X) = maximum([r[1] for r in keys(X)])
# function select(::Val{:sparse}, X, r::Integer, c::Symbol)
#     try
#         X[r,c][1]
#     catch exception
#         exception isa KeyError || throw(exception)
#         missing
#     end
# end
# select(::Val{:sparse}, X, r::AbstractVector{<:Integer}, c::Symbol) = [select(X, s, c) for s in r]
# select(::Val{:sparse}, X, ::Colon, c::Symbol) = [select(X, s, c) for s in 1:nrows(X)]
# selectrows(::Val{:sparse}, X, r::Integer) = X[r:r,:]
# selectrows(::Val{:sparse}, X, r) = X[r,:]
# selectcols(::Val{:sparse}, X, c::Symbol) = select(X, :, c)
# selectcols(::Val{:sparse}, X, c::AbstractVector{Symbol}) = X[:,sort(c)]
# selectcols(::Val{:sparse}, X, ::Colon) = X
# select(::Val{:sparse}, X, r::Integer, c::AbstractVector{Symbol}) = X[r,sort(c)]
# select(::Val{:sparse}, X, r::Integer, ::Colon) = X[r,:]
# select(::Val{:sparse}, X, r, c) = X[r,sort(c)]

## TRANSFORMING BETWEEN CATEGORICAL ELEMENTS AND RAW VALUES

_err_missing_class(c) =  throw(DomainError(
    "Value $c not pool"))

function transform_(pool, x)
    ismissing(x) && return missing
    x in levels(pool) || _err_missing_class(x)
    return pool[get(pool, x)]
end

transform_(pool, X::AbstractArray) = broadcast(x -> transform_(pool, x), X)

"""
    transform(e::Union{CategoricalElement,CategoricalArray,CategoricalPool},  X)

Transform the specified object `X` into a categorical version, using
the pool contained in `e`. Here `X` is a raw value (an element of
`levels(e)`) or an `AbstractArray` of such values.

```julia
v = categorical([:x, :y, :y, :x, :x])
julia> transform(v, :x)
CategoricalValue{Symbol,UInt32} :x

julia> transform(v[1], [:x :x; missing :y])
2×2 CategoricalArray{Union{Missing, Symbol},2,UInt32}:
 :x       :x
 missing  :y

"""
MLJModelInterface.transform(e::Union{CategoricalArray, CategoricalValue},
                            arg) = transform_(CategoricalArrays.pool(e), arg)
MLJModelInterface.transform(e::CategoricalPool, arg) =
    transform_(e, arg)


text:
using CpuId

num_machines() = static(1)
num_sockets() = static(1)

function _get_num_threads()::Int
  (Sys.CPU_THREADS)::Int >> (Sys.ARCH !== :aarch64)
end

_get_num_cores()::Int = clamp(CpuId.cpucores(), 1, (Sys.CPU_THREADS)::Int)

let nc = static(_get_num_cores())
  global num_l1cache() = nc
  global num_cores() = nc
end
let syst = static((Sys.CPU_THREADS)::Int)
  global sys_threads() = syst
  global num_threads() = syst
end
@static if Sys.ARCH === :aarch64
  num_l2cache() = static(1)
  num_l3cache() = static(0)
else
  num_l2cache() = num_l1cache()
  num_l3cache() = static(1)
end
num_l4cache() = static(0)

const PrecompiledCacheSize = CpuId.cachesize()
const PrecompiledCacheInclusive = CpuId.cacheinclusive()
cache_inclusive(_) = False()
@noinline function _eval_cache_size(cachesize)
  for (i, csi) in enumerate(cachesize)
    @eval cache_size(::Union{Val{$i},StaticInt{$i}}) = $(static(csi))
  end
end
@noinline function _eval_cache_inclusive(cacheinclusive)
  for (i, cii) in enumerate(cacheinclusive)
    @eval cache_inclusive(::Union{Val{$i},StaticInt{$i}}) = $(static(cii != 0))
  end
end
_eval_cache_size(PrecompiledCacheSize)
_eval_cache_inclusive(PrecompiledCacheInclusive)
# TODO: implement
cache_associativity(_) = static(0)

cache_type(::Union{Val{1},StaticInt{1}}) = Val{:Data}()
cache_type(_) = Val{:Unified}()
# cache_type(::Union{Val{2},StaticInt{2}}) = Val{:Unified}()
# cache_type(::Union{Val{3},StaticInt{3}}) = Val{:Unified}()
let lnsize = static(CpuId.cachelinesize())
  global cache_linesize(_) = lnsize
end
cache_size(_) = StaticInt{0}()

# cache_size(::Union{Val{3},StaticInt{3}}) = num_cores() * StaticInt{1441792}()
function _extra_init()
  nc = _get_num_cores()
  if (nc != CpuId.cpucores())
    cache_l3_per_core = CpuId.cachesize(3) ÷ max(CpuId.cpucores(), 1)
    @eval cache_size(::Union{Val{3},StaticInt{3}}) = $(static(cache_l3_per_core * nc))
  end
  cs = CpuId.cachesize()
  cs === PrecompiledCacheSize || _eval_cache_size(cs)
  ci = CpuId.cacheinclusive()
  ci === PrecompiledCacheInclusive || _eval_cache_inclusive(ci)
end

ext: jl


text:
# A walk-through example (using StanSample.jl)

using StanSample
using Distributions
using DataFrames
using MonteCarloMeasurements
using Tables

N = 100
df = DataFrame(
    height = rand(Normal(10, 2), N),
    leg_prop = rand(Uniform(0.4, 0.5), N),
)
df.leg_left = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
df.leg_right = df.leg_prop .* df.height + rand(Normal(0, 0.03), N)

stan6_1 = "
data {
  int <lower=1> N;
  vector[N] H;
  vector[N] LL;
  vector[N] LR;
}
parameters {
  real a;
  vector[2] b;
  real <lower=0> sigma;
}
model {
  vector[N] mu;
  mu = a + b[1] * LL + b[2] * LR;
  a ~ normal(10, 100);
  b ~ normal(2, 10);
  sigma ~ exponential(1);
  H ~ normal(mu, sigma);
}
";

m6_1s = SampleModel("m6.1s", stan6_1);
data = (H = df.height, LL = df.leg_left, LR = df.leg_right, N = size(df, 1))
rc6_1s = stan_sample(m6_1s; data);

if success(rc6_1s)
    chns = read_samples(m6_1s, :nesteddataframe)
    chns |> display

    # Obtain  b.1 , b,2 as a matrix

    chns = matrix(chns, :b)

    # Or use read_samples to only use chains 2 and 4 using the chains kwarg.

    chns2 = read_samples(m6_1s, :dataframe; chains=[2, 4])
    mean.(eachcol(chns2)) |> display

end

# Some more options:

init = (a = 2.0, b = [1.0, 2.0], sigma = 1.0)
#rc6_2s = stan_sample(m6_1s; data, init);
rc6_2s = stan_sample(m6_1s; data);

if success(rc6_2s)

    read_summary(m6_1s, true)
    println()

    post6_1s_df = read_samples(m6_1s, :dataframe)

    part6_1s = read_samples(m6_1s, :particles)
    part6_1s |> display
    println()

    nt6_1s = read_samples(m6_1s, :namedtuple)
    mean(nt6_1s.b, dims=2)
end


text:
using VetFun
using Test

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


text:
function topological_sort{T}(outgoing::Dict{T, Vector{T}})
    # kahn topological sort
    # todo: use better data structures
    order = T[]
    incoming = Dict{T,Vector}()
    for (n, children) in outgoing
        for child in children
            if haskey(incoming, child)
                push!(incoming[child], n)
            else
                incoming[child] = [n]
            end
        end
        if !haskey(incoming, n)
            incoming[n] = T[]
        end
    end
    S = T[]

    outgoing_counts = Dict{T,Int}()
    for (child, parents) in outgoing
        if isempty(parents)
            push!(S, child)
        end
        outgoing_counts[child] = length(parents)
    end
    while !isempty(S)
        n = pop!(S)
        push!(order, n)
        children = incoming[n]
        for child in children
            outgoing_counts[child] -= 1
            if outgoing_counts[child] == 0
                push!(S, child)
            end
        end
    end
    if length(order) != length(outgoing)
        error("Cyclic dependencies detected")
    end
    return order
end

function get_dependency_graph()
    deps = [filename=>afile.deps for (filename, afile) in files]
    return deps
end

function _extract_deps(e::Expr, deps::Vector{Dep})
    if e.head == :call
        if e.args[1] == :include
            if isa(e.args[2], AbstractString)
                push!(deps, Dep(false, e.args[2]))
            end
        # elseif e.args[1] == :require
        #     if isa(e.args[2], AbstractString)
        #         push!(deps, Dep(true, e.args[2]))
        #     end
        end
    elseif e.head == :import
        # push!(deps, Dep(true,    string(e.args[1])))
    end
    for arg in e.args
        if isa(arg, Expr)
            _extract_deps(arg, deps)
        end
    end
end

function extract_deps(e::Expr)
    deps = Dep[]
    _extract_deps(e, deps)
    return deps
end


text:
# Model
include_model("hopper")
model_ft = free_time_model(no_slip_model(model))

# Horizon]
T = 51

# Time step
tf = 0.5
h = tf / (T - 1)

# Bounds
_uu = Inf * ones(model_ft.m)
_uu[model_ft.idx_u] = model_ft.uU
_uu[end] = 2.0 * h
_ul = zeros(model_ft.m)
_ul[model_ft.idx_u] = model_ft.uL
_ul[end] = 0.5 * h
ul, uu = control_bounds(model_ft, T, _ul, _uu)

# Initial and final states
z_h = 0.25
q1 = [0.0, 0.5 + z_h, 0.0, 0.25]

xl, xu = state_bounds(model_ft, T,
		[model_ft.qL; model_ft.qL],
		[model_ft.qU; model_ft.qU],
        x1 = [q1; Inf * ones(model.nq)])

# Objective
include_objective("velocity")
obj_tracking = quadratic_time_tracking_objective(
    [Diagonal(zeros(model_ft.n)) for t = 1:T],
    [Diagonal([1.0, 1.0, zeros(model_ft.m - model_ft.nu - 1)..., 0.0]) for t = 1:T-1],
    [zeros(model_ft.n) for t = 1:T],
    [zeros(model_ft.m) for t = 1:T],
    1.0)
obj_contact_penalty = PenaltyObjective(1.0e5, model_ft.m - 1)
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model_ft.nq)) for t = 1:T-1],
    model_ft.nq,
    h = h,
    idx_angle = collect([3]))
obj = MultiObjective([obj_tracking, obj_contact_penalty, obj_velocity])

# Constraints
include_constraints(["free_time", "contact_no_slip", "loop"])
con_free_time = free_time_constraints(T)
con_contact = contact_no_slip_constraints(model_ft, T)
con_loop = loop_constraints(model, (1:model.n), 1, T)
con = multiple_constraints([con_free_time, con_contact, con_loop])

# Problem
prob = trajectory_optimization_problem(model_ft,
               obj,
               T,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con)

# Trajectory initialization
x0 = [[q1; q1] for t = 1:T] # linear interpolation on state
u0 = [[1.0e-5 * rand(model_ft.m-1); h] for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

optimize = true
if optimize
	include_snopt()
	@time z̄ , info = solve(prob, copy(z0),
		nlp = :ipopt,
		tol = 1.0e-5, c_tol = 1.0e-5, mapl = 5,
		time_limit = 60)
	@show check_slack(z̄, prob)
	x̄, ū = unpack(z̄, prob)
	tf, t, h̄ = get_time(ū)

	# projection
	Q = [Diagonal(ones(model_ft.n)) for t = 1:T]
	R = [Diagonal(0.1 * ones(model_ft.m)) for t = 1:T-1]
	x_proj, u_proj = lqr_projection(model_ft, x̄, ū, h̄[1], Q, R)

	@show tf
	@show h̄[1]
	@save joinpath(pwd(), "examples/trajectories/hopper_vertical_gait_no_slip.jld2") x̄ ū h̄ x_proj u_proj
else
	@load joinpath(pwd(), "examples/trajectories/hopper_vertical_gait_no_slip.jld2") x̄ ū h̄ x_proj u_proj
end

using Plots
plot(hcat(ū...)[1:2, :]',
    linetype = :steppost,
    label = "",
    color = :red,
    width = 2.0)

plot!(hcat(u_proj...)[1:2, :]',
    linetype = :steppost,
    label = "", color = :black)

plot(hcat(state_to_configuration(x̄)...)',
    color = :red,
    width = 2.0,
	label = "")
plot!(hcat(state_to_configuration(x_proj)...)',
    color = :black,
	label = "")

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model_ft, state_to_configuration(x_proj), Δt = h̄[1])


text:
data = joinpath("data", "testset1")
result = AnymodResult(data; threaded=false)

@testset "Testset 1" begin
    @testset "Read-in" begin
        @test typeof(result) == AnymodResult
        @test result.scenarios == ["desintegriert", "integriert"]
        @test size(result.summarytable) == (21340, 11)
        @test size(result.costtable) == (6116, 11)
        @test size(result.exchangetable) == (6338, 10)

        summarytable_names = [
            "variable",
            "value",
            "scenario",
            "region_dispatch_1",
            "region_dispatch_2",
            "technology_1",
            "technology_2",
            "technology_3",
            "carrier_1",
            "carrier_2",
            "timestep_superordinate_dispatch_1"
        ]

        @test names(result.summarytable) == summarytable_names

        costtable_names = [
            "variable",
            "value",
            "scenario",
            "region_1",
            "region_2",
            "technology_1",
            "technology_2",
            "technology_3",
            "carrier_1",
            "carrier_2",
            "timestep_superordinate_dispatch_1"
        ]

        @test names(result.costtable) == costtable_names

        exchangetable_names = [
            "value",
            "variable",
            "scenario",
            "region_from_1",
            "region_from_2",
            "region_to_1",
            "region_to_2",
            "carrier_1",
            "carrier_2",
            "timestep_superordinate_dispatch_1"
        ]

        @test names(result.exchangetable) == exchangetable_names
    end

    file_path = joinpath("data", "testset1_results")
    @testset "Tables" begin

        mask1 = Mask(
            Scenario(),
            Technology(1 => ["wind_offshore","wind_onshore","pv"]),
            Region(2 => "DE", mode=:occursin),
            Variable("capaConv"),
            Timestep("2050")
        )

        mask1_result = CSV.read(
            joinpath(file_path, "mask1_result.csv"),
            DataFrame
        )

        @test pivotresult(result, mask1) == mask1_result

        mask2 = Mask(
            Technology(2),
            Region(1 => ["AT", "DE", "CH", "FR", "NL", "DK"]),
            Variable("capaConv", "capaStOut"),
            Scenario("integriert")
        )

        mask2_result = CSV.read(
            joinpath(file_path, "mask2_result.csv"),
            DataFrame
        )

        df2 = pivotresult(result, mask2)
        @test compare_table_values(mask2_result, df2)

        mask3 = Mask(
            Carrier(),
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Scenario("integriert"),
            Region(2 => "DE", mode = :occursin),
            Timestep(["2050"])
        )

        mask3_result = CSV.read(
            joinpath(file_path, "mask3_result.csv"),
            DataFrame
        )

        df3 = pivotresult(result, mask3)
        @test compare_table_values(mask3_result, df3)

        mask_warn = Mask(
            Carrier(4 => ["asdf"]),
            Region(5 => "DE", mode = :occursin)
        )

        @suppress @test_throws ArgumentError pivotresult(result, mask_warn)

        result_mask1 = Mask(
            [Scenario(), Carrier()],
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Region(2 => "DE", mode = :occursin)
        )

        stacked_result1 = CSV.read(
            joinpath(file_path, "stacked_result_1.csv"),
            DataFrame
        )

        df_stacked = stackedresult(result, result_mask1)
        @test compare_table_values(stacked_result1, df_stacked)
    end

    @testset "Plots" begin

        ### Plot 1 ###
        result_mask1 = Mask(
            [Scenario(), Carrier()],
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Region(2 => "DE", mode = :occursin)
        )

        p1 = json(plot_anymod_result(result, result_mask1))
        p1_result = JSON.read(joinpath(file_path, "plot1.json"), String)

        @test p1[1:1000] == p1_result[1:1000]

        ### Plot 2 ###
        result_mask2 = Mask(
            Scenario(),
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Region(2 => "DE", mode = :occursin),
            Carrier("electricity")
        )

        p2 = plot_anymod_result(result, result_mask2) |> json
        p2_result = JSON.read(joinpath(file_path, "plot2.json"), String)
        @test p2[1:1000] == p2_result[1:1000]

        ### Plot 3 ###
        result_mask3 = Mask(
            Scenario(),
            nothing,
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Region(2 => "DE", mode = :occursin),
            Carrier("electricity")
        )

        p3 = plot_anymod_result(result, result_mask3) |> json
        p3_result = JSON.read(joinpath(file_path, "plot3.json"), String)
        @test p3[1:1000] == p3_result[1:1000]

        ### Plot 4 ###
        result_mask4 = Mask(
            [Scenario(), Carrier()],
            nothing,
            Variable(["demand", "use","gen","export","import","trdBuy"]),
            Region(2 => "DE", mode = :occursin)
        )

        p4 = plot_anymod_result(result, result_mask4) |> json
        p4_result = JSON.read(joinpath(file_path, "plot4.json"), String)
        @test p4[1:1000] == p4_result[1:1000]

    end
end


text:
abstract type AnanthModel <: DiabaticModel end

NQCModels.ndofs(::AnanthModel) = 1
NQCModels.nstates(::AnanthModel) = 2

"""
    AnanthModelOne(a=0.01, b=1.6, c=0.005, d=1.0)

Ananth's simple avoided crossing model (similar to Tully's first model) from [J. Chem. Phys. 127, 084114 (2007)](http://dx.doi.org/10.1063/1.2759932).
"""
Parameters.@with_kw struct AnanthModelOne{A,B,C,D} <: AnanthModel
    a::A = 0.01
    b::B = 1.6
    c::C = 0.005
    d::D = 1.0
end

function NQCModels.potential(model::AnanthModelOne, q::Real)
    Parameters.@unpack a, b, c, d = model
    V11 = a * tanh(b*q)
    V22 = -a * tanh(b*q)
    V12 = c * exp(-d*q^2) # sets both off-diagonals as V is Hermitian
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::AnanthModelOne, q::Real)
    Parameters.@unpack a, b, c, d = model
    D11 = a * b * (1 - tanh(b*q)^2)
    D22 = -D11
    D12 = -2 * c * d * q * exp(-d*q^2)
    return Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end

"""
    AnanthModelTwo(a=0.04, b=0.01, c=0.005, d=1.0, e=0.7, f=1.6)

Ananth's asymmetric model from [J. Chem. Phys. 127, 084114 (2007)](http://dx.doi.org/10.1063/1.2759932).
"""
Parameters.@with_kw struct AnanthModelTwo{A,B,C,D,E,F} <: AnanthModel
    a::A = 0.04
    b::B = 0.01
    c::C = 0.005
    d::D = 1.0
    e::E = 0.7
    f::F = 1.6
end

function NQCModels.potential(model::AnanthModelTwo, q::Real)
    Parameters.@unpack a, b, c, d, e, f = model
    V11 = a * tanh(f*q)
    V22 = -b * tanh(f*q)
    V12 = c * exp(-d*(q+e)^2)
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::AnanthModelTwo, q::Real)
    Parameters.@unpack a, b, c, d, e, f = model
    D11 = a * f * (1 - tanh(f*q)^2)
    D22 = -b * f * (1 - tanh(f*q)^2)
    D12 = -2 * c * d * (q+e) * exp(-d*(q+e)^2)

    Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end


text:
# BioCore.jl
# ==========
#
# Core types and methods common to many packages in the BioJulia ecosystem.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

__precompile__()

module BioCore

include("declare.jl")
include("Exceptions.jl")
include("IO.jl")
include("Mem.jl")
include("Ragel.jl")
include("ReaderHelper.jl")
include("RecordHelper.jl")
include("StringFields.jl")
include("Testing.jl")

end # module BioCore


text:
using Distributions
using ParticleMDI

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

## Test data types
# Gaussian
n = 1000
test_data = rand(Normal(0, 1), n, 1)
test_cluster = ParticleMDI.GaussianCluster(test_data)
# Basic check it's been initialised
@test test_cluster.n == 0
# Check the logprob
# ParticleMDI.calc_logprob(test_data[1, :], test_cluster)
for obs in test_data
    ParticleMDI.cluster_add!(test_cluster, [obs], [true])
end

@test test_cluster.n == n
@test isapprox(test_cluster.Σ[1], sum(test_data))
@test isapprox(test_cluster.μ[1], test_cluster.Σ[1] / (n + 0.001))
xbar = test_cluster.Σ[1] / n
s2 = sum((test_data .- xbar) .^ 2)
beta = 0.5 + 0.5 * (s2 + (0.001 * n * xbar ^ 2) / (n + 0.001))
@test isapprox(test_cluster.β[1], beta)
@test isapprox(test_cluster.λ[1], ((0.5 + n * 0.5) * (n + 0.001)) / (test_cluster.β[1] * (n + 1.001)))

xcentred = (test_data[end, :][1] - test_cluster.μ[1]) * sqrt(test_cluster.λ[1])
truelogprob = logpdf(TDist(test_cluster.n + 1), xcentred) + 0.5 * log(test_cluster.λ[1])
estlogprob = ParticleMDI.calc_logprob(test_data[end, :], test_cluster, [true])
@test isapprox(truelogprob, estlogprob)

# Categorical
test_data = rand(1:10, 1000, 1)
test_cluster = ParticleMDI.CategoricalCluster(test_data)
@test test_cluster.n == 0

for obs in test_data
    ParticleMDI.cluster_add!(test_cluster, [obs], [true])
end

@test test_cluster.n == 1000
test_cluster.counts
for x in unique(test_data)
    @test sum(test_data .== x) == test_cluster.counts[x, 1]
end

@test isapprox(ParticleMDI.calc_logprob([1], test_cluster, [true]),
      log((sum(test_data .== 1) + 0.5) / (1005)))

# Test normalising constant calculation
for N = 2:20
    for K = 1:5
        c_combn = Matrix{Int64}(undef, N ^ K, K)
        for k in 1:K
            c_combn[:, K - k + 1] = div.(0:(N ^ K - 1), N ^ (K - k)) .% N .+ 1
        end
        γc = rand(Gamma(1.0 / N, 1), N, K)
        Γ = Matrix{Float64}(undef, N ^ K, K)
        for i in 1:(N ^ K)
            for k in 1:K
                Γ[i, k] = γc[c_combn[i], k]
            end
        end
        if K > 1
            Φ = rand(Gamma(1, 5), binomial(K, 2), 1)
        else
            Φ = zeros(1)
        end
        phi_index = zeros(Int64, K- 1, K)
        num = 1
        for i in 1:(K - 1)
            for j in (i + 1):K
                phi_index[i, j] = num
                num += 1
            end
        end
        Z = 0.0
        for i in 1:(N ^K)
            tmp = prod(Γ[i, :])
            if K > 1
                for k1 in 1:(K - 1)
                    for k2 in (k1 + 1):K
                        tmp *=  (1 + (Φ[phi_index[k1, k2]] * (c_combn[i, k1] == c_combn[i, k2])))
                    end
                end
            end
            Z += tmp
        end
        Φ_index = K > 1 ? Matrix{Bool}(undef, N ^ K, Int64(K * (K - 1) / 2)) : fill(1, (N, 1))
        if K > 1
            i = 1
            for k1 in 1:(K - 1)
                for k2 in (k1 + 1):K
                    Φ_index[:, i] = (c_combn[:, k1] .== c_combn[:, k2])
                    i += 1
                end
            end
        end
        @test isapprox(Z, ParticleMDI.update_Z(Φ, Φ_index, log.(Γ)))
    end
end

# Test cluster alignment
# All datasets strongly agree
K = 5
N = 10
s = rand(1:N, 10000, K)
Φ = [10 for k in 1:binomial(K, 2)]
γ = rand(Gamma(1.0 / N, 1), N, K)
# Perfect agreement up to label permutation
for k in 2:K
    shuf = Distributions.shuffle(1:N)
    s[:, k] .= shuf[s[:, 1]]
    γ[:, k] = γ[indexin(1:N, shuf), 1]
end
# No guarantee that agreement occurs at first attempt
# Run over and over and check that if all the cluster values align
# So do all the gammas
for i = 1:10
    ParticleMDI.align_labels!(s, Φ, γ, N, K)
    @test all(s[:, 2:K] .== s[:, 1]) == all(γ[:, 2:K] .== γ[:, 1])
end
ParticleMDI.align_labels!(s, Φ, γ, N, 0)

@test all(s[:, 2:K] .== s[:, 1])
@test all(γ[:, 2:K] .== γ[:, 1])

# Check the particle IDs are doing what they should
# ie particle IDs
data1 = [rand(Normal(2, 1), 50, 16);
        rand(Normal(-2, 1), 50, 16)]
data2 = [rand(Normal(2, 1), 50, 16);
        rand(Normal(-2, 1), 50, 16)]
data3 = [rand(Normal(2, 1), 50, 16);
        rand(Normal(-2, 1), 50, 16)]
dataFiles = [data1, data2, data3]
dataTypes = [ParticleMDI.GaussianCluster, ParticleMDI.GaussianCluster, ParticleMDI.GaussianCluster]
particle, cluster, counts, n_ops = ParticleMDI.__pmdi(dataFiles, dataTypes, 10, 2, 25 / 100, 1, featureSelect = nothing)
@test all([sum([cluster[k][particle[i, j, k]].n for i in 1:10]) for j in 1:2 for k in 1:3] .== 100)
# Check that cluster counts is working correctly
for k in 1:3
    for i in 1:21
    @test counts[i, k] == count(x -> x == i, particle[:, :, k])
    end
end

particle, cluster, counts, n_ops = ParticleMDI.__pmdi(dataFiles, dataTypes, 10, 1024, 25 / 100, 100, featureSelect = nothing)
@test all([sum([cluster[k][particle[i, j, k]].n for i in 1:10]) for j in 1:1024 for k in 1:3] .== 100)
# Check that cluster counts is working correctly
for k in 1:3
    for i in 1:10241
    @test counts[i, k] == count(x -> x == i, particle[:, :, k])
    end
end


text:
function boot_up(;verbose = false)
    if verbose
        @info "Booting up..."
    end
    img = rand(RGB{N0f8}, (72, 128))
    frame = VsFrame(image = img)
    config = vs_setup(PokemonBattle, init_descriptors = false, set_global = false)
    ctx = VsRecorderBase.initialize(config)
    for T in VsRecorder.Scenes.AvailableScenes
        if verbose
            @info "Compiling for $T..."
        end
        vs_tryparse_scene(T, frame, ctx, force = true)
    end
    if verbose
        @info "Booted."
    end
    nothing
end


text:
# module for drawing moving robots
# one of them move normally and the others sometimes stop
# because of stuck

module AnimeMoveStuck
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function animate_per_time(time_interval, world, stuck_robot,
                            no_stuck_robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(stuck_robot)
    draw!(no_stuck_robot)
  end

  function main(time_interval=0.1; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    init_pose = [0.0, 0.0, 0.0]
    circling = Agent(0.2, 10.0/180*pi)

    sr = DifferentialWheeledRobot(init_pose, 0.2, "gray",
                                  circling, time_interval,
                                  exp_stuck_time=60.0,
                                  exp_escape_time=60.0)

    nsr = DifferentialWheeledRobot(init_pose, 0.2, "red",
                                   circling, time_interval)

    if is_test == false
      anime = @animate for t in 0:time_interval:30
        animate_per_time(t, w, sr, nsr)
      end
      path = "src/model/uncertainty/movement/stuck/anime_move_stuck.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
    else
      for t in 0:time_interval:10
        animate_per_time(t, w, sr, nsr)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
using SpinMC
using LinearAlgebra

a1 = (3/2, sqrt(3)/2)
a2 = (3/2, -sqrt(3)/2)
uc = UnitCell(a1,a2)

b1 = addBasisSite!(uc, (0.0, 0.0))
b2 = addBasisSite!(uc, (1.0, 0.0))

M = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
addInteraction!(uc, b1, b2, M, (0, 0))
addInteraction!(uc, b1, b2, M, (0, -1))
addInteraction!(uc, b1, b2, M, (-1, 0))

L = (16, 16)
lattice = Lattice(uc, L)

thermalizationSweeps = 50000
measurementSweeps = 50000
beta = 10.0
m = MonteCarlo(lattice, beta, thermalizationSweeps, measurementSweeps)
run!(m)

# Fourier transform correlations to compute structure factor.
N = 256
correlation = mean(m.observables.correlation) # The correlation is measured with respect to the first spin, i.e. the i-th entry is the correlation dot(S_1,S_i).
kx = collect(range(-2pi,2pi,length=N))
ky = collect(range(-2pi,2pi,length=N))
structurefactor = zeros(N,N)
for i in 1:N
    for j in 1:N
        z = 0.0
        # Compute Fourier transformation at momentum (kx, ky). The real-space position of the i-th spin is obtained via getSitePosition(lattice,i).
        for k in 1:length(lattice)
            z += cos(dot((kx[i],ky[j]),getSitePosition(lattice,k))) * correlation[k]
        end
        structurefactor[j,i] = z / length(lattice)
    end
end

# Plot result
using Plots
heatmap(kx,ky,structurefactor,aspect_ratio=1,xrange=(-2pi,2pi),yrange=(-2pi,2pi),clims=(0,1),xlabel="kx", ylabel="ky")
xs = [4.0pi/3.0, 2.0pi/3.0, -2.0pi/3.0, -4.0pi/3.0, -2.0pi/3.0, 2.0pi/3.0, 4.0pi/3.0]
ys = [0.0, 2.0pi/sqrt(3.0), 2.0pi/sqrt(3.0), 0.0, -2.0pi/sqrt(3.0), -2.0pi/sqrt(3.0), 0.0]
plot!(xs, ys, label="Extended BZ")


text:
using StatsBase
using JLD2

daily(times::AbstractArray{T} where T, max_days::Integer) = fit(Histogram, times, 0:(max_days) ).weights
daily(times::AbstractVector{T} where T) = daily(times, maximum(times))

struct DailyTrajectories <: Output
  file::JLD2.JLDFile
end

DailyTrajectories(fname::AbstractString) = DailyTrajectories(JLD2.jldopen(fname, "w+", compress=true))
DailyTrajectories(fname::AbstractString, ::Integer) = DailyTrajectories(fname)

function pushtrajectory!(d::DailyTrajectories, trajectory_id::Integer, writelock::Base.AbstractLock, state::MocosSim.SimState, params::MocosSim.SimParams, cb::DetectionCallback)
  try
    lock(writelock)
    trajectory_group = JLD2.Group(d.file, string(trajectory_id))
    save_daily_trajectories(trajectory_group, state, params, cb)
  finally
    unlock(writelock)
  end
  nothing
end

aftertrajectories(d::DailyTrajectories) = close(d.file)

function save_daily_trajectories(dict, state::MocosSim.SimState, params::MocosSim.SimParams, cb::DetectionCallback)
  max_days = MocosSim.time(state) |> floor |> Int
  num_individuals = MocosSim.numindividuals(state)

  infection_times = Vector{OptTimePoint}(missing, num_individuals)
  contact_kinds = Vector{MocosSim.ContactKind}(undef, num_individuals)

  for i in 1:num_individuals
    event = MocosSim.backwardinfection(state, i)
    kind = contactkind(event)
    contact_kinds[i] = kind
    infection_times[i] = ifelse(kind == MocosSim.NoContact, time(event), missing)
  end

  hospitalization_progressions = getproperty.(params.progressions, :severe_symptoms_time)
  death_progressions = getproperty.(params.progressions, :death_time)

  dict["infection_times"] = daily(filter(!ismissing,infection_times), max_days)
  dict["detection_tiems"] = daily(filter(!ismissing, cb.detection_times), max_days)
  dict["daily_deaths"] = daily(filter(!ismissing, infection_times.+death_progressions), max_days)
  dict["daily_hospitalizations"] = daily(filter(!ismissing, infection_times.+hospitalization_progressions), max_days)
  for kind in instances(ContactKind)
    if kind != NoContact
      dict["daily_" * lowercase(string(kind))] = daily(infection_times[contact_kinds.==Int(kind)], max_days)
    end
  end
end

max_stars_repo_name: Novirmor/MocosSimLauncher

max_issues_repo_name: Novirmor/MocosSimLauncher

max_forks_repo_name: Novirmor/MocosSimLauncher


text:
import Base: keys, iterate, length, getindex

"""
    GDTable

Structure representing a grouped `DTable`.
It wraps over a `DTable` object and provides additional information on how the table is grouped.
To represent the grouping a `cols` field is used, which contains the column symbols used for
grouping and an `index`, which allows to effectively lookup the partitions grouped under a single key.
"""
mutable struct GDTable
    dtable::DTable
    cols::Union{Vector{Symbol}, Nothing}
    index::Dict
    grouping_function::Union{Function, Nothing}

    GDTable(dtable, cols, index) = GDTable(dtable, cols, index, nothing)
    GDTable(dtable, cols, index, grouping_function) =
        new(dtable, cols, deepcopy(index), grouping_function)
end

DTable(gd::GDTable) = DTable(gd.dtable.chunks, gd.dtable.tabletype)

fetch(gd::GDTable) = fetch(gd.dtable)
fetch(gd::GDTable, sink) = fetch(gd.dtable, sink)

"""
    grouped_cols(gd::GDTable) -> Vector{Symbol}

Returns the symbols of columns used in the grouping.
In case grouping on a function was performed a `:KEYS` symbol will be returned.
"""
grouped_cols(gd::GDTable) = gd.cols === nothing ? [:KEYS] : gd.cols

"""
    keys(gd::GDTable) -> KeySet

Returns the keys that `gd` is grouped by.
"""
keys(gd::GDTable) = keys(gd.index)

partition(gd::GDTable, key) = partition(gd, gd.index[key])
partition(gd::GDTable, indices::Vector{UInt}) = DTable(VTYPE(getindex.(Ref(gd.dtable.chunks), indices)), gd.dtable.tabletype)

length(gd::GDTable) = length(keys(gd.index))

iterate(gd::GDTable) = _iterate(gd, iterate(gd.index))
iterate(gd::GDTable, index_iter_state) = _iterate(gd, iterate(gd.index, index_iter_state))

function _iterate(gd::GDTable, it)
    if it !== nothing
        ((key, partition_indices), index_iter_state) = it
        return key => partition(gd, partition_indices), index_iter_state
    end
    return nothing
end

"""
    trim!(gd::GDTable) -> GDTable

Removes empty chunks from `gd` and unused keys from its index.
"""
function trim!(gd::GDTable)
    d = gd.dtable
    check_result = [Dagger.@spawn isnonempty(c) for c in d.chunks]
    results = fetch.(check_result)

    ok_indices = filter(x -> results[x], 1:length(results))
    d.chunks = getindex.(Ref(d.chunks), sort(ok_indices))

    offsets = zeros(UInt, length(results))

    counter = zero(UInt)
    for (i, r) in enumerate(results)
        counter = r ? counter : counter + 1
        offsets[i] = counter
    end

    for key in keys(gd.index)
        ind = gd.index[key]
        filter!(x -> results[x], ind)

        if isempty(ind)
            delete!(gd.index, key)
        else
            gd.index[key] = ind .- getindex.(Ref(offsets), ind)
        end
    end
    gd
end

"""
    trim(gd::GDTable) -> GDTable

Returns `gd` with empty chunks and keys removed.
"""
trim(gd::GDTable) = trim!(GDTable(DTable(gd), gd.cols, gd.index))

"""
    tabletype!(gd::GDTable)

Provides the type of the underlying table partition and caches it in `gd`.

In case the tabletype cannot be obtained the default return value is `NamedTuple`.
"""
tabletype!(gd::GDTable) = gd.dtable.tabletype = resolve_tabletype(gd.dtable)

"""
    tabletype(gd::GDTable)

Provides the type of the underlying table partition.
Uses the cached tabletype if available.

In case the tabletype cannot be obtained the default return value is `NamedTuple`.
"""
tabletype(gd::GDTable) = gd.dtable.tabletype === nothing ? resolve_tabletype(gd.dtable) : gd.dtable.tabletype

show(io::IO, gd::GDTable) = show(io, MIME"text/plain"(), gd)

function show(io::IO, ::MIME"text/plain", gd::GDTable)
    tabletype = gd.dtable.tabletype === nothing ? "unknown (use `tabletype!(::GDTable)`)" : gd.dtable.tabletype
    grouped_by_cols = gd.cols === nothing ? string(gd.grouping_function) : grouped_cols(gd)
    println(io, "GDTable with $(length(gd.dtable.chunks)) partitions and $(length(keys(gd.index))) keys")
    println(io, "Tabletype: $tabletype")
    print(io, "Grouped by: $grouped_by_cols")
    nothing
end

"""
    getindex(gdt::GDTable, key) -> DTable

Retrieves a `DTable` from `gdt` with rows belonging to the provided grouping key.
"""
function getindex(gdt::GDTable, key)
    ck = convert(keytype(gdt.index), key)
    ck ∉ keys(gdt) && throw(KeyError(ck))
    partition(gdt, ck)
end


text:
using Wavelets
using Test
using LinearAlgebra
using DelimitedFiles

# modified from Base.Test
function vecnorm_eq(va, vb, Eps, astr="a", bstr="b")
    if length(va) != length(vb)
        #error("lengths of ", astr, " and ", bstr, " do not match: ",
        #      "\n  ", astr, " (length $(length(va))) = ", va,
        #      "\n  ", bstr, " (length $(length(vb))) = ", vb)
        return false
    end

    diff = norm(va - vb)
    if !isnan(Eps) && diff > Eps
        #sdiff = string("|", astr, " - ", bstr, "| <= ", Eps)
        #error("assertion failed: ", sdiff,
        #      "\n  difference = ", diff, " > ", Eps)
        return false
    end
    return true
end
macro vecnorm_eq_eps(a, b, c)
    :(vecnorm_eq($(esc(a)), $(esc(b)), $(esc(c)), $(string(a)), $(string(b))))
end

@testset "Util" begin include("util.jl") end
@testset "Transforms" begin include("transforms.jl") end
@testset "Threshold" begin include("threshold.jl") end
@testset "Plot" begin include("plot.jl") end


text:
# SPDX-License-Identifier: MIT

"""
    abstract type AbstractParameter{T} end

Abstract type for parameters of a member of the exponential family.
"""
abstract type AbstractParameter{T} end

"""
    naturalform(param)

Returns the natural form of the parameter.

See also: [`realform`](@ref).
"""
naturalform(param::AbstractParameter)

"""
    realform(param)

Returns the vector of parameters as stored in `param`. Note that this
function is just an accessor of the internal storage of the parameter,
modifying the returned value should modify the parameter accordingly.

See also: [`naturalform`](@ref).
"""
realform(param::AbstractParameter)

"""
    jacobian(param)

Returns the Jacobian of ξ (the real form) w.r.t. the natural form η.
"""
jacobian(::AbstractParameter)

"""
    reallocate(param, bufferType)

Reallocate (i.e. copy) `param` with all interal buffers stored as
`bufferType`.
"""
reallocate(::AbstractParameter, bufferType::Type)

"""
    todict(param)

Return the object as a dictionary.

See also [`fromdict`](@ref).
"""
todict(param::AbstractParameter)

"""
    fromdict(T, dict)

Return a parameter object of type T from the dictionary provided.

See also [`todict`](@ref).
"""
fromdict(::Type{<:AbstractParameter}, dict::AbstractDict)

#######################################################################
# Default parameter: store the parameters in their natural form.

struct DefaultParameter{T} <: AbstractParameter{T}
    ξ::T
end

naturalform(p::DefaultParameter) = p.ξ
realform(p::DefaultParameter) = p.ξ
jacobian(p::DefaultParameter) = Diagonal(ones(eltype(p.ξ), length(p.ξ)))
reallocate(p::DefaultParameter, bufferType::Type) =
    DefaultParameter(bufferType(realform(p)))
todict(p::DefaultParameter) = Dict(:ξ => p.ξ)
fromdict(T::Type{<:DefaultParameter}, d::AbstractDict) = T(d[:ξ])

max_stars_repo_name: lucasondel/ExpFamilyDistributions-jl

max_issues_repo_name: BUTSpeechFIT/ExpFamilyDistributions

max_forks_repo_name: lucasondel/ExpFamilyDistributions