text:
using BCTRNN
using DiffEqSensitivity
using OrdinaryDiffEq
import DiffEqFlux: FastChain, FastDense
import Flux: ClipValue, ADAM

# Not in Project.toml
using Plots
gr()

include("sine_wave_dataloader.jl")

function train_sine_ss_fc(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32, model_size=5,
  mtkize=false, gen_jac=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToAll()
  sm = BCTRNN.SynsFullyConnected()
  #om = BCTRNN.OutputAll()
  om = BCTRNN.OutputIdxs(collect(1:model_size))
  wiring = BCTRNN.WiringConfig(f_in, model_size, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTCSynStateMTK(wiring, solver, sensealg; T, mtkize, gen_jac, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine_ss_fc(100, Vern7(); model_size=8, mtkize=false, lr=0.01,
)

function train_sine_ss_ncp(epochs, solver=Tsit5();
  sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP(true)),
  T=Float32,
  n_sensory=2, n_inter=5, n_command=5, n_motor=2,
  sensory_inter=2, inter_command=3, command_command=2, command_motor=2,
  mtkize=false, gen_jac=false, lr=0.02, kwargs...)

  train_dl = generate_2d_data(T)

  f_in = 2
  f_out = 1

  im = BCTRNN.InputAllToFirstN(n_sensory)
  sm = BCTRNN.SynsNCP(; n_sensory, n_inter, n_command, n_motor,
    sensory_inter, inter_command, command_command, command_motor
  )
  om = BCTRNN.OutputAll()
  om = BCTRNN.OutputIdxs(collect(1:n_sensory+n_inter+n_command+n_motor))
  wiring = BCTRNN.WiringConfig(f_in, im,sm,om)

  model = FastChain(BCTRNN.Mapper(f_in),
                    BCTRNN.LTCSynStateMTK(wiring, solver, sensealg; T, mtkize, gen_jac, kwargs...),
                    FastDense(wiring.n_out, f_out))

  hs = []
  for (k,v) in wiring.matrices
    push!(hs, heatmap(v, title=k))
  end
  display(plot(hs..., layout=length(hs)))

  cb = BCTRNN.MyCallback(T; ecb=mycb, nepochs=epochs, nsamples=length(train_dl))
  #opt = GalacticOptim.Flux.Optimiser(ClipValue(0.5), ADAM(0.02))
  opt = BCTRNN.ClampBoundOptim(BCTRNN.get_bounds(model,T)..., ClipValue(T(1.0)), ADAM(T(lr)))
  BCTRNN.optimize(model, BCTRNN.loss_seq, cb, opt, train_dl, epochs, T), model
end

@time train_sine_ss_ncp(100, ROCK4(), n_sensory=5, n_inter=10, n_command=10, mtkize=false, lr=0.01,
)

@time train_sine_ss_ncp(100, TRBDF2(linsolve=LinSolveGMRES()), n_sensory=3, n_inter=5, n_command=5, mtkize=false, lr=0.01,
)


text:
# test/runtests.jl
# Test UnitfulCurrency

using UnitfulCurrency

ext: jl


text:
function find_package(f, registry::String="")
    for reg in Pkg.Registry.reachable_registries()
        if isempty(registry) || registry == reg.name
            data = get_registry_file(reg, "Registry.toml")
            for (uuid, pkginfo) in data["packages"]
                f(uuid, pkginfo) && return (;uuid, reg, name=pkginfo["name"], path=pkginfo["path"])
            end
        end
    end
    return
end

function find_package(package::String, registry::String="")
    find_package(registry) do uuid, pkginfo
        pkginfo["name"] == package
    end
end

function find_max_version(package::String, registry::String="")
    info = find_package(package, registry)
    info === nothing && return

    versions = get_registry_file(info.reg, joinpath(info.path, "Versions.toml"))
    max_version = findmax(VersionNumber.(keys(versions)))[1]
    return max_version
end

function get_registry_file(reg, path::String)
    return if isnothing(reg.in_memory_registry)
        TOML.parsefile(joinpath(reg.path, path))
    else
        TOML.parse(reg.in_memory_registry[path])
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 680fdae0-1b3c-11ec-0de6-eb64d535f45f

# ╔═╡ 67f96cb0-1b3c-11ec-2f27-3bf1bcfaca4b

# ╔═╡ 67f945a0-1b3c-11ec-0b73-6bf3aba70083

# ╔═╡ 67cdefde-1b3c-11ec-218e-17e0fe9beb7c

# ╔═╡ Cell order:
# ╠═680fdae0-1b3c-11ec-0de6-eb64d535f45f
# ╠═67f96cb0-1b3c-11ec-2f27-3bf1bcfaca4b
# ╠═67f945a0-1b3c-11ec-0b73-6bf3aba70083
# ╠═67cdefde-1b3c-11ec-218e-17e0fe9beb7c


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

helpers for dense factorizations and linear solves
=#

import LinearAlgebra.BlasReal
import LinearAlgebra.BlasFloat
import LinearAlgebra.BlasInt
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra.LAPACK.liblapack
import LinearAlgebra.copytri!

# helpers for in-place matrix inverses (updates upper triangle only in some cases)

# Cholesky, BlasFloat
function inv_fact!(
    mat::Matrix{R},
    fact::Cholesky{R, Matrix{R}},
    ) where {R <: BlasFloat}
    copyto!(mat, fact.factors)
    LAPACK.potri!(fact.uplo, mat)
    return mat
end

# Cholesky, generic
function inv_fact!(
    mat::Matrix{R},
    fact::Cholesky{R, Matrix{R}},
    ) where {R <: RealOrComplex{<:Real}}
    # this is how Julia computes the inverse, but it could be implemented better
    copyto!(mat, I)
    ldiv!(fact, mat)
    return mat
end

# BunchKaufman, BlasReal
function inv_fact!(
    mat::Matrix{T},
    fact::BunchKaufman{T, Matrix{T}},
    ) where {T <: BlasReal}
    @assert fact.rook
    copyto!(mat, fact.LD)
    LAPACK.sytri_rook!(fact.uplo, mat, fact.ipiv),
    return mat
end

# LU, BlasReal
function inv_fact!(
    mat::Matrix{T},
    fact::LU{T, Matrix{T}},
    ) where {T <: BlasReal}
    copyto!(mat, fact.factors)
    LAPACK.getri!(mat, fact.ipiv)
    return mat
end

# LU, generic
function inv_fact!(
    mat::Matrix{T},
    fact::LU{T, Matrix{T}},
    ) where {T <: Real}
    # this is how Julia computes the inverse, but it could be implemented better
    copyto!(mat, I)
    ldiv!(fact, mat)
    return mat
end

# helpers for updating symmetric/Hermitian eigendecomposition

update_eigen!(X::Matrix{<:BlasFloat}) = LAPACK.syev!('V', 'U', X)[1]

function update_eigen!(X::Matrix{<:RealOrComplex{<:Real}})
    F = eigen(Hermitian(X, :U))
    copyto!(X, F.vectors)
    return F.values
end

# helpers for symmetric outer product (upper triangle only)
# B = alpha * A' * A + beta * B

outer_prod!(
    A::Matrix{T},
    B::Matrix{T},
    alpha::Real,
    beta::Real,
    ) where {T <: LinearAlgebra.BlasReal} =
    BLAS.syrk!('U', 'T', alpha, A, beta, B)

outer_prod!(
    A::AbstractMatrix{Complex{T}},
    B::AbstractMatrix{Complex{T}},
    alpha::Real,
    beta::Real,
    ) where {T <: LinearAlgebra.BlasReal} =
    BLAS.herk!('U', 'C', alpha, A, beta, B)

outer_prod!(
    A::AbstractMatrix{R},
    B::AbstractMatrix{R},
    alpha::Real,
    beta::Real,
    ) where {R <: RealOrComplex} =
    mul!(B, A', A, alpha, beta)

# ensure diagonal terms in square matrix are not too small
function increase_diag!(A::Matrix{T}) where {T <: Real}
    diag_pert = 1 + T(1e-5)
    diag_min = 1000 * eps(T)
    @inbounds for j in 1:size(A, 1)
        A[j, j] = diag_pert * max(A[j, j], diag_min)
    end
    return A
end

# helpers for spectral outer products

function spectral_outer!(
    mat::AbstractMatrix{T},
    vecs::Union{Matrix{T}, Adjoint{T, Matrix{T}}},
    diag::AbstractVector{T},
    temp::Matrix{T},
    ) where {T <: Real}
    mul!(temp, vecs, Diagonal(diag))
    mul!(mat, temp, vecs')
    return mat
end

function spectral_outer!(
    mat::AbstractMatrix{T},
    vecs::Union{Matrix{T}, Adjoint{T, Matrix{T}}},
    symm::Symmetric{T},
    temp::Matrix{T},
    ) where {T <: Real}
    mul!(temp, vecs, symm)
    mul!(mat, temp, vecs')
    return mat
end

#=
nonsymmetric square: LU
=#

function nonsymm_fact_copy!(
    mat2::Matrix{T},
    mat::Matrix{T},
    ) where {T <: Real}
    copyto!(mat2, mat)
    fact = lu!(mat2, check = false)

    if !issuccess(fact)
        copyto!(mat2, mat)
        increase_diag!(mat2)
        fact = lu!(mat2, check = false)
    end

    return fact
end

#=
symmetric indefinite: BunchKaufman (rook pivoting) and LU for generic fallback
NOTE if better fallback becomes available (eg dense LDL), use that
=#

symm_fact!(A::Symmetric{T, Matrix{T}}) where {T <: BlasReal} =
    bunchkaufman!(A, true, check = false)

symm_fact!(A::Symmetric{T, Matrix{T}}) where {T <: Real} =
    lu!(A, check = false)

function symm_fact_copy!(
    mat2::Symmetric{T, Matrix{T}},
    mat::Symmetric{T, Matrix{T}},
    ) where {T <: Real}
    copyto!(mat2, mat)
    fact = symm_fact!(mat2)

    if !issuccess(fact)
        copyto!(mat2, mat)
        increase_diag!(mat2.data)
        fact = symm_fact!(mat2)
    end

    return fact
end

#=
symmetric positive definite: unpivoted Cholesky
NOTE pivoted seems slower than BunchKaufman
=#

posdef_fact!(A::Symmetric{T, Matrix{T}}) where {T <: Real} =
    cholesky!(A, check = false)

function posdef_fact_copy!(
    mat2::Symmetric{T, Matrix{T}},
    mat::Symmetric{T, Matrix{T}},
    try_shift::Bool = true,
    ) where {T <: Real}
    copyto!(mat2, mat)
    fact = posdef_fact!(mat2)

    if !issuccess(fact)
        # try using symmetric factorization instead
        copyto!(mat2, mat)
        fact = symm_fact!(mat2)

        if try_shift && !issuccess(fact)
            copyto!(mat2, mat)
            increase_diag!(mat2.data)
            fact = symm_fact!(mat2)
        end
    end

    return fact
end


text:
# We probably won't use this. Prefer Distances.jl

function hellinger_distance!(p_counts::AbstractDict{T,V}, q_counts::AbstractDict{T,V}) where {T, V}
    p_sum = sum(values(p_counts))
    q_sum = sum(values(q_counts))
    total = float(zero(V))
    for (key, val) in p_counts
        if haskey(q_counts, key)
            total += (sqrt(val / p_sum) - sqrt(q_counts[key] / q_sum))^2
            delete!(q_counts, key)
        else
            total += val / p_sum
        end
    end
    total += sum(x -> x / q_sum, values(q_counts))
    dist = sqrt(total / 2)
    return dist
end

hellinger_distance(p_counts::AbstractDict, q_counts::AbstractDict) =
    hellinger_distance!(copy(p_counts), copy(q_counts))


text:
const hash_mask = typemax(UInt) >>> 0x01
const deletion_mask = hash_mask + 0x01

mutable struct Indices{I} <: AbstractIndices{I}
    # The hash table
    slots::Vector{Int}

    # Hashes and values
    hashes::Vector{UInt} # Deletion marker stored in high bit
    values::Vector{I}

    holes::Int # Number of "vacant" slots in hashes and values
end

_slots(inds::Indices) = getfield(inds, :slots)
_hashes(inds::Indices) = getfield(inds, :hashes)
_values(inds::Indices) = getfield(inds, :values)
_holes(inds::Indices) = getfield(inds, :holes)

Indices(; sizehint = 8) = Indices{Any}(; sizehint = sizehint)

function Indices{I}(; sizehint = 8) where {I}
    newsize = Base._tablesz((3 * sizehint) >> 0x01);
    Indices{I}(fill(0, newsize), Vector{UInt}(), Vector{I}(), 0)
end

"""
    Indices(iter)
    Indices{I}(iter)

Construct a `Indices` with indices from iterable container `iter`.

Note that the elements of `iter` must be distinct/unique. Instead, the `distinct` function
can be used for finding the unique elements.

# Examples

```julia
julia> Indices([1,2,3])
3-element Indices{Int64}
 1
 2
 3

julia> Indices([1,2,3,3])
ERROR: IndexError("Indices are not unique (inputs at positions 3 and 4) - consider using the distinct function")
Stacktrace:
 [1] Indices{Int64}(::Array{Int64,1}) at /home/ferris/.julia/dev/Dictionaries/src/Indices.jl:92
 [2] Indices(::Array{Int64,1}) at /home/ferris/.julia/dev/Dictionaries/src/Indices.jl:53
 [3] top-level scope at REPL[12]:1

julia> distinct([1,2,3,3])
3-element Indices{Int64}
 1
 2
 3
```
"""
function Indices(iter)
    if Base.IteratorEltype(iter) === Base.EltypeUnknown()
        iter = collect(iter)
    end

    return Indices{eltype(iter)}(iter)
end

function Indices{I}(iter) where {I}
    iter_size = Base.IteratorSize(iter)
    if iter_size isa Union{Base.HasLength, Base.HasShape}
        values = Vector{I}(undef, length(iter))
        @inbounds for (i, value) in enumerate(iter)
            values[i] = value
        end
    else
        values = Vector{I}()
        @inbounds for value in iter
            push!(values, value)
        end
    end
    return Indices{I}(values)
end

function Indices{I}(values::Vector{I}) where {I}
    # The input must have unique elements (the constructor is not to be used in place of `distinct`)
    hashes = map(v -> hash(v) & hash_mask, values)

    # Incrementally build the hashmap and throw if duplicates detected
    newsize = Base._tablesz(3*length(values) >> 0x01)
    bit_mask = newsize - 1 # newsize is a power of two
    slots = zeros(Int, newsize)
    @inbounds for index in keys(hashes)
        full_hash = hashes[index]
        trial_slot = reinterpret(Int, full_hash) & bit_mask
        @inbounds while true
            trial_slot = (trial_slot + 1)
            if slots[trial_slot] == 0
                slots[trial_slot] = index
                break
            else
                # TODO make this check optional
                if isequal(values[index], values[slots[trial_slot]])
                    throw(IndexError("Indices are not unique (inputs at positions $(slots[trial_slot]) and $index) - consider using the distinct function"))
                end
            end
            trial_slot = trial_slot & bit_mask
            # This is potentially an infinte loop and care must be taken not to overfill the container
        end
    end
    return Indices{I}(slots, hashes, values, 0)
end

Base.convert(::Type{AbstractIndices{I}}, inds::AbstractIndices) where {I} = convert(Indices{I}, inds) # the default AbstractIndices type
Base.convert(::Type{AbstractIndices{I}}, inds::AbstractIndices{I}) where {I} = inds

Base.convert(::Type{AbstractIndices{I}}, inds::Indices) where {I} = convert(Indices{I}, inds)
Base.convert(::Type{Indices}, inds::AbstractIndices{I}) where {I} = convert(Indices{I}, inds)

Base.convert(::Type{Indices{I}}, inds::Indices{I}) where {I} = inds
function Base.convert(::Type{Indices{I}}, inds::AbstractIndices) where {I}
    # Fast path
    if inds isa Indices && _holes(inds) == 0
        # Note: `convert` doesn't have copy semantics
        return Indices{I}(_slots(inds), _hashes(inds), convert(Vector{I}, _values(inds)), 0)
    end

    # The input is already unique
    values = collect(I, inds)
    hashes = map(v -> hash(v) & hash_mask, values)

    # Incrementally build the hashmap
    newsize = Base._tablesz(3*length(values) >> 0x01)
    bit_mask = newsize - 1 # newsize is a power of two
    slots = zeros(Int, newsize)
    @inbounds for index in keys(hashes)
        full_hash = hashes[index]
        trial_slot = reinterpret(Int, full_hash) & bit_mask
        @inbounds while true
            trial_slot = (trial_slot + 1)
            if slots[trial_slot] == 0
                slots[trial_slot] = index
                break
            else
                # TODO make this check optional
                if isequal(values[index], values[slots[trial_slot]])
                    throw(IndexError("Indices are not unique (inputs at positions $(slots[trial_slot]) and $index)"))
                end
            end
            trial_slot = trial_slot & bit_mask
            # This is potentially an infinte loop and care must be taken not to overfill the container
        end
    end
    return Indices{I}(slots, hashes, values, 0)
end

"""
    copy(inds::AbstractIndices)
    copy(inds::AbstractIndices, ::Type{I})

Create a shallow copy of the indices, optionally changing the element type.

(Note that `copy` on a dictionary does not copy its indices).
"""
function Base.copy(indices::Indices{I}, ::Type{I2}) where {I, I2}
    return Indices{I}(copy(_slots(indices)), copy(_hashes(indices)), collect(I2, _values(indices)), _holes(indices))
end

function Base.copy(indices::ReverseIndices{I,Indices{I}}, ::Type{I2}) where {I, I2}
    p = parent(indices)
    l = length(_values(p)) + 1
    old_slots = _slots(p)
    new_slots = similar(old_slots)
    @inbounds for i in keys(_slots(p))
        index = old_slots[i]
        if index > 0
            new_slots[i] = l - index
        else
            new_slots[i] = index
        end
    end
    return Indices{I2}(new_slots, reverse(_hashes(p)), collect(I2, Iterators.reverse(_values(p))), _holes(p))
end

# private (note that newsize must be power of two)
function rehash!(indices::Indices{I}, newsize::Int, values = (), include_last_values::Bool = true) where {I}
    slots = resize!(_slots(indices), newsize)
    fill!(slots, 0)
    bit_mask = newsize - 1 # newsize is a power of two

    if _holes(indices) == 0
        for (index, full_hash) in enumerate(_hashes(indices))
            trial_slot = reinterpret(Int, full_hash) & bit_mask
            @inbounds while true
                trial_slot = (trial_slot + 1)
                if slots[trial_slot] == 0
                    slots[trial_slot] = index
                    break
                else
                    trial_slot = trial_slot & bit_mask
                end
                # This is potentially an infinte loop and care must be taken not to overfill the container
            end
        end
    else
        # Compactify _values(indices), _hashes(indices) and the values while we are at it
        to_index = Ref(1) # Reassigning to to_index/from_index gives the closure capture boxing issue, so mutate a reference instead
        from_index = Ref(1)
        n_values = length(_values(indices))
        @inbounds while from_index[] <= n_values
            full_hash = _hashes(indices)[from_index[]]
            if full_hash & deletion_mask === zero(UInt)
                trial_slot = reinterpret(Int, full_hash) & bit_mask
                @inbounds while true
                    trial_slot = trial_slot + 1
                    if slots[trial_slot] == 0
                        slots[trial_slot] = to_index[]
                        _hashes(indices)[to_index[]] = _hashes(indices)[from_index[]]
                        _values(indices)[to_index[]] = _values(indices)[from_index[]]
                        if include_last_values || from_index[] < n_values
                            # Note - the last slot might end up with a random value (or
                            # GC'd reference). It's the callers responsibility to ensure the
                            # last slot is written to after this operation.
                            map(values) do (vals)
                                @inbounds vals[to_index[]] = vals[from_index[]]
                            end
                        end
                        to_index[] += 1
                        break
                    else
                        trial_slot = trial_slot & bit_mask
                    end
                end
            end

            from_index[] += 1
        end

        new_size = n_values - _holes(indices)
        resize!(_values(indices), new_size)
        resize!(_hashes(indices), new_size)
        map(values) do (vals)
            resize!(vals, new_size)
        end
        setfield!(indices, :holes, 0)
    end
    return indices
end

Base.length(indices::Indices) = length(_values(indices)) - _holes(indices)

# Token interface
istokenizable(::Indices) = true

tokentype(::Indices) = Int

# Duration iteration the token cannot be used for deletion - we do not worry about the slots
@propagate_inbounds function iteratetoken(indices::Indices)
    if _holes(indices) == 0
        return length(indices) > 0 ? ((0, 1), 1) : nothing
    end
    index = 1
    @inbounds while index <= length(_hashes(indices))
        if _hashes(indices)[index] & deletion_mask === zero(UInt)
            return ((0, index), index)
        end
        index += 1
    end
    return nothing
end

@propagate_inbounds function iteratetoken(indices::Indices, index::Int)
    index += 1
    if _holes(indices) == 0 # apparently this is enough to make it iterate as fast as `Vector`
        return index <= length(_values(indices)) ? ((0, index), index) : nothing
    end
    @inbounds while index <= length(_hashes(indices))
        if _hashes(indices)[index] & deletion_mask === zero(UInt)
            return ((0, index), index)
        end
        index += 1
    end
    return nothing
end

@propagate_inbounds function iteratetoken_reverse(indices::Indices)
    index = length(indices)
    if _holes(indices) == 0
        return index > 0 ? ((0, index), index) : nothing
    end
    @inbounds while index > 0
        if _hashes(indices)[index] & deletion_mask === zero(UInt)
            return ((0, index), index)
        end
        index -= 1
    end
    return nothing
end

@propagate_inbounds function iteratetoken_reverse(indices::Indices, index::Int)
    index -= 1
    if _holes(indices) == 0 # apparently this is enough to make it iterate as fast as `Vector`
        return index > 0 ? ((0, index), index) : nothing
    end
    @inbounds while index > 0
        if _hashes(indices)[index] & deletion_mask === zero(UInt)
            return ((0, index), index)
        end
        index -= 1
    end
    return nothing
end

function gettoken(indices::Indices{I}, i::I) where {I}
    full_hash = hash(i) & hash_mask
    n_slots = length(_slots(indices))
    bit_mask = n_slots - 1 # n_slots is always a power of two

    trial_slot = reinterpret(Int, full_hash) & bit_mask
    @inbounds while true
        trial_slot = (trial_slot + 1)
        trial_index = _slots(indices)[trial_slot]
        if trial_index > 0
            value = _values(indices)[trial_index]
            if i === value || isequal(i, value)
                return (true, (trial_slot, trial_index))
            end
        elseif trial_index === 0
            return (false, (0, 0))
        end

        trial_slot = trial_slot & bit_mask
        # This is potentially an infinte loop and care must be taken upon insertion not
        # to completely fill the container
    end
end

@propagate_inbounds function gettokenvalue(indices::Indices, (_slot, index))
    return _values(indices)[index]
end

@propagate_inbounds function gettokenvalue(indices::Indices, index::Int)
    return _values(indices)[index]
end

# Insertion interface
isinsertable(::Indices) = true

function gettoken!(indices::Indices{I}, i::I, values = ()) where {I}
    full_hash = hash(i) & hash_mask
    n_slots = length(_slots(indices))
    bit_mask = n_slots - 1 # n_slots is always a power of two
    n_values = length(_values(indices))

    trial_slot = reinterpret(Int, full_hash) & bit_mask
    trial_index = 0
    deleted_slot = 0
    @inbounds while true
        trial_slot = (trial_slot + 1)
        trial_index = _slots(indices)[trial_slot]
        if trial_index == 0
            break
        elseif trial_index < 0
            if deleted_slot == 0
                deleted_slot = trial_slot
            end
        else
            value = _values(indices)[trial_index]
            if i === value || isequal(i, value)
                return (true, (trial_slot, trial_index))
            end
        end

        trial_slot = trial_slot & bit_mask
        # This is potentially an infinte loop and care must be taken upon insertion not
        # to completely fill the container
    end

    new_index = n_values + 1
    if deleted_slot == 0
        # Use the trail slot
        _slots(indices)[trial_slot] = new_index
    else
        # Use the deleted slot
        _slots(indices)[deleted_slot] = new_index
    end
    push!(_hashes(indices), full_hash)
    push!(_values(indices), i)
    map(values) do (vals)
        resize!(vals, length(vals) + 1)
    end

    # Expand the hash map when it reaches 2/3rd full
    if 3 * new_index > 2 * n_slots
        # Grow faster for small hash maps than for large ones
        newsize = n_slots > 16000 ? 2 * n_slots : 4 * n_slots
        rehash!(indices, newsize, values, false)

        # The index has changed
        new_index = length(_values(indices))

        # The slot also has changed
        bit_mask = newsize - 1
        trial_slot = reinterpret(Int, full_hash) & bit_mask
        @inbounds while true
            trial_slot = (trial_slot + 1)
            if _slots(indices)[trial_slot] == new_index
                break
            end
            trial_slot = trial_slot & bit_mask
        end
    end

    return (false, (trial_slot, new_index))
end

@propagate_inbounds function deletetoken!(indices::Indices{I}, (slot, index), values = ()) where {I}
    @boundscheck if slot == 0
        error("Cannot use iteration token for deletion")
    end
    _slots(indices)[slot] = -index
    _hashes(indices)[index] = deletion_mask
    isbitstype(I) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), _values(indices), index-1)
    setfield!(indices, :holes, _holes(indices) + 1)

    # Recreate the hash map when 1/3rd of the values are deletions
    n_values = length(_values(indices)) - _holes(indices)
    if 3 * _holes(indices) > n_values
        # Halve if necessary
        n_slots = length(_slots(indices))
        halve = 4 * n_values < n_slots && n_slots > 8
        rehash!(indices, halve ? n_slots >> 0x01 : n_slots, values)
    end

    return indices
end

function Base.empty!(indices::Indices{I}, values = ()) where {I}
    setfield!(indices, :hashes, Vector{UInt}())
    setfield!(indices, :values, Vector{I}())
    setfield!(indices, :slots, fill(0, 8))
    setfield!(indices, :holes, 0)

    foreach(empty!, values)

    return indices
end

# Accelerated filtering

function Base.filter!(pred, indices::Indices)
    _filter!(token -> pred(@inbounds gettokenvalue(indices, token)), indices, ())
end

@inline function _filter!(pred, indices::Indices, values = ())
    indices_values = _values(indices)
    hashes = _hashes(indices)
    n = length(indices_values)
    i = Ref(0)
    j = Ref(0)
    @inbounds while i[] < n
        i[] += 1
        if hashes[i[]] & deletion_mask === zero(UInt) && pred(i[])
            j[] += 1
            indices_values[j[]] = indices_values[i[]]
            hashes[j[]] = hashes[i[]]
            map(vec -> @inbounds(vec[j[]] = vec[i[]]), values)
        end
    end
    newsize = j[]
    resize!(indices_values, newsize)
    resize!(hashes, newsize)
    map(vec -> resize!(vec, newsize), values)
    setfield!(indices, :holes, 0)
    newsize = Base._tablesz(3*length(_values(indices)) >> 0x01)
    rehash!(indices, newsize, values)
end

# Factories

# TODO make this generic... maybe a type-based `empty`?
function _distinct(::Type{Indices}, itr)
    if Base.IteratorEltype(itr) === Base.HasEltype()
        return _distinct(Indices{eltype(itr)}, itr)
    end

    tmp = iterate(itr)
    if tmp === nothing
        return Indices{Base.@default_eltype(itr)}()
    end
    (x, s) = tmp
    indices = Indices{typeof(x)}()
    insert!(indices, x)
    return __distinct!(indices, itr, s, x)
end

# An auto-widening constructor for insertable AbstractIndices
function __distinct!(indices::AbstractIndices, itr, s, x_old)
    T = eltype(indices)
    tmp = iterate(itr, s)
    while tmp !== nothing
        (x, s) = tmp
        if !isequal(x, x_old) # Optimized for repeating elements of `itr`, e.g. if `itr` is sorted
            if !(x isa T) && promote_type(typeof(x), T) != T
                new_indices = copy(indices, promote_type(T, typeof(x)))
                set!(new_indices, x)
                return __distinct!(new_indices, itr, s, x)
            end
            set!(indices, x)
            x_old = x
        end
        tmp = iterate(itr, s)
    end
    return indices
end

function randtoken(rng::Random.AbstractRNG, inds::Indices)
    if inds.holes === 0
        return (0, rand(rng, Base.OneTo(length(inds))))
    end

    # Rejection sampling to handle deleted tokens (which are sparse)
    range = Base.OneTo(length(_hashes(inds)))
    while true
        i = rand(rng, range)
        if inds.hashes[i] !== deletion_mask
            return (0, i)
        end
    end
end


text:
using FileIO
using JLD2
using Pkg.Artifacts
using ClimateMachine.ArtifactWrappers

# Get bomex_edmf dataset folder:
bomex_edmf_dataset = ArtifactWrapper(
    joinpath(@__DIR__, "Artifacts.toml"),
    "bomex_edmf",
    ArtifactFile[ArtifactFile(
        url = "https://caltech.box.com/shared/static/jbhcy6ncc5wh1hg9hcea5f45w6t22kk2.jld2",
        filename = "bomex_edmf.jld2",
    ),],
)
bomex_edmf_dataset_path = get_data_folder(bomex_edmf_dataset)
data_file = joinpath(bomex_edmf_dataset_path, "bomex_edmf.jld2")

updraft_vars(N_up, st, tol, sym...) = Dict(ntuple(N_up) do i
    "turbconv.updraft[$i]." * join(string.(sym), ".") => (st, tol)
end)

vars_to_compare(N_up) = Dict(
    "ρ" => (Prognostic(), 1),
    "ρu[1]" => (Prognostic(), 1),
    "ρu[2]" => (Prognostic(), 1),
    "ρu[3]" => (Prognostic(), 1),
    "ρe" => (Prognostic(), 1e5),
    "moisture.ρq_tot" => (Prognostic(), 1e-2),
    "turbconv.environment.ρatke" => (Prognostic(), 5e-1),
    "turbconv.environment.T" => (Auxiliary(), 300),
    "turbconv.environment.cld_frac" => (Auxiliary(), 1),
    "turbconv.environment.buoyancy" => (Auxiliary(), 1e-2),
    updraft_vars(N_up, Prognostic(), 1 * 0.1, :ρa)...,
    updraft_vars(N_up, Prognostic(), 1 * 0.1, :ρaw)...,
    updraft_vars(N_up, Prognostic(), 300 * 0.1, :ρaθ_liq)...,
    updraft_vars(N_up, Prognostic(), 1e-2 * 0.1, :ρaq_tot)...,
    updraft_vars(N_up, Auxiliary(), 1e-2, :buoyancy)...,
)

compare = Dict()
@testset "Regression Test" begin
    N_up = n_updrafts(solver_config.dg.balance_law.turbconv)
    numerical_data =
        dict_of_nodal_states(solver_config, ["z"], (Prognostic(), Auxiliary()))
    data_to_compare = Dict()
    for (ftc, v) in vars_to_compare(N_up)
        data_to_compare[ftc] = numerical_data[ftc]
    end

    export_new_solution_jld2 = false
    if export_new_solution_jld2
        save("bomex_edmf.jld2", data_to_compare)
    end

    all_data_ref = load(data_file)

    @test all(k in keys(all_data_ref) for k in keys(data_to_compare))
    N_up = n_updrafts(solver_config.dg.balance_law.turbconv)
    comparison_vars = vars_to_compare(N_up)
    for k in keys(all_data_ref)
        data = data_to_compare[k]
        ref_data = all_data_ref[k]
        tol = comparison_vars[k][2]
        s = length(data) / 100

        @test !any(isnan.(ref_data))
        @test !any(isnan.(data))
        absΔdata = abs.(data .- ref_data)
        T1 = isapprox(norm(absΔdata), 0, atol = tol * 0.01 * s) # norm
        T2 = isapprox(maximum(absΔdata), 0, atol = tol * 0.01) # max of local
        compare[k] = (norm(absΔdata), maximum(absΔdata), tol)
        (!T1 || !T2) && @show k, norm(absΔdata), maximum(absΔdata), tol
        @test isapprox(norm(absΔdata), 0, atol = tol * 0.01 * s) # norm
        @test isapprox(maximum(absΔdata), 0, atol = tol * 0.01) # max of local
    end
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
struct Subsequence{VT<:AbstractVector,T}
    c::VT
    fillvalue::T
    winlen::Int
    noverlap::Int
    step::Int
    lpadlen::Int
    rpadlen::Int
end

function getpadlen(winlen)
    if mod(winlen, 2) == 0
        lpadlen = (winlen-1)÷2
        rpadlen = winlen÷2
    else
        lpadlen = rpadlen = winlen÷2
    end
    lpadlen, rpadlen
end

function Subsequence(c::AbstractVector{T}, winlen, noverlap; fillvalue=zero(T)) where T
    winlen > length(c) && throw(ArgumentError("`winlen` has to be smaller than the signal length."))
    step = winlen-noverlap
    lpadlen, rpadlen = getpadlen(winlen)
    Subsequence(c,
                fillvalue,
                winlen,
                noverlap,
                step,
                lpadlen,
                rpadlen)
    # Subsequence(PaddedView(fill, c, (1-lpadlen:length(c)+rpadlen,)), winlen, noverlap, step, lpadlen, rpadlen)
end

function Base.iterate(subseq::Subsequence, state=1)
    lenc = length(subseq.c)
    state > lenc && return nothing
    # return @view(subseq.c[state-subseq.lpadlen:state+subseq.rpadlen]), state+subseq.step
    if state <= subseq.lpadlen
        return vcat(fill(subseq.fillvalue, subseq.lpadlen-(state-1)), @view(subseq.c[1:state+subseq.rpadlen])), state+subseq.step
    elseif state >= lenc-subseq.rpadlen
        return vcat(@view(subseq.c[state-subseq.lpadlen:end]), fill(subseq.fillvalue, subseq.rpadlen-(lenc-state))), state+subseq.step
    else
        return @view(subseq.c[state-subseq.lpadlen:state+subseq.rpadlen]), state+subseq.step
    end
end

Base.length(subseq::Subsequence) = ceil(Int64, length(subseq.c)/subseq.step)

Base.getindex(subseq::Subsequence, i::Number) = iterate(subseq, (1:subseq.step:length(subseq.c))[i])[1]


text:
module QJuliaGaugeUtils

import QJuliaInterface
import QJuliaEnums
import QUDARoutines

debug_constructGaugeField = false

@inline function accumulateConjugateProduct(a::Complex{T}, b::Complex{T}, c::Complex{T}, sign::Float64) where T <: AbstractFloat
   local tmp::Complex{T} = b*c
   a += Complex{T}(sign*real(tmp), -sign*imag(tmp))
   return a
end

# normalize the vector a
function normalize(a::AbstractArray, len::Int)
  local sum::Float64 = 0.0;
  for i in 1:len; sum += abs2(a[i]);  end
  a[:] /= sqrt(sum)
end

# orthogonalize vector b to vector a
function orthogonalize(a::AbstractArray, b::AbstractArray, len::Int)
  local dot::Complex{Float64} = 0.0;
  for i in 1:len; dot  += conj(a[i])*b[i]; end
  b[:] -= (dot*a[:])
end

#Main methods:

function applyGaugeFieldScaling!(gauge::Matrix{Complex{T}}, param::QJuliaInterface.QJuliaGaugeParam_qj) where T <: AbstractFloat

  vol    = param.X[1]*param.X[2]*param.X[3]*param.X[4]
  volh   = Int(vol / 2)
  volh_t = Int(param.X[1] / 2)*param.X[2]*param.X[3]*(param.X[4]-1)

  # Apply spatial scaling factor (u0) to spatial links
  gauge[:, 1:3] /= param.anisotropy;

  # only apply T-boundary at edge nodes (always true for the single device)
  local last_node_in_t = (QUDARoutines.commCoords_qj(3) == QUDARoutines.commDim_qj(3)-1) ? true : false

  # create time direction views:
  even_tlinks = view(gauge, (9volh_t+1):9volh, 4)
  odd_tlinks  = view(gauge, 9(volh+volh_t)+1:9vol, 4)

  # Apply boundary conditions to temporal links
  if param.t_boundary == QJuliaEnums.QJULIA_ANTI_PERIODIC_T && last_node_in_t
    println("Applying antiperiodic BC.")
    even_tlinks[:] *= -1.0
    odd_tlinks[:]  *= -1.0
  end

  if param.gauge_fix == QJuliaEnums.QJULIA_GAUGE_FIXED_YES
    println("Applying gauge fixing.")
    # set all gauge links (except for the last X[1]*X[2]*X[3]/2) to the identity,
    # to simulate fixing to the temporal gauge.
    local iMax = last_node_in_t ? volh_t : volh

    even_gauge = view(gauge, 1:9iMax, 4)
    odd_gauge  = view(gauge, 9volh+1:9(volh+iMax), 4)

    for i in 0:(iMax-1)
      for m in 0:2
	for n in 0:2
	  even_gauge[9i + 3m + n + 1] = (m==n) ? 1.0 : 0.0;
	  odd_gauge[ 9i + 3m + n + 1] = (m==n) ? 1.0 : 0.0;
	end # for n
      end # for m
    end #for i

  end # param.gauge_fix

end #applyGaugeFieldScaling!

function constructUnitGaugeField!(gauge::Matrix{Complex{T}}, param::QJuliaInterface.QJuliaGaugeParam_qj) where T <: AbstractFloat
  vol  = param.X[1]*param.X[2]*param.X[3]*param.X[4]
  volh = Int(vol / 2)

  even_gauge = view(gauge, 1:9volh, :)
  odd_gauge  = view(gauge, 9volh+1:9vol, :)

  for d in 1:4
    for i in 0:(volh-1)
      for m in 0:2
        for n in 0:2
          even_gauge[9i + 3m + n + 1, d] = (m==n) ? 1.0 : 0.0;
          odd_gauge[9i + 3m + n + 1 , d] = (m==n) ? 1.0 : 0.0;
        end # for n
      end # for m
    end # for i
  end # for d

  applyGaugeFieldScaling!(gauge, param)

end # constructUnitGaugeField!

function constructGaugeField!(gauge::Matrix{Complex{T}}, param::QJuliaInterface.QJuliaGaugeParam_qj) where T <: AbstractFloat

  vol  = param.X[1]*param.X[2]*param.X[3]*param.X[4]
  volh = Int(vol / 2)

  println("Gauge field volume: ",  vol)

  even_gauge = view(gauge, 1:9volh, :)
  odd_gauge  = view(gauge, 9volh+1:9vol, :)

  #println(typeof(even_gauge))
  #println(typeof(odd_gauge) )
  #we do always need offset = 1
  for d in 1:4
    for i in 0:(volh-1)
      for m in 0:2
        for n in 0:2
          even_gauge[9i + 3m + n + 1, d] = Complex{T}(rand(), rand())
          odd_gauge[9i + 3m + n + 1 , d] = Complex{T}(rand(), rand())
        end
      end
      #create a view for a given link
      local c = 1;
      gauge_link_col0 = view(even_gauge, (9i       +1):(9i+3c    ), d)
      gauge_link_col1 = view(even_gauge, (9i+3c    +1):(9i+3(c+1)), d)
      gauge_link_col2 = view(even_gauge, (9i+3(c+1)+1):(9i+3(c+2)), d)

      normalize(gauge_link_col1, 3)
      orthogonalize(gauge_link_col1, gauge_link_col2, 3)
      normalize(gauge_link_col2, 3)

      for i in 1:3; gauge_link_col0[i] = 0.0; end

      gauge_link_col0[1] = accumulateConjugateProduct(gauge_link_col0[1], gauge_link_col1[2], gauge_link_col2[3], +1.0)
      gauge_link_col0[1] = accumulateConjugateProduct(gauge_link_col0[1], gauge_link_col1[3], gauge_link_col2[2], -1.0)
      gauge_link_col0[2] = accumulateConjugateProduct(gauge_link_col0[2], gauge_link_col1[3], gauge_link_col2[1], +1.0)
      gauge_link_col0[2] = accumulateConjugateProduct(gauge_link_col0[2], gauge_link_col1[1], gauge_link_col2[3], -1.0)
      gauge_link_col0[3] = accumulateConjugateProduct(gauge_link_col0[3], gauge_link_col1[1], gauge_link_col2[2], +1.0)
      gauge_link_col0[3] = accumulateConjugateProduct(gauge_link_col0[3], gauge_link_col1[2], gauge_link_col2[1], -1.0)

      gauge_link_col0 = view(odd_gauge, (9i       +1):(9i+3c    ), d)
      gauge_link_col1 = view(odd_gauge, (9i+3c+1):(9i+3(c+1)), d)
      gauge_link_col2 = view(odd_gauge, (9i+3(c+1)+1):(9i+3(c+2)), d)

      normalize(gauge_link_col1, 3)
      orthogonalize(gauge_link_col1, gauge_link_col2, 3)
      normalize(gauge_link_col2, 3)

      for i in 1:3; gauge_link_col0[i] = 0.0; end

      gauge_link_col0[1] = accumulateConjugateProduct(gauge_link_col0[1], gauge_link_col1[2], gauge_link_col2[3], +1.0)
      gauge_link_col0[1] = accumulateConjugateProduct(gauge_link_col0[1], gauge_link_col1[3], gauge_link_col2[2], -1.0)
      gauge_link_col0[2] = accumulateConjugateProduct(gauge_link_col0[2], gauge_link_col1[3], gauge_link_col2[1], +1.0)
      gauge_link_col0[2] = accumulateConjugateProduct(gauge_link_col0[2], gauge_link_col1[1], gauge_link_col2[3], -1.0)
      gauge_link_col0[3] = accumulateConjugateProduct(gauge_link_col0[3], gauge_link_col1[1], gauge_link_col2[2], +1.0)
      gauge_link_col0[3] = accumulateConjugateProduct(gauge_link_col0[3], gauge_link_col1[2], gauge_link_col2[1], -1.0)

    end # for i
  end # for d

  if param.gtype == QJuliaEnums.QJULIA_WILSON_LINKS
    println("Applying scaling/BC on the gauge links")
    applyGaugeFieldScaling!(gauge, param)
  end

  if debug_constructGaugeField == true
    println("::> Begin debug info for function constructGaugeField...")
    dir = 1
    for i in 1:16
      println("Cehck value for index ", i, ",dir " , dir, " complex value is = ", gauge[i, dir])
    end
    println("<:: End debug info.")
  end

end #constructGaugeField!

function construct_gauge_field!(gauge::Matrix{Complex{T}}, gtype, param::QJuliaInterface.QJuliaGaugeParam_qj) where T <: AbstractFloat

  println("Working with gauge field type:",  typeof(gauge), " : ", length(gauge))

  if gtype == 0
    println("Construct unit gauge field")
    constructUnitGaugeField!(gauge, param)
  elseif gtype == 1
    println("Construct random gauge field")
    constructGaugeField!(gauge, param)
  else
    println("Apply scaling...")
    applyGaugeFieldScaling!(gauge, param)
  end
end #construct_gauge_field!

end #QJuliaGaugeUtils

max_stars_repo_name: alexstrel/QJulia

max_issues_repo_name: alexstrel/QJulia

max_forks_repo_name: alexstrel/QJulia


text:
# Methods that deal with interpolating the adaptive kernel's kernel parameters
#   or the warp map from samples.

# # One single interpolator for the size of the entire image.
# function getθset{KT}(itp, itp_set, θ::KT, scale_factor::Float64,
#                 ratio::Vector{Float64},
#                 offset::Vector{Float64})::Matrix{AdaptiveKernelType{KT}}
#
#     θ_set = Array{AdaptiveKernelType{KT}}(size(itp_set))
#
#     for i = 1:length(θ_set)
#         θ_set[i] = AdaptiveKernelType(θ, x->scale_factor*itp[(x./ratio+offset)...] )
#     end
#
#     return θ_set
# end
#
# # The version that uses an array of smaller interpolators, with oversized phi patches.
# function getθsetlocalinterpolators2{KT}(itp_sets::Vector,
#                 θ::KT, scale_factor::Float64,
#                 ratio::Vector{Float64},
#                 X_set::Matrix{Matrix{Vector{Float64}}},
#                 overlap::Int)::Matrix{AdaptiveKernelType{KT}}
#
#     N_rows, N_cols = size(itp_sets[1])
#     θ_set = Array{AdaptiveKernelType{KT}}(N_rows,N_cols)
#
#     M = length(itp_sets)
#
#     # Top left. Put into for-loop block so anchor and offset gets its own scope.
#     for i = 1:1
#         anchor = [1; 1]
#         offset = -X_set[1,1][1]./ratio + anchor
#
#         warpfunc_set = Vector{Function}(M)
#         for m = 1:M
#             itp_set = itp_sets[m]
#             warpfunc_set[m] = x->scale_factor*itp_set[1,1][(x./ratio+offset)...]
#         end
#
#         θ_set[1,1] = AdaptiveKernelType(θ,warpfunc_set)
#     end
#
#     # Left column, second row until end.
#     for i = 2:N_rows
#         anchor = [overlap+1; 1]
#         offset = -X_set[i,1][1]./ratio + anchor
#
#         warpfunc_set = Vector{Function}(M)
#         for m = 1:M
#             itp_set = itp_sets[m]
#             warpfunc_set[m] = x->scale_factor*itp_set[i,1][(x./ratio+offset)...]
#         end
#
#         θ_set[i,1] = AdaptiveKernelType(θ,warpfunc_set)
#     end
#
#     # Top row, second column until end.
#     for j = 2:N_cols
#         anchor = [1; overlap+1]
#         offset = -X_set[1,j][1]./ratio + anchor
#
#         warpfunc_set = Vector{Function}(M)
#         for m = 1:M
#             itp_set = itp_sets[m]
#             warpfunc_set[m] = x->scale_factor*itp_set[1,j][(x./ratio+offset)...]
#         end
#
#         θ_set[1,j] = AdaptiveKernelType(θ,warpfunc_set)
#     end
#
#     # Central regions, second row ⊗ second column until end.
#     for i = 2:N_rows
#         for j = 2:N_cols
#             anchor = [overlap+1; overlap+1]
#             offset = -X_set[i,j][1]./ratio + anchor
#
#             warpfunc_set = Vector{Function}(M)
#             for m = 1:M
#                 itp_set = itp_sets[m]
#                 warpfunc_set[m] = x->scale_factor*itp_set[i,j][(x./ratio+offset)...]
#             end
#
#             θ_set[i,j] = AdaptiveKernelType(θ,warpfunc_set)
#         end
#     end
#
#     return θ_set
# end

# options for warp map samples.
# this is the case when the target density is unknown, but realizations are available.
function getwarpmapsamplecustom( y::Array{T,D},
                            ω_set,
                            pass_band_factor) where {T,D}
    #
    N_bands = length(ω_set)

    Y = y

    #### Split-band analysis.
    ϕY, ψY = SignalTools.runsplitbandanalysis(Y, ω_set, SignalTools.getGaussianfilters)
    ηY = SignalTools.runbandpassanalysis(Y, ω_set, pass_band_factor, SignalTools.getGaussianfilters)

    # #### Riesz transform on the different filtered signals.
    # H, ordering = gethigherorderRTfilters(Y,order)
    #
    # 𝓡ϕY = collect( RieszAnalysisLimited(ϕY[s],H) for s = 1:N_bands)
    # 𝓡ψY = collect( RieszAnalysisLimited(ψY[s],H) for s = 1:N_bands)
    # 𝓡ηY = collect( RieszAnalysisLimited(ηY[s],H) for s = 1:N_bands)

    ϕ_set = ηY

    ϕ = reduce(+,ϕ_set)./N_bands

    return ϕY, ψY, ηY
end

function getwarpmaplinear(ϕ::Array{T,D}) where {T,D}

    itp_ϕ = Interpolations.interpolate(ϕ,
                Interpolations.BSpline(Interpolations.Linear()))

    etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())

    return etp_ϕ
end

function getwarpmap(ϕ::Array{T,D}) where {T,D}

    itp_ϕ = Interpolations.interpolate(ϕ,
                Interpolations.BSpline(Interpolations.Cubic(
                    Interpolations.Flat(Interpolations.OnGrid()))))

    etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())
    #etp_ϕ = Interpolations.extrapolate(itp_ϕ, 0)

    return etp_ϕ
end

function getwarpmap(ϕ::Array{T,D}, x_ranges::Vector, amplification_factor::T) where {T,D}

    @assert D == length(x_ranges)
    N_array = collect( length(x_ranges[d]) for d = 1:D )

    itp_ϕ = Interpolations.interpolate(ϕ,
                Interpolations.BSpline(Interpolations.Cubic(
                    Interpolations.Flat(Interpolations.OnGrid()))))

    etp_ϕ = Interpolations.extrapolate(itp_ϕ, Interpolations.Line())
    #etp_ϕ = Interpolations.extrapolate(itp_ϕ, 0)

    st = collect( x_ranges[d][1] for d = 1:D )
    fin = collect( x_ranges[d][end] for d = 1:D )
    f = xx->etp_ϕ(interval2itpindex(xx,
                st,
                fin,
                N_array)...)*amplification_factor

    # chain rule for first derivatives.
    df = xx->( Interpolations.gradient(etp_ϕ,interval2itpindex(xx,
                st,
                fin,
                N_array)...) .* derivativeinterval2itpindex(st,fin,N_array) .*amplification_factor )

    # chain rule for second derivatives.
    d2f = xx->( Interpolations.hessian(etp_ϕ,interval2itpindex(xx,
                st,
                fin,
                N_array)...) .* (derivativeinterval2itpindex(st,fin,N_array).^2) .*amplification_factor )

    return f, df, d2f
end
# # test code.
#
# A = randn(N,N)
# itp_ϕ = Interpolations.interpolate(A,
#             Interpolations.BSpline(Interpolations.Cubic(
#                 Interpolations.Flat(Interpolations.OnGrid()))))
#
# P = [1.1; 2.3]
# ϕ_map_func, d_ϕ_map_func = getwarpmap(A, x_ranges[1:2], amplification_factor)
#
#
# dϕ_ND = xx->Calculus.gradient(ϕ_map_func, xx)
# dϕ_ND(P)
#
# d_ϕ_map_func(P) ./ dϕ_ND(P)

max_stars_repo_name: RoyCCWang/PatchMixtureKriging

max_issues_repo_name: RoyCCWang/PatchMixtureKriging


text:
function get_combined_provider(::Type{TransformProviderUtils}, arg0::TransformProvider, arg1::TransformProvider)
    return jcall(TransformProviderUtils, "getCombinedProvider", TransformProvider, (TransformProvider, TransformProvider), arg0, arg1)
end

function get_reversed_provider(::Type{TransformProviderUtils}, arg0::TransformProvider)
    return jcall(TransformProviderUtils, "getReversedProvider", TransformProvider, (TransformProvider,), arg0)
end


text:
#///////////////////////////////////////
#// File Name: rs_sac_controller_test.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2021/02/26
#// Description: Test code for src/rs_sac_controller.jl
#///////////////////////////////////////

using DataStructures
using LinearAlgebra
using Random
using RobotOS

@testset "RSSAC Controller Test" begin
# Cost parameters
Cep = Matrix(1.0I, 2, 2);
Cu = Matrix(1.0I, 2, 2);
β_pos = 0.6;
β_col = 0.4;
α_col = 100.0;
λ_col = 1.0;
σ_risk = 1.0;
cost_param = CostParameter(Cep, Cu, β_pos, α_col, β_col, λ_col, σ_risk);
# Simulation parameters
conf_file_name = "config.json";
test_data_name = "eth_test.pkl";
test_scene_id = 0;
start_time_idx = 50;
device = py"torch".device("cpu");
incl_robot_name = false;
scene_param = TrajectronSceneParameter(conf_file_name, test_data_name,
                                       test_scene_id, start_time_idx,
                                       incl_robot_name);
scene_loader = TrajectronSceneLoader(scene_param, verbose=false);
num_samples = 50;
prediction_steps = 12;
use_robot_future = false;
deterministic = false;
rng_seed_py = 1;
predictor_param = TrajectronPredictorParameter(prediction_steps, num_samples,
                                               use_robot_future, deterministic,
                                               rng_seed_py);
traj_predictor = TrajectronPredictor(predictor_param,
                                     scene_loader.model_dir,
                                     scene_loader.param.conf_file_name,
                                     device, verbose=false);
initialize_scene_graph!(traj_predictor, scene_loader);
ado_states = fetch_ado_positions!(scene_loader, return_full_state=true);
ado_positions = reduce_to_positions(ado_states);

dtc = 0.01;
sim_param = SimulationParameter(scene_loader, traj_predictor, dtc, cost_param);
# Initial conditions
ep = [-1., 3.5];
ev = [1., 1.];
t_init = Time(1.0);
e_init = RobotState([ep; ev], t_init);
ap_dict_init = convert_nodes_to_str(ado_positions);
t_last_m = Time(0.8);
w_init = WorldState(e_init, ap_dict_init, t_last_m);
# Target Trajectory
u_nominal_base = [0.0, 0.0];
u_array_base = [u_nominal_base for ii = 1:480];
wp1 = WayPoint2D([0.0, 0.0], Time(1.0));
wp2 = WayPoint2D([0.0, 0.0], Time(5, 8e8));
target_trajectory = Trajectory2D([wp1, wp2]);
# ControlParameter
eamax = 5.0;
tcalc = 0.1;
dtexec = [0.05, 0.01, 0.02];
dtr = 0.1;
nominal_search_depth = 2;
constraint_time = nothing;
u_nominal_cand = append!([[0.0, 0.0]],
                         [round.([a*cos(deg2rad(θ)), a*sin(deg2rad(θ))], digits=5)
                          for a = [2., 4.] for θ = 0.:45.:(360. - 45.)])             # nominal control candidate value [ax, ay] [m/s^2]
cnt_param = ControlParameter(eamax, tcalc, dtexec, dtr, u_nominal_base,
                             u_nominal_cand, nominal_search_depth,
                             constraint_time=constraint_time);
# Prediction Dict
prediction_dict = sample_future_ado_positions!(traj_predictor,
                                               ado_states);
num_controls = length(cnt_param.u_nominal_cand)^nominal_search_depth;
for key in keys(prediction_dict)
    prediction_dict[key] = repeat(prediction_dict[key], outer=(num_controls, 1, 1));
end
@test prediction_dict["PEDESTRIAN/25"][1:50, :, :] == prediction_dict["PEDESTRIAN/25"][51:100, :, :]

# # Helper functions test
# convert_to_schedule test
u_schedule = convert_to_schedule(w_init.t, u_array_base, sim_param);

@test length(u_schedule) == 480
@test all(collect(values(u_schedule)) .== [[0.0, 0.0]])
@test haskey(u_schedule, Time(1.0))
@test haskey(u_schedule, Time(1.01))
@test haskey(u_schedule, Time(5.79))

# get_nominal_u_arrays test: get nominal u_arrays from u_schedule and nominal control candidates
#u_nominal_mod_init_time = w_init.t + Duration(cnt_param.tcalc);
#u_nominal_mod_final_time = u_nominal_mod_init_time +
#                           Duration(sim_param.dto) -
#                           Duration(sim_param.dtc);
schedule_before = u_schedule;
u_arrays = get_nominal_u_arrays(u_schedule, sim_param, cnt_param);
@test all(u_arrays[1][1:end] .== [[0.0, 0.0]])
@test all([all(u_arrays[ii][1:10] .== [[0.0, 0.0]]) &&
           all(u_arrays[ii][11:50] .== [u_nominal_cand[div(ii - 1, length(cnt_param.u_nominal_cand)) + 1]]) &&
           all(u_arrays[ii][51:90] .== [u_nominal_cand[mod(ii - 1, length(cnt_param.u_nominal_cand)) + 1]]) &&
           all(u_arrays[ii][91:end] .== [[0.0, 0.0]])
           for ii = 2:length(cnt_param.u_nominal_cand)])
@test u_schedule == schedule_before

# get_robot_present_and_future_test
begin
    e_state = RobotState([-5.0, 0.0, 0.0, 1.0]);
    u_nominal = [[0.0, 0.0] for ii = 1:Int64(round(prediction_steps*sim_param.dto/dtc))];
    u_s = convert_to_schedule(e_state.t, u_nominal, sim_param);
    rpf = get_robot_present_and_future(e_state, u_s, sim_param, cnt_param);
    @test size(rpf) == (17^nominal_search_depth, 13, 6)
    @test all(isapprox.(cumsum(rpf[:, 1:end-1, 5:6], dims=2).*sim_param.dto .+ rpf[:, 1:1, 3:4], rpf[:, 2:end, 3:4], atol=1e-6))

end

# Get simulation result and best nominal u_array
sim_result, u_nominal_array = simulate(w_init, u_arrays, target_trajectory,
                                       prediction_dict, sim_param, cnt_param.constraint_time);

# get_control_coeffs_test
coeff_matrix, coeff_matrix_constraint = get_control_coeffs(sim_result, sim_param, cnt_param);
#@test size(coeff_matrix) == (2, 20);
@test size(coeff_matrix) == (2, 480);
@test isnothing(coeff_matrix_constraint)

test_time_id = 11;
den_test = sum(exp.(σ_risk.*sim_result.sampled_total_costs));
num_test = sum(exp.(σ_risk.*sim_result.sampled_total_costs).*
               [sim_result.e_costate_array[:, test_time_id, jj] for jj = 1:sim_param.num_samples]);

@test transition_control_coeff(sim_result.e_state_array[test_time_id])'*
    (num_test./den_test) ≈ coeff_matrix[:, test_time_id];

# get_control_schedule test
control_schedule_array = get_control_schedule(sim_result, u_nominal_array, coeff_matrix,
                                              sim_param, cnt_param);
#@test length(control_schedule_array) == Int64((tcalc+dtr)/dtc);
@test length(control_schedule_array) == Int64(round(sim_param.dto*sim_param.prediction_steps/dtc, digits=5));
@test maximum(map(s -> norm(vec(s.u)), control_schedule_array)) < eamax ||
      maximum(map(s -> norm(vec(s.u)), control_schedule_array)) ≈ eamax;
@test !any(isnan.(map(s -> s.cost, control_schedule_array)));

# determine_control_time test
control_chosen = determine_control_time(sim_result,
                                        control_schedule_array, sim_param,
                                        cnt_param);
@test t_init + Duration(cnt_param.tcalc) + Duration(maximum(cnt_param.dtexec)) <= control_chosen.t;
#@test control_chosen.t <= t_init + Duration(cnt_param.tcalc) + Duration(cnt_param.dtr)
@test control_chosen.t <= t_init + Duration(sim_param.dto*sim_param.prediction_steps);
@test norm(control_chosen.u) < eamax || norm(control_chosen.u) ≈ eamax;

# sac_control_update_test
tcalc_actual, best_schedule, sim_result =
    sac_control_update(w_init, u_schedule, target_trajectory, prediction_dict,
                       sim_param, cnt_param);
perturbation_times = [k for (k,v) in best_schedule if v==control_chosen.u];
if length(perturbation_times) > 1
    @test best_schedule[control_chosen.t - Duration(dtc)] == control_chosen.u
    @test best_schedule[control_chosen.t] != control_chosen.u
    @test in(to_sec(perturbation_times[end] - perturbation_times[1]) + sim_param.dtc, cnt_param.dtexec)
end

# # Main functions test
u_schedule = convert_to_schedule(w_init.t, u_array_base, sim_param);

controller = RSSACController(traj_predictor, u_schedule, sim_param, cnt_param);
scene_loader = TrajectronSceneLoader(scene_param, verbose=false);
traj_predictor = TrajectronPredictor(predictor_param,
                                     scene_loader.model_dir,
                                     scene_loader.param.conf_file_name,
                                     device, verbose=false);
initialize_scene_graph!(traj_predictor, scene_loader);
ado_states = fetch_ado_positions!(scene_loader, return_full_state=true);
ado_positions = reduce_to_positions(ado_states);

schedule_prediction!(controller, ado_states);
wait(controller.prediction_task);
@test istaskdone(controller.prediction_task);
@test controller.prediction_dict_tmp["PEDESTRIAN/25"][1:50, :, :] ==
      controller.prediction_dict_tmp["PEDESTRIAN/25"][51:100, :, :]

schedule_control_update!(controller, w_init, target_trajectory);
@test !isnothing(controller.prediction_dict);
@test isnothing(controller.prediction_task);
wait(controller.control_update_task);
@test istaskdone(controller.control_update_task);

ado_positions_str = convert_nodes_to_str(ado_positions);
latest_ado_pos_dict = deepcopy(ado_positions_str);
ado_1 = collect(keys(latest_ado_pos_dict))[1];
ado_2 = collect(keys(latest_ado_pos_dict))[2];
prediction_array_ado_2 = copy(controller.prediction_dict[ado_2]);
latest_ado_pos_dict[ado_2] += [0.5, 0.5];
pop!(latest_ado_pos_dict, ado_1);
latest_ado_pos_dict["PEDESTRIAN/999"] = [1.0, 1.0];
adjust_old_prediction!(controller, ado_positions_str, latest_ado_pos_dict);
@test all(controller.prediction_dict[ado_2][:, :, 1] .- prediction_array_ado_2[:, :, 1] .≈ 0.5);
@test all(controller.prediction_dict[ado_2][:, :, 2] .- prediction_array_ado_2[:, :, 2] .≈ 0.5);
@test !in(controller.prediction_dict, ado_1);
@test size(controller.prediction_dict["PEDESTRIAN/999"]) == (num_samples*num_controls, sim_param.prediction_steps, 2);
@test all(controller.prediction_dict["PEDESTRIAN/999"][:, :, 1] .== 1.0)
@test all(controller.prediction_dict["PEDESTRIAN/999"][:, :, 2] .== 1.0)

u_1 = control!(controller, Time(1.0));
@test u_1 == [0.0, 0.0];
@test !isnothing(controller.sim_result);
@test !isnothing(controller.tcalc_actual);
#@test !isnothing(controller.u_init_time);
#@test !isnothing(controller.u_last_time);
#@test !isnothing(controller.u_value);
@test isnothing(controller.control_update_task);

# # constraint_time tests
constraint_time = 0.1;
cnt_param = ControlParameter(eamax, tcalc, dtexec, dtr, u_nominal_base,
                             u_nominal_cand, nominal_search_depth,
                             constraint_time=constraint_time);
# Get simulation result and best nominal u_array
sim_result, u_nominal_array = simulate(w_init, u_arrays, target_trajectory,
                                       prediction_dict, sim_param, cnt_param.constraint_time);
# get_control_coeffs_test
coeff_matrix, coeff_matrix_constraint = get_control_coeffs(sim_result, sim_param, cnt_param);
#@test size(coeff_matrix) == (2, 20);
@test size(coeff_matrix) == (2, 480);
@test !isnothing(coeff_matrix_constraint)
@test size(coeff_matrix_constraint) == (2, 10);

u_array_constraint, cost_array_constraint = solve_multi_qcqp(u_nominal_array,
                                                             coeff_matrix,
                                                             coeff_matrix_constraint,
                                                             sim_param, cnt_param);
@test all([norm(u) <= cnt_param.eamax for u in u_array_constraint])
control_schedule_array_constraint =
    get_control_schedule(sim_result, u_nominal_array, coeff_matrix,
                         sim_param, cnt_param, coeff_matrix_constraint);
@test all([s.u == u for (s, u) in zip(control_schedule_array_constraint[1:10], u_array_constraint)])
@test all([s.cost == cost for (s, cost) in zip(control_schedule_array_constraint[1:10], cost_array_constraint)])
@test all([s.t == s_constraint.t for (s, s_constraint) in zip(control_schedule_array[1:10], control_schedule_array_constraint[1:10])])
@test all([s.u == s_constraint.u for (s, s_constraint) in zip(control_schedule_array[11:20], control_schedule_array_constraint[11:20])])
@test all([s.cost == s_constraint.cost for (s, s_constraint) in zip(control_schedule_array[11:20], control_schedule_array_constraint[11:20])])
@test all([s.t == s_constraint.t for (s, s_constraint) in zip(control_schedule_array[11:20], control_schedule_array_constraint[11:20])])
# sac_control_update_test
tcalc_actual, best_schedule, sim_result =
    sac_control_update(w_init, u_schedule, target_trajectory, prediction_dict,
                       sim_param, cnt_param);
end


text:
#not currently used

#=

function get_data_info_oldversion(at_set, dim)

    # very  minor changes, see differs. This is soley for format conversion from old coefs to new format, which has extra 3body terms.

    n_2body = 5
    n_2body_onsite = 4
    n_2body_S = 7
    n_3body = 5 #differs. is 7 in the main version
    n_3body_same = 5
    n_3body_onsite = 4
    n_3body_onsite_same = 2

    data_info = Dict{Tuple, Array{Int64,1}}()
    orbs = []
    if dim == 2 #2body

        at_list = [i for i in at_set]
#        println(at_list)
        if length(at_list) == 1
            at_list = [at_list[1], at_list[1]]
        end
        sort!(at_list)
#        println(at_list)

        orbs1 = atoms[at_list[1]].orbitals
        orbs2 = atoms[at_list[2]].orbitals

        at1 = at_list[1]
        at2 = at_list[2]

        if at1 == at2
            same_at = true
        else
            same_at = false
        end

#        orbs = []

        #2body part
        function get2bdy(n, symb)
            tot=0
            for o1 in orbs1
                for o2 in orbs2
                    if same_at && ((o2 == :s && o1 == :p) || (o2 == :s && o1 == :d) || (o2 == :p && o1 == :d))
                        continue
                    end
#                    push!(orbs, (o1, o2, symb))
                    if o1 == :s && o2 == :s
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n
                        tot += n

                    elseif (o1 == :s && o2 == :p ) || (o1 == :p && o2 == :s )
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n
                        tot += n
#                        if same_at
#                            data_info[(o2, o1, symb)] = data_info[(o1, o2, symb)]
#                        end

                    elseif (o1 == :p && o2 == :p )
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n*2
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n*2
                        tot += n*2

#                    elseif (o1 == :p && o2 == :p )
#                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n*2
#                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n*2
#                        tot += n*2

                    elseif (o1 == :s && o2 == :d ) || (o1 == :d && o2 == :s )
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n
                        tot += n

                    elseif (o1 == :p && o2 == :d ) || (o1 == :d && o2 == :p )
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n*2
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n*2
                        tot += n*2

                    elseif (o1 == :d && o2 == :d )
                        data_info[(at1, o1, at2, o2, symb)] = tot+1:tot+n*3
                        data_info[(at2, o2, at1, o1, symb)] = tot+1:tot+n*3
                        tot += n*3

                    end
                end
            end
            return tot
        end

        totH = get2bdy(n_2body, :H)
        totS = get2bdy(n_2body_S, :S)
#        println("totH $totH totS $totS")
        #onsite part
        function getonsite(atX,orbsX, tot, n)
            for o1 in orbsX
                for o2 in orbsX
                    if (o2 == :s && o1 == :p) || (o2 == :s && o1 == :d) || (o2 == :p && o1 == :d)
                        continue
                    end

#                    push!(orbs, (o1, o2, :O))
                    if o1 == :s && o2 == :s
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n
#                        println("data_info" , (atX, o1, o2, :O), tot+1:tot+n)
                        tot += n
                    elseif (o1 == :s && o2 == :p )
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n
                        data_info[(atX, o2, o1, :O)] = data_info[(atX, o1, o2, :O)]
                        tot += n
                    elseif (o1 == :p && o2 == :p )
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n*2
                        tot += n*2

                    elseif o1 == :s && o2 == :d
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n
                        data_info[(atX, o2, o1, :O)] = data_info[(atX, o1, o2, :O)]
                        tot += n

                    elseif o1 == :p && o2 == :d
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n
                        data_info[(atX, o2, o1, :O)] = data_info[(atX, o1, o2, :O)]
                        tot += n

                    elseif o1 == :d && o2 == :d
                        data_info[(atX, o1, o2, :O)] = tot+1:tot+n*2
                        data_info[(atX, o2, o1, :O)] = data_info[(atX, o1, o2, :O)]
                        tot += n*2

                    end
                end
            end
            return tot
        end

        if same_at #true onsite terms
           for o in orbs1
#                println("true onsite ", o)
                data_info[(at1, o, :A)] = [totH+1]
                totH += 1
            end
        end

        totHO = getonsite(at1, orbs1, totH, n_2body_onsite)

        if !(same_at) #need reverse if not same atom
            totHO = getonsite(at2, orbs2, totHO, n_2body_onsite)
        end

    elseif dim == 3 #3body

        totS = 0 #no 3body overlap terms

        at_list = [i for i in at_set]
        sort!(at_list)

        if length(at_list) == 1

            #permutations are trivial
            perm_ij = [[at_list[1], at_list[1], at_list[1]]]
            perm_on = [[at_list[1], at_list[1], at_list[1]]]
        elseif length(at_list) == 2

            #unique permutations
            perm_ij = [[at_list[1], at_list[1], at_list[2]] ,
                       [at_list[2], at_list[2], at_list[1]] ,
                       [at_list[1], at_list[2], at_list[1]] ,
                       [at_list[1], at_list[2], at_list[2]] ]

            perm_on = [[at_list[1], at_list[2], at_list[2]] ,
                        [at_list[1], at_list[1], at_list[2]] ,
                        [at_list[2], at_list[1], at_list[1]] ,
                        [at_list[2], at_list[1], at_list[2]] ]

        elseif length(at_list) == 3

            #all permutations exist hij
            perm_ij = [[at_list[1], at_list[2], at_list[3]] ,
                       [at_list[1], at_list[3], at_list[2]] ,
                       [at_list[2], at_list[1], at_list[3]] ,
                       [at_list[2], at_list[3], at_list[1]] ,
                       [at_list[3], at_list[1], at_list[2]] ,
                       [at_list[3], at_list[2], at_list[1]] ]

            #onsite can flip last 2 atoms
            perm_on = [[at_list[1], at_list[2], at_list[3]] ,
                       [at_list[2], at_list[1], at_list[3]] ,
                       [at_list[3], at_list[1], at_list[2]]]

        else
            println("ERROR  get_data_info $dim $at_set $at_list")
        end

        function get3bdy(n, symb, start, at1, at2, at3)
            tot=start

            orbs1 = atoms[at1].orbitals
            orbs2 = atoms[at2].orbitals

            if at1 == at2
                same_at = true
            else
                same_at = false
            end

            for o1 in orbs1
                for o2 in orbs2
                    if same_at && ((o2 == :s && o1 == :p) || (o2 == :s && o1 == :d) || (o2 == :p && o1 == :d))
                        continue
                    end

#                    push!(orbs, (o1, o2, symb))

                    if same_at
                        data_info[(at1, o1, at2, o2, at3,  symb)] = collect(tot+1:tot+n)
                        data_info[(at2, o2, at1, o1, at3,  symb)] = collect(tot+1:tot+n)
                    else                                                  #[1 2 3 4 5 6 7 8  9  10 11 12 13 14 15 16 17 18]
                        data_info[(at1, o1, at2, o2, at3,  symb)] = collect(tot+1:tot+n)
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 4 6 2 5 3 7 10 12 8  11 9  13 16 18 14 17 15]' #switch 2 4 and 3 6
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 4 6 2 5 3 7 10 12 8  11 9  ]' #switch 2 4 and 3 6
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 3 2 4 5 7 6 8 9]' #switch 2 4 and 3 6
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 3 2 4 5 7 6]' #switch 2 4 and 3 6
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 3 2  4 6 5]' #switch 2 4 and 3 6
#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1 3 2 4 6 5  7 9 8 ]' #switch 2 4 and 3 6

