text:
module GeneratorsTest

using Test
using ForneyLab
using LinearAlgebra: Diagonal
using ForneyLab: entropiesSourceCode, energiesSourceCode, freeEnergySourceCode, marginalTableSourceCode, inboundSourceCode, scheduleSourceCode, removePrefix, vagueSourceCode, initializationSourceCode, optimizeSourceCode, algorithmSourceCode, valueSourceCode, posteriorFactorSourceCode

# create custom module
module Foo
    struct Baz
        f::Int64
    end
    struct Bar{T}
        f::T
    end
    struct Bak{T, K}
        f::T
        g::K
    end
end

@testset "removePrefix" begin
    # initial check
    @test removePrefix(ForneyLab.SPGaussianMeanPrecisionOutNPP) == "SPGaussianMeanPrecisionOutNPP"

    # initialize structures
    baz = Foo.Baz(1)
    bar = Foo.Bar(baz)
    bak = Foo.Bak(bar, bar)

    # check whether the removePrefix functions works as desired
    @test removePrefix(typeof(baz)) == "Baz"
    @test removePrefix(typeof(bar)) == "Bar{Baz}"
    @test removePrefix(typeof(bak)) == "Bak{Bar{Baz}, Bar{Baz}}"
end

@testset "valueSourceCode" begin
    @test valueSourceCode(1) == "1"
    @test valueSourceCode([1]) == "[1]"
    @test valueSourceCode(mat(1)) == "mat(1)"
    @test valueSourceCode(reshape([1,2,3], 3, 1)) == "reshape([1, 2, 3], (3, 1))"
    @test valueSourceCode([1 2; 2 1]) == "[1 2; 2 1]"
    @test valueSourceCode(Diagonal([1])) == "Diagonal([1])"
end

f() = 1.0 # Define a function

@testset "inboundSourceCode" begin
    # custom inbound
    inbound = Dict(:keyword => false,
                   :inx     => 1)
    @test inboundSourceCode(inbound) == "1"

    inbound = Dict(:keyword => false,
                   :g       => f)
    @test inboundSourceCode(inbound) == "f"

    inbound = Dict(:g => f)
    @test inboundSourceCode(inbound) == "g=f"

    # value inbound
    g = FactorGraph()
    var = Variable()
    inbound = Clamp(var, 1.0)
    inbound.dist_or_msg = Message
    @test inboundSourceCode(inbound) == "Message(Univariate, PointMass, m=1.0)"

    # buffer inbound
    g = FactorGraph()
    var = Variable()
    inbound = Clamp(var, 1.0)
    inbound.dist_or_msg = Message
    inbound.buffer_id = :x
    inbound.buffer_index = 1
    @test inboundSourceCode(inbound) == "Message(Univariate, PointMass, m=data[:x][1])"

    g = FactorGraph()
    var = Variable()
    inbound = Clamp(var, 1.0)
    inbound.dist_or_msg = Message
    inbound.buffer_id = :x
    inbound.buffer_index = 0
    @test inboundSourceCode(inbound) == "Message(Univariate, PointMass, m=data[:x])"

    # marginal
    inbound = MarginalEntry()
    inbound.marginal_id = :x
    @test inboundSourceCode(inbound) == "marginals[:x]"

    # message
    inbound = ScheduleEntry()
    inbound.schedule_index = 1
    @test inboundSourceCode(inbound) == "messages[1]"

    # nothing
    @test inboundSourceCode(nothing) == "nothing"
end

@testset "vagueSourceCode" begin
    entry = ScheduleEntry()
    entry.family = GaussianMeanPrecision
    entry.dimensionality = ()
    entry_code = vagueSourceCode(entry)
    @test entry_code == "vague(GaussianMeanPrecision)"

    entry = ScheduleEntry()
    entry.family = GaussianMeanPrecision
    entry.dimensionality = (1,)
    entry_code = vagueSourceCode(entry)
    @test entry_code == "vague(GaussianMeanPrecision, (1,))"
end

@testset "marginalTableSourceCode" begin
    inbounds = Vector{ScheduleEntry}(undef, 2)
    inbounds[1] = ScheduleEntry()
    inbounds[1].schedule_index = 1
    inbounds[2] = ScheduleEntry()
    inbounds[2].schedule_index = 2

    marginal_table = Vector{MarginalEntry}(undef, 3)
    marginal_table[1] = MarginalEntry()
    marginal_table[1].marginal_id = :x
    marginal_table[1].marginal_update_rule = Nothing
    marginal_table[1].inbounds = [inbounds[1]]
    marginal_table[2] = MarginalEntry()
    marginal_table[2].marginal_id = :y
    marginal_table[2].marginal_update_rule = ForneyLab.Product
    marginal_table[2].inbounds = inbounds
    marginal_table[3] = MarginalEntry()
    marginal_table[3].marginal_id = :z
    marginal_table[3].marginal_update_rule = ForneyLab.MGaussianMeanPrecisionGGD
    marginal_table[3].inbounds = inbounds

    marginal_table_code = marginalTableSourceCode(marginal_table)

    @test occursin("marginals[:x] = messages[1].dist", marginal_table_code)
    @test occursin("marginals[:y] = messages[1].dist * messages[2].dist", marginal_table_code)
    @test occursin("marginals[:z] = ruleMGaussianMeanPrecisionGGD(messages[1], messages[2])", marginal_table_code)
end

@testset "scheduleSourceCode" begin
    schedule = Vector{ScheduleEntry}(undef, 2)
    schedule[1] = ScheduleEntry()
    schedule[1].schedule_index = 1
    schedule[1].message_update_rule = Nothing
    schedule[1].inbounds = []
    schedule[2] = ScheduleEntry()
    schedule[2].message_update_rule = ForneyLab.SPGaussianMeanPrecisionOutNPP
    schedule[2].schedule_index = 2
    schedule[2].inbounds = [schedule[1]]

    schedule_code = scheduleSourceCode(schedule)
    @test occursin("messages[2] = ruleSPGaussianMeanPrecisionOutNPP(messages[1])", schedule_code)
end

@testset "entropiesSourceCode" begin
    inbounds = Vector{MarginalEntry}(undef, 2)
    inbounds[1] = MarginalEntry()
    inbounds[1].marginal_id = :y_x
    inbounds[2] = MarginalEntry()
    inbounds[2].marginal_id = :x

    entropies_vect = [Dict(:counting_number => -1,
                           :inbound         => inbounds[1]),
                      Dict(:counting_number => 2,
                           :inbound         => inbounds[2])]

    entropies_code = entropiesSourceCode(entropies_vect)

    @test occursin("F -= -1*differentialEntropy(marginals[:y_x])", entropies_code)
    @test occursin("F -= 2*differentialEntropy(marginals[:x])", entropies_code)
end

@testset "energiesSourceCode" begin
    inbound = MarginalEntry()
    inbound.marginal_id = :x
    energies_vect = [Dict(:counting_number => 1,
                          :node            => GaussianMeanPrecision,
                          :inbounds        => [inbound])]
    energies_code = energiesSourceCode(energies_vect)

    @test occursin("F += averageEnergy(GaussianMeanPrecision, marginals[:x])", energies_code)
end

@testset "initializationSourceCode" begin
    fg = FactorGraph()
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg, id=:X)
    algo = InferenceAlgorithm(pfz)
    pf.algorithm_id = algo.id
    pf.initialize = true
    entry = ScheduleEntry()
    entry.schedule_index = 1
    entry.initialize = true
    entry.family = GaussianMeanPrecision
    entry.dimensionality = ()
    pf.schedule = [entry]

    pf_code = initializationSourceCode(pf)
    @test occursin("function initX()", pf_code)
    @test occursin("messages[1] = Message(vague(GaussianMeanPrecision))", pf_code)
end

@testset "optimizeSourceCode" begin
    fg = FactorGraph()
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg, id=:X)
    algo = InferenceAlgorithm(pfz)
    pf.algorithm_id = algo.id
    pf.optimize = true

    pf_code = optimizeSourceCode(pf)
    @test occursin("function optimizeX!(data::Dict, marginals::Dict=Dict(), messages::Vector{Message}=initX())", pf_code)
end

@testset "posteriorFactorSourceCode" begin
    fg = FactorGraph()
    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg, id=:X)
    algo = InferenceAlgorithm(pfz)
    pf.algorithm_id = algo.id
    pf.schedule = []
    pf.marginal_table = []

    pf_code = posteriorFactorSourceCode(pf)
    @test occursin("function stepX!(data::Dict, marginals::Dict=Dict(), messages::Vector{Message}=Array{Message}(undef, 0))", pf_code)
end

@testset "algorithmSourceCode" begin
    algo = InferenceAlgorithm()
    code = algorithmSourceCode(algo)
    @test occursin("begin", code)
end

@testset "freeEnergySourceCode" begin
    algo = InferenceAlgorithm()
    algo.posterior_factorization.free_energy_flag = true
    free_energy_code = freeEnergySourceCode(algo)
    @test occursin("function freeEnergy(data::Dict, marginals::Dict)", free_energy_code)
end

end # module


text:
function datasets()
    dataset_folder = joinpath(@__DIR__, "..", "..", "datasets")
    datasets = Any[]
    for entry in readdir(dataset_folder)
        path = joinpath(dataset_folder, entry)
        isdir(path) || continue
        if entry == "ILSVRC"
            push!(datasets, ImageNet.RawFS(path, :kaggle))
            continue
        elseif entry == "cifar-10-batches-bin"
            push!(datasets, CIFAR10.BinPackedFS(path))
            continue
        end
    end
    datasets
end

"""
Get the appropriate dataset from anywhere we can find.
Available options: ImageNet
"""
function dataset(which)
    if which === ImageNet
        sets = collect(Iterators.filter(x->isa(x, ImageNet.DataSet), datasets()))
        isempty(sets) && error("No ImageNet dataset available. "*
            "See datasets/README.md for download instructions")
        return first(sets)
    elseif which == CIFAR10
        sets = collect(Iterators.filter(x->isa(x, CIFAR10.DataSet), datasets()))
        if isempty(sets)
            download(CIFAR10)
            sets = collect(Iterators.filter(x->isa(x, CIFAR10.DataSet), datasets()))
        end
        return first(sets)
    else
        error("Autodetection not supported for $(which)")
    end
end
valimgs(m::Module) = valimgs(dataset(m))
testimgs(m::Module) = testimgs(dataset(m))
trainimgs(m::Module) = trainimgs(dataset(m))

function download(which)
    if which === ImageNet
        error("ImageNet is not automatiacally downloadable. See instructions in datasets/README.md")
    elseif which == CIFAR10
        local_path = joinpath(@__DIR__, "..", "..", "datasets", "cifar-10-binary.tar.gz")
        dir_path = joinpath(@__DIR__,"..","..","datasets")
        if(!isdir(joinpath(dir_path, "cifar-10-batches-bin")))
            if(!isfile(local_path))
                Base.download("https://www.cs.toronto.edu/~kriz/cifar-10-binary.tar.gz", local_path)
            end
            run(`tar -xzvf $local_path -C $dir_path`)
        end
    else
        error("Download not supported for $(which)")
    end
end

ext: jl


text:
import UUIDs: UUID, uuid1
import .ExpressionExplorer: SymbolsState

Base.@kwdef struct CellOutput
    body::Union{Nothing,String,Vector{UInt8},Dict,Main.PlutoShared.UI}=nothing
    mime::MIME=MIME("text/plain")
    rootassignee::Union{Symbol,Nothing}=nothing

    "Time that the last output was created, used only on the frontend to rerender the output"
    last_run_timestamp::Float64=0

    "Whether `this` inside `<script id=something>` should refer to the previously returned object in HTML output. This is used for fancy animations. true iff a cell runs as a reactive consequence."
    persist_js_state::Bool=false
end

struct CellDependencies{T} # T == Cell, but this has to be parametric to avoid a circular dependency of the structs
    downstream_cells_map::Dict{Symbol,Vector{T}}
    upstream_cells_map::Dict{Symbol,Vector{T}}
    precedence_heuristic::Int
end

"The building block of a `Notebook`. Contains code, output, reactivity data, mitochondria and ribosomes."
Base.@kwdef mutable struct Cell
    "Because Cells can be reordered, they get a UUID. The JavaScript frontend indexes cells using the UUID."
    cell_id::UUID=uuid1()

    code::String=""
    code_folded::Bool=false

    output::CellOutput=CellOutput()
    queued::Bool=false
    running::Bool=false

    errored::Bool=false
    runtime::Union{Nothing,UInt64}=nothing

    # note that this field might be moved somewhere else later. If you are interested in visualizing the cell dependencies, take a look at the cell_dependencies field in the frontend instead.
    cell_dependencies::CellDependencies{Cell}=CellDependencies{Cell}(Dict{Symbol,Vector{Cell}}(), Dict{Symbol,Vector{Cell}}(), 99)

    owner::Union{Nothing,Symbol}=nothing
end

Cell(cell_id, code) = Cell(cell_id=cell_id, code=code)
Cell(code) = Cell(uuid1(), code)

cell_id(cell::Cell) = cell.cell_id

function Base.convert(::Type{Cell}, cell::Dict)
	owner = get(cell, "owner", nothing)
	if owner !== nothing
		owner = Symbol(owner)
	end
	Cell(
        cell_id=UUID(cell["cell_id"]),
        code=cell["code"],
        code_folded=cell["code_folded"],
        owner=owner,
    )
end
function Base.convert(::Type{UUID}, string::String)
    UUID(string)
end


text:
# SAM Writer
# ==========

type SAMWriter{T<:IO} <: Bio.IO.AbstractWriter
    stream::T
end

function Base.close(writer::SAMWriter)
    close(writer.stream)
end

function Base.write(writer::SAMWriter, header::SAMHeader)
    n = 0
    if haskey(header, "HD")
        n += write_headerline(writer.stream, "HD", header["HD"])
    end

    for tag in keys(header)
        if tag == "HD"
            continue
        end
        for h in header[tag]
            n += write_headerline(writer.stream, tag, h)
        end
    end

    return n
end

function write_headerline(io, tag, line)
    n = 0
    n += write(io, '@', tag)
    for (t, val) in line
        checkkeytag(t)
        n += write(io, '\t', t, ':', string(val))
    end
    n += write(io, '\n')
    return n
end

function Base.write(writer::SAMWriter, record::SAMRecord)
    stream = writer.stream

    n = 0
    n += write(
        stream,
        record.name, '\t',
        string(record.flag), '\t',
        record.refname, '\t',
        string(record.pos), '\t',
        string(record.mapq), '\t',
        record.cigar, '\t',
        record.next_refname, '\t',
        string(record.next_pos), '\t',
        string(record.tlen), '\t',
        convert(String, record.seq), '\t',
        [UInt8(q + 33) for q in record.qual])
    n += write_optfields(stream, record.optional_fields)
    n += write(stream, '\n')

    return n
