text:
export MultirateInfinitesimalStep,
    TimeScaledRHS,
    MISRK1,
    MIS2,
    MISRK2a,
    MISRK2b,
    MIS3C,
    MISRK3,
    MIS4,
    MIS4a,
    MISKWRK43,
    TVDMISA,
    TVDMISB,
    getnsubsteps

"""
    TimeScaledRHS(a, b, rhs!)

When evaluate at time `t`, evaluates `rhs!` at time `a + bt`.
"""
mutable struct TimeScaledRHS{N, RT}
    a::RT
    b::RT
    rhs!::Any
    function TimeScaledRHS(a, b, rhs!)
        RT = typeof(a)
        if isa(rhs!, Tuple)
            N = length(rhs!)
        else
            N = 1
        end
        new{N, RT}(a, b, rhs!)
    end
end

function (o::TimeScaledRHS{1, RT} where {RT})(dQ, Q, params, tau; increment)
    o.rhs!(dQ, Q, params, o.a + o.b * tau; increment = increment)
end

function (o::TimeScaledRHS{2, RT} where {RT})(dQ, Q, params, tau, i; increment)
    o.rhs![i](dQ, Q, params, o.a + o.b * tau; increment = increment)
end

mutable struct OffsetRHS{AT}
    offset::AT
    rhs!::Any
    function OffsetRHS(offset, rhs!)
        AT = typeof(offset)
        new{AT}(offset, rhs!)
    end
end

function (o::OffsetRHS{AT} where {AT})(dQ, Q, params, tau; increment)
    o.rhs!(dQ, Q, params, tau; increment = increment)
    dQ .+= o.offset
end

"""
MultirateInfinitesimalStep(slowrhs!, fastrhs!, fastmethod,
                           α, β, γ,
                           Q::AT; dt=0, t0=0) where {AT<:AbstractArray}

This is a time stepping object for explicitly time stepping the partitioned differential
equation given by right-hand-side functions `f_fast` and `f_slow` with the state `Q`, i.e.,

```math
  \\dot{Q} = f_fast(Q, t) + f_slow(Q, t)
```

with the required time step size `dt` and optional initial time `t0`.  This
time stepping object is intended to be passed to the `solve!` command.

The constructor builds a multirate infinitesimal step Runge-Kutta scheme
based on the provided `α`, `β` and `γ` tableaux and `fastmethod` for solving
the fast modes.

### References
 - [KnothWensch2014](@cite)
"""
mutable struct MultirateInfinitesimalStep{
    T,
    RT,
    AT,
    FS,
    Nstages,
    Nstagesm1,
    Nstagesm2,
    Nstages_sq,
} <: AbstractODESolver
    "time step"
    dt::RT
    "time"
    t::RT
    "elapsed time steps"
    steps::Int
    "storage for y_n"
    yn::AT
    "Storage for ``Y_nj - y_n``"
    ΔYnj::NTuple{Nstagesm2, AT}
    "Storage for ``f(Y_nj)``"
    fYnj::NTuple{Nstagesm1, AT}
    "Storage for offset"
    offset::AT
    "slow rhs function"
    slowrhs!::Any
    "RHS for fast solver"
    tsfastrhs!::TimeScaledRHS{N, RT} where {N}
    "fast rhs method"
    fastsolver::FS
    "number of substeps per stage"
    nsubsteps::Int
    α::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
    β::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
    γ::SArray{NTuple{2, Nstages}, RT, 2, Nstages_sq}
    d::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
    c::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
    c̃::SArray{NTuple{1, Nstages}, RT, 1, Nstages}
    function MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q::AT;
        dt = 0,
        t0 = 0,
    ) where {AT <: AbstractArray}

        T = eltype(Q)
        RT = real(T)

        Nstages = size(α, 1)

        yn = similar(Q)
        ΔYnj = ntuple(_ -> similar(Q), Nstages - 2)
        fYnj = ntuple(_ -> similar(Q), Nstages - 1)
        offset = similar(Q)
        tsfastrhs! = TimeScaledRHS(RT(0), RT(0), fastrhs!)
        fastsolver = fastmethod(tsfastrhs!, Q)

        d = sum(β, dims = 2)

        c = similar(d)
        for i in eachindex(c)
            c[i] = d[i]
            if i > 1
                c[i] += sum(j -> (α[i, j] + γ[i, j]) * c[j], 1:(i - 1))
            end

            # When d[i] = 0, we do not perform fast substepping, therefore
            # we do not need to scale the β, γ coefficients
            if !(abs(d[i]) < 1.e-10)
                β[i, :] ./= d[i]
                γ[i, :] ./= d[i]
            end
        end
        c̃ = α * c

        new{
            T,
            RT,
            AT,
            typeof(fastsolver),
            Nstages,
            Nstages - 1,
            Nstages - 2,
            Nstages^2,
        }(
            RT(dt),
            RT(t0),
            0,
            yn,
            ΔYnj,
            fYnj,
            offset,
            slowrhs!,
            tsfastrhs!,
            fastsolver,
            nsubsteps,
            α,
            β,
            γ,
            d,
            c,
            c̃,
        )
    end
end

function MultirateInfinitesimalStep(
    mis,
    op::TimeScaledRHS{2, RT} where {RT},
    fastmethod,
    Q = nothing;
    dt = 0,
    t0 = 0,
    nsubsteps = 1,
) where {AT <: AbstractArray}

    return mis(
        op.rhs![1],
        op.rhs![2],
        fastmethod,
        nsubsteps,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function dostep!(
    Q,
    mis::MultirateInfinitesimalStep,
    p,
    time::Real,
    nsubsteps::Int,
    iStage::Int,
    slow_δ = nothing,
    slow_rv_dQ = nothing,
    slow_scaling = nothing,
)
    if isa(mis.slowrhs!, OffsetRHS{AT} where {AT})
        mis.slowrhs!.offset = slow_rv_dQ
    else
        mis.slowrhs! = OffsetRHS(slow_rv_dQ, mis.slowrhs!)
    end
    for i in 1:nsubsteps
        dostep!(Q, mis, p, time)
        time += mis.fastsolver.dt
    end
end

function dostep!(Q, mis::MultirateInfinitesimalStep, p, time)
    dt = mis.dt
    FT = eltype(dt)
    α = mis.α
    β = mis.β
    γ = mis.γ
    yn = mis.yn
    ΔYnj = mis.ΔYnj
    fYnj = mis.fYnj
    offset = mis.offset
    d = mis.d
    c = mis.c
    c̃ = mis.c̃
    slowrhs! = mis.slowrhs!
    fastsolver = mis.fastsolver
    fastrhs! = mis.tsfastrhs!
    nsubsteps = mis.nsubsteps

    nstages = size(α, 1)

    copyto!(yn, Q) # first stage
    for i in 2:nstages
        slowrhs!(fYnj[i - 1], Q, p, time + c[i - 1] * dt, increment = false)

        groupsize = 256
        event = Event(array_device(Q))
        event = update!(array_device(Q), groupsize)(
            realview(Q),
            realview(offset),
            Val(i),
            realview(yn),
            map(realview, ΔYnj[1:(i - 2)]),
            map(realview, fYnj[1:(i - 1)]),
            α[i, :],
            β[i, :],
            γ[i, :],
            dt;
            ndrange = length(realview(Q)),
            dependencies = (event,),
        )
        wait(array_device(Q), event)

        # When d[i] = 0, we do not perform fast substepping;
        # instead we just update the slow tendency
        if abs(d[i]) < 1.e-10
            Q .+= dt .* offset
        else
            fastrhs!.a = time + c̃[i] * dt
            fastrhs!.b = (c[i] - c̃[i]) / d[i]

            τ = zero(FT)
            nsubstepsLoc = ceil(Int, nsubsteps * d[i])
            dτ = d[i] * dt / nsubstepsLoc
            updatetime!(fastsolver, τ)
            updatedt!(fastsolver, dτ)
            # TODO: we want to be able to write
            #   solve!(Q, fastsolver, p; numberofsteps = mis.nsubsteps)  #(1c)
            # especially if we want to use StormerVerlet, but need some way to pass in `offset`
            dostep!(
                Q,
                fastsolver,
                p,
                τ,
                nsubstepsLoc,
                i,
                FT(1),
                realview(offset),
                nothing,
            )  #(1c)
        end
    end
end

@kernel function update!(
    Q,
    offset,
    ::Val{i},
    yn,
    ΔYnj,
    fYnj,
    αi,
    βi,
    γi,
    dt,
) where {i}
    e = @index(Global, Linear)
    @inbounds begin
        if i > 2
            ΔYnj[i - 2][e] = Q[e] - yn[e] # is 0 for i == 2
        end
        Q[e] = yn[e] # (1a)
        offset[e] = (βi[1]) .* fYnj[1][e] # (1b)
        @unroll for j in 2:(i - 1)
            Q[e] += αi[j] .* ΔYnj[j - 1][e] # (1a cont.)
            offset[e] += (γi[j] / dt) * ΔYnj[j - 1][e] + βi[j] * fYnj[j][e] # (1b cont.)
        end
    end
end

function MISRK1(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)
    α = zeros(2, 2)
    β = beta(MISRK1, RT)
    γ = zeros(2, 2)
    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MISRK1), RT::DataType)
    β = [
        0 0
        1 0
    ]
end

function MIS2(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0
        0 0 0 0
        0 RT(0.536946566710) 0 0
        0 RT(0.480892968551) RT(0.500561163566) 0
    ]

    β = beta(MIS2, RT)

    γ = [
        0 0 0 0
        0 0 0 0
        0 RT(0.652465126004) 0 0
        0 RT(-0.0732769849457) RT(0.144902430420) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MIS2), RT::DataType)
    β = [
        0 0 0 0
        RT(0.126848494553) 0 0 0
        RT(-0.784838278826) RT(1.37442675268) 0 0
        RT(-0.0456727081749) RT(-0.00875082271190) RT(0.524775788629) 0
    ]
end

function MISRK2a(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0
        0 0 0
        0 1 0
    ]

    β = beta(MISRK2a, RT)

    γ = zeros(3, 3)

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MISRK2a), RT::DataType)
    β = [
        0 0 0
        RT(0.5) 0 0
        RT(-0.5) 1 0
    ]
end

function MISRK2b(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0
        0 0 0
        0 1 0
    ]

    β = beta(MISRK2b, RT)

    γ = zeros(3, 3)

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MISRK2b), RT::DataType)
    β = [
        0 0 0
        1 0 0
        RT(-0.5) RT(0.5) 0
    ]
end

function MIS3C(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0
        0 0 0 0
        0 RT(0.589557277145) 0 0
        0 RT(0.544036601551) RT(0.565511042564) 0
    ]

    β = beta(MIS3C, RT)

    γ = [
        0 0 0 0
        0 0 0 0
        0 RT(0.142798786398) 0 0
        0 RT(-0.0428918957402) RT(0.0202720980282) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MIS3C), RT::DataType)
    β = [
        0 0 0 0
        RT(0.397525189225) 0 0 0
        RT(-0.227036463644) RT(0.624528794618) 0 0
        RT(-0.00295238076840) RT(-0.270971764284) RT(0.671323159437) 0
    ]
end

function MISRK3(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)
    α = zeros(4, 4)
    β = beta(MISRK3, RT)
    γ = zeros(4, 4)
    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MISRK3), RT::DataType)
    β = [
        0 0 0 0
        RT(0.3333333333333333) 0 0 0
        0 RT(0.5) 0 0
        0 0 1 0
    ]
end

function MIS4(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0 0
        0 0 0 0 0
        0 RT(0.914092810304) 0 0 0
        0 RT(1.14274417397) RT(-0.295211246188) 0 0
        0 RT(0.112965282231) RT(0.337369411296) RT(0.503747183119) 0
    ]

    β = beta(MIS4, RT)

    γ = [
        0 0 0 0 0
        0 0 0 0 0
        0 RT(0.678951983291) 0 0 0
        0 RT(-1.38974164070) RT(0.503864576302) 0 0
        0 RT(-0.375328608282) RT(0.320925021109) RT(-0.158259688945) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MIS4), RT::DataType)
    β = [
        0 0 0 0 0
        RT(0.136296478423) 0 0 0 0
        RT(0.280462398979) RT(-0.0160351333596) 0 0 0
        RT(0.904713355208) RT(-1.04011183154) RT(0.652337563489) 0 0
        RT(0.0671969845546) RT(-0.365621862610) RT(-0.154861470835) RT(0.970362444469) 0
    ]
end

function MIS4a(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0 0
        0 0 0 0 0
        0 RT(0.52349249922385610) 0 0 0
        0 RT(1.1683374366893629) RT(-0.75762080241712637) 0 0
        0 RT(-0.036477233846797109) RT(0.56936148730740477) RT(0.47746263002599681) 0
    ]

    # β[5,1] in the paper is incorrect
    # the correct value is used below (from authors)
    β = beta(MIS4a, RT)

    γ = [
        0 0 0 0 0
        0 0 0 0 0
        0 RT(0.13145089796226542) 0 0 0
        0 RT(-0.36855857648747881) RT(0.33159232636600550) 0 0
        0 RT(-0.065767130537473045) RT(0.040591093109036858) RT(0.064902111640806712) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MIS4a), RT::DataType)
    β = [
        0 0 0 0 0
        RT(0.38758444641450318) 0 0 0 0
        RT(-0.025318448354142823) RT(0.38668943087310403) 0 0 0
        RT(0.20899983523553325) RT(-0.45856648476371231) RT(0.43423187573425748) 0 0
        RT(-0.10048822195663100) RT(-0.46186171956333327) RT(0.83045062122462809) RT(0.27014914900250392) 0
    ]