#                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1, 3, 2, 4, 5, 7, 6  ] #switch 2 4 and 3 6

                        data_info[(at2, o2, at1, o1, at3,  symb)] = tot .+ [1, 3, 2, 4, 5  ] #differs
                    end

                    tot += n

#                    if same_at
#                        data_info[(o2, o1, symb)] = data_info[(o1, o2, symb)]
#                    end

                end
            end
            return tot
        end

#        if at_list[2] == at_list[3]
#            same_at_on = true
#        else
#            same_at_on = false
 #       end

        function get3bdy_onsite(n, same_at,symb, start, at1, at2, at3)
#            if at2 == at3  #|| at1 == at2 || at1 == at3
#                same_at = true
#            else
#                same_at = false
#            end

            orbs1 = atoms[at1].orbitals

#            println("get3bdy_onsite $at1 $at2 $at3 $n")

            tot=start
            for o1 in orbs1
#                data_info[(at1, o1,at2, at3,  symb)] = collect(tot+1:tot+n)
#                data_info[(at1, o1,at3, at2,  symb)] = collect(tot+1:tot+n)

#                push!(orbs, (at1, o1,at2, at3,  symb))
#                push!(orbs, (at1, o1,at3, at2,  symb))

                if same_at
                    data_info[(at1, o1,at2, at3,  symb)] = collect(tot+1:tot+n)
                    data_info[(at1, o1,at3, at2,  symb)] = collect(tot+1:tot+n)
                else
                    data_info[(at1, o1,at2, at3,  symb)] = collect(tot+1:tot+n)