end

function write_optfields(io, optfields)
    n = 0
    for (tag, val) in optfields
        checkkeytag(tag)
        n += write(io, '\t', tag, ':', auxtypechar[typeof(val)], ':', string(val))
    end
    return n
end


text:
@testset "reshape" begin
    test_rrule(reshape, rand(4, 5), (2, 10) ⊢ nothing)
    test_rrule(reshape, rand(4, 5), 2 ⊢ nothing, 10 ⊢ nothing)
end

@testset "hcat" begin
    A = randn(3, 2)
    B = randn(3)
    C = randn(3, 3)
    test_rrule(hcat, A, B, C; check_inferred=false)
end

@testset "reduce hcat" begin
    A = randn(3, 2)
    B = randn(3, 1)
    C = randn(3, 3)
    test_rrule(reduce, hcat ⊢ nothing, [A, B, C])
end

@testset "vcat" begin
    A = randn(2, 4)
    B = randn(1, 4)
    C = randn(3, 4)
    test_rrule(vcat, A, B, C; check_inferred=false)
end

@testset "reduce vcat" begin
    A = randn(2, 4)
    B = randn(1, 4)
    C = randn(3, 4)
    test_rrule(reduce, vcat ⊢ nothing, [A, B, C])
end

@testset "fill" begin
    test_rrule(fill, 44.0, 4 ⊢ nothing; check_inferred=false)
    test_rrule(fill, 2.0, (3, 3, 3) ⊢ nothing)
end


text:
using ChainRulesCore
using Documenter

makedocs(modules=[ChainRulesCore],
         sitename="ChainRulesCore",
         authors="Jarrett Revels and other contributors",
         pages=["Introduction" => "index.md",
                "Getting Started" => "getting_started.md",
                "ChainRulesCore API Documentation" => "api.md"])

deploydocs(repo="github.com/JuliaDiff/ChainRulesCore.jl.git")


text:
"""
    converting(; fromfield, tofield) :: Lens

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> halve(x) = x / 2;

julia> double(x) = 2x;

julia> l = (@lens _.y[2]) ∘ converting(fromfield = halve, tofield = double)
(@lens _.y[2]) ∘ (←double|halve→)

julia> obj = (x=0, y=(1, 2, 3));

julia> @assert get(obj, l) == 1.0 == 2/2

julia> set(obj, l, 0.5)
(x = 0, y = (1, 1.0, 3))
```
"""
converting

"""
    setting(xf::TransformVariables.AbstractTransform) :: Lens

Lens to set value transformed by `xf` (and get value via the inverse
transformation).

# Examples
```jldoctest
julia> using Setfield, Kaleido, TransformVariables

julia> l = (@lens _.y[2]) ∘ setting(as𝕀)
(@lens _.y[2]) ∘ (←|as𝕀→)

julia> obj = (x=0, y=(1, 0.5, 3));

julia> get(obj, l)
0.0

julia> @assert set(obj, l, Inf).y[2] ≈ 1

julia> @assert set(obj, l, -Inf).y[2] ≈ 0.0
```
"""
function setting end

"""
    getting(xf::TransformVariables.AbstractTransform) :: Lens

Lens to get value transformed by `xf` (and set value via the inverse
transformation).
"""
function getting end

abstract type Bijection end

struct FunctionPair{TO, TI} <: Bijection
    fromfield::TO
    tofield::TI
end

struct XFBijection{INV, T} <: Bijection
    xf::T
end

Base.inv(bijection::FunctionPair) =
    FunctionPair(bijection.tofield, bijection.fromfield)

tofield(b::FunctionPair, x) = b.tofield(x)
fromfield(b::FunctionPair, y) = b.fromfield(y)

struct BijectionLens{TB <: Bijection} <: KaleidoLens
    bijection::TB
end

Setfield.get(obj, l::BijectionLens) = fromfield(l.bijection, obj)
Setfield.set(::Any, l::BijectionLens, x) = tofield(l.bijection, x)

BijectionLens(fromfield, tofield) = BijectionLens(FunctionPair(fromfield, tofield))
# TODO: remove `BijectionLens(fromfield, tofield)`

converting(; fromfield, tofield) =
    BijectionLens(FunctionPair(
        prefer_singleton_callable(fromfield),
        prefer_singleton_callable(tofield),
    ))

_setting(thing) = BijectionLens(Bijection(thing))
_getting(thing) = BijectionLens(inv(Bijection(thing)))

Base.show(io::IO, lens::BijectionLens{<:FunctionPair}) =
    print_apply(io, typeof(lens), _getfields(lens.bijection))

# Taken from TransformVariables:
logistic(x::Real) = inv(one(x) + exp(-x))
logit(x::Real) = log(x / (one(x) - x))

logneg(x) = log(-x)
negexp(x) = -exp(x)

"""
    settingasℝ₊ :: BijectionLens

This is a stripped-down version of `setting(asℝ₊)` that works without
TransformVariables.jl.

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> l = (@lens _.y[2]) ∘ settingasℝ₊
(@lens _.y[2]) ∘ (←exp|log→)

julia> obj = (x=0, y=(0, 1, 2));

julia> @assert get(obj, l) == 0.0 == log(obj.y[2])

julia> @assert set(obj, l, -1) == (x=0, y=(0, exp(-1), 2))
```
"""
const settingasℝ₊ = converting(fromfield=log, tofield=exp)
const gettingasℝ₊ = converting(fromfield=exp, tofield=log)

"""
    settingasℝ₋ :: BijectionLens

This is a stripped-down version of `setting(asℝ₋)` that works without
TransformVariables.jl.

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> l = (@lens _.y[2]) ∘ settingasℝ₋
(@lens _.y[2]) ∘ (←negexp|logneg→)

julia> obj = (x=0, y=(0, -1, 2));

julia> @assert get(obj, l) == 0.0 == log(-obj.y[2])

julia> @assert set(obj, l, 1) == (x=0, y=(0, -exp(1), 2))
```
"""
const settingasℝ₋ = converting(fromfield=logneg, tofield=negexp)
const gettingasℝ₋ = converting(fromfield=negexp, tofield=logneg)

"""
    settingas𝕀 :: BijectionLens

This is a stripped-down version of `setting(as𝕀)` that works without
TransformVariables.jl.

# Examples
```jldoctest
julia> using Setfield, Kaleido

julia> l = (@lens _.y[2]) ∘ settingas𝕀
(@lens _.y[2]) ∘ (←logistic|logit→)

julia> obj = (x=0, y=(0, 0.5, 2));

julia> get(obj, l)
0.0

julia> @assert set(obj, l, Inf).y[2] ≈ 1

julia> @assert set(obj, l, -Inf).y[2] ≈ 0
```
"""
const settingas𝕀 = converting(fromfield=logit, tofield=logistic)
const gettingas𝕀 = converting(fromfield=logistic, tofield=logit)


text:
using Test
using FSM
using CSV
using DataFrames

data_force = CSV.File("../data/met_CdP_0506.csv") |> DataFrame

input = Input{Float64}(
    data_force.year,
    data_force.month,
    data_force.day,
    data_force.hour,
    data_force.SW,
    data_force.LW,
    data_force.Sf,
    data_force.Rf,
    data_force.Ta,
    data_force.RH,
    data_force.Ua,
    data_force.Ps,
    )

@testset "complete sim" begin

    path_ref = "../fortran/output/"
    files_ref = readdir(path_ref)

    for file_ref in files_ref

        # Initilize model

        ebm = EBM{Float64}(
            am=parse(Int, file_ref[14]),
            cm=parse(Int, file_ref[15]),
            dm=parse(Int, file_ref[16]),
            em=parse(Int, file_ref[17]),
            hm=parse(Int, file_ref[18]),
            zT=1.5,
            zvar=false,
            Tsoil=[282.98, 284.17, 284.70, 284.70]
        )

        cn = Constants{Float64}()

        snowdepth = similar(input.Ta)
        SWE = similar(input.Ta)
        Tsurf = similar(input.Ta)

        # Run model

        run!(ebm, cn, snowdepth, SWE, Tsurf, input)

        # Read reference data

        data_ref = CSV.File(joinpath(path_ref, file_ref), header=["year", "month", "day", "hour", "alb", "Roff", "snowdepth", "SWE", "Tsurf", "Tsoil"], delim=" ", ignorerepeated=true) |> DataFrame

        # Compute error

        println(file_ref)

        err_snowdepth = snowdepth - data_ref.snowdepth
        @test maximum(abs.(err_snowdepth)) < 0.01
        @info maximum(abs.(err_snowdepth))

        err_swe = SWE - data_ref.SWE
        @test maximum(abs.(err_swe)) < 0.01
        @info maximum(abs.(err_swe))

        err_Tsurf = Tsurf - (data_ref.Tsurf .+ 273.15)
        @test maximum(abs.(err_Tsurf)) < 0.01
        @info maximum(abs.(err_Tsurf))

    end

end


text:
# CUDA errors

const driver_error_descriptions = Dict{Int,ASCIIString}(
    0 => "Success",
    1 => "Invalid value",
    2 => "Out of memory",
    3 => "Driver not initialized",
    4 => "Driver being shutdown",
    5 => "Profiler disabled",
    6 => "Profiler not initialized",
    7 => "Profiler already started",
    8 => "Profiler already stopped",
    100 => "No CUDA-capable device",
    101 => "Invalid device ordinal",
    200 => "Invalid kernel image",
    201 => "Invalid context",
    202 => "Context already current",
    205 => "Map operation failed",
    206 => "Unmap operation failed",
    207 => "Array mapped",
    208 => "Resource already mapped",
    209 => "No kernel image available/suitable for GPU",
    210 => "Resource already acquired",
    211 => "Resource not mapped",
    212 => "Resource not mapped as array",
    213 => "Resource not mapped as pointer",
    214 => "Uncorrectable ECC error detected",
    215 => "Unsupported limit",
    216 => "Context already in use",
    217 => "Peer access not supported",
    218 => "PTX JIT compilation failed",
    300 => "Invalid kernel source",
    301 => "File not found",
    302 => "Shared object symbol not found",
    303 => "Shared object initialization failed",
    304 => "OS call failed",
    400 => "Invalid handle",
    500 => "Named symbol not found",
    600 => "Not ready",
    700 => "Kernel launch failed",
    701 => "Launch out of resources",
    702 => "Launch timeout",
    703 => "Incompatible texturing mode",
    704 => "Peer access already enabled",
    705 => "Peer access not enabled",
    708 => "Primary context already active",
    709 => "Context destroyed",
    710 => "Assertion triggered failure",
    711 => "Too many peers",
    712 => "Host memory already registered",
    713 => "Host memory not registered",
    714 => "Device stack error",
    715 => "Device illegal instruction",
    716 => "Device load or store instruction on unaligned memory address",
    717 => "Device instruction on address not belonging to an allowed space",
    718 => "Device program counter wrapped its address space",
    719 => "Device exception while executing a kernel",
    800 => "Operation not permitted",
    801 => "Operation not supported",
    999 => "Unknown error"
)

immutable CuDriverError
    code::Int
end

description(err::CuDriverError) = driver_error_descriptions[err.code]


text:
module Halite3

import Base: +, -, abs, getindex
import JSON
import Logging
import Random

export end_turn, spawn, move, stay_still, get_ships, is_occupied, init_game,
start_game, update_frame!, is_full, is_empty, random_direction, make_dropoff,
get_ship, random_ship, get_dropoff, get_dropoffs, has_ship, has_structure,
structure_type, mark_unsafe!, normalize, calculate_distance, get_target_direction,
get_unsafe_moves, naive_navigate, get_all_cardinals, get_surrounding_cardinals

include("commands.jl")
include("constants.jl")
include("positionals.jl")
include("game_map.jl")
include("networking.jl")

end # module

max_stars_repo_name: johnnykwwang/Halite-III

max_issues_repo_name: johnnykwwang/Halite-III

max_forks_repo_name: johnnykwwang/Halite-III


text:
"""
    grm(vcffile::AbstractString, method = :GRM, minmaf=0.01, colinds=nothing)

Computes the genetic relationship matrix of a VCF file.

# Inputs
- `vcffile`: VCF file path

# Optional Inputs
- `method`: Method to compute GRM. Can be `:Robust` (default), `:GRM`, or `:MoM`
- `minmaf`: columns (SNPs) with MAF less than `minmaf` are excluded; default 0.01.
- `cinds`: indices or mask of columns to be used for calculating GRM (default
    `nothing`)
- `t`: Float type for calculating GRM; default `Float64`.
- `scale_missing`: If `false`, all missing data is imputed to the mean. If `true`,
    missing data is handled according to this paper:
    `https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6220858/` (default `false`)
"""
function grm(
    vcffile::AbstractString;
    method::Symbol = :Robust,
    minmaf::Real = 0.01,
    cinds::Union{Nothing, AbstractVector{<:Integer}} = nothing,
    t::Type{T} = Float64,
    scale_missing::Bool = false
    ) where T <: Real

    # load genotype into memory. Each row is a sample
    x = convert_gt(t, vcffile, model=:additive, msg="importing genotypes")
    mis, maf, wts, cts = summarize(x)

    colinds = something(cinds, maf .≥ minmaf)

    return @views grm!(x[:, colinds], method, maf[colinds], wts[colinds],
        cts[colinds], mis, scale_missing)
end

"""
    grm!(G, method, maf, wts, cts)

Computes GRM of `G` where `Gᵢ ∈ {0, 1, 2, missing}`. `G` is modified.

If `scale_missing = true`, missing data is handled according to:
`https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6220858/`
"""
function grm!(
    G::AbstractMatrix,
    method::Symbol,
    maf::AbstractVector,
    wts::AbstractVector,
    cts::AbstractVector,
    mis::AbstractVector,
    scale_missing::Bool = false
    )
    m, n = size(G)
    if method == :GRM
        _impute_center_scale!(G, wts, cts, impute=true, center=true, scale=true)
        Φ = G * transpose(G)
        Φ ./= 2n
    elseif method == :MoM
        _impute_center_scale!(G, wts, cts, impute=true)
        G .-= 1
        Φ = G * transpose(G)
        c = sum(x -> abs2(x) + abs2(1 - x), maf)
        shft, scal = n / 2 - c, 1 / (n - c)
        @inbounds @simd for i in eachindex(Φ)
            Φ[i] = (Φ[i] / 2 + shft) * scal
        end
    elseif method == :Robust
        _impute_center_scale!(G, wts, cts, impute=true, center=true)
        scal = sum(x -> 4x * (1 - x), maf)
        Φ = G * transpose(G)
        Φ ./= scal
    else
        throw(ArgumentError("method should be :GRM, :MoM, or :Robust"))
    end

    # scale missing data
    if scale_missing
        @inbounds for i in 1:m
            @simd for j in 1:m
                Φ[j, i] /= ((1 - mis[i]) * (1 - mis[j]))
            end
        end
    end

    return Φ