end

function MISKWRK43(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0 0
        0 0 0 0 0
        0 1 0 0 0
        0 0 1 0 0
        0 0 0 1 0
    ]

    β = beta(MISKWRK43, RT)

    γ = zeros(5, 5)

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(MISKWRK43), RT::DataType)
    β = [
        0 0 0 0 0
        0.5 0 0 0 0
        -RT(2 // 3) RT(2 // 3) 0 0 0
        RT(0.5) -1 1 0 0
        -RT(1 // 6) RT(2 // 3) -RT(2 // 3) RT(1 // 6) 0
    ]
end

function TVDMISA(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0
        0 0 0 0
        0 RT(0.1946360605647457) 0 0
        0 RT(0.3971200136786614) RT(0.2609434606211801) 0
    ]

    β = beta(TVDMISA, RT)

    γ = [
        0 0 0 0
        0 0 0 0
        0 RT(0.5624048933209129) 0 0
        0 RT(0.4408467475713277) RT(-0.2459300561692391) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(TVDMISA), RT::DataType)
    β = [
        0 0 0 0
        RT(2 // 3) 0 0 0
        RT(-0.28247174703488398) RT(4 // 9) 0 0
        RT(-0.31198081960042401) RT(0.18082737579913699) RT(9 // 16) 0
    ]
end

function TVDMISB(
    slowrhs!,
    fastrhs!,
    fastmethod,
    nsubsteps,
    Q::AT;
    dt = 0,
    t0 = 0,
) where {AT <: AbstractArray}
    FT = eltype(Q)
    RT = real(FT)

    α = [
        0 0 0 0
        0 0 0 0
        0 RT(0.42668232863311001) 0 0
        0 RT(0.26570779016173801) RT(0.41489966891866698) 0
    ]

    β = beta(TVDMISB, RT)

    γ = [
        0 0 0 0
        0 0 0 0
        0 RT(0.28904389120139701) 0 0
        0 RT(0.45113560071334202) RT(-0.25006656847591002) 0
    ]

    MultirateInfinitesimalStep(
        slowrhs!,
        fastrhs!,
        fastmethod,
        nsubsteps,
        α,
        β,
        γ,
        Q;
        dt = dt,
        t0 = t0,
    )
end

function beta(::typeof(TVDMISB), RT::DataType)
    β = [
        0 0 0 0
        RT(2 // 3) 0 0 0
        RT(-0.25492859100078202) RT(4 // 9) 0 0
        RT(-0.26452517179288798) RT(0.11424084424766399) RT(9 // 16) 0
    ]
end

function getnsubsteps(mis, ns::Int, RT::DataType)
    d = sum(beta(mis, RT), dims = 2)
    return d ./ ns
end

max_forks_repo_name: jm-c/CLIMA


text:
function test_hailstone_nexthailstone()
    @test_throws DomainError nexthailstone(0)
    @test nexthailstone(1) == 4
    @test nexthailstone(2) == 1
end

function test_hailstone_hailstonelength()
    @test_throws DomainError hailstonelength(0)
    @test hailstonelength(1) == 1
    @test hailstonelength(13) == 10
end

function test_hailstone_allhailstone()
    @test_throws DomainError allhailstone(0)
    @test collect(allhailstone(1)) == [1]
    @test collect(allhailstone(13)) == [13, 40, 20, 10, 5, 16, 8, 4, 2, 1]
    @test eltype(allhailstone(1)) == typeof(1)
    @test eltype(allhailstone(13)) == typeof(13)
end

function test_hailstone_all()
    print(rpad("Math.NumberTheory.Hailstone...", 50, ' '))

    test_hailstone_nexthailstone()
    test_hailstone_hailstonelength()
    test_hailstone_allhailstone()

    println("PASS")
end

ext: jl


text:
using IndirectLikelihood

using IndirectLikelihood:
    # imported for testing
    vec_parameters

import IndirectLikelihood:
    # problem API
    simulate_data, common_random, common_random!, MLE, loglikelihood

using Base.Test

using ContinuousTransformations
using Distributions: logpdf, Normal, MvNormal
using Optim
using StatsBase: Weights
using Suppressor

const RNG = Base.Random.GLOBAL_RNG

"Extract first line."
firstline(lines) = split(lines, '\n')[1]

include("utilities.jl")
include("mvnormal.jl")
include("ols.jl")
include("general.jl")
include("../docs/make.jl")


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    abstract type BAT.AbstractModeEstimator

Abstract type for BAT optimization algorithms.

A typical application for optimization in BAT is mode estimation
(see [`bat_findmode`](@ref)),
"""
abstract type AbstractModeEstimator end

"""
    bat_findmode(
        [rng::AbstractRNG],
        target::BAT.AnySampleable,
        [algorithm::BAT.AbstractModeEstimator]
    )::DensitySampleVector

Estimate the global mode of `target`.

Returns a NamedTuple of the shape

```julia
(result = X::DensitySampleVector, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_findmode`, add methods to
    `bat_findmode_impl` instead.
"""
function bat_findmode end
export bat_findmode

function bat_findmode_impl end

function bat_findmode(
    rng::AbstractRNG,
    target::AnySampleable,
    algorithm = bat_default_withdebug(bat_findmode, Val(:algorithm), target);
    kwargs...
)
    r = bat_findmode_impl(rng, target, algorithm; kwargs...)
    result_with_args(r, (algorithm = algorithm,))
end

function bat_findmode(
    target::AnySampleable,
    algorithm = bat_default_withdebug(bat_findmode, Val(:algorithm), target);
    kwargs...
)
    rng = bat_default_withinfo(bat_findmode, Val(:rng), target)
    r = bat_findmode_impl(rng, target, algorithm; kwargs...)
    result_with_args(r, (algorithm = algorithm,))
end

function argchoice_msg(::typeof(bat_findmode), ::Val{:rng}, x::AbstractRNG)
    "Initializing new RNG of type $(typeof(x))"
end

function argchoice_msg(::typeof(bat_findmode), ::Val{:algorithm}, x::AbstractModeEstimator)
    "Using mode estimator $x"
end

"""
    bat_marginalmode(
        samples::DensitySampleVector;
        nbins::Union{Integer, Symbol} = 200
    )::DensitySampleVector

*Experimental feature, not part of stable public API.*

Estimates a marginal mode of `samples` by finding the maximum of marginalized posterior for each dimension.

Returns a NamedTuple of the shape

```julia
(result = X::DensitySampleVector, ...)
```

`nbins` specifies the number of bins that are used for marginalization. The default value is `nbins=200`. The optimal number of bins can be estimated using  the following keywords:

* `:sqrt`  — Square-root choice

* `:sturges` — Sturges' formula

* `:rice` — Rice Rule

* `:scott` — Scott's normal reference rule

* `:fd` —  Freedman–Diaconis rule

Returns a NamedTuple of the shape

```julia
(result = X::DensitySampleVector, ...)
```

Result properties not listed here are algorithm-specific and are not part
of the stable public API.

!!! note

    Do not add add methods to `bat_marginalmode`, add methods to
    `bat_marginalmode_impl` instead.
"""
function bat_marginalmode end
export bat_marginalmode

function bat_marginalmode_impl end

function bat_marginalmode(samples::DensitySampleVector; kwargs...)
    r = bat_marginalmode_impl(samples::DensitySampleVector; kwargs...)
    result_with_args(r, NamedTuple(), kwargs)
end


text:
const keywordlist = [
    "ELEMENT",
    "FUNCTION",
    "FUNCT",
    "TYPE_DEFINITION",
    "PHASE",
    "CONSTITUENT",
    "CONST",
    "PARAMETER",
    "PAR",
]

function read_tdb(io::IO)
    elems = Vector{Element}()
    funcs = Vector{GFunction}()
    phass = Vector{Phase}()
    types = Vector{TypeDefinition}()

    # Parse each elements
    while !eof(io)
        block = readline(io, keep = true)

        block = replace(block, '\t' => "") # remove tabs

        isempty(block) && continue # blank line
        !occursin(r"\S+", block) && continue # blank line
        occursin(r"\$", block) && continue # comment line

        i = findfirst('!', block)
        if isnothing(i)
            block *= readuntil(io, '!')
        else
            block = block[1:i-1] # remove '!'
        end

        block = replace(block, '\r' => "")
        block = replace(block, '\n' => ' ')
        strs = split(block)
        isempty(strs) && continue

        keyword = uppercase(strs[1])

        !(keyword in keywordlist) && continue
        try
            if keyword == "ELEMENT"
                elem = parse_element(block)
                push!(elems, elem)
            elseif keyword in ["FUNCTION", "FUNCT"]
                func = parse_gfunction(block)
                push!(funcs, func)
            elseif keyword == "TYPE_DEFINITION"
                type = parse_typedefinition(block)
                push!(types, type)
            elseif keyword == "PHASE"
                phas = parse_phase(block)
                push!(phass, phas)
            elseif keyword in ["CONSTITUENT", "CONST"]
                phasname, cons = parse_constituent(block)
                phas = filter(phas -> uppercase(phas.name) == uppercase(phasname), phass)
                @assert length(phas) == 1
                phas[1].cons = cons
            elseif keyword in ["PARAMETER", "PAR"]
                phasname, param = parse_parameter(block)
                isnothing(param) && continue
                phas = filter(phas -> uppercase(phas.name) == uppercase(phasname), phass)
                @assert length(phas) == 1
                push!(phas[1].params, param)
            end
        catch
            error("Unrecognized description in the tdb file:\n$block")
        end
    end

    return Database(elems, funcs, phass, types)
end

function read_tdb(tdb::AbstractString)
    open(tdb, "r") do io
        read_tdb(io)
    end
end

function parse_element(text::AbstractString)
    strs = split(text)
    @assert length(strs) == 6
    @assert match(r"ELEMENT"i, strs[1]) ≠ ""
    name = format_element(strs[2])
    refstate = strs[3]
    mass, H298, S298 = parse.(Float64, strs[4:6])
    return Element(name, refstate, mass, H298, S298)
end

function parse_gfunction(text::AbstractString)
    strs = split(text)
    @assert match(r"FUNCTION"i, strs[1]) ≠ ""
    name = strs[2]
    temp, funcstr = parse_function(name, strs[3:end])
    return GFunction(name, temp, funcstr)
end

function parse_function(name::AbstractString, strs::Vector{SubString{String}})
    Ts = Vector{AbstractString}()
    funcstrs = Vector{AbstractString}()

    push!(Ts, strs[1]) # ex. 298.15
    k = 2

    while true
        push!(funcstrs, replace(strs[k], "#" => ""))
        k += 1
        while isnothing(findfirst(';', funcstrs[end]))
            funcstrs[end] = funcstrs[end] * strs[k]
            k += 1
        end
        push!(Ts, strs[k])
        k += 1
        if strs[k] == "N"
            k += 1
            break
        elseif strs[k] == "Y"
            k += 1
            continue
        else
            error("unexpected format in the tdb file")
        end
    end

    N = length(funcstrs)
    str = "function $(name)(T)::typeof(T)\n"
    for i in 1:N
        funcstrs[i] = strip(funcstrs[i], ';')

        str *= i == 1 ? "\tif " : "\telseif "
        str *= "$(Ts[i]) ≤ T ≤ $(Ts[i+1])\n"
        str *= "\t\t" * funcstrs[i] * "\n"
    end
    str *= "\tend\nend\n"

    Tl, Tu = tryparse.(Int, [Ts[1], Ts[end]])
    if nothing in [Tl, Tu]
        Tl, Tu = parse.(Float64, [Ts[1], Ts[end]])
    end

    return (Tl, Tu), str
end

function parse_typedefinition(text::AbstractString)
    strs = split(text)
    @assert match(r"TYPE_DEFINITION"i, strs[1]) ≠ ""
    code = strs[2]
    @assert length(code) == 1
    code = code[1]
    key = strs[3]
    if !(key in ["SEQ", "GES"])
        @warn "Unsupported keyword, $key, for TYPE_DEFINITION"
    end
    args = join(strs[4:end], ' ')
    return TypeDefinition(code, key, args)
end

function parse_phase(text::AbstractString)
    strs = split(text)
    @assert match(r"PHASE"i, strs[1]) ≠ ""
    model = strs[3] # "%", "%&", "%Q+", etc.
    K = parse(Int, strs[4])
    i = 5 + K - 1
    sites = tryparse.(Int, strs[5:i])
    if nothing in sites
        sites = parse.(Float64, strs[5:i])
    end
    text = replace(strs[2], " " => "")
    strs = split(text, ':')
    @assert length(strs) ≤ 2
    name = strs[1]
    state = length(strs) == 1 ? 'S' : strs[2][1]
    return Phase(name, state[1], model, sites)
end

function parse_constituent(text::AbstractString)
    strs = split(text)
    @assert match(r"CONSTITUENT"i, strs[1]) ≠ ""
    phas = split(strs[2], ':')[1] # "Liquid", "Bcc", "Fcc", etc.
    cons_text = join(strs[3:end]) # ":Cu,Zn", ":Cu,Zn:Cu,Zn:", etc.
    cons_text = format_constitution(cons_text)
    cons_strs = split(cons_text, ':')
    cons_strs = filter(e -> e ≠ "", cons_strs)
    cons = map(text -> split(text, ','), cons_strs)
    return phas, cons
end

function parse_parameter(text::AbstractString)
    # Parameters
    strs = split(text)
    @assert match(r"PAR"i, strs[1]) ≠ ""

    text = join(strs[2:end], ' ') # ex. G(BCC2,...

    k = findfirst('(', text)
    l = findfirst(')', text)
    name = replace(text[1:l], " " => "") # G(BCC_B2,Cu:Cu,Zn;0)
    name = format_constitution(name)
    funcname = getfuncname(name) # GBCC_B2CuCuZn0
    symbol = split(name, '(')[1] # G

    # Skip MatCalc specific parameters
    symbol == "HMVA" && return nothing, nothing
    symbol == "SE" && return nothing, nothing

    symbol = symbol[1]

    strs = split(text[l+1:end]) # 298.15...

    text = text[k+1:l-1] # BCC_B2,Cu,Zn:Cu,Zn;0
    text = format_constitution(text)

    k = findfirst(',', text)
    phas = replace(text[1:k-1], " " => "") # BCC_B2

    text = text[k+1:end] # Cu:Cu,Zn;0
    name = symbol * '(' * text * ')'
    comb_text, order_text = split(text, ';') # Cu:Cu,Zn, 0
    comb = parse_combination(comb_text)
    order = parse(Int, order_text)

    temp, funcstr = parse_function(funcname, strs)

    return phas, Parameter(name, symbol, comb, order, temp, funcname, funcstr)
end

function parse_combination(str::AbstractString)
    comb_strs = split(str, ':') # ["Cu", "Cu,Zn"]
    comb = map(text -> split(text, ','), comb_strs) # [[Cu], [Cu,Zn]]
end

function getfuncname(str::AbstractString)
    filter(c -> !in(c, ['(', ')', ',', ':', ';']), str)
end

function localname(param::Parameter)
    str = ""
    S = length(param.comb)
    for s in 1:S
    latt = param.comb[s]
        for cons in latt
            str *= subscript(cons)
            if cons ≠ latt[end]
                str *= ','
            end
        end
        if s ≠ S
            str *= ':'
        end
    end
    str *= ';' * string(param.order)
end

function localname(func::GFunction)
    return func.name
end

function format_constitution(str::AbstractString)
    # ex. str = "G(FCC,Cu:Cu,Zn,P;0)"
    replace(str, r"(?<=,|;|:).+?(?=,|;|:)" => format_specie)
end

function format_specie(str::AbstractString)
    # ex. str = "Cu2S1", "CU2O"
    str = replace(str, "%" => "")
    str = replace(str, r"(?<=\D)1(?=\D|$)" => "")
    replace(str, r"\D{1,2}(?=\d{1,}|$)" => format_element)
end

function format_element(str::AbstractString)
    N = length(str)
    N > 2 && error("Invalid element name: $str")
    N == 2 && return uppercase(str[1]) * lowercase(str[2])
    N == 1 && return uppercase(str)
end


text:
"""
    sfgintensity(ω, β, χ2_eff; <keyword arguments>)

Compute the intentensity of the sum frequency
in the reflected direction according to

``I_\\mathrm{SF} = \\frac{8\\pi^3\\omega^2\\sec^2\\beta}{c^3 n n_1 n_2}
\\left| \\chi_\\mathrm{eff}^{(2)} \\right|^2 I_1 I_2``.

# Arguments
* `ω`:  sum frequency
* `β`:  reflection angle of the sum frequency field
* `χ2_eff`: effective susceptibility

# Keyword Arguments
* `n`:  refractive index of the bulk medium for `ω`
* `n1`: refractive index of the bulk medium for `ω1`
* `n2`: refractive index of the bulk medium for `ω2`
* `I1`: intensity of `ω1`
* `I2`: intensity of `ω2`

https://doi.org/10.1103/PhysRevB.59.12632
"""
function sfgintensity(ω, β, χ2_eff; n=1.0, n1=1.0, n2=1.0, I1=1.0, I2=1.0)
    c = 299792458.0 # speed of light in vacuum
    8π^3 * ω^2 * sec(β)^2 / (c^3 * n * n1 * n2) * abs2(χ2_eff) * I1 * I2
end

"""
    sfgintensity(s::Setup, χ2_eff::Number; <keyword arguments>)

Same as `sfgintensity(ω, β, χ2_eff; <keyword arguments>)`.

# Arguments
* `s`:  `Setup` object
* `χ2_eff`: effective susceptibility

# Keyword Arguments
* `I1`: intensity of `ω1`
* `I2`: intensity of `ω2`

https://doi.org/10.1103/PhysRevB.59.12632
"""
function sfgintensity(s, χ2_eff; I1=1.0, I2=1.0)
    n  = s.Ω.n1
    n1 = s.Ω1.n1
    n2 = s.Ω2.n1
    sfgintensity(s.Ω.ω, s.Ω.β, χ2_eff; n=n, n1=n1, n2=n2, I1=I1, I2=I2)
end

"""
    sfangle(ω, ω1, ω2, β1, β2)

Compute the angle of the reflected sum frequency field.
In a counter-propagating geometry one of the `β`s gets negative.

# Arguments
* `ω`:  sum frequency
* `ω1`: frequency of wave 1
* `ω2`: frequency of wave 2
* `β1`:  incident angle for field with frequency `ω1`
* `β2`:  incident angle for field with frequency `ω2`
"""
function sfangle(ω, ω1, ω2, β1, β2)
    asin(ω1 / ω * sin(β1) + ω2 / ω * sin(β2))
end

"""
    wl2freq(λ)

Compute the frequency of a field with wavelength `λ` in nm.
"""
wl2freq(λ) = 299792458.0 / (λ * 1e-9)

"""
    freq2wl(ω)

Compute the wavelength of a field with frequency `ω` in 1/s.
The wavelength will be given in nm.
"""
freq2wl(ω) = 299792458.0 / ω * 1e9


text:
# test the vm2d_particle.advance() function
include("vm2d.jl")
timestep = 0.01

# first test the working case
# create vortex ring geometry
radius1 = 0.1
radius2 = 0.1
spacing = 0.05
points = [[radius1,0.0,0.0],[-radius1,0.0,0.0],[radius2,0.0,spacing],[-radius2,0.0,spacing]]

# set vectorial vortex strengths
gamma = 0.1
gamma_upper = [0,gamma,0]
gamma_lower = [0,-gamma,0]
gammas = [gamma_upper,gamma_lower,gamma_upper,gamma_lower]

# set initial velocities
v0 = 0.0
vs = [[v0,0.0,0.0],[v0,0.0,0.0],[0.0,0.0,0.0],[0.0,0.0,0.0]]

particles = place(points,gammas,vs)

# now test the trouble case

# create circle geometry
center = [0.0,0.0,0.0]
radius = 0.1
refinement = 3
gamma = -0.1
velocity=[0.0,0.0,0.0]
spacing = 0.05

ring1 = vortexring(radius,center,refinement,gamma;velocity=velocity)
ring2 = vortexring(radius,center.+[0.0,0.0,spacing],refinement,gamma;velocity=velocity)
particles2 = vcat(ring1,ring2)

println("===== test_advance() =====")
println("Manual case particles velocity BEFORE advance(): ")
for particle in particles
    println(particle.v)
end
println("")
println("Trouble case particles velocity BEFORE advance():")
for particle in particles2
    println(particle.v)
end
println("")
println("")

advance(particles,timestep,Uinf)
advance(particles2,timestep,Uinf)

println("Manual case particles velocity after advance(): ")
for particle in particles
    println(particle.v)
end
println("")
println("Trouble case particles velocity after advance():")
for particle in particles2
    println(particle.v)
end
println("")
println("")
println("===== END TEST =====")


text:
# Julia wrapper for header: /usr/local/include/lz4frame.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0
# docstrings copied from /usr/local/include/lz4frame.h

###########################
 # PRIVATE DEFINITIONS : Do Not Export
 # Do not use these definitions.
 # They are exposed to allow static allocation of `LZ4_streamHC_t`.
 # Using these definitions makes the code vulnerable to potential API break when upgrading LZ4
############################
const LZ4HC_DICTIONARY_LOGSIZE = 17
const LZ4HC_MAXD = 1 << LZ4HC_DICTIONARY_LOGSIZE
const LZ4HC_MAXD_MASK = LZ4HC_MAXD - 1
const LZ4HC_HASH_LOG = 15
const LZ4HC_HASHTABLESIZE = 1 << LZ4HC_HASH_LOG
const LZ4HC_HASH_MASK = LZ4HC_HASHTABLESIZE - 1
const LZ4_STREAMHCSIZE = 4LZ4HC_HASHTABLESIZE + 2LZ4HC_MAXD + 56
const LZ4_STREAMHCSIZE_SIZET = floor(Int, LZ4_STREAMHCSIZE / sizeof(Csize_t))

struct LZ4_streamHC_t
    table::NTuple{LZ4_STREAMHCSIZE_SIZET, Csize_t}
end

const LZ4_MEMORY_USAGE = 14
const LZ4_MAX_INPUT_SIZE = 0x7e000000

const LZ4_HASHLOG = LZ4_MEMORY_USAGE - 2
const LZ4_HASHTABLESIZE = 1 << LZ4_MEMORY_USAGE
const LZ4_HASH_SIZE_U32 = 1 << LZ4_HASHLOG

const LZ4_STREAMDECODESIZE_U64 = 4
const LZ4_STREAMSIZE_U64 =((1 << (LZ4_MEMORY_USAGE-3)) + 4)

struct LZ4_stream_t
    table::NTuple{LZ4_STREAMSIZE_U64, Culonglong}
end

struct LZ4_streamDecode_t
    table::NTuple{LZ4_STREAMDECODESIZE_U64, Culonglong}
end

# Constants
const LZ4F_VERSION = 100
const LZ4F_HEADER_SIZE_MAX = 19

# Block Size
const LZ4F_default = (UInt32)(0)
const LZ4F_max64KB = (UInt32)(4)
const LZ4F_max256KB = (UInt32)(5)
const LZ4F_max1MB = (UInt32)(6)
const LZ4F_max4MB = (UInt32)(7)

# Block Mode
const LZ4F_blockLinked = (UInt32)(0)
const LZ4F_blockIndependent = (UInt32)(1)

# Content Checksum Flag
const LZ4F_noContentChecksum = (UInt32)(0)
const LZ4F_contentChecksumEnabled = (UInt32)(1)

# Block Checksum Flag
const LZ4F_noBlockChecksum = (UInt32)(0)
const LZ4F_blockChecksumEnabled = (UInt32)(1)

# Frame Type
const LZ4F_frame = (UInt32)(0)
const LZ4F_skippableFrame = (UInt32)(1)

@enum BlockSizeID::Cuint default_size = 0 max64KB = 4 max256KB = 5 max1MB = 6 max4MB = 7
@enum BlockMode::Cuint block_linked = 0 block_independent = 1
@enum FrameType::Cuint normal_frame = 0 skippable_frame = 1

struct LZ4F_compressOptions_t
    stableSrc::Cuint
    reserved::NTuple{3, Cuint}
end

struct LZ4F_decompressOptions_t
    stableDst::Cuint
    reserved::NTuple{3, Cuint}
end

struct LZ4F_frameInfo_t
    blockSizeID::Cuint
    blockMode::Cuint
    contentChecksumFlag::Cuint
    frameType::Cuint
    contentSize::Culonglong
    dictID::Cuint
    blockChecksumFlag::Cuint
end

function LZ4F_frameInfo_t(;
    blocksizeid::BlockSizeID=default_size,
    blockmode::BlockMode=block_linked,
    contentchecksum::Bool=false,
    frametype::FrameType=normal_frame,
    contentsize::Integer=0,
    blockchecksum::Bool=false,
)
    LZ4F_frameInfo_t(Cuint(blocksizeid), Cuint(blockmode), Cuint(contentchecksum), Cuint(frametype), Culonglong(contentsize), Cuint(0), Cuint(blockchecksum))
end

struct LZ4F_preferences_t
    frameInfo::LZ4F_frameInfo_t
    compressionLevel::Cint
    autoFlush::Cuint
    reserved::NTuple{4, Cuint}
end

function LZ4F_preferences_t(frame_info::LZ4F_frameInfo_t; compressionlevel::Integer=0, autoflush::Bool=false)
    LZ4F_preferences_t(frame_info, Cint(compressionlevel), Cuint(autoflush), (0,0,0,0))
end

struct LZ4F_CDict
    dictContent::Ptr{Cvoid}
    fastCtx::Ptr{LZ4_stream_t}
    HCCtx::Ptr{LZ4_streamHC_t}
end

struct XXH32_state_t
   total_len_32::Cuint
   large_len::Cuint
   v1::Cuint
   v2::Cuint
   v3::Cuint
   v4::Cuint
   mem32::NTuple{4,UInt32}
   memsize::Cuint
   reserved::Cuint
end

struct LZ4F_cctx
    prefs::LZ4F_preferences_t
    version::UInt32
    cStage::UInt32
    cdict::Ptr{LZ4F_CDict}
    maxBlockSize::Csize_t
    maxBufferSize::Csize_t
    tmpBuff::Ptr{Cuchar}
    tmpIn::Ptr{Cuchar}
    tmpInSize::Csize_t
    totalInSize::UInt64
    xxh::XXH32_state_t
    lz4CtxPtr::Ptr{Cvoid}
    lz4CtxLevel::UInt32
end

struct LZ4F_dctx
    frameInfo::LZ4F_frameInfo_t
    version::UInt32
    dStage::UInt32
    frameRemainingSize::UInt64
    maxBlockSize::Csize_t
    maxBufferSize::Csize_t
    tmpIn::Ptr{Cuchar}
    tmpInSize::Csize_t
    tmpInTarget::Csize_t
    tmpOutBuffer::Ptr{Cuchar}
    dict::Ptr{Cuchar}
    dictSize::Csize_t
    tmpOut::Ptr{Cuchar}
    tmpOutSize::Csize_t
    tmpOutStart::Csize_t
    xxh::XXH32_state_t
    blockChecksum::XXH32_state_t
    header::NTuple{LZ4F_HEADER_SIZE_MAX, Cuchar}
end

function check_context_initialized(ctx::Ptr{LZ4F_cctx})
    if ctx == Ptr{LZ4F_cctx}(C_NULL)
        throw(LZ4Exception("LZ4F_cctx", "Uninitialized compression context"))
    end
end

function check_context_initialized(ctx::Ptr{LZ4F_dctx})
    if ctx == Ptr{LZ4F_dctx}(C_NULL)
        throw(LZ4Exception("LZ4F_dctx", "Uninitialized decompression context"))
    end
end

function LZ4F_isError(code::Csize_t)
    err = ccall((:LZ4F_isError, liblz4), UInt32, (Csize_t,), code)
    convert(Bool, err)
end

function LZ4F_getErrorName(code::Csize_t)
    str = ccall((:LZ4F_getErrorName, liblz4), Cstring, (Csize_t,), code)
    unsafe_string(str)
end

"""
Gets the current LZ4F version.
"""
function LZ4F_getVersion()
    ccall((:LZ4F_getVersion, liblz4), UInt32, ())
end

"""
The first thing to do is to create a compressionContext object, which will be used in all compression operations.
This is achieved using `LZ4F_createCompressionContext()`, which takes as argument a version.
The version provided MUST be the current version. It is intended to track potential version mismatch, notably when using DLL.
The function will provide a pointer to a fully allocated `LZ4F_cctx` object.
Will throw an error if there was an error during context creation.
"""
function LZ4F_createCompressionContext(cctxPtr::Ref{Ptr{LZ4F_cctx}}, version::UInt32)
    ret = ccall((:LZ4F_createCompressionContext, liblz4), Csize_t, (Ref{Ptr{LZ4F_cctx}}, UInt32), cctxPtr, version)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_createCompressionContext", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Releases the memory of a `LZ4F_cctx`.
"""
function LZ4F_freeCompressionContext(cctx::Ptr{LZ4F_cctx})
    ccall((:LZ4F_freeCompressionContext, liblz4), Csize_t, (Ptr{LZ4F_cctx},), cctx)
end

"""
Will write the frame header into `dstBuffer`.

`dstCapacity` must be >= `LZ4F_HEADER_SIZE_MAX` bytes.
`prefsPtr` is optional : you can provide `C_NULL` as argument, all preferences will then be set to default.
Returns the number of bytes written into `dstBuffer` for the header or throws an error.
 """
function LZ4F_compressBegin(cctx::Ptr{LZ4F_cctx}, dstBuffer, dstCapacity::Csize_t, prefsPtr::Ref{LZ4F_preferences_t})
    check_context_initialized(cctx)
    ret = ccall((:LZ4F_compressBegin, liblz4), Csize_t, (Ptr{LZ4F_cctx}, Ptr{Cvoid}, Csize_t, Ref{LZ4F_preferences_t}), cctx, dstBuffer, dstCapacity, prefsPtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_compressBegin" ,LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Provides minimum `dstCapacity` for a given `srcSize` to guarantee operation success in worst case situations.

`prefsPtr` is optional : when `C_NULL` is provided, preferences will be set to cover worst case scenario.
Result is always the same for a `srcSize` and `prefsPtr`, so it can be trusted to size reusable buffers.
When `srcSize==0`, `LZ4F_compressBound()` provides an upper bound for `LZ4F_flush()` and `LZ4F_compressEnd()` operations.
"""
function LZ4F_compressBound(srcSize::Csize_t, prefsPtr::Ref{LZ4F_preferences_t})
    ccall((:LZ4F_compressBound, liblz4), Csize_t, (Csize_t, Ref{LZ4F_preferences_t}), srcSize, prefsPtr)
end

"""
Can be called repetitively to compress as much data as necessary.

An important rule is that `dstCapacity` MUST be large enough to ensure operation success even in worst case situations.
This value is provided by `LZ4F_compressBound()`.
If this condition is not respected, `LZ4F_compress()` will fail.
LZ4F_compressUpdate() doesn't guarantee error recovery. When an error occurs, compression context must be freed or resized.
`cOptPtr` is optional : `C_NULL` can be provided, in which case all options are set to default.
Returns the number of bytes written into `dstBuffer` (it can be zero, meaning input data was just buffered).
or throws an error if it fails.
 """
function LZ4F_compressUpdate(cctx::Ptr{LZ4F_cctx}, dstBuffer, dstCapacity::Csize_t, srcBuffer, srcSize::Csize_t, cOptPtr)
    check_context_initialized(cctx)
    ret = ccall((:LZ4F_compressUpdate, liblz4), Csize_t, (Ptr{LZ4F_cctx}, Ptr{Cvoid}, Csize_t, Ptr{Cvoid}, Csize_t, Ptr{LZ4F_compressOptions_t}), cctx, dstBuffer, dstCapacity, srcBuffer, srcSize, cOptPtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_compressUpdate", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
When data must be generated and sent immediately, without waiting for a block to be completely filled,
it's possible to call `LZ4_flush()`. It will immediately compress any data buffered within `cctx`.
`dstCapacity` must be large enough to ensure the operation will be successful.
`cOptPtr` is optional : it's possible to provide `C_NULL`, all options will be set to default.
Returns the number of bytes written into `dstBuffer` (it can be zero, which means there was no data stored within `cctx`)
or throws an error if it fails.
"""
function LZ4F_flush(cctx::Ptr{LZ4F_cctx}, dstBuffer, dstCapacity::Csize_t, cOptPtr)
    check_context_initialized(cctx)
    ret = ccall((:LZ4F_flush, liblz4), Csize_t, (Ptr{LZ4F_cctx}, Ptr{Cvoid}, Csize_t, Ptr{LZ4F_compressOptions_t}), cctx, dstBuffer, dstCapacity, cOptPtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_flush", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Invoke to properly finish an LZ4 frame.

It will flush whatever data remained within `cctx` (like `LZ4_flush()`)
and properly finalize the frame, with an endMark and a checksum.
`cOptPtr` is optional : `C_NULL` can be provided, in which case all options will be set to default.
Returns the number of bytes written into `dstBuffer` (necessarily >= 4 (endMark), or 8 if optional frame checksum is enabled)
or throws an error if it fails (which can be tested using `LZ4F_isError()`)
A successful call to `LZ4F_compressEnd()` makes `cctx` available again for another compression task.
"""
function LZ4F_compressEnd(cctx::Ptr{LZ4F_cctx}, dstBuffer, dstCapacity::Csize_t, cOptPtr)
    check_context_initialized(cctx)
    ret = ccall((:LZ4F_compressEnd, liblz4), Csize_t, (Ptr{LZ4F_cctx}, Ptr{Cvoid}, Csize_t, Ptr{LZ4F_compressOptions_t}), cctx, dstBuffer, dstCapacity, cOptPtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_compressEnd", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Create an `LZ4F_dctx object`, to track all decompression operations.

The version provided MUST be the current LZ4F version.
The function provides a pointer to an allocated and initialized `LZ4F_dctx` object.
The the function throws an error if the `LZ4F_dctx` object cannot be initialized.
The `dctx` memory can be released using `LZ4F_freeDecompressionContext()`.
"""
function LZ4F_createDecompressionContext(dctxPtr::Ref{Ptr{LZ4F_dctx}}, version::UInt32)
    ret = ccall((:LZ4F_createDecompressionContext, liblz4), Csize_t, (Ref{Ptr{LZ4F_dctx}}, UInt32), dctxPtr, version)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_createDecompressionContext", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Frees the decompressionContext.

The result of `LZ4F_freeDecompressionContext()` is indicative of the current state of decompressionContext when being released.
That is, it should be == 0 if decompression has been completed fully and correctly.
"""
function LZ4F_freeDecompressionContext(dctx::Ptr{LZ4F_dctx})
    ccall((:LZ4F_freeDecompressionContext, liblz4), Csize_t, (Ptr{LZ4F_dctx},), dctx)
end

"""
Extracts frame parameters (max blockSize, dictID, etc.).

Its usage is optional.
Extracted information is typically useful for allocation and dictionary.
This function works in 2 situations :
 - At the beginning of a new frame, in which case
   it will decode information from `srcBuffer`, starting the decoding process.
   Input size must be large enough to successfully decode the entire frame header.
   Frame header size is variable, but is guaranteed to be <= `LZ4F_HEADER_SIZE_MAX` bytes.
   It's allowed to provide more input data than this minimum.
 - After decoding has been started.
   In which case, no input is read, frame parameters are extracted from dctx.
 - If decoding has barely started, but not yet extracted information from header,
   `LZ4F_getFrameInfo()` will fail.
The number of bytes consumed from srcBuffer will be updated within `srcSizePtr` (necessarily <= original value).
Decompression must resume from (`srcBuffer` + `srcSizePtr`).
Returns an hint about how many `srcSize` bytes `LZ4F_decompress()` expects for next call or throws an error.

Note 1 : In case of error, dctx is not modified. Decoding operation can resume from beginning safely.
Note 2 : Frame parameters are *copied into* an already allocated `LZ4F_frameInfo_t` structure.
"""
function LZ4F_getFrameInfo(dctx::Ptr{LZ4F_dctx}, frameInfoPtr::Ref{LZ4F_frameInfo_t}, srcBuffer, srcSizePtr)
    check_context_initialized(dctx)
    ret = ccall((:LZ4F_getFrameInfo, liblz4), Csize_t, (Ptr{LZ4F_dctx}, Ref{LZ4F_frameInfo_t}, Ptr{Cvoid}, Ref{Csize_t}), dctx, frameInfoPtr, srcBuffer, srcSizePtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_getFrameInfo", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Call this function repetitively to regenerate compressed data from `srcBuffer`.

The function will read up to `srcSizePtr` bytes from `srcBuffer`,
and decompress data into `dstBuffer`, of capacity `dstSizePtr`.
The number of bytes consumed from `srcBuffer` will be written into `srcSizePtr` (necessarily <= original value).
The number of bytes decompressed into `dstBuffer` will be written into `dstSizePtr` (necessarily <= original value).
The function does not necessarily read all input bytes, so always check value in `srcSizePtr`.
Unconsumed source data must be presented again in subsequent invocations.

`dstBuffer` can freely change between each consecutive function invocation.
`dstBuffer` content will be overwritten.

Returns an hint of how many `srcSize` bytes `LZ4F_decompress()` expects for next call.
Schematically, it's the size of the current (or remaining) compressed block + header of next block.
Respecting the hint provides some small speed benefit, because it skips intermediate buffers.
This is just a hint though, it's always possible to provide any srcSize.

When a frame is fully decoded, returns 0 (no more data expected).
When provided with more bytes than necessary to decode a frame,
`LZ4F_decompress()` will stop reading exactly at end of current frame, and return 0.

If decompression failed, an error is thrown.
After a decompression error, the `dctx` context is not resumable.
Use `LZ4F_resetDecompressionContext()` to return to clean state.
After a frame is fully decoded, dctx can be used again to decompress another frame.
"""
function LZ4F_decompress(dctx::Ptr{LZ4F_dctx}, dstBuffer, dstSizePtr::Ref{Csize_t}, srcBuffer, srcSizePtr::Ref{Csize_t}, dOptPtr)
    check_context_initialized(dctx)
    ret = ccall((:LZ4F_decompress, liblz4), Csize_t, (Ptr{LZ4F_dctx}, Ptr{Cvoid}, Ref{Csize_t}, Ptr{Cvoid}, Ref{Csize_t}, Ptr{LZ4F_decompressOptions_t}), dctx, dstBuffer, dstSizePtr, srcBuffer, srcSizePtr, dOptPtr)
    if LZ4F_isError(ret)
        throw(LZ4Exception("LZ4F_decompress", LZ4F_getErrorName(ret)))
    end
    ret
end

"""
Re-initializes decompression context

In case of an error, the context is left in "undefined" state.
In which case, it's necessary to reset it, before re-using it.
This method can also be used to abruptly stop any unfinished decompression,
and start a new one using same context resources.
 """
function LZ4F_resetDecompressionContext(dctx::Ptr{LZ4F_dctx})
    check_context_initialized(dctx)
    ccall((:LZ4F_resetDecompressionContext, liblz4), Cvoid, (Ptr{LZ4F_dctx},), dctx)
end


text:
module test_h3_api_indexing

using Test
using H3.API # geoToH3 h3ToGeo h3ToGeoBoundary GeoCoord GeoBoundary ≈

h = 0x85283473fffffff
@test h3GetResolution(h) == 5

location = GeoCoord(0.6518070561696664, -2.128889370371519)
@test geoToH3(location, 1) == 0x081283ffffffffff
@test geoToH3(location, 5) == h
@test geoToH3(location, 10) == 0x08a2834700007fff
@test h3ToGeo(h) == location

verts = h3ToGeoBoundary(0x85283473fffffff)
@test length(verts) == 6
@test verts ≈ [GeoCoord(0.6505078765569766, -2.1278195595404963),
               GeoCoord(0.6519490051151717, -2.126897030193998),
               GeoCoord(0.6532477872571462, -2.1279673831553825),
               GeoCoord(0.6531044519454293, -2.1299602868027208),
               GeoCoord(0.6516632883200013, -2.130879969983952),
               GeoCoord(0.6503654944795706, -2.129809601095088)]

# https://github.com/uber/h3/blob/master/examples/index.c
location = GeoCoord(deg2rad(40.689167), deg2rad(-74.044444))
indexed = geoToH3(location, 10)
verts = h3ToGeoBoundary(indexed)
@test map(vert -> rad2deg.((vert.lat, vert.lon)), verts) ≈ [(40.69005860095358, -74.04415176176158),
                                                               (40.68990769452519, -74.0450617923963),
                                                               (40.68927093604355, -74.04534141750702),
                                                               (40.68878509072403, -74.04471103053613),
                                                               (40.688935992642726, -74.04380102076254),
                                                               (40.689572744390524, -74.04352137709904)]
center = h3ToGeo(indexed)
@test (rad2deg(center.lat), rad2deg(center.lon)) == (40.68942184369929, -74.04443139990863)

@test geoToH3(GeoCoord(0, 0), 0) == 0x08075fffffffffff
@test geoToH3(GeoCoord(0, 0), 5) == 0x085754e67fffffff

# https://github.com/wookay/H3.jl/issues/2
@testset "Compact issue" begin
    hids = Array{UInt64, 1}()

    open(normpath(@__DIR__, "tmp.h3index"), "r") do fp
        while !eof(fp)
            push!(hids, read(fp, H3Index))
        end
    end

    l = length([h for h in compact(hids) if h3IsValid(h)])
    for i in 1:1000
        @test length([h for h in compact(hids) if h3IsValid(h)]) == l
    end
end

end # module test_h3_api_indexing


text:
using LinearAlgebra
@testset "Graph export and import" begin
    fnames = Vector{String}()
    tests = Vector{Any}(undef, 10)
    tests[1:5] = randn(ComplexF64, 5)
    for i = 1:5
        tests[5+i] = randn(ComplexF64, 20, 20)
    end

    den_coeffs = 10 * randn(7)
    num_coeffs = 10 * randn(7)

    (graph, cref) = graph_rational(den_coeffs, num_coeffs, graph_ps)
    fname = string(tempname(), ".cgr")
    export_compgraph(graph, fname)
    graph2 = import_compgraph(fname)
    for i = 1:10
        @test eval_graph(graph, tests[i]) ≈ eval_graph(graph2, tests[i])
    end
    rm(fname)

    graph = Compgraph(ComplexF64, graph)
    x = get_coeffs(graph)
    x = x + 0.1 * randn(ComplexF64, length(x))
    set_coeffs!(graph, x)
    graph1 = deepcopy(graph)

    add_output!(graph1, cref[1][1])
    add_output!(graph1, cref[6][1])
    add_output!(graph1, cref[12][1])

    fname = string(tempname(), ".cgr")
    order = get_topo_order(graph)[1]
    export_compgraph(
        graph1,
        fname,
        fun = "TEST",
        dom = "[-2,0.1]",
        err = "???",
        genby = basename(@__FILE__),
        main_output = :pade,
        order = order,
    )
    graph2 = import_compgraph(fname)
    rm(fname)
    for i = 1:10
        @test eval_graph(graph, tests[i]) != eval_graph(graph1, tests[i])
        @test eval_graph(graph, tests[i]) ≈ eval_graph(graph2, tests[i])
    end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

include("hdf5.jl")
include("bat_io.jl")

# Included via require: "hdf5_specific.jl"


text:
"""
    StokesOutput
Return type for Stokes computations.
Fields are `nothing` on return unless requested.
See individual FMM/direct computation function
documentation for specifics.
"""
mutable struct StokesOutput <: FMMVals
    pot
    pre
    grad

    pottarg
    pretarg
    gradtarg

    ifppregtarg
    ifppreg
end
function StokesOutput(ifppregtarg, ifppreg=nothing)
    StokesOutput(nothing,nothing,nothing,
                 nothing,nothing,nothing,
                 ifppregtarg, ifppreg)
end


text:
# Observer
abstract Observer

# Lazy Object
abstract LazyObject <: Observer

# Process
abstract StochasticProcess
abstract StochasticProcess1D <: StochasticProcess
abstract AbstractBlackScholesProcess <: StochasticProcess1D
abstract AbstractHestonProcess <: StochasticProcess1D
abstract AbstractDiscretization
abstract AbstractHestonDiscretization <: AbstractDiscretization
abstract BlackScholesType

# Methods
abstract Lattice
abstract TreeLattice <: Lattice
abstract AbstractTree
abstract AbstractBinomialTree <: AbstractTree
abstract BinomialTreeType
abstract EqualProbabilitiesBinomialTreeType <: BinomialTreeType
abstract EqualJumpsBinomialTreeType <: BinomialTreeType
abstract FdmSchemeDescType
abstract FdScheme
abstract FdmMesher
abstract Fdm1DMesher
abstract StepCondition
abstract CurveWrapper
abstract StepType
abstract FdmInnerValueCalculator
abstract FdmLinearOpComposite
abstract NinePointLinearOp
abstract BoundaryConditionType
abstract BC_Side
# Monte Carlo
abstract AbstractMonteCarloModel
abstract MCTrait
abstract AbstractPathPricer
abstract EarlyExercisePathPricer <: AbstractPathPricer
abstract LsmBasisSystemPolynomType
abstract LSMBasisSystemFunction
# abstract BoundaryCondition

# Exercise
abstract Exercise
abstract EarlyExercise <: Exercise

# Instruments
abstract Instrument <: LazyObject
abstract PositionType
abstract Bond <: Instrument
abstract AbstractCallableBond <: Bond
abstract AbstractRate <: Instrument
abstract Swap <: Instrument
abstract AbstractClaim
abstract SettlementType
abstract AbstractPayoff
abstract StrikedTypePayoff <: AbstractPayoff
abstract Option{E} <: Instrument
abstract OneAssetOption{E} <: Option{E}
abstract OptionType
abstract SwapType
abstract CDSProtectionSide
abstract Results
abstract CallType
abstract AbstractForward <: Instrument

# Term Structures
abstract TermStructure <: LazyObject
abstract YieldTermStructure <: TermStructure
abstract CreditTermStructure <: TermStructure
# Curves
abstract Curve <: YieldTermStructure
abstract InterpolatedCurve{P, T} <: Curve
abstract AbstractDefaultProbabilityTermStructure <: CreditTermStructure
abstract AbstractDefaultProbabilityCurve <: AbstractDefaultProbabilityTermStructure
abstract InterpolatedDefaultProbabilityCurve{P} <: AbstractDefaultProbabilityCurve

abstract VolatilityTermStructure <: TermStructure
abstract OptionletVolatilityStructure <: VolatilityTermStructure
abstract SwaptionVolatilityStructure <: VolatilityTermStructure
abstract BlackVolTermStructure <: VolatilityTermStructure
abstract LocalVolTermStructure <: VolatilityTermStructure
abstract CallableBondVolatilityStructure <: TermStructure
abstract AbstractSmileSection
abstract VolatilityType
abstract BootstrapTrait
abstract Bootstrap
abstract FittingMethod
abstract BootstrapHelper <: LazyObject
abstract BondHelper <: BootstrapHelper
abstract RateHelper <: BootstrapHelper
abstract AbstractCDSHelper <: BootstrapHelper

# Pricing Engines
abstract PricingEngine
abstract DiscretizedAsset
abstract DiscretizedOption <: DiscretizedAsset
abstract LatticeShortRateModelEngine{S} <: PricingEngine
abstract AbstractHestonEngine <: PricingEngine
abstract AbstractFDVanillaEngine <: PricingEngine
abstract FDMultiPeriodEngine <: AbstractFDVanillaEngine
abstract FDStepConditionEngine <: AbstractFDVanillaEngine
abstract AbstractVanillaEngine <: PricingEngine
abstract MCVanillaEngine <: PricingEngine
abstract MCLongstaffSchwartzEngine <: PricingEngine

# Cash Flows
abstract CashFlows
abstract CashFlow
abstract Coupon <: CashFlow
abstract Duration
abstract CouponPricer
abstract IborCouponPricer <: CouponPricer

# Indexes
abstract InterestRateIndex

# Models
abstract Parameter
abstract CalibrationErrorType
abstract CalibrationHelper <: LazyObject
abstract CalibrationBasketType
abstract ModelType
abstract Model{T <: ModelType} <: LazyObject
abstract ShortRateModel{T} <: Model{T}
abstract OneFactorModel{T} <: ShortRateModel{T}
abstract Gaussian1DModel{T} <: OneFactorModel{T}
abstract TwoFactorModel{T} <: ShortRateModel{T}
abstract ComplexLogFormula
abstract HestonIntegration
abstract ShortRateDynamics
abstract ShortRateTree <: TreeLattice
abstract AbstractMarketModel
abstract AbstractMarketModelEvolver
abstract ExerciseStrategy
abstract MarketModelMultiProduct
abstract MarketModelBasisSystem
abstract MarketModelPathwiseMultiProduct
abstract MultiProductMultiStep <: MarketModelMultiProduct
abstract MarketModelExerciseValue
abstract PiecewiseConstantCorrelation
abstract BrownianGeneratorFactory
abstract BrownianGenerator
abstract SobolOrdering
abstract CurveState

# Currencies
abstract AbstractCurrency


text:
"
Returns (LT, ER, (im*k^2*inv_w).*invV*conj(w_vec)), which connect the effective and average wave through α = LT*A + (im*k^2*inv_w).*invV*conj(w_vec).
The matching region is X[L:end].
"
function match_arrays(ω::T, wave_effs::Vector{EffectiveWave{T}}, L::Int, X::AbstractVector{T},
        medium::Medium{T}, species::Vector{Specie{T}};
        a12k::T = 1.005*T(2)*real(specie.r*ω/medium.c),
        scheme::Symbol = :trapezoidal, θin::T = 0.0) where T<:Number

    J = length(X) - 1
    k = ω/medium.c
    hos = union(w.hankel_order for w in wave_effs)

    if length(hos) > 1
        @warn("Expected all effective waves to have the same hankel order!")
    end
    if abs(1 - (X[2]-X[1])*J/X[J+1]) > eps(T)
        error("Expected a uniform mesh X starting at 0.")
    end

    ho = minimum(hos)
    S = length(species)

    t_vecs = t_vectors(ω, medium, species; hankel_order = ho)

    σ = integration_scheme(X[1:L]; scheme=scheme) # integration scheme: trapezoidal

    w_vec = T(2) .*
        [
            sum(
                exp(im*m*(θin - w.θ_eff) + im*X[L]*(w.k_eff*cos(w.θ_eff) - cos(θin))) *
                species[s].num_density *  w.amplitudes[m+ho+1,s]
            for m = -ho:ho, s = 1:S) / (cos(θin)*(w.k_eff*cos(w.θ_eff) - cos(θin)))
        for w in wave_effs]

    G_arr = [
        (j > L) ? zero(Complex{T}) :
            T(2) * (-im)^T(m-1) * species[s].num_density * exp(im*m*θin - im*X[j]*cos(θin)) * σ[j] / cos(θin)
    for j = 1:(J+1), m = -ho:ho, s = 1:S]

    avg_wave_effs = [AverageWave(X, wave) for wave in wave_effs]
    vs = [
        [w.amplitudes[j,n+ho+1,1] for w in avg_wave_effs]
    for j = L:(J+1), n = -ho:ho]

    invV = inv(sum(conj(vs)[inds] * transpose(vs[inds])  for inds in eachindex(vs)))
    inv_w = one(T)/(transpose(w_vec)*invV*conj(w_vec))

    YTs = [
            invV * [(j < L) ? zero(Complex{T}) : conj(w.amplitudes[j,n+ho+1,1]) for w in avg_wave_effs]
        for j = 1:(J+1), n = -ho:ho]
    YT = hcat(YTs...)

    LT_mat = YT + invV * (conj(w_vec).*inv_w) * (transpose(G_arr[:]) - transpose(w_vec)*YT)

    S_mat = OffsetArray{Complex{T}}(undef, 0:J, -2ho:2ho);
    for j = 0:J, m = -2ho:2ho
        S_mat[j,m] = integrate_S(m, X[j+1]; θin = θin)
    end
    Es = [
        [
            sum(
                (species[s].num_density/(k^2)) * t_vecs[s][m+ho+1] * im^T(n+1) * S_mat[J-l,n-m] *
                exp(im*X[J+1]*wave_effs[p].k_eff*cos(wave_effs[p].θ_eff) - im*n*wave_effs[p].θ_eff) *
                 wave_effs[p].amplitudes[n+ho+1] / (wave_effs[p].k_eff*cos(wave_effs[p].θ_eff) + cos(θin))
            for n = -ho:ho, s = 1:S)
        for l = 0:J, p in eachindex(wave_effs)]
    for m = -ho:ho]

    q = min(Int(floor(a12k/(X[2]-X[1]))), J)
    B_mat = OffsetArray{Complex{T}}(undef, 0:q, -2ho:2ho);
    for j = 0:q, m = -2ho:2ho
        B_mat[j,m] = integrate_B(m, X[j+1], sqrt(abs(a12k^2 -X[j+1]^2)); θin = θin)
    end
    XR = OffsetArray((J:(J+q))*(X[2]-X[1]), J:(J+q));
    # the integration scheme changes with the domain
    σs = OffsetArray(
        [OffsetArray(integration_scheme(XR[J:l+q]; scheme=scheme), J:(l+q)) for l = (J-q+1):J]
    , (J-q+1):J)

    Rs = [
        [
        (l+q <= J) ?
            zero(Complex{T}) :
            sum(
                (species[s].num_density/(k^2)) * t_vecs[s][m+ho+1] * im^T(n) * wave_effs[p].amplitudes[n+ho+1] *
                exp(im*XR[j]*wave_effs[p].k_eff*cos(wave_effs[p].θ_eff) - im*n*wave_effs[p].θ_eff) *
                (B_mat[j-l,n-m] - S_mat[j-l,n-m]) * σs[l][j]
            for j = J:(l+q), n = -ho:ho, s = 1:S)
        for l = 0:J, p in eachindex(wave_effs)]
    for m = -ho:ho]

    ER_mat = vcat((Es+Rs)...)

    return (LT_mat, ER_mat, (im*k^2*inv_w).*invV*conj(w_vec))
end


text:
using Test
using Zarr
using JSON
using Pkg
using PyCall

@testset "Zarr" begin

@testset "ZArray" begin
    @testset "fields" begin
        z = zzeros(Int, 2, 3)
        @test z isa ZArray{Int, 2, Zarr.BloscCompressor,
            Zarr.DictStore}

        @test z.storage.name === "data"
        @test length(z.storage.a) === 3
        @test length(z.storage.a["0.0"]) === 64
        @test eltype(z.storage.a["0.0"]) === UInt8
        @test z.metadata.shape[] === (2, 3)
        @test z.metadata.order === 'C'
        @test z.metadata.chunks === (2, 3)
        @test z.metadata.fill_value === nothing
        @test z.metadata.compressor isa Zarr.BloscCompressor
        @test z.metadata.compressor.blocksize === 0
        @test z.metadata.compressor.clevel === 5
        @test z.metadata.compressor.cname === "lz4"
        @test z.metadata.compressor.shuffle === true
        @test z.attrs == Dict{Any, Any}()
        @test z.writeable === true
    end

    @testset "methods" begin
        z = zzeros(Int, 2, 3)
        @test z isa ZArray{Int, 2, Zarr.BloscCompressor,
            Zarr.DictStore}

        @test eltype(z) === Int
        @test ndims(z) === 2
        @test size(z) === (2, 3)
        @test size(z, 2) === 3
        @test length(z) === 2 * 3
        @test lastindex(z, 2) === 3
        @test Zarr.zname(z) === "data"
    end

    @testset "NoCompressor DirectoryStore" begin
        mktempdir(@__DIR__) do dir
            name = "nocompressor"
            z = zzeros(Int, 2, 3, path="$dir/$name",
                compressor=Zarr.NoCompressor())

            @test z.metadata.compressor === Zarr.NoCompressor()
            @test z.storage === Zarr.DirectoryStore("$dir/$name")
            @test isdir("$dir/$name")
            @test ispath("$dir/$name/.zarray")
            @test ispath("$dir/$name/.zattrs")
            @test ispath("$dir/$name/0.0")
            @test JSON.parsefile("$dir/$name/.zattrs") == Dict{String, Any}()
            @test JSON.parsefile("$dir/$name/.zarray") == Dict{String, Any}(
                "dtype" => "<i8",
                "filters" => nothing,
                "shape" => [3, 2],
                "order" => "C",
                "zarr_format" => 2,
                "chunks" => [3, 2],
                "fill_value" => nothing,
                "compressor" => nothing)
            # call gc to avoid unlink: operation not permitted (EPERM) on Windows
            # might be because files are left open
            # from https://github.com/JuliaLang/julia/blob/f6344d32d3ebb307e2b54a77e042559f42d2ebf6/stdlib/SharedArrays/test/runtests.jl#L146
            GC.gc()
        end
    end
end

@testset "Metadata" begin
    @testset "Data type encoding" begin
        @test Zarr.typestr(Bool) === "<b1"
        @test Zarr.typestr(Int8) === "<i1"
        @test Zarr.typestr(Int64) === "<i8"
        @test Zarr.typestr(UInt32) === "<u4"
        @test Zarr.typestr(UInt128) === "<u16"
        @test Zarr.typestr(Complex{Float32}) === "<c8"
        @test Zarr.typestr(Complex{Float64}) === "<c16"
        @test Zarr.typestr(Float16) === "<f2"
        @test Zarr.typestr(Float64) === "<f8"
        @test Zarr.typestr(Zarr.MaxLengthString{5,UInt8}) === "<S5"
        @test Zarr.typestr(Zarr.MaxLengthString{9,UInt32}) === "<U9"
    end

    @testset "Metadata struct and JSON representation" begin
        A = fill(1.0, 30, 20)
        chunks = (5,10)
        metadata = Zarr.Metadata(A, chunks; fill_value=-1.5)
        @test metadata isa Zarr.Metadata
        @test metadata.zarr_format === 2
        @test metadata.shape[] === size(A)
        @test metadata.chunks === chunks
        @test metadata.dtype === "<f8"
        @test metadata.compressor === Zarr.BloscCompressor(0, 5, "lz4", true)
        @test metadata.fill_value === -1.5
        @test metadata.order === 'C'
        @test metadata.filters === nothing

        jsonstr = json(metadata)
        metadata_cycled = Zarr.Metadata(jsonstr)
        @test metadata == metadata_cycled
    end

    @testset "Fill value" begin
        @test Zarr.fill_value_encoding(Inf) === "Infinity"
        @test Zarr.fill_value_encoding(-Inf) === "-Infinity"
        @test Zarr.fill_value_encoding(NaN) === "NaN"
        @test Zarr.fill_value_encoding(nothing) === nothing
        @test Zarr.fill_value_encoding("-") === "-"

        @test Zarr.fill_value_decoding("Infinity", Float64) === Inf
        @test Zarr.fill_value_decoding("-Infinity", Float64) === -Inf
        @test Zarr.fill_value_decoding("NaN", Float32) === NaN32
        @test Zarr.fill_value_decoding("3.4", Float64) === 3.4
        @test Zarr.fill_value_decoding("3", Int) === 3
        @test Zarr.fill_value_decoding(nothing, Int) === nothing
        @test Zarr.fill_value_decoding("-", String) === "-"
        @test Zarr.fill_value_decoding("", Zarr.ASCIIChar) === nothing
        @test Zarr.fill_value_decoding("", Zarr.MaxLengthString{6,UInt8}) === Zarr.MaxLengthString{6,UInt8}("")
    end
end

@testset "getindex/setindex" begin
  a = zzeros(Int64, 10, 10, chunks = (5,2))
  a[2,:] .= 5
  a[:,3] .= 6
  a[9:10,9:10] .= 2
  a[5,5] = 1

  @test a[2,:] == [5, 5, 6, 5, 5, 5, 5, 5, 5, 5]
  @test a[:,3] == fill(6,10)
  @test a[4,4] == 0
  @test a[5:6,5:6] == [1 0; 0 0]
  @test a[9:10,9:10] == fill(2,2,2)
  # Now with FillValue
  amiss = zzeros(Int64, 10,10,chunks=(5,2), fill_value=-1)
  amiss[:,1] = 1:10
  amiss[:,2] .= missing
  amiss[1:3,4] = [1,missing,3]
  amiss[1,10] = 5
  amiss[1:5,9:10] .= missing

  @test amiss[:,1] == 1:10
  @test all(ismissing,amiss[:,2])
  @test all(i->isequal(i...),zip(amiss[1:3,4],[1,missing,3]))
  # Test that chunk containing only missings is not initialized
  @test !Zarr.isinitialized(amiss.storage,Zarr.citostring(CartesianIndex((1,5))))

end

@testset "resize" begin
  a = zzeros(Int64, 10, 10, chunks = (5,2), fill_value=-1)
  resize!(a,5,4)
  @test size(a)==(5,4)
  resize!(a,10,10)
  @test size(a)==(10,10)
  @test all(ismissing,a[6:end,:])
  xapp = rand(1:10,10,20)
  append!(a,xapp)
  @test size(a)==(10,30)
  @test a[:,11:30] == xapp
  singlevec = rand(1:10,10)
  append!(a,singlevec)
  @test size(a)==(10,31)
  @test a[:,31]==singlevec
  singlerow = rand(1:10,31)
  append!(a,singlerow,dims=1)
  @test size(a)==(11,31)
  @test a[11,:]==singlerow
  append!(a,vcat(singlerow', singlerow'), dims=1)
  @test size(a)==(13,31)
  @test a[12:13,:]==vcat(singlerow', singlerow')
end

@testset "string array getindex/setindex" begin
  using Zarr: MaxLengthString
  aa = ["this", "is", "all ", "ascii"]
  bb = ["And" "Unicode"; "ματριξ" missing]
  a = ZArray(aa)
  b = ZArray(bb, fill_value = MaxLengthString{7,UInt32}(""))
  @test eltype(a) == MaxLengthString{5,UInt8}
  @test eltype(b) == Union{MaxLengthString{7,UInt32},Missing}
  @test a[:] == ["this", "is", "all ", "ascii"]
  @test all(isequal.(b[:,:],["And" "Unicode"; "ματριξ" missing]))
end

include("storage.jl")

include("python.jl")

end  # @testset "Zarr"


text:
function MultivariateFunctionPenaltyAdapter(arg0::MultivariateFunction, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::jdouble, arg4::Vector{jdouble})
    return MultivariateFunctionPenaltyAdapter((MultivariateFunction, Vector{jdouble}, Vector{jdouble}, jdouble, Vector{jdouble}), arg0, arg1, arg2, arg3, arg4)
end

function value(obj::MultivariateFunctionPenaltyAdapter, arg0::Vector{jdouble})
    return jcall(obj, "value", jdouble, (Vector{jdouble},), arg0)
end


text:
using GeodesicLearning

## energy function
f(x) = sum( 3*(x .- 5).^2 )

pop = initialize_pop(f, 20, dim = 3) # initialize_pop with 20 points
##
x = GeodesicLearning.least_squares(pop)
entropy_invert(pop, log(20)/100, x)

max_stars_repo_name: chelate/GeodesicLearning

max_issues_repo_name: chelate/GeodesicLearning

max_forks_repo_name: chelate/GeodesicLearning


text:
# Some misc functions related to quadratic optimization problems
#
import LinearAlgebra: Symmetric, UniformScaling

@doc raw"""
```julia
check_first_order(∇f::AbstractVector{<:Real},
                  xstar::AbstractVector{<:Real},
                  bc::BoundConstraints{<:Real,1})

check_first_order(Q::Symmetric{<:Real},
                  q::AbstractVector{<:Real},
                  xstar::AbstractVector{<:Real},
                  bc::BoundConstraints{<:Real,1})
```

Check First-Order Conditions
(see [Bound Constrained Optimization slides](https://wiki.mcs.anl.gov/leyffer/images/0/01/07-bndCons.pdf))

If ``x^\star=\arg\min f(x), x\in[l,u]`` then:

```math
\partial_i f(x^\star) = \left\{\begin{array}{ll}
\ge 0, & \text{if } x^\star[i] = l[i] \\
= 0, & \text{if } l[i] \le x^\star[i] \le u[i] \\
\le 0, & \text{if } x^\star[i] = u[i] \\
\end{array}
\right.
```

This is equivalent to:
```math
x^\star = P_{[l,u]}(x^\star-\nabla f(x^\star))
```

According to the previous result, this function returns:
```math
\max \mid x^\star - P_{[l,u]}(x^\star-(Q.x^\star+q)) \mid
```

For a local stationary point this quantity must be null

The second function is a wrapper that computes ``∇f=Q.x^\star+q``
"""
function check_first_order(∇f::AbstractVector{<:Real},
                           xstar::AbstractVector{<:Real},
                           bc::BoundConstraints{<:Real,1})

    # the condition assumes feasible constraints
    @assert xstar ∈ bc

    @assert size(∇f)==size(xstar)

    # TODO: refactoring to avoid mem alloc
    # (-> inline computation of max(...sequence...)
    v = xstar-∇f
    v = project!(v,bc)

    maximum(abs.(xstar .- v))
end

function check_first_order(Q::Symmetric{<:Real},
                           q::AbstractVector{<:Real},
                           xstar::AbstractVector{<:Real},
                           bc::BoundConstraints{<:Real,1})
    # TODO: for BlasFloat use BLAS... to avoid unnecessary memory
    # allocs...
    ∇f=Q*xstar+q
    check_first_order(∇f,xstar,bc)
end


text:
mutable struct SeqRing{T <: FieldElem} <: Ring
    base_ring::PolyRing{T}

    function SeqRing{T}(R::PolyRing, cached::Bool = true) where T <: FieldElem
       if cached && haskey(SequenceID, R)
          return SequenceID[R]::SeqRing{T}
       else
          z = new{T}(R)
          if cached
            SequenceID[R] = z
          end
          return z
       end
    end
end

const SequenceID = Dict{Ring, Ring}()

function SequenceRing(R::Ring, s::AbstractString; cached::Bool = true)
    S, v = PolynomialRing(R, s)
    parent_obj = SeqRing{elem_type(base_ring(S))}(S, cached)
    return parent_obj, v
end

SequenceRing(R::PolyRing; cached::Bool = true) = SeqRing{elem_type(base_ring(R))}(R, cached)

function change_coeff_field(F::Field, a::SeqRing{T}) where T <: FieldElem
    R = parent(change_base_ring(F, zero(base_ring(a))))
    SeqRing{elem_type(F)}(R)
end

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

 mutable struct Seq{T <: FieldElem} <: RingElem
    terms::Vector{HyperTerm{T}}
    parent::SeqRing{T}

    Seq{T}() where T <: FieldElem = new{T}(Array{HyperTerm{T}}(undef, 0))

    function Seq{T}(b::Vector{HyperTerm{T}}) where T <: FieldElem
       z = new{T}(b)
       return z
    end
 end

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

parent_type(::Type{Seq{T}}) where T <: FieldElem = SeqRing{T}

elem_type(::Type{SeqRing{T}}) where T <: FieldElem = Seq{T}

base_ring(R::SeqRing{T}) where T <: FieldElem = R.base_ring
base_ring(a::Seq) = base_ring(parent(a))

parent(a::Seq) = a.parent

isdomain_type(::Type{Seq}) = false
isexact_type(::Type{Seq}) = true

Base.iszero(a::Seq) = iszero(length(a.terms))

Base.one(R::SeqRing) = R(1)
Base.zero(R::SeqRing) = R()

Base.one(s::Seq) = one(parent(s))
Base.zero(s::Seq) = zero(parent(s))

function Base.hash(a::Seq, h::UInt)
    b = 0x53dd43cd511044d1%UInt
    for t in a.terms
        b = xor(b, xor(hash(t, h), h))
        b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
    end
    return b
end

function Base.isequal(x::Seq{T}, y::Seq{T}) where {T <: FieldElem}
    parent(x) != parent(y) && return false
    length(x.terms) != length(y.terms) && return false
    for i in 1:length(x.terms)
       x.terms[i] != y.terms[i] && return false
    end
    return true
end

Base.:(==)(x::Seq{T}, y::Seq{T}) where {T <: FieldElem} = isequal(x, y)

canonical_unit(x::Seq) = x

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

function (a::SeqRing{T})() where {T <: FieldElem}
    S = elem_type(base_ring(base_ring(a)))
    z = Seq{S}()
    z.parent = a
    return z
end

function (a::SeqRing{T})(b::Seq{T}) where {T <: FieldElem}
    z = Seq{T}(deepcopy(b.terms))
    z.parent = a
    return z
end

function (a::SeqRing{T})(b::Union{Integer, Rational}) where {T <: FieldElem}
    F = FractionField(base_ring(a))
    return a(F(b))
end

function (a::SeqRing{T})(b::T) where {T <: FieldElem}
    parent(b) != base_ring(base_ring(a)) && error("Unable to coerce to sequence")
    F = FractionField(base_ring(a))
    return a(F(b))
end

function (a::SeqRing{T})(b::PolyElem{T}) where {T <: FieldElem}
    parent(b) != base_ring(a) && error("Unable to coerce to sequence")
    F = FractionField(parent(b))
    return a(F(b))
end

function (a::SeqRing{T})(b::FracElem{<:PolyElem{T}}) where {T <: FieldElem}
    base_ring(b) != base_ring(a) && error("Unable to coerce to sequence")
    # S = elem_type(base_ring(base_ring(a)))
    t = HyperTerm{T}(b)
    z = Seq{T}([t])
    z.parent = a
    return z
end

function (a::SeqRing{T})(b::Union{Integer, Rational}, c::PolyElem{T}) where {T <: FieldElem}
    F = FractionField(base_ring(a))
    return a(F(b), c)
end

function (a::SeqRing{T})(b::T, c::PolyElem{T}) where {T <: FieldElem}
    parent(b) != base_ring(base_ring(a)) && error("Unable to coerce to sequence")
    F = FractionField(base_ring(a))
    return a(F(b), c)
end

function (a::SeqRing{T})(b::PolyElem{T}, c::PolyElem{T}) where {T <: FieldElem}
    parent(b) != base_ring(a) && error("Unable to coerce to sequence")
    F = FractionField(parent(b))
    return a(F(b), c)
end

function (a::SeqRing{T})(b::FracElem{<:PolyElem{T}}, c::PolyElem{T}) where {T <: FieldElem}
    base_ring(b) != base_ring(a) && error("Unable to coerce to sequence")
    S = elem_type(base_ring(base_ring(a)))
    t = HyperTerm{S}(b, c)
    z = Seq{S}([t])
    z.parent = a
    return z
end

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

terms(a::Seq{T}) where {T <: FieldElem} = a.terms
nterms(a::Seq{T}) where {T <: FieldElem} = length(a.terms)

coeffs(a::Seq{T}) where {T <: FieldElem} = iszero(a) ? [FractionField(base_ring(a))(0)] : map(coeff, terms(a))

geometric_terms(a::Seq{T}) where {T <: FieldElem} = map(geometric_term, terms(a))
factorial_terms(a::Seq{T}) where {T <: FieldElem} = map(factorial_term, terms(a))

ishypergeometric(a::Seq{T}) where {T <: FieldElem} = nterms(a) == 1
isrational(a::Seq{T}) where {T <: FieldElem} = iszero(a) || (ishypergeometric(a) && isrational(first(terms(a))))

function Base.gcd(a::Seq{T}, b::Seq{T}) where {T <: FieldElem}
    check_parent(a, b)
    iszero(a) && iszero(b) && return one(a)
    r = reduce(gcd, [a.terms; b.terms])
    z = Seq{T}([r])
    z.parent = a.parent
    return z
end

function Base.div(a::Seq{}, b::Seq{T}) where {T <: FieldElem}
    check_parent(a, b)
    iszero(b) && error("Division by zero")
    !ishypergeometric(b) && error("Division only for hypergeometric sequences supported")
    z = Seq{T}([div(t, b.terms[1]) for t in a.terms])
    z.parent = a.parent
    return z
end

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

needs_parentheses(::Seq) = true
displayed_with_minus_in_front(::Seq) = false

function Base.show(io::IO, a::Seq)
    if iszero(a)
        show(io, 0)
    else
        show(io, a.terms[1])
        for t in a.terms[2:end]
            print(io, " + ")
            show(io, t)
        end
    end
end

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

function check_parent(a::Seq, b::Seq, throw::Bool = true)
    c = parent(a) != parent(b)
    c && throw && error("Incompatible polynomial rings in sequence operation")
    return !c
 end

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

function Base.:-(f::Seq{T}) where {T <: FieldElem}
    z = Seq{T}([-t for t in f.terms])
    z.parent = f.parent
    return z
end

function Base.:+(f::Seq{T}, g::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    terms = deepcopy(f.terms)
    for t in g.terms
        idx = findfirst(x->(issummable(x, t)), terms)
        if idx === nothing
            push!(terms, t)
        else
            terms[idx] += t
        end

    end
    filter!(!iszero, terms)
    isempty(terms) && return zero(f)
    z = Seq{T}(terms)
    z.parent = f.parent
    return z
end

function Base.:-(f::Seq{T}, g::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    terms = deepcopy(f.terms)
    for t in g.terms
        idx = findfirst(x->(issummable(x, t)), terms)
        if idx === nothing
            push!(terms, -t)
        else
            terms[idx] -= t
        end
    end
    filter!(!iszero, terms)
    isempty(terms) && return zero(f)
    z = Seq{T}(terms)
    z.parent = f.parent
    return z
end

function Base.:*(f::Seq{T}, g::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    seqs = Seq{T}[]
    for a in f.terms
        for b in g.terms
            z = Seq{T}([a*b])
            z.parent = f.parent
            push!(seqs, z)
        end
    end
    reduce(+, seqs, init=zero(f))
end

function Base.:^(f::Seq{T}, i::Integer) where T <: FieldElem
    isnegative(i) && error("Only nonnegative powers allowed")
    iszero(i) && return one(f)
    prod(f for _ in 1:i)
end

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

function zero!(f::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    f.terms = Array{HyperTerm{T}}(undef, 0)
    f
end

function add!(f::Seq{T}, g::Seq{T}, h::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    f.terms = (g+h).terms
    f
end

function mul!(f::Seq{T}, g::Seq{T}, h::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    f.terms = (g*h).terms
    f
end

function addeq!(f::Seq{T}, g::Seq{T}) where T <: FieldElem
    check_parent(f, g)
    f.terms = (f+g).terms
    f
end

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

function (a::Seq{T})(b::Integer) where T <: FieldElem
    sum(t(b) for t in a.terms)
end

function (a::Seq{T})(b::Nemo.fmpq_poly) where T <: FieldElem
    # base_ring(a) != parent(b) && error("Polynomial rings do not match")
    (Nemo.degree(b) != 1 || lead(b) != 1 || denominator(Nemo.coeff(b, 0)) != 1) && error("Argument must be of form n+c where c is an integer")
    z = Seq{T}([t(b) for t in a.terms])
    z.parent = a.parent
    return z
end

function change_coeff_field(F::Field, a::Seq{T}) where T <: FieldElem
    ts = [change_coeff_field(F, t) for t in terms(a)]
    z = Seq{elem_type(F)}(ts)
    z.parent = change_coeff_field(F, a.parent)
    return z
end


text:
module SpecialMatrices

using Compat
using Compat.LinearAlgebra
using FFTW

import Compat.LinearAlgebra: Matrix, inv
if VERSION >= v"0.7.0"
    import Compat.LinearAlgebra: mul!
end

import Base: getindex, isassigned, size, *
import Base.\

include("cauchy.jl") #Cauchy matrix
include("companion.jl") #Companion matrix
include("frobenius.jl") #Frobenius matrix
include("hankel.jl") #Hankel matrix
include("hilbert.jl") #Hilbert matrix
include("kahan.jl") #Kahan matrix
include("riemann.jl") #Riemann matrix
include("strang.jl") #Strang matrix
include("toeplitz.jl") #Toeplitz matrix, Circulant matrix
include("vandermonde.jl") #Vandermonde matrix

end # module


text:
using DashBootstrapComponents

list_group = dbc_listgroup([
    dbc_listgroupitem([
        dbc_listgroupitemheading("This item has a heading"),
        dbc_listgroupitemtext("And some text underneath"),
    ]),
    dbc_listgroupitem([
        dbc_listgroupitemheading("This item also has a heading"),
        dbc_listgroupitemtext("And some more text underneath too"),
    ]),
]);

max_stars_repo_name: harisbal/dash-bootstrap-components

max_issues_repo_name: harisbal/dash-bootstrap-components

max_forks_repo_name: harisbal/dash-bootstrap-components


text:
@testset "Core" begin
    A = reshape([0], 1, 1, 1, 1, 1)
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedAdjIndex(1, (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1,)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,), (1.0,)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1,0)), ndims(A))
    @test @inferred(A[wis...]) === 0
    wis = ntuple(d->Interpolations.WeightedArbIndex((1,1), (1.0,0.0)), ndims(A))
    @test @inferred(A[wis...]) === 0.0
end


text:
using BlockMaps
using Test

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


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    RandomFolding(k, shuffle=true)

A method for creating `k` random folds from a spatial object.
Optionally `shuffle` the object before creating the folds.
"""
struct RandomFolding <: FoldingMethod
  k::Int
  shuffle::Bool
end

RandomFolding(k::Int) = RandomFolding(k, true)

function folds(object, method::RandomFolding)
  # retrieve parameters
  k, shuffle = method.k, method.shuffle

  # partition object
  p = partition(object, RandomPartition(k, shuffle))
  s = subsets(p)
  n = length(p)

  function pair(i)
    # source and target subsets
    source = [1:i-1; i+1:n]
    target = [i]

    # indices within subsets
    sinds = reduce(vcat, s[source])
    tinds = reduce(vcat, s[target])

    sinds, tinds
  end

  (pair(i) for i in 1:n)
end


text:
value(x::ReverseDiff.TrackedReal)  = x.value
value(x::ReverseDiff.TrackedArray) = x.value

# Support adaptive with non-tracked time
@inline function ODE_DEFAULT_NORM(u::ReverseDiff.TrackedArray,t) where {N}
  sqrt(sum(abs2,value(u)) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ReverseDiff.TrackedReal,N},t) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::Array{<:ReverseDiff.TrackedReal,N},t) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip((value(x) for x in u),Iterators.repeated(t))) / length(u))
end
@inline ODE_DEFAULT_NORM(u::ReverseDiff.TrackedReal,t) = abs(value(u))

# Support TrackedReal time, don't drop tracking on the adaptivity there
@inline function ODE_DEFAULT_NORM(u::ReverseDiff.TrackedArray,t::ReverseDiff.TrackedReal) where {N}
  sqrt(sum(abs2,u) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::AbstractArray{<:ReverseDiff.TrackedReal,N},t::ReverseDiff.TrackedReal) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
end
@inline function ODE_DEFAULT_NORM(u::Array{<:ReverseDiff.TrackedReal,N},t::ReverseDiff.TrackedReal) where {N}
  sqrt(sum(x->ODE_DEFAULT_NORM(x[1],x[2]),zip(u,Iterators.repeated(t))) / length(u))
end
@inline ODE_DEFAULT_NORM(u::ReverseDiff.TrackedReal,t::ReverseDiff.TrackedReal) = abs(u)

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0::ReverseDiff.TrackedArray,p::ReverseDiff.TrackedArray,args...;kwargs...)
  ReverseDiff.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0,p::ReverseDiff.TrackedArray,args...;kwargs...)
  ReverseDiff.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

function solve_up(prob::DiffEqBase.DEProblem,sensealg::Union{AbstractSensitivityAlgorithm,Nothing},u0::ReverseDiff.TrackedArray,p,args...;kwargs...)
  ReverseDiff.track(solve_up,prob,sensealg,u0,p,args...;kwargs...)
end

ReverseDiff.@grad function solve_up(prob,sensealg,u0,p,args...;kwargs...)
  out = _solve_adjoint(prob,sensealg,ReverseDiff.value(u0),ReverseDiff.value(p),args...;kwargs...)
  Array(out[1]),out[2]
end


text:
using PrimaryRecursive
using Test

include("arithmetic.jl")
include("relation.jl")


text:
#=
  Package: FishEBM.jl
  File: utilities.jl
  Devin Rose
  Utilities for use of FishEBM, plots, memory allocation, etc. during runtime.
  Created: March 2016
=#

"""
  Description: Allocates enough memory for a population density matrix and
    sets all elements in the array to zero.

  Returns: Array(Int64, size of habitat)

  Last update: May 2016
"""
function initPopulationDensity(enviro_a::EnvironmentAssumptions)
  popDensity = Array(Int64, size(enviro_a.habitat)[1], size(enviro_a.habitat)[2])

  for i = 1:size(enviro_a.habitat)[1]
    for j = 1:size(enviro_a.habitat)[2]
      popDensity[i, j] = 0
    end
  end

  return popDensity
end

"""
  Description: Updates a population density matrix in the simulate function.
    The population density is used for visualizing the movement of fish during
    the simulation to monitor the behaviour.

  Precondition: A population density matrix is already allocated prior to
    calling this function.

  Returns: operates directly on pop_density

  Last update: June 2016
"""
function updatePopulationDensity!(agent_db::Vector, pop_density::Array, week_total::Int64)
  totalPop = 0

  # get adult age vector
  ageVec = fill(0, length(agent_db[1].alive))
  for j = 1:length(ageVec)
    ageVec[j] = getAge(week_total, agent_db[1].weekNum[j])
  end

  for k = 1:length(agent_db)
    pop_density[agent_db[k].locationID] = 1

    if isEmpty(agent_db[k]) == false
      for m = 1:length(ageVec)
        if ageVec[m] != 0
          pop_density[agent_db[k].locationID] += agent_db[k].alive[m]
        end
        totalPop += agent_db[k].alive[m]
      end
    end
  end

  return totalPop
end


text:
export SWAP, FLIP
export ROT, IROT
export INC, DEC, NEG, INV, AddConst, SubConst
export HADAMARD

"""
    NoGrad{T} <: IWrapper{T}
    NoGrad(x)

A `NoGrad(x)` is equivalent to `GVar^{-1}(x)`, which cancels the `GVar` wrapper.
"""
@pure_wrapper NoGrad

const NullType{T} = Union{NoGrad{T}, Partial{T}}

# TODO: deprecate
@selfdual -

NEG(a!) = -(a!)
@selfdual NEG

INV(a!) = inv(a!)
@selfdual INV

@inline FLIP(b::Bool) = !b
@selfdual FLIP

"""
    INC(a!) -> a! + 1
"""
@inline function INC(a!::Number)
    a! + one(a!)
end

"""
    DEC(a!) -> a! - 1
"""
@inline function DEC(a!::Number)
    a! - one(a!)
end
@dual INC DEC

"""
    SWAP(a!, b!) -> b!, a!
"""
@inline function SWAP(a!::T, b!::T) where T
    b!, a!
end
@selfdual SWAP

"""
    ROT(a!, b!, θ) -> a!', b!', θ

```math
\\begin{align}
    {\\rm ROT}(a!, b!, \\theta)  = \\begin{bmatrix}
        \\cos(\\theta) & - \\sin(\\theta)\\\\
        \\sin(\\theta)  & \\cos(\\theta)
    \\end{bmatrix}
    \\begin{bmatrix}
        a!\\\\
        b!
    \\end{bmatrix},
\\end{align}
```
"""
@inline function ROT(i::Real, j::Real, θ::Real)
    a, b = rot(i, j, θ)
    a, b, θ
end

"""
    IROT(a!, b!, θ) -> ROT(a!, b!, -θ)
"""
@inline function IROT(i::Real, j::Real, θ::Real)
    i, j, _ = ROT(i, j, -θ)
    i, j, θ
end
@dual ROT IROT

"""
    HADAMARD(x::Real, y::Real)

Hadamard transformation that returns `(x + y)/√2, (x - y)/√2`
"""
function HADAMARD(x::Real, y::Real)
    sqrt(0.5) * (x + y), sqrt(0.5) * (x - y)
end

@selfdual HADAMARD

# more data views
for (DT, OP, NOP) in [(:AddConst, :+, :-), (:SubConst, :-, :+)]
    @eval struct $DT{T}
        x::T
    end

    @eval function (f::$DT)(y::Real)
        $OP(y, f.x)
    end

    @eval NiLangCore.chfield(x::T, ac::$DT, xval::T) where T<:Real = $NOP(xval, ac.x)
end

for F1 in [:(Base.:-), :NEG, :(ac::AddConst), :(sc::SubConst)]
    @eval @inline function $F1(a!::NullType)
        @instr $F1(a! |> value)
        a!
    end
end

for F2 in [:SWAP, :HADAMARD, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    @eval @inline function $F2(a::NullType, b::Real)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b)
        a, b
    end
    @eval @inline function $F2(a::NullType, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a |> value, b |> value)
        a, b
    end
    @eval @inline function $F2(a::Real, b::NullType)
        @instr $(NiLangCore.get_argname(F2))(a, b |> value)
        a, b
    end
end

function type_except(::Type{TT}, ::Type{T2}) where {TT, T2}
    N = length(TT.parameters)
    setdiff(Base.Iterators.product(zip(TT.parameters, repeat([T2], N))...), [ntuple(x->T2, N)])
end

for F3 in [:ROT, :IROT, :((inf::PlusEq)), :((inf::MinusEq)), :((inf::XorEq))]
    PS = (:a, :b, :c)
    for PTS in type_except(Tuple{NullType, NullType, NullType}, Real)
        params = map((P,PT)->PT <: NullType ? :($P |> value) : P, PS, PTS)
        params_ts = map((P,PT)->:($P::$PT), PS, PTS)
        @eval @inline function $F3($(params_ts...))
            @instr $F3($(params...))
            ($(PS...),)
        end
    end
end

# patch for fixed point numbers
function (f::PlusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!+T(x)/y, x, y
end

function (f::MinusEq{typeof(/)})(out!::T, x::Integer, y::Integer) where T<:Fixed
    out!-T(x)/y, x, y
end

for F in [:exp, :log, :sin, :sinh, :asin, :cos, :cosh, :acos, :tan, :tanh, :atan]
    @eval Base.$F(x::Fixed43) = Fixed43($F(Float64(x)))
    @eval (f::PlusEq{typeof($F)})(out!::Fixed43, x::Real) = out! + Fixed43($F(x)), x
    @eval (f::MinusEq{typeof($F)})(out!::Fixed43, x::Real) = out! - Fixed43($F(x)), x
end

Base.:^(x::Integer, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::Fixed43, y::Fixed43) = Fixed43(x^(Float64(y)))
Base.:^(x::T, y::Fixed43) where T<:AbstractFloat = x^(T(y))

function (::PlusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! + convert(T, y), y
end

function (::MinusEq{typeof(convert)})(out!::T, y) where T<:Real
    out! - convert(T, y), y
end

Base.:~(ac::AddConst) = SubConst(ac.x)
Base.:~(ac::SubConst) = AddConst(ac.x)
@dualtype AddConst SubConst

for F in [:INV, :NEG, :FLIP, :INC, :DEC]
    @eval NiLangCore.chfield(x::T, ::typeof($F), xval::T) where T<:Real = (~$F)(xval)
end


text:
using Quiqbox: checkBSList

@testset "Library.jl tests" begin
    @test try checkBSList(); true catch; false end
end


text:
using Documenter, Downloads, ArviZ
using MCMCChains: MCMCChains # make `from_mcmcchains` available for API docs
using SampleChains: SampleChains # make `from_samplechains` available for API docs
using PlutoStaticHTML: PlutoStaticHTML

const DOCS_SRC_PATH = joinpath(@__DIR__, "src")

# generate markdown from Pluto notebooks
output_format = PlutoStaticHTML.documenter_output
build_opts = PlutoStaticHTML.BuildOptions(
    DOCS_SRC_PATH;
    previous_dir=DOCS_SRC_PATH,
    output_format=output_format,
    add_documenter_css=false,
)
PlutoStaticHTML.build_notebooks(build_opts)

const ASSETS_DIR = joinpath(@__DIR__, "src", "assets")
const ARVIZ_ASSETS_URL = "https://raw.githubusercontent.com/arviz-devs/arviz-project/main/arviz_logos"

function download_asset(remote_fn, fn=remote_fn)
    mkpath(ASSETS_DIR)
    return Downloads.download(
        joinpath(ARVIZ_ASSETS_URL, remote_fn), joinpath(ASSETS_DIR, fn); verbose=true
    )
end

# download arviz-devs org logo assets
download_asset("ArviZ.png", "logo.png")
download_asset("ArviZ_white.png", "logo-dark.png")
download_asset("favicon.ico")

makedocs(;
    modules=[ArviZ],
    sitename="ArviZ.jl",
    pages=[
        "Home" => "index.md",
        "Quickstart" => "quickstart.md",
        "API" => "api.md",
        "Example Gallery" =>
            ["Matplotlib" => "mpl_examples.md", "Bokeh" => "bokeh_examples.md"],
        "Reference" => "reference.md",
    ],
    checkdocs=:exports,
    format=Documenter.HTML(;
        prettyurls=haskey(ENV, "CI"),
        assets=["assets/favicon.ico", "assets/custom.css"],
        sidebar_sitename=false,
        canonical="stable",
    ),
    linkcheck=true,
    analytics="G-W1G68W77YV",
)

deploydocs(; repo="github.com/arviz-devs/ArviZ.jl.git", devbranch="main", push_preview=true)


text:
using Flux, Statistics
using Flux.Data: DataLoader
using Flux: onehotbatch, onecold, logitcrossentropy, throttle, @epochs
using Base.Iterators: repeated
using Parameters: @with_kw
using CUDAapi
using MLDatasets
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw mutable struct Args
    η::Float64 = 3e-4       # learning rate
    batchsize::Int = 1024   # batch size
    epochs::Int = 10        # number of epochs
    device::Function = gpu  # set as gpu, if gpu available
end

function getdata(args)
    # Loading Dataset
    xtrain, ytrain = MLDatasets.MNIST.traindata(Float32)
    xtest, ytest = MLDatasets.MNIST.testdata(Float32)

    # Reshape Data for flatten the each image into linear array
    xtrain = Flux.flatten(xtrain)
    xtest = Flux.flatten(xtest)

    # One-hot-encode the labels
    ytrain, ytest = onehotbatch(ytrain, 0:9), onehotbatch(ytest, 0:9)

    # Batching
    train_data = DataLoader(xtrain, ytrain, batchsize=args.batchsize, shuffle=true)
    test_data = DataLoader(xtest, ytest, batchsize=args.batchsize)

    return train_data, test_data
end

function build_model(; imgsize=(28,28,1), nclasses=10)
    return Chain(
 	    Dense(prod(imgsize), 32, relu),
            Dense(32, nclasses))
end

function loss_all(dataloader, model)
    l = 0f0
    for (x,y) in dataloader
        l += logitcrossentropy(model(x), y)
    end
    l/length(dataloader)
end

function accuracy(data_loader, model)
    acc = 0
    for (x,y) in data_loader
        acc += sum(onecold(cpu(model(x))) .== onecold(cpu(y)))*1 / size(x,2)
    end
    acc/length(data_loader)
end

function train(; kws...)
    # Initializing Model parameters
    args = Args(; kws...)

    # Load Data
    train_data,test_data = getdata(args)

    # Construct model
    m = build_model()
    train_data = args.device.(train_data)
    test_data = args.device.(test_data)
    m = args.device(m)
    loss(x,y) = logitcrossentropy(m(x), y)

    ## Training
    evalcb = () -> @show(loss_all(train_data, m))
    opt = ADAM(args.η)

    @epochs args.epochs Flux.train!(loss, params(m), train_data, opt, cb = evalcb)

    @show accuracy(train_data, m)

    @show accuracy(test_data, m)
end

cd(@__DIR__)
train()

max_stars_repo_name: natema/model-zoo

max_issues_repo_name: natema/model-zoo

max_forks_repo_name: natema/model-zoo