#                    data_info[(at1, o1,at3, at2,  symb)] = collect(tot+1:tot+n)
                    data_info[(at1, o1,at3, at2,  symb)] = tot .+ [1, 3, 2, 4]

#                    data_info[(at1, o1,at2, at3,  symb)] = collect(tot+1:tot+n)
#                    data_info[(at1, o1,at3, at2,  symb)] = tot .+ [1 3 2 4]'
                end
                tot += n                               #       1 2 3 4 5 6 7 8

            end
            return tot
        end

        tot_size = 0
        for p in perm_ij
            if p[1] == p[2]
                tot_size = get3bdy(n_3body_same, :H, tot_size, p[1], p[2], p[3])
            else
                tot_size = get3bdy(n_3body, :H, tot_size, p[1], p[2], p[3])
            end
        end

        for p in perm_on
            if  (p[1] == p[2] ||  p[2] == p[3] || p[1] == p[3])
                tot_size = get3bdy_onsite(n_3body_onsite_same,true, :O, tot_size, p[1], p[2], p[3]) #all diff
            else
                tot_size = get3bdy_onsite(n_3body_onsite,false, :O, tot_size, p[1], p[2], p[3]) #
            end
        end

        totHO = tot_size

    else
        println("error, only 2 or 3 body terms, you gave me : ", at_list)
    end
    return totHO ,totS, data_info, orbs