end

function _impute_center_scale!(
    x::AbstractMatrix,
    wts::AbstractVector,
    cts::AbstractVector;
    impute::Bool=false,
    center::Bool=false,
    scale::Bool=false
    )
    m, n = size(x)
    @inbounds for i in 1:n
        ct = cts[i]
        wt = wts[i]
        for j in 1:m
            impute && ismissing(x[j, i]) && (x[j, i] = ct)
            center && (x[j, i] -= ct)
            scale && wt > 0 && (x[j, i] *= wt)
        end
    end
    return nothing
end

"""
    summarize(x::AbstractMatrix)

Given genotype matrix `x`, calculates missing proportion for each sample
and minor allele frequency, inverse std, and mean for each SNP.

# Inputs
- `X`: Raw genotype matrix, each row is a sample. `X[i, j] ∈ [0, 1, 2, missing]`
"""
function summarize(x::AbstractMatrix)
    n, p = size(x)
    mis = zeros(n) # proportion of missingness SNPs in each sample
    maf = zeros(p) # minor allele frequency
    wts = zeros(p) # inverse standard deviation
    cts = zeros(p) # mean
    @inbounds for i in 1:p
        n00 = n01 = n11 = nmissing = 0
        for j in 1:n
            if ismissing(x[j, i])
                nmissing += 1
                mis[j] += 1
            elseif x[j, i] == 0
                n00 += 1
            elseif x[j, i] == 1
                n01 += 1
            elseif x[j, i] == 2
                n11 += 1
            else
                error("genotype is not 0, 1, 2, or missing!")
            end
        end
        n0 = 2n00 + n01
        n1 = 2n11 + n01
        altfreq = n1 / (n0 + n1)
        reffreq = n0 / (n0 + n1)
        maf[i] = altfreq < reffreq ? altfreq : reffreq
        cts[i] = 2altfreq
        wts[i] = altfreq == 0 ? 1.0 : 1.0 / √(2altfreq * (1 - altfreq))
    end

    mis ./= p

    return mis, maf, wts, cts
end

"""
    robust_GRM_skipmissing(X::AbstractMatrix, maf::AbstractVector)

Calculates robust GRM estimator for genotype matrix `X`, but the inner product
between 2 samples are only performed on SNPs present in both samples.

# Inputs
- `X`: Raw genotype matrix, each row is a sample. `X[i, j] ∈ [0, 1, 2, missing]`
"""
function robust_GRM_skipmissing(X::AbstractMatrix)
    n, p = size(X)
    mis, maf, wts, cts = VCFTools.summarize(X)
    Φ = Matrix{eltype(X)}(undef, n, n)
    @inbounds for i in 1:n, j in 1:n
        Φij = 0.0
        scal = 0.0
        for k in 1:p
            if !ismissing(X[i, k]) && !ismissing(X[j, k])
                scal += 4 * maf[k] * (1 - maf[k])
                Φij += (X[i, k] - cts[k]) * (X[j, k] - cts[k])
            end
        end
        Φ[i, j] = Φij / scal
    end
    return Φ
end


text:
# kelvin helmholtz (wakatani with no alpha terms)
using LinearAlgebra, SparseArrays, ProgressMeter, ForwardDiff, RecipesBase, Plots, Distributed

include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/ImmersedMirrors/mirror.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/ImmersedMirrors/grid.jl")
include("C:/Users/lucas/OneDrive/Documents/GitHub/GDB.jl/src/ImmersedMirrors/extrapolator.jl")

# code taken from test_schemes.jl
function derivative_1d(N, n, m)

    # make stencil
    ic = Int(ceil(m/2.0))
    M = zeros((m, m))
    for i=1:m
        for j=1:m
            M[i, j] = (i - ic)^(j-1) / factorial(j-1)
        end
    end
    stencil = inv(M)[n+1,:]

    D = sparse(zeros(Float64, (N, N)))
    for i=1:m
        k = i - ic
        vec = stencil[i] * ones(N - abs(k))
        D += spdiagm(k => vec)
    end

    return D
end

function mixed_derivative(nx, ny, m, grd)

    operator_to_grid(grd) do
        Dx = derivative_1d(grd._Nx, nx, m)
        Dy = derivative_1d(grd._Ny, ny, m)
        kron(Dy, Dx) / grd.dr^(nx + ny)
    end
end

function Laplacian(m, grd)
    Dxx = mixed_derivative(2, 0, m, grd)
    Dyy = mixed_derivative(0, 2, m, grd)

    return Dxx + Dyy
end

function hyperdiffusion(n, mu, grd)

    operator_to_grid(grd) do
        Dx = derivative_1d(grd._Nx, 2*n, 2*n+1)
        Dy = derivative_1d(grd._Ny, 2*n, 2*n+1)
        Ix = sparse(I, grd._Nx, grd._Nx)
        Iy = sparse(I, grd._Ny, grd._Ny)

        I + (-1)^n*mu*(kron(Iy, Dx) + kron(Dy, Ix))
    end

end

function direct_solve(A, b, bval, extr)

    Anew, bnew = extr(A, b, bval)
    bnew = transpose(A)*b
    Anew = transpose(A)*A
    return extr(Anew \ bnew, bval)
end

function kelvinhelmholtz_update!(w, phi, n, Dx, Dy, w_ex, phi_ex, n_ex, L, Dh, dt)
    # goes through and solves each value twice - once using first two columns, then purely with second column
    bval = zeros(size(w)[1])

    phi_x = Dx*phi[:,2]
    phi_y = Dy*phi[:,2]
    n_x = Dx*n[:,2]
    n_y = Dy*n[:,2]
    w_x = Dx*w[:,2]
    w_y = Dy*w[:,2]

    w_t = -(phi_x.*w_y - phi_y.*w_x)
    # w_t = zeros(size(phi_x))
    n_t = -(phi_x.*n_y - phi_y.*n_x)
    # n_t = zeros(size(phi_x))

    w[:,3] = (w[:,1] + w[:,2])/2 + dt*w_t
    n[:,3] = (n[:,1] + n[:,2])/2 + dt*n_t

    w[:,3] = direct_solve(Dh, w[:,3], bval, w_ex)
    n[:,3] = direct_solve(Dh, n[:,3], bval, n_ex)
    phi[:,3] = direct_solve(L, w[:,3], bval, phi_ex)

    phi_x = Dx*phi[:,3]
    phi_y = Dy*phi[:,3]
    n_x = Dx*n[:,3]
    n_y = Dy*n[:,3]
    w_x = Dx*w[:,3]
    w_y = Dy*w[:,3]

    w_t = -(phi_x.*w_y - phi_y.*w_x)
    n_t = -(phi_x.*n_y - phi_y.*n_x)

    w[:,3] = w[:,2] + dt*w_t
    n[:,3] = n[:,2] + dt*n_t

    w[:,3] = direct_solve(Dh, w[:,3], bval, w_ex)
    n[:,3] = direct_solve(Dh, n[:,3], bval, n_ex)
    phi[:,3] = direct_solve(L, w[:,3], bval, phi_ex)

    w[:,1:2] = w[:,2:3]
    n[:,1:2] = n[:,2:3]
    phi[:,1:2] = phi[:,2:3]

    return w, n, phi
end

#=
Adding parallel derivative of current
Current density given by eq 19 in Ben's GDB paper
jparallel = c/4pi * perpendicular laplacian of psi
How to calculate flux?
- input function
=#

function kelvinhelmholtzcurrent_update!(w, phi, n, psi, j, Dx, Dy, w_ex, phi_ex, n_ex, psi_ex, j_ex, L, Dh, dt, eta)
    # goes through and solves each value twice - once using first two columns, then purely with second column
    bval = zeros(size(w)[1])

    phi_x = Dx*phi[:,2]
    phi_y = Dy*phi[:,2]
    n_x = Dx*n[:,2]
    n_y = Dy*n[:,2]
    w_x = Dx*w[:,2]
    w_y = Dy*w[:,2]
    psi_x = Dx*psi[:,2]
    psi_y = Dy*psi[:,2]
    j_x = Dx*j[:,2]
    j_y = Dy*j[:,2]

    w_t = -(phi_x.*w_y - phi_y.*w_x) + (psi_x.*j_y - psi_y.*j_x)
    n_t = -(phi_x.*n_y - phi_y.*n_x)
    psi_t = -(phi_x.*psi_y - phi_y.*psi_x) + eta*j[:,2]

    w[:,3] = (w[:,1] + w[:,2])/2 + dt*w_t
    n[:,3] = (n[:,1] + n[:,2])/2 + dt*n_t
    psi[:,3] = (psi[:,1] + psi[:,2])/2 + dt*psi_t

    w[:,3] = direct_solve(Dh, w[:,3], bval, w_ex)
    n[:,3] = direct_solve(Dh, n[:,3], bval, n_ex)
    psi[:,3] = direct_solve(Dh, psi[:,3], bval, psi_ex)
    phi[:,3] = direct_solve(L, w[:,3], bval, phi_ex)
    j[:,3] = L*psi[:,3]

    phi_x = Dx*phi[:,3]
    phi_y = Dy*phi[:,3]
    n_x = Dx*n[:,3]
    n_y = Dy*n[:,3]
    w_x = Dx*w[:,3]
    w_y = Dy*w[:,3]
    psi_x = Dx*psi[:,3]
    psi_y = Dy*psi[:,3]
    j_x = Dx*j[:,3]
    j_y = Dy*j[:,3]

    w_t = -(phi_x.*w_y - phi_y.*w_x) + (psi_x.*j_y - psi_y.*j_x)
    n_t = -(phi_x.*n_y - phi_y.*n_x)
    psi_t = -(phi_x.*psi_y - phi_y.*psi_x) + eta*j[:,3]

    w[:,3] = w[:,2] + dt*w_t
    n[:,3] = n[:,2] + dt*n_t
    psi[:,3] = psi[:,2] + dt*psi_t

    w[:,3] = direct_solve(Dh, w[:,3], bval, w_ex)
    n[:,3] = direct_solve(Dh, n[:,3], bval, n_ex)
    psi[:,3] = direct_solve(Dh, psi[:,3], bval, psi_ex)
    phi[:,3] = direct_solve(L, w[:,3], bval, phi_ex)
    j[:,3] = L*psi[:,3]

    w[:,1:2] = w[:,2:3]
    n[:,1:2] = n[:,2:3]
    phi[:,1:2] = phi[:,2:3]
    psi[:,1:2] = psi[:,2:3]
    j[:,1:2] = j[:,2:3]

    return w, n, phi, psi, j
end

# MAIN
m = 5
dt = .01
N = 900
mu = .3
eta = .001

thetas = LinRange(0, 2*pi, 17)[1:16]
c1 = zeros(2, 16)
c2 = zeros(2, 16)
for i=1:16
  c1[:,i] = [cos(thetas[i]), sin(thetas[i])]
  c2[:,i] = .5*[cos(thetas[i]), -sin(thetas[i])]
end
M = Mirror([c1, c2])

r0 = [-1.1, -1.1]
r1 = [1.1, 1.1]
h = .1
grid = Grid(M, h, r0, r1, 8)

n(x,y) = sqrt(x^2+y^2) - 1/3*x^3
n0 = function_to_grid(n, grid)
zros = zeros(grid.Nk)
n0 = hcat(n0, n0, n0)
phi(x,y) = 1/(8*pi^2)*cos(2*pi*x)*sin(2*pi*y)
phi0 = function_to_grid(phi, grid)
phi0 = hcat(phi0, phi0, phi0)

# vorticity determined by initial conditions - choose an easy one (parallel laplacian of phi equals vorticity)
w(x,y) = -cos(2*pi*x)*sin(2*pi*y)
w0 = function_to_grid(w, grid)
w0 = hcat(w0, w0, w0)

j(x,y) = 1/(2*pi)*(x^2+y^2-2)*exp(-1/2*(x^2+y^2))
j0 = function_to_grid(j, grid)
j0 = hcat(j0, j0, j0)
psi(x,y) = 1/(2*pi)*exp(-1/2*(x^2+y^2))
psi0 = function_to_grid(psi, grid)
psi0 = hcat(psi0, psi0, psi0)

println(grid.Nk)

Dx = mixed_derivative(1, 0, m, grid)
Dy = mixed_derivative(0, 1, m, grid)

extrall = Extrapolator(M, grid)

Dh = hyperdiffusion(2, mu, grid)
L = Laplacian(m, grid)

wfin, nfin, phifin, psifin, jfin = kelvinhelmholtzcurrent_update!(w0, phi0, n0, psi0, j0, Dx, Dy, extrall, extrall, extrall, extrall, extrall, L, Dh, dt, eta)

anim = @animate for i=1:N
    wfin, nfin, phifin, psifin, jfin = kelvinhelmholtzcurrent_update!(wfin, phifin, nfin, psifin, jfin, Dx, Dy, extrall, extrall, extrall, extrall, extrall, L, Dh, dt, eta)
    heatmap(nfin[:,1], grid)
    println(i)
end

gif(anim, fps=100)


text:
using Test, MLToolkit
using MLToolkit: _lbeta, _beta

@testset "Special" begin
    n_tests = 5
    sz = (10, 5)

    @testset "lbeta" begin
        for _ = 1:n_tests
            x = randn(FT, sz...).^2
            y = randn(FT, sz...).^2
            @test _lbeta.(AT(x), AT(y)) ≈ lbeta.(x, y)
        end
    end

    @testset "beta" begin
        for _ = 1:n_tests
            x = randn(FT, sz...).^2
            y = randn(FT, sz...).^2
            @test _beta.(AT(x), AT(y)) ≈ beta.(x, y)
        end
    end

    @warn "`logit` is not tested."
    @warn "`logsumexp` is not tested."
end


text:
"""
    MPS(lb, bulk, rb)
constructs an `MPS` where the left-most tensor is the rank-2 tensor lb,
the right-most tensor is the rank-2 tensor rb and in-between are the rank-4 tensors
in the vector `bulk`.
The leg-numbering should be

    (1)                 (2)                (2)              (2)
     ↑         …         ↑                  ↑          …     ↑
    [lb] →(2)   (1)→[bulk[i]]→(3) (1)→[bulk[i+1]]→(3)   (1)→[rb]

The tensors can then be accessed by indexing into the mps.
"""
mutable struct MPS{L, T, TE <: AbstractTensor{T,2}, TB <: AbstractTensor{T,3}} <: AbstractMPS{L,T,TE,TB}
    sites::Vector{Union{TE,TB}}
    MPS(a::TE, b::Vector{TB}, c::TE) where {T, TE <: AbstractTensor{T,2},
        TB <: AbstractTensor{T,3}} = new{length(b)+2,T,TE,TB}([a; b; c])