end

function fix_format_change(datH, totHnew, dim, at_list, data_info)

    #this converts from old data_info to new data_info, setting the extra terms to zero.

    totH,totS, data_info_old, orbs = get_data_info_oldversion(at_list, dim)

    datHnew = zeros(totHnew)

    for k in keys(data_info)
        dnew = data_info[k]
        dold = data_info_old[k]

        n=length(dold)
        datHnew[dnew[1:n] ] = datH[dold]  #the extra terms are left as zero
    end

    return datHnew

end
=#


text:
using NonlinearEigenproblems
using SparseArrays

export RKNEP
export MatrixAndFunction
export LowRankMatrixAndFunction

export get_rk_nep

import Base.size
import ..NEPCore.compute_Mder
import ..NEPCore.compute_Mlincomb
import ..NEPTypes.get_Av
import ..NEPTypes.get_fv
import ..NEPTypes.LowRankFactorizedNEP

"NEP instance supplemented with various structures used for rational Krylov problems."
struct RKNEP{S<:AbstractMatrix{<:Number}, T<:Number}
    nep::NEP            # Original NEP problem
    spmf::Bool          # Whether this NEP is a sum a products of matrices and functions
    p::Int              # Order of polynomial part
    q::Int              # Number of nonlinear matrices and functions
    BBCC::S             # The polynomial and nonlinear matrices concatenated vertically
    is_low_rank::Bool   # Whether the nonlinear matrices of this NEP have an attached low rank LU factorization
    r::Int              # Sum of ranks of low rank nonlinear matrices
    iL::Vector{Int}     # Vector with indices of L-factors
    iLr::Vector{Int}    # Vector with row indices of L-factors that have any entries present
    L::Vector{S}        # Vector of L factors of low rank LU factorization of nonlinear part
    LL::Vector{SparseVector{T,Int}} # Rows of L factors concatenated into vectors (only rows with nnz > 0 are stored)
    UU::S               # The U factors concatenated vertically
end

RKNEP(::Type{T}, nep::NEP) where T<:Number =
    RKNEP(nep, false, 0, 0, Matrix{T}(undef, 0, 0), false, 0, Vector{Int}(), Vector{Int}(), Vector{Matrix{T}}(), Vector{SparseVector{T,Int}}(), Matrix{T}(undef, 0, 0))

RKNEP(nep::NEP, p, q, BBCC::AbstractMatrix{T}) where T<:Number =
    RKNEP(nep, true, p, q, BBCC, false, 0, Vector{Int}(), Vector{Int}(), Vector{typeof(BBCC)}(), Vector{SparseVector{T,Int}}(), similar(BBCC, 0, 0))

RKNEP(nep::NEP, p, q, BBCC, r, iL, iLr, L, LL, UU) =
    RKNEP(nep, true, p, q, BBCC, true, r, iL, iLr, L, LL, UU)

struct LowRankMatrixAndFunction{S<:AbstractMatrix{<:Number}}
    A::S
    L::S        # L factor of LU-factorized A
    U::S        # U factor of LU-factorized A
    f::Function
end

"Create low rank LU factorization of A."
function LowRankMatrixAndFunction(A::AbstractMatrix{<:Number}, f::Function)
    L, U = low_rank_lu_factors(A)
    LowRankMatrixAndFunction(A, L, U, f)