end

Base.getindex(mps::AbstractMPS{L}, i) where L = mps.sites[i]
Base.setindex!(mps::AbstractMPS{L}, tensor, i) where {L} = setindex!(mps.sites, tensor, i)
"""
    length(mps::AbstractMPS)
return the number of sites in `mps`
"""
Base.length(::AbstractMPS{L}) where L = L
"""
    adjoint(mps::AbstractMPS)
conjugate each tensor in the mps
"""
Base.adjoint(mps::MPS{L,T,TE,TB}) where {L,T,TE,TB} =
    MPS(mps[1]', adjoint.(convert(Array{TB},mps[2:L-1])), mps[L]')

"""
    CanonicalMPS(lb, bulk, rb, cl, zs)
returns a canonical MPS with tensors `lb`, `bulk` and `rb` (see `?MPS`) and
a rank-2 tensor zs (zerosite) between tensor `cl` and `cl+1`.
The index labeling is

    (1)                 (2)                               (2)                 (2)
     ↑         …         ↑                                 ↑          …        ↑
    [lb] →(2)   (1)→[bulk[cl]]→(3)  (1)→[zs]→(2)  (1)→[bulk[cl+1]]→(3)   (1)→[rb]
"""
mutable struct CanonicalMPS{L, T, TE <: AbstractTensor{T,2}, TB <: AbstractTensor{T,3}} <: AbstractMPS{L,T,TE,TB}
    #= CMPS with centerlink = i
      (1)                 (2)                               (2)                 (2)
       ↑         …         ↑                                 ↑          …        ↑
    [ledge] →(2)   (1)→[bulk(i)]→(3)  (1)→[zs]→(2)  (1)→[bulk(i+1)]→(3)   (1)→[redge]
    =#
    sites::Vector{Union{TB,TE}}
    centerlink::Int
    zerosite::TE
    CanonicalMPS(a::TE, b::Vector{TB}, c::TE, cl::Int, zs::TE) where {T, TE <: AbstractTensor{T,2},
        TB <: AbstractTensor{T,3}} = new{length(b)+2,T,TE,TB}([a; b; c], cl, zs)
end