end

export LowRankFactorizedNEP

# Additional constructor for particle example
function LowRankFactorizedNEP(Amf::AbstractVector{LowRankMatrixAndFunction{S}}) where {T<:Number, S<:AbstractMatrix{T}}
    # if A is not specified, create it from LU factors
    A = [isempty(M.A) ? M.L * M.U' : M.A for M in Amf]
    L = getfield.(Amf, :L)
    U = getfield.(Amf, :U)
    f = getfield.(Amf, :f)
    rank = mapreduce(M -> size(M.U, 2), +, Amf)
    return LowRankFactorizedNEP(SPMF_NEP(A, f, align_sparsity_patterns=true), rank, L, U)
end

function low_rank_lu_factors(A::SparseMatrixCSC{<:Number,Int})
    n = size(A, 1)
    idx = findall(!iszero, A)
    r = extrema(getindex.(idx, 1))
    c = extrema(getindex.(idx, 2))
    B = A[r[1]:r[2], c[1]:c[2]]
    L, U = lu(Matrix(B); check = false)
    Lc, Uc = compactlu(sparse(L), sparse(U))
    Lca = spzeros(n, size(Lc, 2))
    Lca[r[1]:r[2], :] = Lc
    Uca = spzeros(size(Uc, 1), n)
    Uca[:, c[1]:c[2]] = Uc
    return Lca, sparse(Uca')

    # TODO use this; however we then need to support permutation and scaling
    #F = lu(B)
    #Lcf,Ucf = compactlu(sparse(F.L),sparse(F.U))
    #Lcaf = spzeros(n, size(Lcf, 2))
    #Lcaf[r[1]:r[2], :] = Lcf
    #Ucaf = spzeros(size(Ucf, 1), n)
    #Ucaf[:, c[1]:c[2]] = Ucf
    #Ucaf = Ucaf'
end

function compactlu(L, U)
    n = size(L, 1)
    select = map(i -> nnz(L[i:n, i]) > 1 || nnz(U[i, i:n]) > 0, 1:n)
    return L[:,select], U[select,:]
end

"Create RKNEP instance, exploiting the type of the input NEP as much as possible"
function get_rk_nep(::Type{T}, nep::NEP) where T<:Real
    # Most generic case: No coefficient matrices, all we have is M(λ)
    if !isa(nep, AbstractSPMF)
        return RKNEP(T, nep)
    end

    Av = get_Av(nep)
    BBCC = vcat(Av...)::eltype(Av)

    # Polynomial eigenvalue problem
    if isa(nep, PEP)
        return RKNEP(nep, length(Av) - 1, 0, BBCC)
    end

    # If we can't separate the problem into a PEP + SPMF, consider it purely SPMF
    if !isa(nep, SPMFSumNEP{PEP,S} where S<:AbstractSPMF)
        return RKNEP(nep, -1, length(Av), BBCC)
    end

    p = length(get_Av(nep.nep1)) - 1
    q = length(get_Av(nep.nep2))

    # Case when there is no low rank structure to exploit
    if q == 0 || !isa(nep.nep2, LowRankFactorizedNEP{S} where S<:Any)
        return RKNEP(nep, p, q, BBCC)
    end

    # L and U factors of the low rank nonlinear part
    L = nep.nep2.L
    UU = hcat(nep.nep2.U...)::eltype(nep.nep2.U)
    r = nep.nep2.r
    iL = zeros(Int, r)
    c = 0
    for ii = 1:q
        ri = size(L[ii], 2)
        iL[c+1:c+ri] .= ii
        c += ri
    end

    # Store L factors in a compact format to speed up system solves later on
    LL = Vector{SparseVector{eltype(L[1]),Int}}()
    iLr = Vector{Int}()
    for ri = 1:size(nep, 1)
        row = reduce(vcat, [L[i][ri,:] for i=1:length(L)])
        if nnz(row) > 0
            push!(LL, row)
            push!(iLr, ri)
        end
    end

    return RKNEP(nep, p, q, BBCC, r, iL, iLr, L, LL, UU)
end


text:
using GLMakie, GeometryBasics, RPRMakie, RadeonProRender
using Colors, FileIO
using Colors: N0f8
RPR = RadeonProRender

context = RPR.Context()
matsys = RPR.MaterialSystem(context, 0)

materials = [
    RPR.DiffuseMaterial(matsys) RPR.MicrofacetMaterial(matsys);
    RPR.ReflectionMaterial(matsys) RPR.RefractionMaterial(matsys);
    RPR.EmissiveMaterial(matsys) RPR.UberMaterial(matsys);
]

cat = load(GLMakie.assetpath("cat.obj"))

fig = Figure(resolution=(1000, 1000))
ax = LScene(fig[1, 1], scenekw=(show_axis=false,))
palette = reshape(Makie.default_palettes.color[][1:6], size(materials))
for i in CartesianIndices(materials)
    x, y = Tuple(i)
    catmesh = mesh!(ax, cat, material=materials[i], color=palette[i])
    translate!(catmesh, x, y, 0)
end
# materials[3, 1].color = Vec4(200)
display(fig)
context, task = RPRMakie.replace_scene_rpr!(ax.scene, context, matsys)
# fetch(task)

volmat = materials[end, end]

volmat.scattering = Vec3(0, 0, 0)
volmat.absorption = RGB(0.01, 0.01, 0.01)
volmat.multiscatter = true
# volmat.emission = RPR.RPR_MATERIAL_INPUT_EMISSION,
# volmat.scatter_direction = RPR.RPR_MATERIAL_INPUT_G,


text:
module DF

using Base.Iterators
import IterTools: nth

import Base: show, first, last, iterate, size, length, ndims,
              summary, getindex, setindex!, ∘, iszero, isone,
              isless
#import Base.∘

export Segment, DirectProduct, BooleanCube
export DiscreteFunction, ResidueFunction, ExtendedResidueFunction, BooleanFunction
export FunProduct, FunTupling
export tablegen
export domain, codomain
export BooleanGenerator, ResidueGenerator, ExtendedResidueGenerator
export FunProductGenerator, FunTuplingGenerator

# short names
export DProd, BCube
export DFun, ResFun, ExtFun, BoolFun, FProd, FTuple
export BoolGen, ResGen, ExtGen, FTGen, FPGen

include("finite_sets.jl")
include("types.jl")
include("helper_utils.jl")
include("generators.jl")
include("boolean/boolean.jl")
include("compositions.jl")
end # module

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
using StatsFuns: logistic
using Base: setindex!, getindex

"""
    Nonparametric discrete distribution with parameters in logit space.

This represents a discrete distribution with p.m.f of

```math
P(K = k) = (1 - \\rho_k) \\prod_{i = 1}^{k - 1} \\rho_i.
```
"""
struct LogitNPD{T<:Real,A<:AbstractVector{T}} <: Distributions.DiscreteUnivariateDistribution
    logitρ  ::  A
    l_init  ::  T   # initial value for each logitρᵢ
end

LogitNPD(; l_init=zero(Double)) = LogitNPD(0; l_init=l_init)
LogitNPD(alpha::AbstractFloat; l_init=zero(Double)) = LogitNPD(ceil(Int, alpha); l_init=l_init)
LogitNPD(k_init::Int; l_init=zero(Double)) = LogitNPD(ones(Double, k_init) * l_init, l_init)

function getlogitρ(lnpd::LogitNPD{T,A}, k::Int) where {T<:Real,A<:AbstractVector{T}}
    l = length(lnpd.logitρ)
    if k > l
        append!(lnpd.logitρ, ones(T, k - l) * lnpd.l_init...)
    end
    return lnpd.logitρ[k]
end

function getlogitρ(lnpd::LogitNPD{T,A}, k1::Int, k2::Int) where {T<:Real,A<:AbstractVector{T}}
    l = length(lnpd.logitρ)
    if k2 > l
        append!(lnpd.logitρ, ones(T, k2 - l) * lnpd.l_init...)
    end
    return lnpd.logitρ[k1:k2]
end

getρ(lnpd::LogitNPD, k...) = logistic.(getlogitρ(lnpd, k...))

function getlogρ(lnpd::LogitNPD, k...)
    l = getlogitρ(lnpd, k...)
    return l .- log1pexp.(l)
end

function pdf(lnpd::LogitNPD, k::Int)
    if k <= 0
        return 0
    elseif k == 1
        return 1 - getρ(lnpd, 1)
    else
        return prod(getρ(lnpd, 1, k - 1)) * (1 - getρ(lnpd, k))
    end
end

function logpdf(lnpd::LogitNPD, k::Int)
    if k <= 0
        return -Inf
    elseif k == 1
        return log(1 - getρ(lnpd, 1))
    else
        return sum(getlogrho(lnpd, 1, k - 1)) + log(1 - getρ(lnpd, k))
    end
end

function ccdf(lnpd::LogitNPD{T,A}, k::Int) where {T<:Real,A<:AbstractVector{T}}
    if k <= 0
        return one(T)
    else
        return prod(getρ(lnpd, 1, k))
    end
end

function cdf(lnpd::LogitNPD, k::Int)
    return 1 - ccdf(lnpd, k)
end

function rand(lnpd::LogitNPD{T,A}) where {T<:Real,A<:AbstractVector{T}}
    i = 1
    while true
        u = rand(T)
        if u > getρ(lnpd, i)
            return i
        end
        i = i + 1
    end
end

function rand(lnpd::LogitNPD{T,A}, n::Int) where {T<:Real,A<:AbstractVector{T}}
    return Int[rand(lnpd) for _ = 1:n]
end

function mode(lnpd::LogitNPD)
    p = pdf(lnpd, 1)
    p_acc = p_max = p
    k = k_max = 1
    while p_max < 1 - p_acc
        p_next = pdf(lnpd, k + 1)
        p_acc += p_next
        if p_next > p_max
            p_max = p_next
            k_max = k
        end
        k = k + 1
    end
    return k_max
end

function invlogcdf(lnpd::LogitNPD, lc::AbstractFloat)
    p_target = exp(lc)
    k = 1
    p = pdf(lnpd, 1)
    p_acc = p
    while p_acc < p_target
        k += 1
        p_acc += pdf(lnpd, k)
    end
    return k
end


text:
abstract type Event end
abstract type ObjEvent <: Event end
abstract type BaseEvent <: ObjEvent end # Basic || influencial
abstract type HighLvlEvent <: ObjEvent end #
abstract type UniverseEvent <: Event end

# Obj Events

#- Base
#-- Disappear
struct Disappear <: BaseEvent
end

#-- Appear
struct Appear <: BaseEvent
end
#-- Move
struct Move <: BaseEvent
    dr::Vector{Int64}
end

#-- ReSize
struct ReSize <: BaseEvent
    newSize::Int64
end

#- High level
struct ReDirection <: HighLvlEvent
    newDirection::Vector{Int64}
end

#-- Swallow
struct Swallow <: HighLvlEvent
end

#-- Slip
struct Slip <: HighLvlEvent
    upAndDown::Int64
end

#-- ReFre
struct ReFres <: HighLvlEvent
    newFres::Vector{Float64}
end

# Universe Event
struct IntroduceRandomDefects <: UniverseEvent
    numDefects::Int64
    maxSize::Int64
end

# Event Container
mutable struct EventContainer{T<:Event}
    fre::Float64
    fres::Vector{Float64}
    events::Vector{T}
end

max_stars_repo_name: permissionx/mc-ke

max_issues_repo_name: permissionx/mcke

max_forks_repo_name: permissionx/mcke


text:
"""
Constant that represents document term vector (DTV) models used in text embedding.
"""
const DTVModel{T,S} = Union{
    StringAnalysis.RPModel{S,T,<:AbstractMatrix{T},<:Integer},
    StringAnalysis.LSAModel{S,T,<:AbstractMatrix{T},<:Integer}
}

"""
Structure for document embedding using DTV's.
"""
struct DTVEmbedder{T,S} <: AbstractEmbedder{T,S}
    model::DTVModel{T,S}
    config::NamedTuple
end

DTVEmbedder(mtype::Type{<:DTVModel}, dtm, config; kwargs...) = DTVEmbedder(mtype(dtm; kwargs...), config)

# Document to vector embedding function
function __document2vec(embedder::DTVEmbedder{T,S},
                        document::Vector{String};  # a vector of sentences
                        isregex::Bool=false,
                        kwargs...  # for the unused arguments
                       )::SparseVector{T, Int} where {T,S}
    dtv_function = ifelse(isregex, dtv_regex, dtv)
    words = Vector{String}()
    for sentence in document
        for word in tokenize(sentence, method=:stringanalysis)
            push!(words, word)
        end
    end
    vocab_hash = embedder.model.vocab_hash
    model = embedder.model
    v = dtv_function(words, vocab_hash, T;
                     ngram_complexity=embedder.config.ngram_complexity,
                     tokenizer=DEFAULT_TOKENIZER,
                     lex_is_row_indices=true)
    embedded_document = StringAnalysis.embed_document(model, v)
    return embedded_document
end

function document2vec(embedder::DTVEmbedder{T,S},
                      document::Vector{String};  # a vector of sentences
                      isregex::Bool=false,
                      kwargs...  # for the unused arguments
                     )::Tuple{SparseVector{T, Int}, Bool} where {T,S}
    embedded_document = __document2vec(embedder,
                                       document;
                                       isregex=isregex,
                                       kwargs...)
    is_embedded = !iszero(embedded_document)
    # Check for OOV (out-of-vocabulary) policy
    if embedder.config.oov_policy == :large_vector && !is_embedded
        embedded_document .= T(DEFAULT_OOV_VAL)
    end
    return embedded_document, is_embedded
end

# Dimensionality function
function dimensionality(embedder::DTVEmbedder)
    # Return output dimensionality (second dim)
    return size(embedder.model)[2]
end


text:
# 1.1 UE-BMI under benchmark scenario
# NOTE: it is a 4-subfigure figure
   # 1.1.1 second, define a temp expr to save code lines
   tmpexpr = :(
      xlabel("Year"); ylabel("Percentage (%)"); grid(true);
      xlim([ idx_year2plot[1] - 1, idx_year2plot[end] + 1 ]);
      plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* dict_Demog_base["AgingPopuRatio"][idx_plot], "--b" );
      plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* (1.0 .- dict_Demog_base["WorkPopuRatio"] )[idx_plot], "-.r" );
   )
   # 1.1.2 then, do plotting
   figure( figsize = (13,8) )
      subplot(2,2,1)  # gap / pooling account's expenditure
         plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* (DatPkg_base.Dt[:LI]./DatPkg_base.Dt[:AggPoolExp])[idx_plot] )
         eval(tmpexpr)
         legend( ("Pool Gap/Pool Benefits","Aging Population Share (65+)",L"\rho"), loc = "best")
      subplot(2,2,2)  # gap / pooling account's income
         plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* (DatPkg_base.Dt[:LI]./DatPkg_base.Dt[:AggPoolIn])[idx_plot] )
         eval(tmpexpr)
         legend( ("Pool Gap/Pool Incomes","Aging Population Share (65+)",L"\rho"), loc = "best")
      subplot(2,2,3)  # gap / GDP
         plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* (DatPkg_base.Dt[:LI]./DatPkg_base.Dt[:Y])[idx_plot] )
         eval(tmpexpr)
         legend( ("Pool Gap/GDP","Aging Population Share (65+)",L"\rho"), loc = "best")
      subplot(2,2,4)  # gap / fiscal incomes
         plot( DatPkg_base.Dt[:Year][idx_plot], 100 .* (DatPkg_base.Dt[:LI]./(DatPkg_base.Dt[:TRw] .+ DatPkg_base.Dt[:TRc]))[idx_plot] )
         eval(tmpexpr)
         legend( ("Pool Gap/Tax Revenues","Aging Population Share (65+)",L"\rho"), loc = "best")
      tight_layout()  # tight layout of the figure
   # 1.1.3 finally, save the figure as a pdf file
   savefig( "./output/BenchProfile.pdf", format = "pdf" )
# -----------------------------------------
# 1.2 simulation v.s. accounting
   # 1.2.1 do underlying plotting & read in accounting data
   EasyPlot.Plot_Calibrate( DatPkg_base.Dt, DatPkg_base.Dst, DatPkg_base.Pt, DatPkg_base.Ps, DatPkg_base.Pc, env,
      YearRange = ( 2010, 2050 ), LineWidth = 1.0,   outpdf = nothing, picsize = (12,9),
      tmpLayout = (2,1)
   )
   tmpDat = EasyIO.readcsv("./data/Calib_统筹账户收支核算结果v3_190403.csv")
   tmpEndTime = 40 + 2
   # 1.2.2 decoration (NOTE: do not new a figure GUI, just decorate the current one)
   subplot(2,1,1)
      plot( tmpDat[2:tmpEndTime,1] , 100.0 .* tmpDat[2:tmpEndTime,4] ./ tmpDat[2:tmpEndTime,2] , "-.r" )
      legend(["Baseline simulation","Pooling gap / Pooling account expenditure"],
         fontsize = 14)
   subplot(2,1,2)
      plot( tmpDat[2:tmpEndTime,1] , 100.0 .* tmpDat[2:tmpEndTime,4] ./ tmpDat[2:tmpEndTime,3] , "-.r" )
      legend(["Baseline simulation","Pooling gap / Pooling account revenues"],
         fontsize = 14)
   tight_layout()
   # 1.2.3 save figure
   savefig( "./output/BenchmarkCpAccount.pdf", format = "pdf" )

#


text:
using Caesar, Caesar.ZmqCaesar
using Distributions
using LinearAlgebra

# Distribution: AliasingScalarSampler
packed = Packed_AliasingScalarSampler(collect(1:10), collect(11:20), 0, "AliasingScalarSampler")
json = JSON.json(packed)
sampler = convert(IncrementalInference.AliasingScalarSampler, JSON.parse(json))
packedCompare = convert(Packed_AliasingScalarSampler, sampler)
jsonCompare = JSON.json(packedCompare)

# Distribution: Normal
normal = Normal(5.0, 10.0)
packed = convert(Dict{String, Any}, normal)
back = convert(Normal, JSON.parse(JSON.json(packed)))

# Factor: PartialPriorRollPitchZ
a = MvNormal([0, 0, 0, 0], Matrix(LinearAlgebra.I, 4, 4))
b = Normal(0, 1)
partPrior = PartialPriorRollPitchZ(a, b)
j = JSON.json(convert(Dict{String, Any}, partPrior))
# Test the deserializer
back = convert(PartialPriorRollPitchZ, JSON.parse(j))
jback = JSON.json(convert(Dict{String, Any}, back))
jback = j