"""
    centerlink(mps::CanonicalMPS)
returns the centerlink `cl` of the canonical MPS `mps` such that the zerosite-tensor is located
between site `cl` and `cl+1`
"""
@inline centerlink(mps::CanonicalMPS) = mps.centerlink
"""
    zerosite(mps::CanonicalMPS)
return the rank-2 zerosite tensor of the canonical MPS `mps`.
Its location can be accessed with `centerlink(mps)`.
"""
@inline zerosite(mps::CanonicalMPS) = mps.zerosite
Base.adjoint(cmps::CanonicalMPS{L,T,TE,TB}) where {L,T,TE,TB} =
    CanonicalMPS(cmps[1]', adjoint.(convert(Array{TB},cmps[2:L-1])), cmps[L]', centerlink(cmps), zerosite(cmps)')

"""
    MPO(lb, bulk, rb)
return an MPO with left-most site rank-3 tensor `lb`, right-most site rank-3 tensor `rb`
and rank-4 bulk tensors in `bulk`.
The labeling is

    (3)          (4)         (3)
     ↑            ↑           ↑
    [lb]→(2) (2)→[b]→(3) (2)→[rb]
     ↑            ↑           ↑
    (1)          (1)         (1)

"""
mutable struct MPO{L, T, TE <: AbstractTensor{T,3}, TB <: AbstractTensor{T,4}} <: AbstractMPO{L,T,TE,TB}
    sites::Vector{Union{TE,TB}}
    MPO{L}(a::TE, b::TB, c::TE) where {T, TE <: AbstractTensor{T,3},
        TB <: AbstractTensor{T,4},L} = new{L,T,TE,TB}([a; fill(b ,L-2); c])
end

@inline Base.getindex(mpo::AbstractMPO{L},i) where L = mpo.sites[i]
@inline Base.setindex!(mpo::AbstractMPO{L}, tensor, i) where L =
    setindex!(mpo.sites, tensor, i)

"""
    todense(mps::AbstractMPS)
apply `todense` to all tensors in mps
"""
function TensorNetworkTensors.todense(mps::MPS{L,T,TE,TB}) where {L,T,TE,TB}
    le = mps[1]
    re = mps[L]
    blk = mps[2:L-1]
    return MPS( todense(le), map(todense,blk), todense(re))
end

"""
    todense(mpo::AbstractMPO)
apply `todense` to all tensors in mpo
"""
function TensorNetworkTensors.todense(mps::MPO{L,T,TE,TB}) where {L,T,TE,TB}
    le = mps[1]
    re = mps[L]
    blk = mps[2]
    return MPO{L}( todense(le), todense(blk), todense(re))
end

"""
    charge(mps::AbstractMPS)
returns the charge of the mps by calculating the charge of each tensor and adding them up.
"""
TensorNetworkTensors.charge(mps::MPS{L}) where L = reduce(⊕,charge.(mps[1:L]))
TensorNetworkTensors.charge(mps::CanonicalMPS{L}) where L =
    reduce(⊕, charge.(mps[1:L])) ⊕ charge(zerosite(mps))

"""
    randmps(N, T, (d, χ))
return an mps with `N` sites of dense tensors with element-type `T`, virtual dimension `χ`
and physical dimension `d`.
"""
randmps(N, T::Type, (d, χ)::Tuple{Int,Int}) =  randmps(N, DTensor{T}((χ,d,χ)))

"""
    randmps(N, T, (sym, (pchs,vchs), (pds, vds))
return an mps with `N` sites of symmetric tensors with symmetries `sym`,
tensors with element-type `T`, physical charges `pchs` with dimensions `pds` and
virtual charges `vchs` with dimensions `vds`.
The tensors legs are oriented as

    1→[b]→3
       ↑
       2
"""
function randmps(N, T::Type, (sym, (pchs, vchs), (pds, vds)))
    a = DASTensor{T,3}(sym, (vchs, pchs, vchs), (vds, pds, vds), InOut(-1,-1,1))
    randmps(N, a)
end

"""
    randmps(N, b::AbstractTensor)
build a random mps where the bulk-tensors have the same structure as `b` and boundaries.
"""
function randmps(N, A::AbstractTensor{T,3}) where T
    if A isa DASTensor
        in_out(A) == InOut(-1,-1,1) || throw(
            ArgumentError("DASTensor for MPS needs InOut(-1,-1,1)"))
    end
    bulk = [similar(A) for i in 2:N-1]
    lb = checked_similar_from_indices(nothing, T, (2,3) , (), A)
    rb = checked_similar_from_indices(nothing, T, (1,2) , (), A)
    initwithrand!.(bulk)
    initwithrand!(lb)
    initwithrand!(rb)
    return MPS(lb, bulk, rb)
end

#= mpscontract to scalar =#
function contractsites(A::T, B::T) where T <: AbstractTensor{<:Any,N} where N
    tensorcontract(A, 1:N, B, 1:N)
end

function contractsites((B,C)::Tuple{AbstractTensor{T,2},AbstractTensor{T,2}},
    A::AbstractTensor{T,2}) where T
    @tensor res = A[1,2] * B[-1,1] * C[-1, 2]
end

function contractsites(A::AbstractTensor{T,2},
    (B,C)::Tuple{AbstractTensor{T,2},AbstractTensor{T,2}}) where T
    @tensor res = A[1,2] * B[1,-1] * C[2, -1]
end

function contractsites((B,C,D)::Tuple{AbstractTensor{T,2},AbstractTensor{T,3},AbstractTensor{T,2}},
    A::AbstractTensor{T,3}) where T
    @tensor res = A[1,2,3] * B[-1,1] * C[-1,2,-2] * D[-2,3]
end

function contractsites(A::AbstractTensor{T,3},
    (B,C,D)::Tuple{AbstractTensor{T,2},AbstractTensor{T,3},AbstractTensor{T,2}}) where T
    @tensor res = A[1,2,3] * B[1,-1] * C[-1,2,-2] * D[3,-2]
end

#= mpscontract to grow edges (rank == 2 || rank == 3) =#
# leftedge rank 3
function contractsites(A::AbstractTensor{T,3},
     (C,D,E)::Tuple{AbstractTensor{T,3}, AbstractTensor{T,4}, AbstractTensor{T,3}}) where T
    @tensor res[1,2,3] := A[-1,-2,-3] * C[-1,-4,1] * D[-4,-2,2,-5] * E[-3,-5,3]
    return res
end

# leftedge rank 2
function contractsites(A::AbstractTensor{T,2},
     (C,D)::Tuple{AbstractTensor{T,3}, AbstractTensor{T,3}}) where T
    @tensor res[1,2] := A[-1,-2] * C[-1,-3,1] * D[-2,-3,2]
    return res
end

#rightedge rank 3
function contractsites((C,D,E)::Tuple{AbstractTensor{T,3}, AbstractTensor{T,4}, AbstractTensor{T,3}},
     A::AbstractTensor{T,3}) where T
    @tensor res[1,2,3] := A[-1,-2,-3] * C[1,-4, -1] * D[-4,2,-2,-5] * E[3,-5,-3]
    return res
end

#rightedge rank 2
function contractsites((C,D)::Tuple{AbstractTensor{T,3}, AbstractTensor{T,3}},
     A::AbstractTensor{T,2}) where T
    @tensor res[1,2] := A[-1,-2] * C[1,-3,-1] * D[2,-3,-2]
    return res
end

#= mpscontract to initialize edges (rank == 2 || rank == 3) =#

#left edge, rank 2
function contractsites(::Tuple{},(A,B)::Tuple{AbstractTensor{T,2},AbstractTensor{T,2}}) where T
    @tensor res[1,2] := A[-1,1] * B[-1,2]
    return res
end
#left edge, rank 3
function contractsites(::Tuple{},(A,B,C)::Tuple{AbstractTensor{T,2},AbstractTensor{T,3},AbstractTensor{T,2}}) where T
    @tensor res[1,2,3] := A[-1,1] * B[-1,2,-2] * C[-2,3]
end
#right edge, rank 2
function contractsites((A,B)::Tuple{AbstractTensor{T,2},AbstractTensor{T,2}},::Tuple{}) where T
    @tensor res[1,2] := A[1,-1] * B[2,-1]
end
#right edge, rank 3
function contractsites((A,B,C)::Tuple{AbstractTensor{T,2},AbstractTensor{T,3},AbstractTensor{T,2}},::Tuple{}) where T
    @tensor res[1,2,3] := A[1,-1] * B[-1,2,-2] * C[3,-2]
end

#= mpscontract to contract neighboring sites =#
#rank 3 with rank 2 on the right
contractsite(A::AbstractTensor{T,3},Bs::AbstractTensor{T,2}...) where T =
    contractsite(A, contractsite(Bs...))

function contractsite(A::AbstractTensor{T,3},B::AbstractTensor{T,2}) where T
    @tensor res[1,2,3] := A[1,2,-1] * B[-1,3]
end

#rank 3 with rank 2 on the left
contractsite(A::AbstractTensor{T,2}, B::AbstractTensor{T,2}, C::AbstractTensor{T,3}) where
    T = contractsite(contractsite(A,B),C)

function contractsite(A::AbstractTensor{T,2}, B::AbstractTensor{T,3}) where T
    @tensor res[1,2,3] := A[1,-1] * B[-1,2,3]
end

#rank 2 with rank 2 on the left
contractsite(A::AbstractTensor{T,2}, Bs::AbstractTensor{T,2}...) where T =
    foldl(contractsite, Bs, init=A)

function contractsite(A::AbstractTensor{T,2},B::AbstractTensor{T,2}) where T
    @tensor res[1,2] := A[1,-1] * B[-1,2]
end

"""
    inner(mps)
return inner(mps,mps)
"""
inner(A) = inner(A,A)

"""
    inner(a, b)
return the inner product between MPS `a` and `b`.
"""
inner(A::MPS{L,T,TE,TB}, B::MPS{L,T,TE,TB}) where {L,T,TE,TB} =
    foldl(contractsites, zip(A[:], B[:]'), init=())

inner(A::MPS{L,T,TE,TB}, O::MPO{L,T}, B::MPS{L,T,TE,TB}) where
    {L,T,TE,TB}  = foldl(contractsites, zip(A[:], O[1:L], B[:]'), init=())
inner(A::CanonicalMPS{L,T,TE,TB}, O::MPO{L,T}, B::MPS{L,T,TE,TB}) where
    {L,T,TE,TB}  = foldl(contractsites, zip(convert(MPS,A)[:], O[1:L], B[:]'), init=())
inner(A::CanonicalMPS{L,T,TE,TB}, O::MPO{L,T}, B::CanonicalMPS{L,T,TE,TB}) where
    {L,T,TE,TB}  = foldl(contractsites, zip(convert(MPS,A)[:], O[1:L], convert(MPS,B)[:]'), init=())
inner(A::MPS{L,T,TE,TB}, O::MPO{L,T}, B::CanonicalMPS{L,T,TE,TB}) where
    {L,T,TE,TB}  = foldl(contractsites, zip(A[:], O[1:L], convert(MPS,B)[:]'), init=())

function Base.convert(::Type{MPS}, A::CanonicalMPS{L,T,TE,TB}, left = true) where {L,T,TE,TB}
    l::TE = A[1]
    b::Vector{TB} = A[2:L-1]
    r::TE = A[L]
    mps = MPS(l, b, r)
    cl = centerlink(A)
    z  = zerosite(A)
    if left
        mps[cl]   = contractsite(mps[cl], z)
    else
        mps[cl+1] = contractsite(z, mps[cl+1])
    end
    return mps
end

inner(A::CanonicalMPS{L,T}, B::MPS{L,T}) where {L,T} = inner(convert(MPS, A), B)
inner(A::MPS{L,T}, B::CanonicalMPS{L,T}) where {L,T} = inner(B, A)

inner(A::CanonicalMPS{L,T}, B::MPO{L,T}, C::MPS{L,T}) where {L,T} = inner(convert(MPS, A), B, C)
inner(A::MPS{L,T}, B::MPO{L,T}, C::CanonicalMPS{L,T}) where {L,T} = inner(A, B, convert(MPS, C))

function inner(A::CanonicalMPS{L,T}, B::CanonicalMPS{L,T}) where {L,T}
    if A == B
        zs = zerosite(A)
        @tensor res[] := zs[1, 2] * zs'[1, 2]
        return scalar(res)
    else
        return inner(convert(MPS, A), convert(MPS, B))
    end
end

function onesitesvd2right(A::AbstractTensor{T,3}) where T
    #    ↑          ↑
    #  →[A]→  ⇒ ( →[U]→ , →[S]→, →[V]→ )
    return tensorsvd(A,((1,2),(3,)))
end

onesitesvd2right(A::AbstractTensor{<:Any,2}) = tensorsvd(A)

onesitesvd2left(A::AbstractTensor{<:Any,3}) = tensorsvd(A, ((1,),(2,3)))

onesitesvd2left(A::AbstractTensor{T,2}) where T = tensorsvd(A)

"""
    leftcanonicalize!(mps, site::Int)
return an `MPS` where all sites to the left of `site` are canonicalized, e.g.

    IN:
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    o-o-o-o-o-o-o, e.g uptolink = 3

    OUT:
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    >->->-ø-o-o-o where > are leftcanonical
                        ø is changed
"""
function leftcanonicalize!(mps::MPS{L,T}, uptolink::Int) where {L,T}
    uptolink < 1 && return mps
    uptolink > L-1 && throw(ArgumentError(
    "can't bring in to leftcanonical-form up to link $uptolink with only $(L-1) links"))

    for i in 1:uptolink
        mps[i], Sm, Vm = onesitesvd2right(mps[i])
        mps[i+1] = contractsite(Sm,Vm,mps[i+1])
    end

    return mps
end

"""
    rightcanonicalize!(mps, site::Int)
return an `MPS` where all sites to the right of `site` are canonicalized, e.g.
    IN:
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    o-o-o-o-o-o-o, e.g downtolink = 3

    OUT:
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    o-o-ø-<-<-<-<, where < is rightcanonicalize
                         ø is changed

"""
function rightcanonicalize!(mps::AbstractMPS{L,T}, downtolink::Int) where {L,T}
    downtolink >= L && return mps
    downtolink < 1 && throw( ArgumentError(
    "can't bring in to rightcanonical-form down to link $downtolink < 1 "))

    for i in reverse(downtolink+1:L)
        Um, Sm, mps[i] = onesitesvd2left(mps[i])
        mps[i-1] = contractsite(mps[i-1], Um, Sm)
    end

    return mps
end

"""
    canonicalize(mps [, link=1])
return a `CanonicalMPS` with centerlink `link`, e.g.

    IN:     ↑ ↑ ↑ ↑ ↑ ↑ ↑
    mps =   o-o-o-o-o-o-o
            link = 3

    Out:    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    mps =   >->->x<-<-<-<

where `x` is the zerosite.
"""
function canonicalize(mps::MPS{L,T,TE,TB}, link = 1) where {L,T,TE,TB}
    mps = deepcopy(mps)
    1 <= link <= (L-1) || throw(
        ArgumentError("link=$link must be between 1 and $(L-1)"))

    leftcanonicalize!(mps, link)
    #=    ↑(link)
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    >->->-ø-o-o-o
    =#
    rightcanonicalize!(mps, link)
    #=    ↑(link)
    ↑ ↑ ↑ ↑ ↑ ↑ ↑
    >->->-ø-<-<-<
    =#
    mps[link], Sm, Vm = onesitesvd2right(mps[link])
    Sm = contractsite(Sm, Vm)
    local l::TE = mps[1]
    local r::TE = mps[L]
    local b::Vector{TB} = mps[2:L-1]

    return CanonicalMPS(l, b, r, link, Sm)::CanonicalMPS{L}
end

"""
    canonicalize!(cmps::CanonicalMPS, link)
move the centerlink/zerosite to `link`.
"""
function canonicalize!(mps::CanonicalMPS{L}, link::Int) where {L}
    1 <= link <= (L-1) || throw(
        ArgumentError("link must be between 1 and $(L-1), is $link"))
    clink = centerlink(mps)
    clink == link && return mps
    while centerlink(mps) != link
        movecenterlink!(mps, link > centerlink(mps))
    end
    return mps::CanonicalMPS{L}
end

canonicalize(mps::CanonicalMPS{L,T}, link) where {L,T} = canonicalize!(deepcopy(mps), link)

function movecenterlink!(mps::CanonicalMPS{L,T}, up::Bool) where {L,T}
    clink = centerlink(mps)
    up && clink == L-1 && throw(
        ArgumentError("cannot move centerlink off the right boundary") )
    !up && clink == 1 && throw(
        ArgumentError("cannot move centerlink off the left boundary") )
    S = zerosite(mps)
    if up
        mps[clink+1], S, V = onesitesvd2right(contractsite(S, mps[clink+1]))
        mps.centerlink += 1
        mps.zerosite = contractsite(S, V)
    else
        U, S, mps[clink] = onesitesvd2left(contractsite(mps[clink], S))
        mps.centerlink -= 1
        mps.zerosite = contractsite(U, S)
    end
    return mps
end

"""
    normalize!(cmps::CanonicalMPS)
normalize `cmps` by normalising its zerosite.
"""
function LinearAlgebra.normalize!(mps::CanonicalMPS)
    fac = sqrt(inner(mps))
    apply!(zerosite(mps), x -> x .= x ./ fac)
    return mps
end

"""
    normalize!(mps::MPS)
normalize `mps` by first canonicalizing it and then normalize the `CanonicalMPS`.
"""
LinearAlgebra.normalize(mps::MPS) = normalize!(canonicalize(mps))
LinearAlgebra.normalize(mps::CanonicalMPS) = normalize!(deepcopy(mps))

"""
    onesitexp(mps, op, site)
return the one-site expectationvalue of the rank-2 tensor `op` applied at site `site`.
"""
function onesiteexp(mps::AbstractMPS{L,T,TE}, op::TE, site) where {T,L,TE}
    mps2 = applyonesite(mps, op, site)
    return inner(mps,mps2)
end

"""
    onesitevar(mps, op, site)
return the variance of the operator `op` at `site`, i.e.

    <mps|op^2|mps> - <mps|op|mps>^2
"""
function onesitevar(mps::AbstractMPS{L,T,TE}, op::TE, site) where {T,L,TE}
    mps1 = applyonesite(mps, op, site)
    mps2 = applyonesite(mps1, op, site)
    return inner(mps2,mps) - inner(mps1,mps)^2
end

"""
    twositecorr(mps, (op1,op2), (site1,sit2))
return the two-site correlator with `op1` applied at `site1`
and `op2` applied at `site2`.
"""
function twositecorr(mps::AbstractMPS{L,T,TE},
        (op1,op2)::NTuple{2,TE},
        (site1,site2)::NTuple{2,Int}) where {L,T,TE}
     mps1 = applyonesite(mps, op1, site1)
     applyonesite!(mps1, op2, site2)
     return inner(mps1, mps)
end

function applyonesite!(mps::AbstractMPS{L,T,TE}, op::TE, site::Int) where {L,T,TE}
    a = mps[site]
    if site == 1
        @tensor a[1,2] := a[-1,2] * op[-1,1]
    elseif site == L
        @tensor a[1,2] := a[1,-1] * op[-1,2]
    else
        @tensor a[1,2,3] := a[1,-1,3] * op[-1,2]
    end
    mps[site] = a
    return mps
end

applyonesite(mps, op, site) = applyonesite!(deepcopy(mps), op, site)

"""
    entropy(mps::AbstractMPS, link)
return the entropy of entanglement accross link `link` connecting site
`link` to site `link+1`.
"""
function entropy(mps::AbstractMPS, link::Int)
    zs = zerosite(canonicalize(mps, link))
    Σ = tensorsvd(zs)[2]
    schmidtvals = abs.(diag(Σ))
    sv2 = schmidtvals.^2
    return -(sv2' * log.(sv2))
end


text:
using BenchmarkTools, McCormick, Random
using Statistics: std

# Note: Run time is still random so results are't strictly reproducable...
Random.seed!(1234)

ex_softplus(x) = log(1 + exp(x))
ex_maxsig(x) =  max(x, 1/(1 + exp(-(x))))
ex_maxtanh(x) = max(x,tanh(x))
ex_softsign(x) = (x)/(1 + x)
ex_sigmoid(x) = 1/(1 + exp(-(x)))
ex_silu(x) = (x)/(1 + exp(-(x)))
ex_gelu(x) = (x)*(1 + erf((x)/sqrt(2)))/2

const TRIAL_REPEAT = 10000

function loop(func, x)
    foreach(i -> func(x), 1:TRIAL_REPEAT)
end
function summarize!(n, f_trial, ex_trial)
    env_min = median(f_trial.times)/TRIAL_REPEAT
    expr_min = median(ex_trial.times)/TRIAL_REPEAT
    ratio_percent = 100*(env_min/expr_min - 1)
    println("$n | env = $(env_min) | exp = $(expr_min) | ratio = $(ratio_percent)")
end
#start_mc() = MC{1,NS}(0.0, Interval(-3,3), 1)
start_mc() = MC{1,NS}(6*rand()-3, Interval(-3,3), 1)

f_trial = @benchmark loop(softplus,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_softplus,x) setup=(x=start_mc())
summarize!("Softplus", f_trial, ex_trial)

f_trial = @benchmark loop(maxsig,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_maxsig,x) setup=(x=start_mc())
summarize!("Maxsig", f_trial, ex_trial)

f_trial = @benchmark loop(maxtanh,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_maxtanh,x) setup=(x=start_mc())
summarize!("MaxTanh", f_trial, ex_trial)

f_trial = @benchmark loop(softsign,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_softsign,x) setup=(x=start_mc())
summarize!("Softsign", f_trial, ex_trial)

f_trial = @benchmark loop(sigmoid,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_sigmoid,x) setup=(x=start_mc())
summarize!("Sigmoid", f_trial, ex_trial)

f_trial = @benchmark loop(swish,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_silu,x) setup=(x=start_mc())
summarize!("Silu", f_trial, ex_trial)

f_trial = @benchmark loop(gelu,x) setup=(x=start_mc())
ex_trial = @benchmark loop(ex_gelu,x) setup=(x=start_mc())
summarize!("Gelu", f_trial, ex_trial)

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
dir = @__DIR__
"""
genplain(s)

Small helper function to run some code and redirect the output (stdout) to a file.
"""
function genplain(s::String)
    open(joinpath(dir, "output", "$(splitext(s)[1]).out"), "w") do outf
        redirect_stdout(outf) do
            include(joinpath(dir, s))
        end
    end
end
# run `script1.jl` and redirect what it prints to `output/script1.out`

genplain("mandP.jl")
include("vismandP.jl")


text:
"""
    ReplaySaver()

ReplaySaver is used to save/load SLAM result for the visualizer.
If added to the SlamManager instead of the visualizer, it accumulates
poses needed for the visualization/replay and can be later used
to replay the results.

# Usage:

To use it, pass it to `visualizer` keyword argument, when creating `SlamManager`.
Do some work, and in the end, save it to the directory.

```julia
saver = ReplaySaver()
SlamManager(params, camera; visualizer=saver)
# ...run slam...
SLAM.save(saver, "save_dir")
```

Then you can load it from the save directory and use it to replay results.

```julia
SLAM.load!(saver, "save_dir")
```
"""
mutable struct ReplaySaver <: SLAMIO
    ids::Dict{Int64, Int64} # frame id -> position id
    positions::Vector{Point3f0}

    positions_lock::ReentrantLock
end

ReplaySaver() = ReplaySaver(Dict{Int64, Int64}(), Point3f0[], ReentrantLock())

"""
    set_frame_wc!(saver::ReplaySaver, frame_id, wc)

Add new pose to store.
"""
function set_frame_wc!(saver::ReplaySaver, frame_id, wc)
    lock(saver.positions_lock) do
        base_position = SVector{4, Float64}(0, 0, 0, 1)
        position = (wc * base_position)[[1, 3, 2]]

        pid = get(saver.ids, frame_id, -1)
        if pid == -1
            push!(saver.positions, position)
            saver.ids[frame_id] = length(saver.positions)
        else
            saver.positions[pid] = position
        end
    end
end

"""
    save(saver::ReplaySaver, save_dir)

Save results into the given directory.
"""
function save(saver::ReplaySaver, save_dir)
    isdir(save_dir) || mkdir(save_dir)
    positions_file = joinpath(save_dir, "positions.bson")
    ids_file = joinpath(save_dir, "ids.bson")

    positions = saver.positions
    @save positions_file positions

    ids = saver.ids
    @save ids_file ids
end

"""
    load!(saver::ReplaySaver, save_dir)

Load results from a given directory.
"""
function load!(saver::ReplaySaver, save_dir)
    isdir(save_dir) || error("Directory `$save_dir` does not exist.")

    positions_file = joinpath(save_dir, "positions.bson")
    isfile(positions_file) || error("Positions file `$positions_file` not found.")

    ids_file = joinpath(save_dir, "ids.bson")
    isfile(ids_file) || error("Ids file `$ids_file` not found.")

    @load positions_file positions
    @load ids_file ids

    saver.ids = ids
    saver.positions = positions
end


text:
addprocs(2)
# Visualize SONAR with Director

using TransformUtils
using KernelDensityEstimate
using IncrementalInference
using RoME
using Caesar
using DrakeVisualizer, GeometryTypes
using ColorTypes: RGBA
using MeshIO, FileIO
using CoordinateTransformations
using Rotations

# load necessary data
rovt = loadmodel(:rov)

# true position
gt = Dict{Symbol, Vector{Float64}}()
gt[:l1] = [4.0;0.0;-0.7]

# start visualizer and models
vc = startdefaultvisualization()
rovt(vc)

# and some default object in the world
boxdata = GeometryData(HyperRectangle(Vec(4.0,0.0,-0.7), Vec(5.0,5.0,1.4)))
boxdata.color = RGBA(0.5,0.1,0.0,0.5)
model = Visualizer(boxdata,100)

# l1p = [4.0; 0; 0.7]
# featuredata = GeometryData(HyperSphere(Point(l1p...), 0.15) )
# featuredata.color = RGBA(0.0, 1.0, 0.0, 1.0)
# model = Visualizer(featuredata,101)

# l1p = [4.0; 0; -0.7]
# featuredata = GeometryData(HyperSphere(Point(l1p...), 0.15) )
# featuredata.color = RGBA(0.0, 0.0, 1.0, 1.0)
# model = Visualizer(featuredata,102)

@async begin
  rho = 5.0
  pnt = (4.5,0.0)
  for t in linspace(0,1,100)
    the = t*(pi/2+pi/6)
    x = -rho*cos(the)+pnt[1]
    y = -rho*sin(the)
    rovt(vc, Translation(x,y,0) ∘ LinearMap(Rotations.AngleAxis(the,0,0,1.0)) )
    # draw(modelrov, [Translation(x, y, 0.25) ∘ LinearMap(Rotations.AngleAxis(pi+the,0,0,1.0))])
    sleep(0.05)
  end
  for t in linspace(0,1,25)
    the = (pi/2+pi/6)
    themax = pi/2+pi/6
    xmax = -rho*cos(themax)+pnt[1]
    ymax = -rho*sin(themax)
    z = t
    rovt(vc, Translation(xmax,ymax,-z) ∘ LinearMap(Rotations.AngleAxis(themax,0,0,1.0)) )
    # draw(modelrov, [Translation(xmax, ymax, 0.25-z) ∘ LinearMap(Rotations.AngleAxis(pi+themax,0,0,1.0))])
    sleep(0.05)
  end
  for t in linspace(1,0.25,75)
    the = t*(pi/2+pi/6)
    themax = pi/2+pi/6
    x = -rho*cos(the)+pnt[1]
    y = -rho*sin(the)
    xmax = -rho*cos(themax)+pnt[1]
    ymax = -rho*sin(themax)
    rovt(vc, Translation(x,y,-1.0) ∘ LinearMap(Rotations.AngleAxis(the,0,0,1.0)) )
    # draw(modelrov, [Translation(x, y, -0.75) ∘ LinearMap(Rotations.AngleAxis(pi+the,0,0,1.0))])
    sleep(0.05)
  end
end

# covariance parameters
initCov = 0.01*eye(6); [initCov[i,i] = 0.001 for i in 4:6];
odoCov, rangecov, bearingcov = deepcopy(initCov), 3e-4, 2e-4

# okay build the graph
fg = identitypose6fg(initCov=initCov)

println("Adding :l1 LinearRangeBearingElevation to graph...")
addLinearArrayConstraint(fg, (4.12, 0.0), :x1, :l1, rangecov=rangecov,bearingcov=bearingcov)

addOdoFG!(fg, Pose3Pose3(MvNormal([2.0;-4.0;0.0;0.0;0.0;pi/3], odoCov) )
addOdoFG!(fg, Pose3Pose3(MvNormal([2.0;-4.0;0.0;0.0;0.0;pi/3], odoCov) )
visualizeallposes!(vc, fg, drawlandms=false)

rho2, bearing2, elev2 = projectrbe(fg,:x2,gt[:l1])
addLinearArrayConstraint(fg, (rho2, bearing2), :x2, :l1, rangecov=rangecov,bearingcov=bearingcov)

rho3, bearing3, elev3 = projectrbe(fg,:x3,gt[:l1])
addLinearArrayConstraint(fg, (rho3, bearing3), :x3, :l1, rangecov=rangecov,bearingcov=bearingcov)

visualizeDensityMesh!(vc, fg, :l1, meshid=2)

addOdoFG!(fg, Pose3Pose3(MvNormal([1.0;2.0;-1.0;0.0;0.0;-pi/6], odoCov) )
visualizeallposes!(vc, fg, drawlandms=false)

rho4, bearing4, elev4 = projectrbe(fg,:x4,gt[:l1])
bearing4 = 0.0
addLinearArrayConstraint(fg, (rho4, bearing4), :x4, :l1, rangecov=rangecov,bearingcov=bearingcov)

addOdoFG!(fg, Pose3Pose3(MvNormal([2.0;4.0;0.0;0.0;0.0;-pi/3], odoCov) )

visualizeallposes!(vc, fg, drawlandms=false)
rho5, bearing5, elev5 = projectrbe(fg,:x5,:l1)
bearing5 = 0.0
addLinearArrayConstraint(fg, (rho5, bearing5), :x5, :l1, rangecov=rangecov,bearingcov=bearingcov)

batchSolve(fg)

# visualize(fg, vc, densitymeshes=[:l1;:x5])
visualizeallposes!(vc, fg, drawlandms=false)

visualizeDensityMesh!(vc, fg, :x4, meshid=3)

#
# rovdata = GeometryData(rov )
# rovdata.color = RGBA(0., 1.0, 0.5, 0.3)
# modelrov = Visualizer(rovdata, 99)
# rho = 5.0
# pnt = (4.5,0.0)
# for t in linspace(0,1,100)
#   the = t*(pi/2+pi/6)
#   x = -rho*cos(the)+pnt[1]
#   y = -rho*sin(the)
#   draw(modelrov, [Translation(x, y, 0.25) ∘ LinearMap(Rotations.AngleAxis(pi+the,0,0,1.0))])
#   sleep(0.05)
# end

@show

#


text:
struct Put <: OptionType end
struct Call <: OptionType end

value(::Put) = -1
value(::Call) = 1

struct PlainVanillaPayoff{OT <: OptionType} <: StrikedTypePayoff
  optionType::OT
  strike::Float64
end

(payoff::PlainVanillaPayoff)(price::Float64) = _get_payoff(payoff.optionType, payoff, price)

_get_payoff(::Call, payoff::PlainVanillaPayoff, price::Float64) = max(price - payoff.strike, 0.0)
_get_payoff(::Put, payoff::PlainVanillaPayoff, price::Float64) = max(payoff.strike - price, 0.0)

struct ForwardTypePayoff{P <: PositionType} <: AbstractPayoff
  position::P
  strike::Float64
end


text:
@testset "Join" begin

    left = GorFile("left.gor")
    right = GorFile("right.gor")

    j = GOR.join(left, right)

    @test left |> GOR.join(right) |> collect |> length == 100

    @test eltype(j) == NamedTuple{(:Chrom, :Pos, :Val, :Chromx, :Posx, :Valx),Tuple{String,Int64,String,String,Int64,String}}

    using DataFrames, Tables

    l = Tables.namedtupleiterator(DataFrame(
        Chrom = repeat(["chr1"], 11),
        bpStart = collect(0:10:100),
        bpEnd = collect(10:10:110)))

    r = Tables.namedtupleiterator(DataFrame(
        Chrom = ["chr1", "chr1", "chr1"],
        bpStart = [15, 55, 75],
        bpEnd = [45, 85, 130]))

    @test r |> GOR.join(l; leftjoin=false, kind = :segsnp) |>
        DataFrame |> size == (9, 6)

    @test r |> GOR.join(l; leftjoin = false, kind = :segseg) |>
        DataFrame |> size == (12, 6)

    @test l |> GOR.join(r; leftjoin = false, kind = :segseg) |>
        DataFrame |> size == (12, 6)

    @test l |> GOR.join(r) |> eltype ==
        NamedTuple{(:Chrom, :bpStart, :bpEnd, :Chromx, :bpStartx, :bpEndx),Tuple{String,Int64,Int64,String,Int64,Int64}}

    @test l |> GOR.join(r; leftjoin = true) |> eltype ==
        NamedTuple{(:Chrom, :bpStart, :bpEnd, :Chromx, :bpStartx, :bpEndx),Tuple{String,Int64,Int64,Union{Missing, String},Union{Missing, Int64},Union{Missing, Int64}}}

end


text:
@doc raw"""
`LDAE`: Lagrangian Differential Algebraic Equation

Defines a Lagrangian differential algebraic initial value problem, that is
a special implicit initial value problem
```math
\begin{aligned}
\dot{q} (t) &= v(t) + \lambda(t), &
q(t_{0}) &= q_{0} , \\
\dot{p} (t) &= f(t, q(t), v(t)) + g(t, q(t), \lambda(t)) + \bar{g} (t, q(t), \mu(t)) , &
p(t_{0}) &= p_{0} , \\
p(t) &= ϑ(t, q(t), v(t)) , \\
0 &= \phi (t, q(t), p(t)) , \\
0 &= \psi (t, q(t), p(t), \dot{q}(t), \dot{p}(t)) ,
\end{aligned}
```
with vector field ``f``, the momentum defined by ``p``, initial conditions ``(q_{0}, p_{0})`` and the solution
``(q,p)`` taking values in ``\mathbb{R}^{d} \times \mathbb{R}^{d}`` and
the algebraic variables ``(v, \lambda, \mu)`` taking values in
``\mathbb{R}^{d} \times \mathbb{R}^{m} \times \mathbb{R}^{m}``.
This is a special case of a differential algebraic equation with dynamical
variables ``(q,p)`` and algebraic variables ``v``, ``\lambda`` and ``\mu``.

### Parameters

* `DT <: Number`: data type
* `TT <: Real`: time step type
* `AT <: AbstractArray{DT}`: array type
* `ϑType <: Function`: type of `ϑ`
* `fType <: Function`: type of `f`
* `uType <: Function`: type of `u`
* `gType <: Function`: type of `g`
* `ϕType <: Function`: type of `ϕ`
* `ūType <: Function`: type of `ū`
* `ḡType <: Function`: type of `ḡ`
* `ψType <: Function`: type of `ψ`
* `ωType <: Function`: type of `ω`
* `v̄Type <: Function`: type of `v̄`
* `f̄Type <: Function`: type of `f̄`
* `lagType <: Function`: Lagrangian type
* `invType <: OptionalNamedTuple`: invariants type
* `parType <: OptionalNamedTuple`: parameters type
* `perType <: OptionalArray{AT}`: periodicity type

### Fields

* `d`: dimension of dynamical variables ``q`` and ``p`` as well as the vector fields ``f`` and ``p``
* `ϑ`: function determining the momentum
* `f`: function computing the vector field
* `u`: function computing the projection for ``q``, for a degenerate system given by ``\lambda``
* `g`: function computing the projection for ``p``, for a degenerate system given by ``\nabla \vartheta (q) \cdot \lambda``
* `ϕ`: primary constraints, for a degenerate system given by ``p - \vartheta (q)``
* `ū`: function computing the secondary projection field ``\bar{u}``, for a degenerate system given by ``\lambda`` (optional)
* `ḡ`: function computing the secondary projection field ``\bar{g}``, for a degenerate system given by ``\lambda \cdot \nabla \vartheta (q)`` (optional)
* `ψ`: secondary constraints, for a degenerate system given by ``\dot{p} - \dot{q} \cdot \nabla \vartheta (q)`` (optional)
* `ω`: function computing the symplectic matrix
* `v̄`: function computing an initial guess for the velocity field ``v`` (optional)
* `f̄`: function computing an initial guess for the force field ``f`` (optional)
* `t₀`: initial time (optional)
* `q₀`: initial condition for dynamical variable `q`
* `p₀`: initial condition for dynamical variable `p`
* `λ₀`: initial condition for algebraic variable `λ`
* `μ₀`: initial condition for algebraic variable `μ` (optional)
* `lagrangian`: function computing the Lagrangian ``L``
* `invariants`: either a `NamedTuple` containing the equation's invariants or `nothing`
* `parameters`: either a `NamedTuple` containing the equation's parameters or `nothing`
* `periodicity`: determines the periodicity of the state vector `q` for cutting periodic solutions

The functions `ϑ` and `f` must have the interface
```julia
    function ϑ(t, q, v, p)
        p[1] = ...
        p[2] = ...
        ...
    end
```
and
```julia
    function f(t, q, v, f)
        f[1] = ...
        f[2] = ...
        ...
    end
```
where `t` is the current time, `q` is the current solution vector, `v` is the
current velocity and `f` and `p` are the vectors which hold the result of
evaluating the functions ``f`` and ``ϑ`` on `t`, `q` and `v`.
The funtions `g`, `v̄` and `f̄` are specified by
```julia
    function g(t, q, λ, g)
        g[1] = ...
        g[2] = ...
        ...
    end

    function v̄(t, q, v)
        v[1] = ...
        v[2] = ...
        ...
    end

    function f̄(t, q, v, f)
        f[1] = ...
        f[2] = ...
        ...
    end
```

### Constructors

```julia
LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, t₀, q₀, p₀, λ₀, μ₀, lagrangian, invariants, parameters, periodicity)

LDAE(ϑ, f, u, g, ϕ, l, ω, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, l, ω, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, l, ω, t₀, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, l, ω, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)

LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...)
```

### Keyword arguments

* `v̄ = (t,q,v) -> nothing`
* `f̄ = f`
* `invariants = nothing`
* `parameters = nothing`
* `periodicity = nothing`

"""
struct LDAE{ϑType <: Callable, fType <: Callable,
            uType <: Callable, gType <: Callable, ϕType <: Callable,
            ūType <: OptionalCallable, ḡType <: OptionalCallable, ψType <: OptionalCallable,
            ωType <: Callable, v̄Type <: Callable, f̄Type <: Callable,
            lagType <: Callable,
            invType <: OptionalInvariants,
            parType <: OptionalParameters,
            perType <: OptionalPeriodicity} <: AbstractEquationPDAE{invType,parType,perType,ψType}

    ϑ::ϑType
    f::fType
    u::uType
    g::gType
    ϕ::ϕType
    ū::ūType
    ḡ::ḡType
    ψ::ψType
    ω::ωType

    v̄::v̄Type
    f̄::f̄Type

    lagrangian::lagType
    invariants::invType
    parameters::parType
    periodicity::perType

    function LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, lagrangian, invariants, parameters, periodicity)
        @assert !isempty(methods(ϑ))
        @assert !isempty(methods(f))
        @assert !isempty(methods(u))
        @assert !isempty(methods(g))
        @assert !isempty(methods(ϕ))
        @assert !isempty(methods(ū)) || ū === nothing
        @assert !isempty(methods(ḡ)) || ḡ === nothing
        @assert !isempty(methods(ψ)) || ψ === nothing
        @assert !isempty(methods(ω))
        @assert !isempty(methods(v̄))
        @assert !isempty(methods(f̄))
        @assert !isempty(methods(lagrangian))

        new{typeof(ϑ), typeof(f),
            typeof(u), typeof(g), typeof(ϕ),
            typeof(ū), typeof(ḡ), typeof(ψ),
            typeof(ω), typeof(v̄), typeof(f̄),
            typeof(lagrangian), typeof(invariants), typeof(parameters), typeof(periodicity)}(
                ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, lagrangian, invariants, parameters, periodicity)
    end
end

_ldae_default_v̄(t,q,v,params) = nothing

LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, lagrangian; v̄=_ldae_default_v̄, f̄=f, invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity()) = LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, lagrangian, invariants, parameters, periodicity)
LDAE(ϑ, f, u, g, ϕ, ω, lagrangian; kwargs...) = LDAE(ϑ, f, u, g, ϕ, nothing, nothing, nothing, ω, lagrangian; kwargs...)

GeometricBase.invariants(equation::LDAE) = equation.invariants
GeometricBase.parameters(equation::LDAE) = equation.parameters
GeometricBase.periodicity(equation::LDAE) = equation.periodicity

hasvectorfield(::LDAE) = true
haslagrangian(::LDAE) = true

function check_initial_conditions(equ::LDAE, ics::NamedTuple)
    haskey(ics, :q) || return false
    haskey(ics, :v) || haskey(ics, :p) || return false
    haskey(ics, :λ) || return false
    if hassecondary(equ)
        haskey(ics, :μ) || return false
        eltype(ics.λ) == eltype(ics.μ) || return false
        typeof(ics.λ) == typeof(ics.μ) || return false
        axes(ics.λ) == axes(ics.μ) || return false
    end
    return true
end

function check_methods(equ::LDAE, tspan, ics::NamedTuple, params)
    applicable(equ.ϑ, tspan[begin], ics.q, zero(ics.q), zero(ics.p), params) || return false
    applicable(equ.f, tspan[begin], ics.q, zero(ics.q), zero(ics.p), params) || return false
    applicable(equ.ϕ, tspan[begin], ics.q, ics.p, ics.λ, zero(ics.λ), params) || return false
    # TODO add missing methods
    return true
end

function datatype(equ::LDAE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return eltype(ics.q)
end

function arrtype(equ::LDAE, ics::NamedTuple)
    @assert check_initial_conditions(equ, ics)
    return typeof(ics.q)
end

_get_ϑ(equ::LDAE, params) = (t,q,v,ϑ)     -> equ.ϑ(t, q, v, ϑ, params)
_get_f(equ::LDAE, params) = (t,q,v,f)     -> equ.f(t, q, v, f, params)
_get_u(equ::LDAE, params) = (t,q,p,λ,u)   -> equ.u(t, q, p, λ, u, params)
_get_g(equ::LDAE, params) = (t,q,p,λ,g)   -> equ.g(t, q, p, λ, g, params)
_get_ϕ(equ::LDAE, params) = (t,q,p,ϕ)     -> equ.ϕ(t, q, p, ϕ, params)
_get_ū(equ::LDAE, params) = (t,q,p,λ,u)   -> equ.ū(t, q, p, λ, u, params)
_get_ḡ(equ::LDAE, params) = (t,q,p,λ,g)   -> equ.ḡ(t, q, p, λ, g, params)
_get_ψ(equ::LDAE, params) = (t,q,p,v,f,ψ) -> equ.ψ(t, q, p, v, f, ψ, params)
_get_v̄(equ::LDAE, params) = (t,q,p,v)     -> equ.v̄(t, q, p, v, params)
_get_f̄(equ::LDAE, params) = (t,q,p,f)     -> equ.f̄(t, q, p, f, params)
_get_ω(equ::LDAE, params) = (t,q,v,ω)     -> equ.ω(t, q, v, ω, params)
_get_l(equ::LDAE, params) = (t,q,v)       -> equ.lagrangian(t, q, v, params)

function _functions(equ::LDAE)
    if hassecondary(equ)
        (ϑ = equ.ϑ, f = equ.f, u = equ.u, g = equ.g, ϕ = equ.ϕ, ū = equ.ū, ḡ = equ.ḡ, ψ = equ.ψ, v̄ = equ.v̄, f̄ = equ.f̄, ω = equ.ω, l = equ.lagrangian)
    else
        (ϑ = equ.ϑ, f = equ.f, u = equ.u, g = equ.g, ϕ = equ.ϕ, v̄ = equ.v̄, f̄ = equ.f̄, ω = equ.ω, l = equ.lagrangian)
    end
end

function _functions(equ::LDAE, params::OptionalParameters)
    if hassecondary(equ)
        (ϑ = _get_ϑ(equ, params), f = _get_f(equ, params), u = _get_u(equ, params), g = _get_g(equ, params), ϕ = _get_ϕ(equ, params), ū = _get_ū(equ, params), ḡ = _get_ḡ(equ, params), ψ = _get_ψ(equ, params), v̄ = _get_v̄(equ, params), f̄ = _get_f̄(equ, params), ω = _get_ω(equ, params), l = _get_l(equ, params))
    else
        (ϑ = _get_ϑ(equ, params), f = _get_f(equ, params), u = _get_u(equ, params), g = _get_g(equ, params), ϕ = _get_ϕ(equ, params), v̄ = _get_v̄(equ, params), f̄ = _get_f̄(equ, params), ω = _get_ω(equ, params), l = _get_l(equ, params))
    end
end

# struct LDAE{dType <: Number, tType <: Real, arrayType <: AbstractArray{dType},
#             ϑType <: Function, fType <: Function,
#             uType <: Function, gType <: Function, ϕType <: Function,
#             ūType <: OptionalFunction, ḡType <: OptionalFunction, ψType <: OptionalFunction,
#             v̄Type <: Function, f̄Type <: Function,

#             invType <: OptionalInvariants,
#             parType <: OptionalParameters,
#             perType <: OptionalPeriodicity} <: AbstractEquationPDAE

#     d::Int
#     m::Int

#     ϑ::ϑType
#     f::fType
#     u::uType
#     g::gType
#     ϕ::ϕType
#     ū::ūType
#     ḡ::ḡType
#     ψ::ψType
#     ω::ωType
#     v̄::v̄Type
#     f̄::f̄Type

#     t₀::tType
#     q₀::Vector{arrayType}
#     p₀::Vector{arrayType}
#     λ₀::Vector{arrayType}
#     μ₀::Vector{arrayType}

#     ω::ωType
#     lagrangian::lagType
#     invariants::invType
#     parameters::parType
#     periodicity::perType

#     function LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄,
#                   t₀::tType,
#                   q₀::Vector{arrayType}, p₀::Vector{arrayType},
#                   λ₀::Vector{arrayType}, μ₀::Vector{arrayType},
#                   lagrangian, invariants, parameters, periodicity) where {
#                         dType <: Number, tType <: Real, arrayType <: AbstractArray{dType}}

#         d = length(q₀[begin])
#         m = length(λ₀[begin])

#         @assert 2d ≥ m

#         @assert length(q₀) == length(p₀) == length(λ₀) == length(μ₀)

#         @assert all(length(q) == d for q in q₀)
#         @assert all(length(p) == d for p in p₀)
#         @assert all(length(λ) == m for λ in λ₀)
#         @assert all(length(μ) == m for μ in μ₀)

#         @assert all([ndims(q) == ndims(p) == ndims(λ) == ndims(μ) for (q,p,λ,μ) in zip(q₀,p₀,λ₀,μ₀)])

#         new{dType, tType, arrayType,
#             typeof(ϑ), typeof(f),
#             typeof(u), typeof(g), typeof(ϕ),
#             typeof(ū), typeof(ḡ), typeof(ψ),
#             typeof(ω), typeof(v̄), typeof(f̄),
#             typeof(lagrangian), typeof(invariants), typeof(parameters), typeof(periodicity)}(
#                 d, m, ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, t₀, q₀, p₀, λ₀, μ₀, lagrangian, invariants, parameters, periodicity)
#     end
# end

# _LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, lagrangian, ω, t₀, q₀, p₀, λ₀, μ₀; invariants=NullInvariants(), parameters=NullParameters(), periodicity=NullPeriodicity(), v̄=(parameters === nothing ? (t,q,v)->nothing : (t,q,v,params)->nothing), f̄=f) = LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, ω, v̄, f̄, t₀, q₀, p₀, λ₀, μ₀, lagrangian, invariants, parameters, periodicity)

# LDAE(ϑ, f, u, g, ϕ, l, ω, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = _LDAE(ϑ, f, u, g, ϕ, nothing, nothing, nothing, l, ω, t₀, q₀, p₀, λ₀, μ₀; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, l, ω, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, l, ω, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, l, ω, t₀, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, l, ω, t₀, [q₀], [p₀], [λ₀], [μ₀]; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, l, ω, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, l, ω, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)

# LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = _LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, q₀, p₀, λ₀, μ₀; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, t₀, [q₀], [p₀], [λ₀], [μ₀]; kwargs...)
# LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, q₀::State, p₀::State, λ₀::State, μ₀::State=zero(λ₀); kwargs...) = LDAE(ϑ, f, u, g, ϕ, ū, ḡ, ψ, l, ω, 0.0, q₀, p₀, λ₀, μ₀; kwargs...)

# const LDAEpsiType{psiT,DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,ΩT,V̄T,F̄T,lagT,invT,parT,perT} = LDAE{DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,parT,perT} # type alias for dispatch on secondary constraint type parameter
# const LDAEinvType{invT,DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,parT,perT} = LDAE{DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,parT,perT} # type alias for dispatch on invariants type parameter
# const LDAEparType{parT,DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,perT} = LDAE{DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,parT,perT} # type alias for dispatch on parameters type parameter
# const LDAEperType{perT,DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,parT} = LDAE{DT,TT,AT,PT,FT,UT,GT,ΦT,ŪT,ḠT,psiT,ΩT,V̄T,F̄T,lagT,invT,parT,perT} # type alias for dispatch on periodicity type parameter

# Base.hash(dae::LDAE, h::UInt) = hash(dae.d, hash(dae.m,
#           hash(dae.ϑ, hash(dae.f,
#           hash(dae.u, hash(dae.g, hash(dae.ϕ,
#           hash(dae.ū, hash(dae.ḡ, hash(dae.ψ,
#           hash(dae.ω, hash(dae.v̄, hash(dae.f̄,
#           hash(dae.t₀, hash(dae.q₀, hash(dae.p₀, hash(dae.λ₀, hash(dae.μ₀,
#           hash(dae.invariants, hash(dae.parameters, hash(dae.periodicity, h)))))))))))))))))))))

# Base.:(==)(dae1::LDAE, dae2::LDAE) = (
#                                 dae1.d == dae2.d
#                              && dae1.m == dae2.m
#                              && dae1.ϑ == dae2.ϑ
#                              && dae1.f == dae2.f
#                              && dae1.u == dae2.u
#                              && dae1.g == dae2.g
#                              && dae1.ϕ == dae2.ϕ
#                              && dae1.ū == dae2.ū
#                              && dae1.ḡ == dae2.ḡ
#                              && dae1.ψ == dae2.ψ
#                              && dae1.ω == dae2.ω
#                              && dae1.v̄ == dae2.v̄
#                              && dae1.f̄ == dae2.f̄
#                              && dae1.t₀ == dae2.t₀
#                              && dae1.q₀ == dae2.q₀
#                              && dae1.p₀ == dae2.p₀
#                              && dae1.λ₀ == dae2.λ₀
#                              && dae1.μ₀ == dae2.μ₀
#                              && dae1.lagrangian == dae2.lagrangian
#                              && dae1.invariants  == dae2.invariants
#                              && dae1.parameters  == dae2.parameters
#                              && dae1.periodicity == dae2.periodicity)

# function Base.similar(equ::LDAE, t₀::Real, q₀::StateVector, p₀::StateVector, λ₀::StateVector, μ₀::StateVector=initial_multiplier(q₀, equ.μ₀); parameters=equ.parameters)
#     @assert all([length(q) == ndims(equ) for q in q₀])
#     @assert all([length(p) == ndims(equ) for p in p₀])
#     @assert all([length(λ) == ndims(equ) for λ in λ₀])
#     @assert all([length(μ) == ndims(equ) for μ in μ₀])
#     _LDAE(equ.ϑ, equ.f, equ.u, equ.g, equ.ϕ, equ.ū, equ.ḡ, equ.ψ, equ.lagrangian, equ.ω, t₀, q₀, p₀, λ₀, μ₀;
#          v̄=equ.v̄, f̄=equ.f̄, invariants=equ.invariants, parameters=parameters, periodicity=equ.periodicity)
# end

# Base.similar(equ::LDAE, q₀, p₀, λ₀=initial_multiplier(q₀, equ.λ₀), μ₀=initial_multiplier(q₀, equ.μ₀); kwargs...) = similar(equ, equ.t₀, q₀, p₀, λ₀, μ₀; kwargs...)
# Base.similar(equ::LDAE, t₀::Real, q₀::State, p₀::State, λ₀::State=initial_multiplier(q₀, equ.λ₀), μ₀=initial_multiplier(q₀, equ.μ₀); kwargs...) = similar(equ, t₀, [q₀], [p₀], [λ₀], [μ₀]; kwargs...)

# hassecondary(::LDAEpsiType{<:Nothing}) = false
# hassecondary(::LDAEpsiType{<:Function}) = true

# hasinvariants(::LDAEinvType{<:NullInvariants}) = false
# hasinvariants(::LDAEinvType{<:NamedTuple}) = true

# hasparameters(::LDAEparType{<:NullParameters}) = false
# hasparameters(::LDAEparType{<:NamedTuple}) = true

# hasperiodicity(::LDAEperType{<:Nothing}) = false
# hasperiodicity(::LDAEperType{<:AbstractArray}) = true

# @inline Base.axes(equation::LDAE) = axes(equation.q₀[begin])
# @inline Base.ndims(equation::LDAE) = equation.d
# @inline GeometricBase.nsamples(equ::LDAE) = length(eachindex(equ.q₀))
# @inline GeometricBase.nconstraints(equation::LDAE) = equation.m

# @inline GeometricBase.periodicity(equation::LDAE) = hasperiodicity(equation) ? equation.periodicity : zero(equation.q₀[begin])
# @inline initial_conditions(equation::LDAE) = (equation.t₀, equation.q₀, equation.p₀, equation.λ₀, equation.μ₀)

# _get_ϑ(equ::LDAE) = hasparameters(equ) ? (t,q,v,ϑ)     -> equ.ϑ(t, q, v, ϑ, equ.parameters) : equ.ϑ
# _get_f(equ::LDAE) = hasparameters(equ) ? (t,q,v,f)     -> equ.f(t, q, v, f, equ.parameters) : equ.f
# _get_u(equ::LDAE) = hasparameters(equ) ? (t,q,v,u)     -> equ.u(t, q, v, u, equ.parameters) : equ.u
# _get_g(equ::LDAE) = hasparameters(equ) ? (t,q,v,g)     -> equ.g(t, q, v, g, equ.parameters) : equ.g
# _get_ϕ(equ::LDAE) = hasparameters(equ) ? (t,q,v,ϕ)     -> equ.ϕ(t, q, v, ϕ, equ.parameters) : equ.ϕ
# _get_ū(equ::LDAE) = hasparameters(equ) ? (t,q,λ,ū)     -> equ.ū(t, q, λ, ū, equ.parameters) : equ.ū
# _get_ḡ(equ::LDAE) = hasparameters(equ) ? (t,q,λ,ḡ)     -> equ.ḡ(t, q, λ, ḡ, equ.parameters) : equ.ḡ
# _get_ψ(equ::LDAE) = hasparameters(equ) ? (t,q,p,v,f,ψ) -> equ.ψ(t, q, p, v, f, ψ, equ.parameters) : equ.ψ
# _get_ω(equ::LDAE) = hasparameters(equ) ? (t,q,v,ω)     -> equ.ω(t, q, v, ω, equ.parameters) : equ.ω
# _get_v̄(equ::LDAE) = hasparameters(equ) ? (t,q,v)       -> equ.v̄(t, q, v, equ.parameters) : equ.v̄
# _get_f̄(equ::LDAE) = hasparameters(equ) ? (t,q,v,f)     -> equ.f̄(t, q, v, f, equ.parameters) : equ.f̄
# _get_l(equ::LDAE) = hasparameters(equ) ? (t,q,v)       -> equ.lagrangian(t, q, v, equ.parameters) : equ.lagrangian

# function functions(equ::LDAE)
#     names = (:ϑ, :f, :u, :g, :ϕ)
#     equs  = (_get_ϑ(equ), _get_f(equ), _get_u(equ), _get_g(equ), _get_ϕ(equ))

#     if hassecondary(equ)
#         names = (names..., :ū, :ḡ, :ψ)
#         equs  = (equs..., _get_ū(equ), _get_ḡ(equ), _get_ψ(equ))
#     end

#     names = (names..., :l, :ω, :v̄, :f̄)
#     equs  = (equs..., _get_l(equ), _get_ω(equ), _get_v̄(equ), _get_f̄(equ))

#     NamedTuple{names}(equs)
# end


text:
# Autogenerated wrapper script for XGBoost_jll for i686-w64-mingw32-cxx11
export libxgboost, xgboost

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("XGBoost")
JLLWrappers.@declare_library_product(libxgboost, "xgboost.dll")
JLLWrappers.@declare_executable_product(xgboost)
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libxgboost,
        "bin\\xgboost.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        xgboost,
        "bin\\xgboost.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module CommunalHelperCassetteMoveBlock

using ..Ahorn, Maple
using Ahorn.CommunalHelper

@mapdef Entity "CommunalHelper/CassetteMoveBlock" CassetteMoveBlock(
    x::Integer,
    y::Integer,
    width::Integer=Maple.defaultBlockWidth,
    height::Integer=Maple.defaultBlockHeight,
    direction::String="Right",
    moveSpeed::Number=60.0,
    index::Integer=0,
    tempo::Number=1.0,
    customColor="",
)

const placements = Ahorn.PlacementDict(
    "Cassette Move Block ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteMoveBlock,
        "rectangle",
        Dict{String,Any}(
            "index" => index,
        ),
    ) for (color, index) in cassetteColorNames
)

Ahorn.editingOptions(entity::CassetteMoveBlock) = Dict{String,Any}(
    "index"     => cassetteColorNames,
    "direction" => Maple.move_block_directions,
    "moveSpeed" => Dict{String,Number}(
        "Slow" => 60.0,
        "Fast" => 75.0,
    ),
)
Ahorn.minimumSize(entity::CassetteMoveBlock) = 16, 16
Ahorn.resizable(entity::CassetteMoveBlock) = true, true

Ahorn.selection(entity::CassetteMoveBlock) = Ahorn.getEntityRectangle(entity)

const arrows = Dict{String,String}(
    "up"    => "objects/CommunalHelper/cassetteMoveBlock/arrow02",
    "left"  => "objects/CommunalHelper/cassetteMoveBlock/arrow04",
    "right" => "objects/CommunalHelper/cassetteMoveBlock/arrow00",
    "down"  => "objects/CommunalHelper/cassetteMoveBlock/arrow06",
)

const block = "objects/cassetteblock/solid"

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CassetteMoveBlock)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    index = Int(get(entity.data, "index", 0))
    color = getCassetteColor(index)

    hexColor = String(get(entity.data, "customColor", ""))
    if hexColor != "" && length(hexColor) == 6
        color = hexToRGBA(hexColor)
    end

    renderCassetteBlock(ctx, 0, 0, width, height, index, color)

    direction = lowercase(get(entity.data, "direction", "up"))
    arrowSprite = Ahorn.getSprite(arrows[lowercase(direction)], "Gameplay")

    Ahorn.drawImage(ctx, arrowSprite, div(width - arrowSprite.width, 2), div(height - arrowSprite.height, 2), tint=color)
end

end

max_stars_repo_name: vitellaryjr/CommunalHelper

max_issues_repo_name: vitellaryjr/CommunalHelper

max_forks_repo_name: vitellaryjr/CommunalHelper


text:
using ProgressMeter
using Printf
using MATLAB
using BSON: @save, @load

include("Aref_LD_grad_function.jl")

@everywhere begin
    include("Aref_LD_grad_function.jl")
    push!(LOAD_PATH, "/home/brandon_behring/Desktop/Leap_Frog_Code")
    using ProgressMeter
    using Printf
    using MATLAB
    using .Aref_LD_grad_function
    using BSON: @save, @load
    # H=range(.11,stop=0.145,length=107)
    # Energy=0.112012
    # Energy=.125789
    Energy=.129
    # Energy=.12501
    # Energy=.12501
    # Energy=.1249
    # Energy=.125
    t_typical=Aref_LD_grad_function.final_T(Energy)
    # file_name=location*h*".fig"
    #defining colors for PSS
    n_iter_Q=1000;#50
    # Q_start=.20
    Q_end=.05
    Q_start=0.0
    # Q_end=.7
    n_iter_P=1000
    # P_start=
    # P_end=.2
    P_end=3.5*1e-3
    P_start=-P_end

    t_end=100*t_typical
    N=n_iter_Q*n_iter_P
    t_end_mesh = t_end * ones(n_iter_Q,n_iter_P)
    Energy_mesh = Energy* ones(n_iter_Q,n_iter_P)
    ArrP=range(P_start,stop=P_end,length=n_iter_P)
    ArrQ=range(Q_start,stop=Q_end,length=n_iter_Q)
    mesh = [(P, Q) for Q in ArrQ, P in ArrP]
    mesh_list = reshape(mesh, 1, :)
end
@everywhere location="/media/brandon_behring/Extra_Space/MATLAB_FIGURES/"
@everywhere h=replace(@sprintf("%.15f",Energy),"."=>"_")
@everywhere h=replace(@sprintf("%.15f",Energy),"."=>"_")
@everywhere nQ = @sprintf("_%d",n_iter_Q)
@everywhere tend= @sprintf("_%d",t_end)
@everywhere tend_T= @sprintf("%4f",t_end)
@everywhere h_BSON="lagrangian_descriptors_data"*h*nQ*tend*".bson"
@everywhere file_name=location*"MATLAB_LD_"*h*nQ*tend*".fig"
@everywhere h_title = @sprintf("h= %.6f",Energy)*"with tend="*tend_T

LD= @showprogress pmap(Aref_LD_grad_function.gradM,mesh, Energy_mesh, t_end_mesh)
SAVE_DATA=Dict("LD"=>LD, "mesh"=>mesh, "ArrP"=>ArrP, "ArrQ"=>ArrQ, "Energy"=>Energy)
@save h_BSON SAVE_DATA
mat"figure();set(gcf, 'Position',  [0, 0, 1500, 1500]); hold on;"
mat"title($h_title)"
biggest=maximum(filter(!isnan,LD))
smallest=minimum(filter(!isnan,LD))

# mat"plot($P_0,$Q_0 ,'b.','MarkerSize',3)"
# mat"plot($P_1,$Q_1 ,'r.','MarkerSize',3)"
# mat"plot($P_2,$Q_2 ,'g.','MarkerSize',3)"
# mat"plot($P_4,$Q_4 ,'c.','MarkerSize',3)"
# mat"plot($P_5,$Q_5 ,'y.','MarkerSize',3)"

# mat"imagesc([0,$P_end],[0,$Q_end ],$LD, $clims)"
# mat"imagesc([0,-$P_end],[0,$Q_end ],$LD, $clims)"
# mat"imagesc([0,$P_end],[0,-$Q_end ],$LD, $clims)"
# mat"imagesc([0,-$P_end],[0,-$Q_end ],$LD, $clims)"
# mat"imagesc([$P_start,$P_end],[$Q_start,$Q_end ],$LD,$clims)"
# mat"imagesc([$P_start,$P_end],[$Q_start,$Q_end ],$LD)"

# mat"imagesc([0, $P_end ],[$Q_start, $Q_end],$LD)"
# mat"imagesc([0,-$P_end ],[$Q_start, $Q_end],$LD)"

clims=[smallest biggest*.005]
mat"imagesc([$P_start,$P_end ],[0,$Q_end],$LD, $clims)"
mat"imagesc([$P_end,$P_start ],[0,-$Q_end],$LD, $clims)"

mat"imagesc([$P_start,$P_end ],[0,$Q_end],$LD)"
mat"imagesc([$P_end,$P_start ],[0,-$Q_end],$LD)"

# mat"imagesc([0,$P_end],[0,$Q_end ],$LD)"
# mat"imagesc([0,-$P_end],[0,$Q_end ],$LD)"
# mat"imagesc([0,$P_end],[0,-$Q_end ],$LD)"
# mat"imagesc([0,-$P_end],[0,-$Q_end ],$LD)"

mat"colorbar"
mat"axis([-$P_end,$P_end,-$Q_end,$Q_end ])"
# mat"axis([ $P_start,$P_end,$Q_start,$Q_end ])"
mat"savefig($file_name)"


text:
using MedicalTutorials
using Test

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


text:
"""
Functionality for dealing with HTTP cookies.
"""
module Cookies

using HTTP, HttpCommon, Nullables
using Genie, Genie.Encryption

"""
    get(res::Response, key::Union{String,Symbol}) :: Nullable{String}

Retrieves a value stored on the cookie as `key` from the `Respose` object.
"""
function get(res::HTTP.Response, key::Union{String,Symbol}) :: Nullable{String}
  if haskey(Dict(res.headers), "Set-Cookie")
    cookies = to_dict(res)
    if haskey(cookies, string(key))
      return Nullable{String}(cookies[string(key)] |> Genie.Encryption.decrypt)
    end
  end

  Nullable{String}()
end

"""
    get(req::Request, key::Union{String,Symbol}) :: Nullable{String}

Retrieves a value stored on the cookie as `key` from the `Request` object.
"""
function get(req::HTTP.Request, key::Union{String,Symbol}) :: Nullable{String}
  if haskey(Dict(req.headers), "Cookie")
    cookies = to_dict(req)
    if haskey(cookies, string(key))
      return Nullable{String}(cookies[string(key)] |> Genie.Encryption.decrypt)
    end
  end

  Nullable{String}()
end

"""
    set!(res::Response, key::Union{String,Symbol}, value::Any, attributes::Dict) :: Dict{String,HttpCommon.Cookie}

Sets `value` under the `key` label on the `Cookie`.
"""
function set!(res::HTTP.Response, key::Union{String,Symbol}, value::Any, attributes::Dict = Dict()) :: Vector{HTTP.Cookies.Cookie}
  normalized_attrs = Dict{Symbol,Any}()
  for (k,v) in attributes
    normalized_attrs[Symbol(lowercase(string(k)))] = v
  end
  HTTP.Cookie(string(key), string(value) |> Genie.Encryption.encrypt; normalized_attrs...)

  HTTP.Cookies.readsetcookies("", "")
end

"""
    to_dict(req::Request) :: Dict{String,String}

Extracts the `Cookie` and `Set-Cookie` data from the `Request` and `Response` objects and converts it into a dict.
"""
function to_dict(r::Union{HTTP.Request,HTTP.Response}) :: Dict{String,String}
  d = Dict{String,String}()
  headers = Dict(r.headers)

  if haskey(headers, "Cookie")
    for cookie in split(headers["Cookie"], ";")
      cookie_parts = split(cookie, "=")
      d[strip(cookie_parts[1])] = cookie_parts[2]
    end
  end

  if haskey(headers, "Set-Cookie")
    for cookie in split(headers["Set-Cookie"], ";")
      cookie_parts = split(cookie, "=")
      d[strip(cookie_parts[1])] = cookie_parts[2]
    end
  end

  d
end

end


text:
@sumProductRule(:node_type     => Nonlinear{Sampling},
                :outbound_type => Message{SampleList},
                :inbound_types => (Nothing, Message),
                :name          => SPNonlinearSOutNM)

@sumProductRule(:node_type     => Nonlinear{Sampling},
                :outbound_type => Message{Function},
                :inbound_types => (Message, Nothing),
                :name          => SPNonlinearSIn1MN)

mutable struct SPNonlinearSOutNGX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSOutNGX}) = Message{SampleList}
function isApplicable(::Type{SPNonlinearSOutNGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    for input_type in input_types[2:end]
        (input_type << Message{Gaussian}) || return false
    end

    return true
end

mutable struct SPNonlinearSInGX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSInGX}) = Message{GaussianWeightedMeanPrecision}
function isApplicable(::Type{SPNonlinearSInGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) && return false # Require any message on out

    nothing_inputs = 0
    gaussian_inputs = 0
    for input_type in input_types[2:end]
        if input_type == Nothing
            nothing_inputs += 1
        elseif input_type << Message{Gaussian}
            gaussian_inputs += 1
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs == total_inputs - 2)
end

mutable struct SPNonlinearSOutNMX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSOutNMX}) = Message{SampleList}
function isApplicable(::Type{SPNonlinearSOutNMX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false

    gaussian_inputs = 0
    for input_type in input_types[2:end]
        (input_type << Message) || return false
        if input_type << Message{Gaussian}
            gaussian_inputs += 1
        end
    end
    (gaussian_inputs == total_inputs - 1) && return false # Rule does not apply if all inputs are Gaussian

    return true
end

mutable struct SPNonlinearSInMX <: SumProductRule{Nonlinear{Sampling}} end
outboundType(::Type{SPNonlinearSInMX}) = Message{Function}
function isApplicable(::Type{SPNonlinearSInMX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) && return false

    nothing_inputs = 0
    gaussian_inputs = 0
    for input_type in input_types[2:end]
        if input_type == Nothing
            nothing_inputs += 1
        elseif input_type << Message{Gaussian}
            gaussian_inputs += 1
        end
    end

    return (nothing_inputs == 1) && (gaussian_inputs != total_inputs - 2) # Rule does not apply if all inbounds are Gaussian
end

mutable struct MNonlinearSInMGX <: MarginalRule{Nonlinear{Sampling}} end
function isApplicable(::Type{MNonlinearSInMGX}, input_types::Vector{<:Type})
    total_inputs = length(input_types)
    (total_inputs > 2) || return false
    (input_types[1] == Nothing) || return false # Indicates marginalization over outbound variable

    for input_type in input_types[2:end]
        (input_type << Message{Gaussian}) || return false
    end

    return true
end