text:
module ReactiveMPMathTest

using Test
using ReactiveMP
using Random

@testset "Math" begin

    @testset "tiny/huge" begin

        @test typeof(tiny) === TinyNumber
        @test typeof(huge) === HugeNumber

        @test convert(Float32, tiny)  == 1f-6
        @test convert(Float64, tiny)  == 1e-12
        @test convert(BigFloat, tiny) == big"1e-24"

        @test convert(Float32, huge)  == 1f+6
        @test convert(Float64, huge)  == 1e+12
        @test convert(BigFloat, huge) == big"1e+24"

        @test @inferred clamp(1f0, tiny, huge)  == 1f0
        @test @inferred clamp(0f0, tiny, huge)  == 1f-6
        @test @inferred clamp(1f13, tiny, huge) == 1f+6

        @test @inferred clamp(1.0, tiny, huge)  == 1.0
        @test @inferred clamp(0.0, tiny, huge)  == 1e-12
        @test @inferred clamp(1e13, tiny, huge) == 1e12

        @test @inferred clamp(big"1.0", tiny, huge)  == big"1.0"
        @test @inferred clamp(big"0.0", tiny, huge)  == big"1e-24"
        @test @inferred clamp(big"1e25", tiny, huge) == big"1e+24"

        for a in (1, 1.0, 0, 0.0, 1f0, 0f0, Int32(0), Int32(1), big"1", big"1.0", big"0", big"0.0")
            T = typeof(a)
            for v in [ tiny, huge ]
                V = typeof(v)

                for op in [ +, -, *, /, >, >=, <, <= ]
                    @test @inferred op(a, v) == op(a, convert(promote_type(T, V), v))
                    @test @inferred op(v, a) == op(convert(promote_type(T, V), v), a)

                    @test @inferred op(a, v) == op(a, v)
                    @test @inferred op(v, a) == op(v, a)
                end

                @test v <= (@inferred clamp(a, v, Inf)) <= Inf
                @test zero(a) <= (@inferred clamp(a, zero(a), v)) <= v

                for size in  [ 3, 5 ]
                    for array in [ fill(a, (size, )), fill(a, (size, size)) ]

                        for op in [ +, -, *, /, >, >=, <, <= ]
                            @test @inferred op.(array, v) == op.(array, convert(promote_type(T, V), v))
                            @test @inferred op.(v, array) == op.(convert(promote_type(T, V), v), array)

                            @test @inferred op.(array, v) == op.(array, v)
                            @test @inferred op.(v, array) == op.(v, array)
                        end

                        @test @inferred clamp.(array, v, Inf)         == clamp.(array, v, Inf)
                        @test @inferred clamp.(array, zero(array), v) == clamp.(array, zero(array), v)
                    end
                end
            end

        end

    end

end

end


text:
"""

For managing a pool of, potentially straggling, workers that communicate with the coordinator using MPI.
"""
module MPIAsyncPools

using MPI

export MPIAsyncPool, waitall!

"""

    MPIAsyncPool(ranks::Vector{<:Integer}; epoch0::Integer=0, nwait::Integer=length(ranks))

Used to manage a pool of, potentially straggling, workers that communicate with the coordinator
using MPI. `nwait` is the default number of workers to wait for in calls to `asyncmap!` and
`epoch0` is the epoch of the first iteration.

```julia
MPIAsyncPool(n)             # create a pool composed of n workers with ranks 1:n
MPIAsyncPool([1, 4, 5])     # create a pool composed of the 3 workers with ranks [1, 4, 5]
```
"""
mutable struct MPIAsyncPool
    ranks::Vector{Int} # MPI ranks of the workers managed by the pool
    responded::Vector{Bool} # for each worker, indicates if it has responded in the current epoch
    sreqs::Vector{MPI.Request} # MPI send requests for each worker
    rreqs::Vector{MPI.Request} # MPI receive requests for each worker
    sepochs::Vector{Int} # epoch that transmission was initiated at for each worker
    repochs::Vector{Int} # send epoch corresponding to the most recently received result for each worker
    active::Vector{Bool} # false for workers with no current computing task
    stimestamps::Vector{Int} # time when transmission was initiated for each worker
    latency::Vector{Float64} # latency of individual workers
    nwait::Int # default number of workers to wait for
    epoch::Int # current epoch
    function MPIAsyncPool(ranks::Vector{<:Integer}; epoch0::Integer=0, nwait::Integer=length(ranks))
        n = length(ranks)
        new(copy(ranks), zeros(Bool, n),
            Vector{MPI.Request}(undef, n), Vector{MPI.Request}(undef, n),
            Vector{Int}(undef, n), fill(epoch0, n),
            zeros(Bool, n),
            zeros(UInt64, n), zeros(Float64, n),
            nwait, epoch0)
    end
end

MPIAsyncPool(n::Integer, args...; kwargs...) = MPIAsyncPool(collect(1:n), args...; kwargs...)

"""
    asyncmap!(pool::MPIAsyncPool, sendbuf::AbstractArray, recvbuf::AbstractArray, isendbuf::AbstractArray, irecvbuf::AbstractArray, comm::MPI.Comm; nwait::Union{<:Integer,Function}=pool.nwait, epoch::Integer=pool.epoch+1, tag::Integer=0)

Send the data in `sendbuf` asynchronously (via `MPI.Isend`) to all workers and wait for some of
them to respond (via a corresponding `MPI.Isend`). If `nwait` is an integer, this function returns
when at least `nwait` workers have responded. On the other hand, if `nwait` is a function, this
function returns once `nwait(epoch, repochs)` evaluates to `true`. The second argument `repochs` is
a vector of length equal to the number of workers in the pool, where the `i`-th element is the
`epoch` that transmission to the corresponding worker was initiated at.

Returns the `repochs` vector. Similarly to the `MPI.Gather!` function, the elements of `recvbuf`
are partitioned into a number of chunks equal to the number of workers (hence, the length of
`recvbuf` must be divisible by the number of workers) and the data recieved from the `j`-th worker
is stored in the `j`-th partition.

The buffers `isendbuf` and `irecvbuf` are for internal use by this function only and should never
be changed or accessed outside of it. The length of `isendbuf` must be equal to the length of
`sendbuf` multiplied by the number of workers, and `irecvbuf` must have length equal to that of
`recvbuf`.
"""
function Base.asyncmap!(pool::MPIAsyncPool, sendbuf::AbstractArray, recvbuf::AbstractArray, isendbuf::AbstractArray, irecvbuf::AbstractArray, comm::MPI.Comm; nwait::Union{<:Integer,Function}=pool.nwait, epoch::Integer=pool.epoch+1, tag::Integer=0, recvf::Union{Nothing,Function}=nothing, latency_tol::Real=1e-2)
    comm_size = length(pool.ranks)
    if typeof(nwait) <: Integer
        0 <= nwait <= comm_size || throw(ArgumentError("nwait must be in the range [0, length(pool.ranks)], but is $nwait"))
    end
    isbitstype(eltype(sendbuf)) || throw(ArgumentError("The eltype of sendbuf must be isbits, but is $(eltype(sendbuf))"))
    isbitstype(eltype(recvbuf)) || throw(ArgumentError("The eltype of sendbuf must be isbits, but is $(eltype(recvbuf))"))
    sizeof(isendbuf) == comm_size*sizeof(sendbuf) || throw(DimensionMismatch("sendbuf is of size $(sizeof(sendbuf)) bytes, but isendbuf is of size $(sizeof(isendbuf)) bytes when $(comm_size*sizeof(sendbuf)) bytes are needed"))
    sizeof(recvbuf) == sizeof(irecvbuf) || throw(DimensionMismatch("recvbuf is of size $(sizeof(recvbuf)) bytes, but irecvbuf is of size $(sizeof(irecvbuf)) bytes"))
    mod(length(recvbuf), comm_size) == 0 || throw(DimensionMismatch("The length of recvbuf and irecvbuf must be a multiple of the number of workers"))
    0 <= latency_tol || throw(ArgumentError("latency_tol must be non-negative, but is $latency_tol"))
    pool.responded .= false

    # send/receive buffers for each worker
    sl = sizeof(sendbuf)
    rl = div(sizeof(irecvbuf), comm_size)
    isendbufs = [view(reinterpret(UInt8, isendbuf), (i-1)*sl+1:i*sl) for i in 1:comm_size]
    irecvbufs = [view(reinterpret(UInt8, irecvbuf), (i-1)*rl+1:i*rl) for i in 1:comm_size]
    recvbufs = [view(reinterpret(UInt8, recvbuf), (i-1)*rl+1:i*rl) for i in 1:comm_size]

    # each call to asyncmap! is the start of a new epoch
    pool.epoch = epoch

    # start nonblocking transmission of sendbuf to all inactive workers
    # all workers are active after this loop
    for i in 1:length(pool.ranks)

        # do nothing with active workers
        if pool.active[i]
            continue
        end

        # worker is inactive
        pool.active[i] = true
        rank = pool.ranks[i]

        # copy send data to a worker-specific buffer to ensure the data isn't changed while sending
        isendbufs[i] .= view(reinterpret(UInt8, sendbuf), :)

        # remember at what epoch data was sent
        pool.sepochs[i] = pool.epoch

        # non-blocking send and receive
        pool.stimestamps[i] = time_ns()
        pool.sreqs[i] = MPI.Isend(isendbufs[i], rank, tag, comm)
        pool.rreqs[i] = MPI.Irecv!(irecvbufs[i], rank, tag, comm)
    end

    # track the number of active workers
    nactive = comm_size

    # record the start time to later compute the overall latency
    time0 = time_ns()

    # continually send and receive until,
    # 1) if nwait is an integer, we've received nwait results from this epoch, or
    # 2) if nwait is a function, when nwait(epoch, repochs) evaluates to true.
    nfresh = 0
    while 0 < nactive

        if typeof(nwait) <: Integer
            if nfresh >= nwait
                break
            end
        elseif typeof(nwait) <: Function
            if nwait(pool.epoch, pool.repochs)::Bool
                break
            end
        else
            error("nwait must be either an Integer or a Function, but is a $(typeof(nwait))")
        end

        # block until we've received a response
        indices, _ = MPI.Waitsome!(pool.rreqs)
        timestamp = time_ns()
        for i in indices

            # record latency
            pool.latency[i] = (timestamp - pool.stimestamps[i]) / 1e9

            # process previous response before overwriting it if we've already received from this worker
            if recvf isa Function && pool.responded[i]
                recvf(i, pool.epoch, pool.repochs[i], recvbufs[i])
            end

            # store the received data and set the epoch
            recvbufs[i] .= irecvbufs[i]
            pool.repochs[i] = pool.sepochs[i]

            # mark the worker as having responded, so that we later can process the reponse
            pool.responded[i] = true

            # cleanup the corresponding send request; should return immediately
            MPI.Wait!(pool.sreqs[i])

            # only receives initiated this epoch count towards completion
            if pool.repochs[i] == pool.epoch
                nfresh += 1
                nactive -= 1
                pool.active[i] = false
            else # otherwise send the worker a new task
                isendbufs[i] .= view(reinterpret(UInt8, sendbuf), :)
                pool.sepochs[i] = pool.epoch
                rank = pool.ranks[i]
                pool.stimestamps[i] = time_ns()
                pool.sreqs[i] = MPI.Isend(isendbufs[i], rank, tag, comm)
                pool.rreqs[i] = MPI.Irecv!(irecvbufs[i], rank, tag, comm)
            end
        end
    end

    # ns between starting to wait for workers until the exit condition evaluated to true
    latency = float(time_ns() - time0)
    extra_latency = latency

    # process any responses received within (latency * (1 + latency_tol)) since the exit condition evaluated to true
    while 0 < nactive && extra_latency <= latency * (1 + latency_tol)
        extra_latency = float(time_ns() - time0)
        indices, _ = MPI.Testsome!(pool.rreqs)
        timestamp = time_ns()
        for i in indices

            # record latency
            pool.latency[i] = (timestamp - pool.stimestamps[i]) / 1e9

            # process previous response before overwriting it if we've already received from this worker
            if recvf isa Function && pool.responded[i]
                recvf(i, pool.epoch, pool.repochs[i], recvbufs[i])
            end

            # store the received data and set the epoch
            recvbufs[i] .= irecvbufs[i]
            pool.repochs[i] = pool.sepochs[i]

            # mark the worker as having responded, so that we later can process the reponse
            pool.responded[i] = true

            # cleanup the corresponding send request; should return immediately
            MPI.Wait!(pool.sreqs[i])

            # mark as no longer active
            pool.active[i] = false
            nactive -= 1
        end
    end

    # process results
    if recvf isa Function
        for i in 1:length(pool.ranks)
            if pool.responded[i]
                recvf(i, pool.epoch, pool.repochs[i], recvbufs[i])
            end
        end
    end

    pool.repochs
end

"""
    waitall!(pool::MPIAsyncPool, recvbuf::AbstractArray, irecvbuf::AbstractArray)

Wait for all workers to respond. All workers are inactive when this function has returned.
"""
function waitall!(pool::MPIAsyncPool, recvbuf::AbstractArray, irecvbuf::AbstractArray; recvf::Union{Nothing,Function}=nothing)
    comm_size = length(pool.ranks)
    isbitstype(eltype(recvbuf)) || throw(ArgumentError("The eltype of sendbuf must be isbits, but is $(eltype(recvbuf))"))
    sizeof(recvbuf) == sizeof(irecvbuf) || throw(DimensionMismatch("recvbuf is of size $(sizeof(recvbuf)) bytes, but irecvbuf is of size $(sizeof(irecvbuf)) bytes"))
    mod(length(recvbuf), comm_size) == 0 || throw(DimensionMismatch("The length of recvbuf and irecvbuf must be a multiple of the number of workers"))
    pool.responded .= false

    nactive = sum(pool.active)
    if nactive == 0
        return pool.repochs
    end

    # receive buffers for each worker
    rl = div(sizeof(irecvbuf), comm_size)
    irecvbufs = [view(reinterpret(UInt8, irecvbuf), (i-1)*rl+1:i*rl) for i in 1:comm_size]
    recvbufs = [view(reinterpret(UInt8, recvbuf), (i-1)*rl+1:i*rl) for i in 1:comm_size]

    # receive from all active workers
    MPI.Waitall!(pool.rreqs)
    for i in 1:length(pool.ranks)
        if pool.active[i]
            pool.latency[i] = (time_ns() - pool.stimestamps[i]) / 1e9
            recvbufs[i] .= irecvbufs[i]
            pool.repochs[i] = pool.sepochs[i]
            pool.active[i] = false
            pool.responded[i] = true
            MPI.Wait!(pool.sreqs[i])
        end
    end

    # process results
    if recvf isa Function
        for i in 1:length(pool.ranks)
            if pool.responded[i]
                recvf(i, pool.epoch, pool.repochs[i], recvbufs[i])
            end
        end
    end

    pool.repochs
end

end

ext: jl


text:
using LogNormals
using Test, Distributions, StatsBase, Missings, MissingStrategies, Random
using OffsetArrays, RecursiveArrayTools, LinearAlgebra
using DistributionVectors # cormatri_for_acf
using Unitful

@testset "sem_cor" begin
    acf0 = [1,0.4,0.1]
    Sigma = cormatrix_for_acf(100, acf0);
    dmn = MvNormal(ones(100), Symmetric(Sigma));
    Random.seed!(1234)
    a = rand(dmn);
    am = allowmissing(a); am[3:4] .= missing;
    az = copy(a); az[3:4] .= 0.0;
    amu = am .* u"m"; # unitful missing vector of meters
    @testset "count_forlags" begin
        pred(x_i, x_iplusk) = ismissing(x_i) || ismissing(x_iplusk)
        x = [1,2,missing,missing,5,6]
        lags = 0:5
        nmiss = @inferred count_forlags(pred, x, lags)
        @test nmiss == [2,3,4,2,0,0]
        x = [1,2,missing,missing,5]
        lags = 0:4
        nmiss = @inferred count_forlags(pred, x, lags)
        @test nmiss == [2,3,3,1,0]
        nmissu = @inferred count_forlags(pred, x.*u"m", lags)
        @test nmissu == nmiss
        # test for custom indices
        nmiss2 = count_forlags(pred, OffsetArray(x,-1), lags)
        @test nmiss2 == nmiss
        # views still work and larger lags return 0
        nmiss3 = @inferred count_forlags(pred, view(x,1:2:5), lags)
        @test nmiss3 == [1,2,0,0,0]
    end;
    @testset "autocorr with missings" begin
        # for comparability center all around nomissing mean
        meana = mean(a)
        lags = 0:8
        acfe1 = @inferred autocor(a, lags) # StatsBase still works
        acfe2 = @inferred autocor(a, lags, PassMissing()) # can supply Missing clause
        @test acfe1 == acfe2
        # get only missing if not explicitly requesting SkipMissing or ExactMissing
        ismissing(@inferred(Vector{Float64},autocor(am, lags, PassMissing(); demean=false)))
        ismissing(@inferred(Vector{Float64},autocor(am, lags)))
        ismissing(@inferred(Vector{Float64},autocor(am, lags; demean=false)))
        acfe = @inferred autocor(a.-meana, lags; demean=false);
        acfemf = @inferred autocor(am.-meana, lags, SkipMissing(); demean=false);
        acfem = @inferred(autocor(am.-meana, lags, ExactMissing(); demean=false))
        acfes = @inferred autocor((a.-meana)[4:end], lags; demean=false) # different nobs
        azdemean = copy(a.-meana); azdemean[ismissing.(am)] .= 0.0
        acfez = @inferred autocor(azdemean, lags; demean=false);
        hcat(acfe, acfem, acfez, acfemf)[1:4,:]
        ## with SkipMissing() rather than ExactMissing(): missing same as vector with zeros
        @test acfemf == acfez
        @test acfem[1] == acfez[1] == 1
        # with excactmissing: correct var(x) by two missings (100-2)
        # lag1: three missings
        @test acfem[2] ≈ acfez[2]*98/97
        # lag2: four missings
        @test acfem[3] ≈ acfez[3]*98/96
        # lag3: three missings
        @test acfem[4] ≈ acfez[4]*98/97
        # lag4: two missings as in z
        @test acfem[5:end] ≈ acfez[5:end]
        # TODO after StatsBase.autocor can deal with unitful
        #acfemu = @inferred autocor((am.-meana).*u"m", lags; demean=false);
        #
        # variant without lags
        ismissing(@inferred autocor(am))
        ismissing(@inferred autocor(am, PassMissing()))
        @test acfemf ≈ (@inferred autocor(am.-meana, SkipMissing(); demean=false))[1:9]
        @test acfem ≈ (@inferred(autocor(am.-meana, ExactMissing(); demean=false)))[1:9]
        #SimpleTraits.istrait(IsEltypeSuperOfMissing{typeof(am.-meana)})
        #
        # with demean
        autocor(hcat(a,a), lags) # original variant
        tmp = @inferred autocor(am, lags, SkipMissing());
        tmp2 = @inferred autocor(hcat(am,am), lags, SkipMissing());
        @test tmp2[:,1] == tmp2[:,2] == tmp
        tmp = @inferred autocor(am, lags, ExactMissing());
        tmp2 = @inferred autocor(hcat(am,am), lags, ExactMissing());
        @test tmp2[:,1] == tmp2[:,2] == tmp
    end;
    @testset "autocor_effective" begin
        effa = @inferred autocor_effective(a, acf0)
        @test effa == acf0
        effa = @inferred autocor_effective(a, [acf0..., -0.001])
        @test effa == acf0
    end;
    @testset "effective_n_cor" begin
        neff = @inferred effective_n_cor(a, acf0)
        @test neff < length(a)
        @test neff ≈ 50.30 atol=0.01 # regression test
        # call with MissingStrategy but nonmissing type
        neffs = @inferred effective_n_cor(a, acf0, ExactMissing())
        @test neffs == neff
        #
        @test ismissing(@inferred(Missing,effective_n_cor(am, acf0, PassMissing())))
        @test ismissing(@inferred(Missing,effective_n_cor(am, acf0)))
        neffm = @inferred effective_n_cor(am, acf0, ExactMissing())
        @test neffm < neff
        @test neffm ≈ 49.61 atol=0.01 # regression test
        # n smaller than length(acf):
        neff2 = @inferred effective_n_cor(a[1:2], [acf0..., 0.05, 0.05])
        @test neff2 ≈ 1.43 atol=0.01 # regression test
        # no correlation estimate for lag 3
        neff3 = @inferred effective_n_cor(am[1:4], [acf0..., 0.05, 0.05], ExactMissing())
        @test neff3 ≈ 1.43 atol=0.01 # regression test
        #
        # without specifying acf
        ismissing(@inferred(Float64,effective_n_cor(am)))
        neff = @inferred effective_n_cor(am, ExactMissing())
        @test neff == effective_n_cor(am, autocor(am, ExactMissing()), ExactMissing())
    end;
    @testset "var_cor" begin
        va1 = @inferred var_cor(a, acf0, ExactMissing())
        va2 = @inferred var_cor(a, acf0, PassMissing())
        va = @inferred var_cor(a, acf0)
        # slight differences between var(x) and var(skipmissing(x))
        @test va ≈ va1 ≈ va2
        # need to specify handling of missing otherwise get missing back
        @test ismissing(@inferred(Float64,var_cor(am, acf0)))
        vam = @inferred var_cor(am, acf0, ExactMissing())
        #@code_warntype var_cor(am, acf0, ms=ExactMissing(); neff=nothing)
        @test isnan(var_cor(am[[3]],acf0, ExactMissing()))
        @test isnan(var_cor(am[2:3],acf0, ExactMissing())) # only one non-missing
        vamf = @inferred var_cor(am, acf0, SkipMissing())
        # larger uncertainty with missings
        @test vam > va
        # not correcting neff for missings overestimates neff
        # and hence underestimates uncertainty
        @test vamf < vam
        # calling without acf
        @test @inferred var_cor(a) > 0
        @test ismissing(@inferred(var_cor(am)))
        @test ismissing(@inferred(var_cor(am, PassMissing())))
        @test @inferred(var_cor(am, ExactMissing())) > 0
    end;
    @testset "semcor with same effective acf" begin
        se_a = @inferred sem_cor(a, acf0)
        ismissing(@inferred(Float64, sem_cor(am, acf0, PassMissing())))
        ismissing(@inferred(Float64, sem_cor(am, acf0)))
        se_am = @inferred sem_cor(am, acf0, ExactMissing())
        @test isnan(sem_cor(am[[3]],acf0, ExactMissing()))
        se_az = @inferred sem_cor(az, acf0)
        ar = filter(x -> !ismissing(x), am)
        se_ar = @inferred Missing sem_cor(ar, acf0)
        [se_a, se_am, se_az, se_ar]
        # the few missing do not lead to large deviations
        @test se_am ≈ se_a atol=0.02
        # skipping missing values has larger stderror, because shorter series
        @test se_ar > se_a
        # caring for missing pairs in correlation also has larger stderror
        # because some of larger relative error in sum_normals
        @test se_am > se_a
        # the SkipMissing() case may lead to lower estimates
        # of unceertainty than without missings?
    end;
    @testset "semcor calling without acf" begin
        @test @inferred sem_cor(a, ExactMissing()) ≈ sem_cor(a, autocor_effective(a))
        @test @inferred sem_cor(a, SkipMissing()) ≈ sem_cor(a, autocor_effective(a))
        @test @inferred sem_cor(a, PassMissing()) == sem_cor(a, autocor_effective(a))
        @test @inferred sem_cor(am, ExactMissing()) == sem_cor(am,
            autocor_effective(am, ExactMissing()), ExactMissing())
        @test @inferred sem_cor(am, SkipMissing()) == sem_cor(am,
            autocor_effective(am, SkipMissing()), SkipMissing())
    end;
    @testset "semcor with many random missings" begin
        b = allowmissing(a);
        Random.seed!(1234)
        b[sample(axes(a,1),60, replace=false)] .= missing;
        se_a = @inferred sem_cor(a, acf0)
        se_b = @inferred Missing sem_cor(b, acf0, ExactMissing())
        se_be = @inferred Missing sem_cor(b, ExactMissing())
        # uncertainty increases due to missings (fewer obs)
        @test se_b > se_a
        #@test se_be > se_a
        # regression tests with random seed 1234
        @test se_b ≈ 0.18 atol = 0.01
    end;
   @testset "semcor short series" begin
        @test isnan(@inferred(var_cor([1.0])))
        @test isnan(@inferred(sem_cor([1.0])))
        @test @inferred sem_cor(1.0:2) == 0.5
        @test @inferred sem_cor(1.0:3) == sqrt(var(1:3)/3)
        @test isnan(@inferred var_cor(fill(1.0,4))) # autocorrelation NaN
        #@test @inferred var_cor(fill(1.0,4)) == 0.0
        @test isnan(@inferred sem_cor(fill(1.0,4))) # autocorrelation NaN
        #@test @inferred sem_cor(fill(1.0,4)) == 0.0
        @test ismissing(@inferred(Float64, sem_cor([1.0,2,missing])))
        @test ismissing(@inferred(Float64, var_cor([1.0,2,missing])))
        @test @inferred var_cor([1.0,2,missing], ExactMissing()) == 0.5
        @test @inferred sem_cor([1.0,2,missing], ExactMissing()) == sqrt(0.5/2)
        @test isnan(@inferred sem_cor([1.0,missing,missing], ExactMissing()))
   end;
end;

@testset "autocorr with many missings" begin
    # ovrriding a from sd testset
    nsum = 1_000
    acf0 = [1,0.4,0.1]
    Sigma = cormatrix_for_acf(nsum, acf0);
    dmn = MvNormal(ones(nsum), Symmetric(Sigma));
    a = rand(dmn);
    probgap = 0.4
    clustersize = 4
    # need to average over some clusters to get reliable estimate
    nboot = 30
    resboot = map(1:nboot) do iboot
        igap0 = sample(0:(nsum-clustersize),trunc(Int,nsum*probgap/clustersize), replace=true);
        b = allowmissing(a);
        # outer product see https://stackoverflow.com/a/44592419
        #igap0.+Base.OneTo(clustersize)'
        b[igap0.+Base.OneTo(clustersize)'] .= missing;
        ae = autocor(b, 1:2, ExactMissing())
        af = autocor(b, 1:2, SkipMissing())
        (ae, af)
    end;
    aeb, afb = VectorOfArray.(vectuptotupvec(resboot));
    ae = mean(aeb; dims=2)
    af = mean(afb; dims=2)
    aorig = autocor(a, 1:2)
    @test ae ≈ aorig atol=0.02
    @test all((aorig .- af)./aorig .> 0.05) # >5% low bias
end;


text:
using JQ
using Test
using JSON

@testset "JQ" begin
    maths_test_dataset_file = joinpath(@__DIR__, "sample_math_test_results.json")
    maths_test_dataset_url = "https://data.cityofnewyork.us/api/views/7yig-nj52/rows.json?accessType=DOWNLOAD"
    # "https://data.cityofnewyork.us/api/views/7vy4-ats6/rows.json?accessType=DOWNLOAD"
    isfile(maths_test_dataset_file) || download(maths_test_dataset_url, maths_test_dataset_file)
    maths_test_dataset = JSON.parsefile(maths_test_dataset_file)

    attr_jqr = jqrunner(".meta.view.attributionLink")
    colname_jqr = jqrunner("[.meta.view.columns[].name]")
    yearselect_jqr = jqrunner("[.data[] | select(.[10]==\"2008\") | .[10]]")
    yearcount_jqr = jqrunner("length")

    @test jqstr(attr_jqr(maths_test_dataset)) == "\"http://schools.nyc.gov/NR/rdonlyres/3E439D63-A14E-4B6E-B019-B634A8915D71/0/DistrictMathResults20062012Public.xlsx\""

    column_names = jqparse(colname_jqr(maths_test_dataset))
    @test column_names == [
        "sid",
        "id",
        "position",
        "created_at",
        "created_meta",
        "updated_at",
        "updated_meta",
        "meta",
        "District",
        "Grade",
        "Year",
        "Demographic",
        "Number Tested",
        "Mean Scale Score",
        "Num Level 1",
        "Pct Level 1",
        "Num Level 2",
        "Pct Level 2",
        "Num Level 3",
        "Pct Level 3",
        "Num Level 4",
        "Pct Level 4",
        "Num Level 3 and 4",
        "Pct Level 3 and 4"
    ]

    @test Set(jqparse(yearselect_jqr(maths_test_dataset))) == Set(["2008"])

    @test parse(Int, jqstr(yearcount_jqr(column_names))) == length(column_names)
end


text:
@kwdef mutable struct SDecorrelator
    stepval::Int = 16
    rngrange::Int = 255
    buffer::CircularBuffer{Bool} = CircularBuffer{Bool}(2)
    counter::Int = 0
    rng::MersenneTwister = MersenneTwister(rand(UInt))

    function SDecorrelator(stepval, rngrange, buffer, counter, rng)
        decorr = new(stepval, rngrange, buffer, counter, rng)
        fill!(decorr.buffer, false)

        return decorr
    end
end
function (op::SDecorrelator)(x::Bool)
    # Increment counters
    op.counter += op.stepval * x

    # Decide output
    z = pop!(op.buffer)
    r = rand(op.rng, 0:op.rngrange)
    if r <= op.counter
        push!(op.buffer, true)
    else
        push!(op.buffer, false)
    end

    # Decrement counters
    op.counter = max(op.counter - op.stepval * z, 0)

    return z
end
function Base.show(io::IO, x::SDecorrelator)
    print(io, "SDecorrelator(")
    join(io, map((:stepval, :rngrange, :counter)) do field
        string(getproperty(x, field))
    end, ", ")
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SDecorrelator)
    print(io, "SDecorrelator(")
    join(io, map((:stepval, :rngrange, :counter)) do field
        string(field) * " = " * string(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

"""
    SSignedDecorrelator

A stochastic bitstream decorrelator.
"""
@kwdef struct SSignedDecorrelator
    pdecorr::SDecorrelator = SDecorrelator()
    ndecorr::SDecorrelator = SDecorrelator()
end
(op::SSignedDecorrelator)(x::SBit) = SBit((op.pdecorr(x.pos), op.ndecorr(x.neg)))
Base.show(io::IO, ::SSignedDecorrelator) = print(io, "SSignedDecorrelator(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedDecorrelator)
    print(io, "SSignedDecorrelator(")
    show(io, x.pdecorr)
    print(io, ", ")
    show(io, x.ndecorr)
    print(io, ")")
end

@kwdef mutable struct SAdder
    counter::Int = 0
end
function (op::SAdder)(x::Bool, y::Bool)
    # Increment counter
    op.counter += x + y

    # Decide output
    z = (op.counter >= 1)

    # Decrement counter
    op.counter = max(op.counter - z, 0)

    return z
end
function Base.show(io::IO, x::SAdder)
    print(io, "SAdder(")
    print(io, string(x.counter))
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SAdder)
    print(io, "SAdder(counter = ")
    print(io, string(x.counter))
    print(io, ")")
end

"""
    SSignedAdder

A signed stochastic bitstream add operator.
"""
@kwdef struct SSignedAdder
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
end
function (op::SSignedAdder)(x::SBit, y::SBit)
    pbit = op.padder(x.pos, y.pos)
    nbit = op.nadder(x.neg, y.neg)

    return SBit((pbit, nbit))
end
Base.show(io::IO, ::SSignedAdder) = print(io, "SSignedAdder(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedAdder)
    print(io, "SSignedAdder(")
    show(io, x.padder)
    print(io, ", ")
    show(io, x.nadder)
    print(io, ")")
end

@kwdef mutable struct SSaturatingSubtractor
    counter::Int = 0
end
function (op::SSaturatingSubtractor)(x::Bool, y::Bool)
    # # Increment counter
    # op.counter += (xor(x, y) && y) - (xor(x, y) && x)

    # # Saturate counters
    # op.counter = max(op.counter, 0)

    # # Decide output
    # z = (op.counter == 0 && x && !y)

    op.counter += x - y

    z = (op.counter >= 1)

    op.counter = max(op.counter - z, 0)

    return z
end
function Base.show(io::IO, x::SSaturatingSubtractor)
    print(io, "SSaturatingSubtractor(")
    print(io, string(x.counter))
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SSaturatingSubtractor)
    print(io, "SSaturatingSubtractor(counter = ")
    print(io, string(x.counter))
    print(io, ")")
end

"""
    SSignedSubtractor

A signed stochastic bitstream subtract operator.
"""
@kwdef struct SSignedSubtractor
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
    ppsub::SSaturatingSubtractor = SSaturatingSubtractor()
    pnsub::SSaturatingSubtractor = SSaturatingSubtractor()
    npsub::SSaturatingSubtractor = SSaturatingSubtractor()
    nnsub::SSaturatingSubtractor = SSaturatingSubtractor()
end
function (op::SSignedSubtractor)(x::SBit, y::SBit)
    pp = op.ppsub(x.pos, y.pos)
    pn = op.pnsub(y.pos, x.pos)
    np = op.npsub(y.neg, x.neg)
    nn = op.nnsub(x.neg, y.neg)
    z = (op.padder(pp, np), op.nadder(pn, nn))

    return SBit(z)
end
Base.show(io::IO, ::SSignedSubtractor) = print(io, "SSignedSubtractor(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedSubtractor)
    print(io, "SSignedSubtractor(")
    join(io, map(fieldnames(SSignedSubtractor)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

struct SMultiplier end
(op::SMultiplier)(x::Bool, y::Bool) = x & y

"""
    SSignedMultiplier

A signed stochastic bitstream multiply operator.
"""
@kwdef struct SSignedMultiplier
    ppmult::SMultiplier = SMultiplier()
    pnmult::SMultiplier = SMultiplier()
    npmult::SMultiplier = SMultiplier()
    nnmult::SMultiplier = SMultiplier()
    sub11::SSaturatingSubtractor = SSaturatingSubtractor()
    sub12::SSaturatingSubtractor = SSaturatingSubtractor()
    sub13::SSaturatingSubtractor = SSaturatingSubtractor()
    sub14::SSaturatingSubtractor = SSaturatingSubtractor()
    padder::SAdder = SAdder()
    nadder::SAdder = SAdder()
end
function (op::SSignedMultiplier)(x::SBit, y::SBit)
    pp = op.ppmult(x.pos, y.pos)
    pn = op.pnmult(x.pos, y.neg)
    np = op.npmult(x.neg, y.pos)
    nn = op.nnmult(x.neg, y.neg)

    s₁₁ = op.sub11(pp, pn)
    s₁₂ = op.sub12(pn, pp)
    s₁₃ = op.sub13(np, nn)
    s₁₄ = op.sub14(nn, np)

    z = (op.padder(s₁₁, s₁₄), op.nadder(s₁₂, s₁₃))

    return SBit(z)
end
Base.show(io::IO, ::SSignedMultiplier) = print(io, "SSignedMultiplier(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedMultiplier)
    print(io, "SSignedMultiplier(")
    join(io, map(fieldnames(SSignedMultiplier)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

@kwdef mutable struct SDivider
    counter::Int = 0
    zand::Bool = false
    rng::MersenneTwister = MersenneTwister(rand(UInt))
end
function (op::SDivider)(x::Bool, y::Bool)
    # Update counter
    op.counter = max(op.counter + 2 * x - 2 * op.zand, -100)

    # Decide output
    r = rand(op.rng, 0:64)
    z = (op.counter > r)

    # Update zand
    op.zand = z && y

    return z
end
function Base.show(io::IO, x::SDivider)
    print(io, "SDivider(")
    print(io, string(x.counter))
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SDivider)
    print(io, "SDivider(counter = ")
    print(io, string(x.counter))
    print(io, ")")
end

"""
    SSignedDivider

A signed stochastic bitstream divide operator.
"""
@kwdef struct SSignedDivider
    pdiv::SDivider = SDivider()
    ndiv::SDivider = SDivider()
    psub::SSaturatingSubtractor = SSaturatingSubtractor()
    nsub::SSaturatingSubtractor = SSaturatingSubtractor()
end
function (op::SSignedDivider)(x::SBit, y::SBit)
    pp = op.pdiv(x.pos, y.pos)
    np = op.ndiv(x.neg, y.pos)

    z = (op.psub(pp, np), op.nsub(np, pp))

    return SBit(z)
end
Base.show(io::IO, ::SSignedDivider) = print(io, "SSignedDivider(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedDivider)
    print(io, "SSignedDivider(")
    join(io, map(fieldnames(SSignedDivider)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

@kwdef mutable struct SFixedGainDivider
    counter::Int = 0
end
function (op::SFixedGainDivider)(x::Bool, y::Real)
    # Update counter
    op.counter += 255 * x

    # Decide output
    z = (op.counter >= round(255 * y))

    # Decrement counter
    op.counter -= z * round(255 * y)

    return z
end
function Base.show(io::IO, x::SFixedGainDivider)
    print(io, "SFixedGainDivider(")
    print(io, string(x.counter))
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SFixedGainDivider)
    print(io, "SFixedGainDivider(counter = ")
    print(io, string(x.counter))
    print(io, ")")
end

"""
    SSignedFixedGainDivider

A stochastic bitstream fixed gain divide operator.
"""
@kwdef struct SSignedFixedGainDivider
    pdiv::SFixedGainDivider = SFixedGainDivider()
    ndiv::SFixedGainDivider = SFixedGainDivider()
end
(op::SSignedFixedGainDivider)(x::SBit, y::Real) = SBit((op.pdiv(x.pos, y), op.ndiv(x.neg, y)))
Base.show(io::IO, ::SSignedFixedGainDivider) = print(io, "SSignedFixedGainDivider(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedFixedGainDivider)
    print(io, "SSignedFixedGainDivider(")
    join(io, map(fieldnames(SSignedFixedGainDivider)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

"""
    SSquareRoot

A stochastic bitstream square root operator.
"""
@kwdef mutable struct SSquareRoot
    counter::Int = 0
    zand::Bool = false
    buffer::CircularBuffer{Bool} = CircularBuffer{Bool}(1)
    rng::MersenneTwister = MersenneTwister(rand(UInt))

    function SSquareRoot(counter, zand, buffer, rng)
        op = new(counter, zand, buffer, rng)
        fill!(op.buffer, false)

        return op
    end
end
function (op::SSquareRoot)(x::SBit)
    # Update counter
    op.counter = max(op.counter + 4 * x.pos - 4 * op.zand, -100)

    # Decide output
    r = rand(op.rng, 0:511)
    z = (op.counter >= r)

    # Update zand
    op.zand = z && pop!(op.buffer)
    push!(op.buffer, z)

    return SBit((z, false))
end
function Base.show(io::IO, x::SSquareRoot)
    print(io, "SSquareRoot(")
    print(io, string(x.counter))
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SSquareRoot)
    print(io, "SSquareRoot(counter = ")
    print(io, string(x.counter))
    print(io, ")")
end

struct SMatMultiplier
    counters::Array{Int, 2}
end
SMatMultiplier(nrows, ncols) = SMatMultiplier(zeros(nrows, ncols))
function (op::SMatMultiplier)(x::VecOrMat{Bool}, y::VecOrMat{Bool})
    # Increment counters
    op.counters .+= x * y

    # Decide output
    z = (op.counters .>= 1)

    # Decrement counters
    op.counters .-= z

    return z
end
(op::SMatMultiplier)(x::BitArray, y::BitArray) =
    op(convert(Array{Bool, ndims(x)}, x), convert(Array{Bool, ndims(y)}, y))
function Base.show(io::IO, x::SMatMultiplier)
    print(io, "SMatMultiplier(")
    show(io, x.counter)
    print(io, ")")
end
function Base.show(io::IO, ::MIME"text/plain", x::SMatMultiplier)
    println(io, "SMatMultiplier(counter = ")
    show(io, MIME("text/plain"), x.counter)
    print(io, "\n)")
end

"""
    SSignedMatMultiplier

A stochastic bitstream matrix multiply operator.
"""
struct SSignedMatMultiplier
    ppmult::SMatMultiplier
    pnmult::SMatMultiplier
    npmult::SMatMultiplier
    nnmult::SMatMultiplier
    sub1::Array{SSaturatingSubtractor, 2}
    sub2::Array{SSaturatingSubtractor, 2}
    sub3::Array{SSaturatingSubtractor, 2}
    sub4::Array{SSaturatingSubtractor, 2}
    padder::Array{SAdder, 2}
    nadder::Array{SAdder, 2}
end
SSignedMatMultiplier(nrows, ncols) =
    SSignedMatMultiplier(
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        SMatMultiplier(nrows, ncols),
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SSaturatingSubtractor() for i in 1:nrows, j in 1:ncols],
        [SAdder() for i in 1:nrows, j in 1:ncols],
        [SAdder() for i in 1:nrows, j in 1:ncols]
    )
function (op::SSignedMatMultiplier)(x::VecOrMat{SBit}, y::VecOrMat{SBit})
    pp = op.ppmult(pos.(x), pos.(y))
    pn = op.pnmult(pos.(x), neg.(y))
    np = op.npmult(neg.(x), pos.(y))
    nn = op.nnmult(neg.(x), neg.(y))

    s₁ = map.(op.sub1, pp, pn)
    s₂ = map.(op.sub2, pn, pp)
    s₃ = map.(op.sub3, np, nn)
    s₄ = map.(op.sub4, nn, np)

    z = zip(map.(op.padder, s₁, s₄), map.(op.nadder, s₂, s₃)) |> collect

    return (size(z, 2) == 1) ? SBit.(dropdims(z; dims = 2)) : SBit.(z)
end
Base.show(io::IO, ::SSignedMatMultiplier) = print(io, "SSignedMatMultiplier(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedMatMultiplier)
    print(io, "SSignedMatMultiplier(")
    join(io, map(fieldnames(SSignedMatMultiplier)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

"""
    SL2Normer

A stochastic bitstream L2-norm operator.
"""
@kwdef struct SL2Normer
    dot::SSignedMatMultiplier = SSignedMatMultiplier(1, 1)
    buffer::CircularBuffer{Matrix{SBit}} = CircularBuffer{Matrix{SBit}}(1)
    root::SSquareRoot = SSquareRoot()
end
function (op::SL2Normer)(x::Vector{SBit})
    # get row vector
    dummybit = fill(SBit((false, false)), size(permutedims(x)))
    xt = isempty(op.buffer) ? dummybit : pop!(op.buffer)
    push!(op.buffer, permutedims(x))

    # compute inner product
    xtx = op.dot(xt, x)[1]

    # compute root
    z = op.root(xtx)

    return z
end
Base.show(io::IO, ::SL2Normer) = print(io, "SL2Normer(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SL2Normer)
    print(io, "SL2Normer(")
    join(io, map((:dot, :root)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

@kwdef mutable struct SSignedMaxer
    sub::SSignedSubtractor = SSignedSubtractor()
    counter::Int = 0
end
function (op::SSignedMaxer)(x::SBit, y::SBit)
    z = op.sub(x, y)
    op.counter += pos(z)
    op.counter -= neg(z)

    return (op.counter >= 0) ? x : y
end
Base.show(io::IO, ::SSignedMaxer) = print(io, "SSignedMaxer(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedMaxer)
    print(io, "SSignedMaxer(")
    join(io, map(fieldnames(SSignedMaxer)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end

struct SSignedNMaxer{N}
    maxers::Vector{SSignedMaxer}

    function SSignedNMaxer{N}(maxers::Vector{SSignedMaxer}) where N
        @assert length(maxers) == N - 1 "Incorrect number of `SSignedMaxer`s ($(length(maxers))) for `SSignedNMaxer{$N}`"

        new{N}(maxers)
    end
end
SSignedNMaxer(N) = SSignedNMaxer{N}([SSignedMaxer() for _ in 1:(N - 1)])
function (op::SSignedNMaxer{N})(xs::Vararg{SBit, N}) where N
    z = foldl(zip(op.maxers, Base.tail(xs)); init = first(xs)) do prev, (maxer, curr)
        maxer(prev, curr)
    end

    return z
end
Base.show(io::IO, ::SSignedNMaxer) = print(io, "SSignedNMaxer(...)")
function Base.show(io::IO, ::MIME"text/plain", x::SSignedNMaxer)
    print(io, "SSignedNMaxer(")
    join(io, map(fieldnames(SSignedNMaxer)) do field
        repr(getproperty(x, field))
    end, ", ")
    print(io, ")")
end


text:
"""
    struct TensorLagInterp{N,Td,T}

Generic Lagrange interpolation over an `N`-dimensional tensor grid. The
implementation uses a multidimensional generalization of the barycentric
formula.

The main constructor takes an `SVector{N,Vector{Td}}` containig the `N`
one-dimensional `nodes` and an `Array{N,T}` of the function `vals` at the tensor
product grid formed by the one-dimensional `nodes`.

# Examples:
```julia
nx = 10
ny = 12
x   = [0.5+0.5cos((2k-1)*π/2nx) for k in 1:nx] # Chebyshev nodes
y   = [0.5+0.5cos((2k-1)*π/2ny) for k in 1:ny] # Chebyshev nodes
f   = (x) -> cos(x[1]*x[2])
vals = [f((x,y)) for x in x, y in y]
p   = TensorLagInterp(SVector(x,y),vals)
p((0.1,0.2)) ≈ f((0.1,0.2))
```
"""
mutable struct TensorLagInterp{N,Td,T}
    vals::Array{T,N}
    nodes1d::NTuple{N,Vector{Td}}
    weights1d::NTuple{N,Vector{Td}}
end

interpolation_nodes(p::TensorLagInterp)                 = Iterators.product(p.nodes1d...)
interpolation_nodes(p::TensorLagInterp,d::Int,i::Int)   = p.nodes1d[d][i]
function interpolation_nodes(p::TensorLagInterp{N},I::CartesianIndex{N}) where {N}
    svector(N) do d
        i = I[d]
        p.nodes1d[d][i]
    end
end

weights1d(p::TensorLagInterp) = p.weights1d
vals(p::TensorLagInterp)      = p.vals

return_type(::TensorLagInterp{_,__,T}) where {_,__,T} = T

ambient_dimension(::TensorLagInterp{N}) where {N} = N

Base.eachindex(p::TensorLagInterp)  = eachindex(vals(p))
function Base.CartesianIndices(p::TensorLagInterp)
    sz = length.(p.nodes1d)
    CartesianIndices(sz)
end

# if passed the values and nodes, construct weights automatically
function TensorLagInterp(vals::Array{T,N},nodes::NTuple{N,Vector{Td}}) where {N,Td,T}
    weights = map(barycentric_lagrange_weights,nodes)
    TensorLagInterp{N,Td,T}(vals,nodes,weights)
end
TensorLagInterp(vals::Vector,nodes::Vector) = TensorLagInterp(vals,(nodes,))
TensorLagInterp(vals::Array,nodes::NTuple)  = TensorLagInterp(vals,SVector(nodes))

# barycentric lagrange interpolation
function (p::TensorLagInterp{N,Td,T})(x::SVector) where {N,Td,T}
    return bclag_interp(x, p.vals, p.nodes1d, p.weights1d, Val{N}(), 1, length(p.vals))
end
(p::TensorLagInterp{N,Td})(x) where {N,Td} =  p(SVector{N,Td}(x))
(p::TensorLagInterp{1})(x::Number) = p((x,))

# low-level function for multidimensional barycentric lagrange interpolation
# (adapted from FastChebInterp package)
@inline function bclag_interp(x::SVector{N,Td}, vals::Array{<:Any,N},
                            nodes::NTuple{N}, weights::NTuple{N},
                            ::Val{dim}, i1, len) where {N,Td,dim}
    T   = eltype(vals)
    n   = size(vals,dim)
    @inbounds xd = x[dim]
    @inbounds W  = weights[dim]
    @inbounds X  = nodes[dim]
    num = zero(T)
    den = zero(Td)
    if dim == 1
        for i in 1:n
            ci     = vals[i1+(i-1)]
            wi     = W[i]
            x_m_xi = xd - X[i]
            if x_m_xi == 0
                return ci
            end
            num += (wi/x_m_xi)*ci
            den += (wi/x_m_xi)
        end
        return num/den
    else
        Δi = len ÷ n # column-major stride of current dimension
        # recurse down on dimension
        dim′ = Val{dim-1}()
        for i in 1:n
            ci     = bclag_interp(x, vals, nodes, weights, dim′, i1+(i-1)*Δi, Δi)
            wi     = W[i]
            x_m_xi = xd - X[i]
            if x_m_xi == 0
                return ci
            end
            num += wi/x_m_xi*ci
            den += wi/x_m_xi
        end
        return num/den
    end
end

# algorithm on https://people.maths.ox.ac.uk/trefethen/barycentric.pdf
function barycentric_lagrange_weights(x::Vector)
    n = length(x) - 1
    w = similar(x)
    w[1] = 1.0
    for j in 1:n
        for k in 0:j-1
            w[k+1] = (x[k+1] - x[j+1]) * w[k+1]
        end
        w[j+1] = prod(0:j-1) do k
            x[j+1] - x[k+1]
        end
    end
    for j in 0:n
        w[j+1] = 1/w[j+1]
    end
    return w
end

# special nodes distributions and the corresponding weights

"""
    cheb1nodes(n,a,b)

Return the `n` Chebyshev points of the first kind on the interval `[a,b]`.
"""
function cheb1nodes(n::NTuple{N},lc,uc) where {N}
    iter = cheb1nodes_iter(n,lc,uc)
    return map(x->SVector{N,Float64}(x),iter)
end
function cheb1nodes(n::NTuple{N},rec::AbstractHyperRectangle{N}) where {N}
    cheb1nodes(n,low_corner(rec),high_corner(rec))
end
function cheb1nodes(n::Integer,a::Number,b::Number)
    xcheb::Vector{Float64} = cheb1nodes(n)
    c0 = (a+b)/2
    c1 = (b-a)/2
    return @. c0 + c1*xcheb
end
function cheb1nodes(n)
    return [-cos( (2*i-1)*π / (2*n) ) for i in 1:n]
end

function cheb1nodes_iter(n::NTuple{N},lc,uc) where {N}
    nodes1d = ntuple(d->cheb1nodes(n[d],lc[d],uc[d]),N)
    return iter = Iterators.product(nodes1d...)
end
function cheb1nodes_iter(n::NTuple{N},rec::AbstractHyperRectangle{N}) where {N}
    cheb1nodes_iter(n,low_corner(rec),high_corner(rec))
end

function cheb1weights(n::NTuple{N}) where {N}
    weights1d = ntuple(d->cheb1weights(n[d]),N)
    iter = Iterators.product(weights1d...)
    return map(x->prod(x),iter)
end
function cheb1weights(i,n)
    sgn = (-1)^(i)
    val = sin((2*i-1)*π/(2*n))
    return sgn*val
end
function cheb1weights(n)
    return [(-1)^(i)*sin((2*i-1)*π/(2*n)) for i in 1:n]
end

"""
    cheb2nodes(n,a,b)

Return the `n` Chebyshev points of the second kind on the interval `[a,b]`. The
nodes are nested in the following sense: `cheb2nodes(n,a,b) ==
cheb2nodes(2n-1,a,b)[1:2:end]`.
"""
function cheb2nodes(n::NTuple{N},lc,uc) where {N}
    iter = cheb2nodes_iter(n,lc,uc)
    return map(x->SVector{N,Float64}(x),iter)
end
function cheb2nodes(n::NTuple{N},rec::AbstractHyperRectangle{N}) where {N}
    cheb2nodes(n,low_corner(rec),high_corner(rec))
end
function cheb2nodes(n::Integer,a::Number,b::Number)
    xcheb::Vector{Float64} = cheb2nodes(n)
    c0 = (a+b)/2
    c1 = (b-a)/2
    return @. c0 + c1*xcheb
end
function cheb2nodes(n)
    return [cos((i-1)*π /(n-1)) for i in 1:n]
end

function cheb2nodes_iter(n::NTuple{N},lc,uc) where {N}
    nodes1d = ntuple(d->cheb2nodes(n[d],lc[d],uc[d]),N)
    iter = Iterators.product(nodes1d...)
    return iter
end
function cheb2nodes_iter(n::NTuple{N},rec::AbstractHyperRectangle{N}) where {N}
    return cheb2nodes_iter(n,low_corner(rec),high_corner(rec))
end

function cheb2weights(n::NTuple{N}) where {N}
    weights1d = ntuple(d->cheb2weights(n[d]),N)
    iter = Iterators.product(weights1d...)
    return map(x->prod(x),iter)
end

function cheb2weights(i,n)
    sgn = (-1)^(i)
    val = i==n || i == 1 ? 0.5 : 1.0
    return sgn*val
end

function cheb2weights(n)
    return [cheb2weights(i,n) for i in 1:n]
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
"""
	$(SIGNATURES)

Given a vector of mutable structs, compute a function of all numeric fields.
Immutable structs are skipped.

Example use case: Compute the mean and standard deviation of simulated model statistics.

`reduceFct` takes a vector of numeric inputs and returns a single numeric output that is of compatible size and type so it can be written back into an instance of the objects contained in `oVecV`.

To handle `Nan` or `missing`, make sure that `reduceFct` appropriately handles or skips them.

Optional input `fieldTypes` indicates which non-numeric field types to use as well.

This relies on `propertynames` to identify the fields to operate on. It uses `getproperty` and `setproperty!` to set those fields. These need to be appropriately overloaded for user defined types that do not store the properties "directly" (as accessed by `getfield`).

Limitation: Because the output object is of the same type as the input objects in `oVecV`, it is not possible to apply `reduceFct` that returns a different type from its inputs. For example, `mean` of integers does not work.

Change: exclude immutable fields, but keep scalars (which are immutable but can be changed!) +++
"""
function reduce_object_vector(oVecV :: Vector{T1},  reduceFct :: Function;
    fieldTypes :: Vector = []) where T1

    # Allocate output object
    oOut = deepcopy(oVecV[1]);
    # Loop over numeric fields
    pnV = propertynames(oOut);
    for pn in pnV
        #if !isimmutable(getproperty(oOut, pn))
            xM = reduce_one_field(oVecV, pn, reduceFct; fieldTypes = fieldTypes);
            if !isnothing(xM)
                setproperty!(oOut, pn, xM);
            end
        #end
    end
    return oOut
end

# Apply `reduceFct` to one field in a vector of objects.
# Returns `nothing` if reduction not possible.
function reduce_one_field(oVecV :: Vector{T1},  pn,  reduceFct :: Function;
    fieldTypes :: Vector = []) where T1

    n = length(oVecV);
    # Use getproperty instead of getfield (b/c getproperty can be overloaded)
    # @assert hasproperty(oVecV[1], pn)  "$(oVecV[1]) does not have property $pn"
    fieldM = getproperty(oVecV[1], pn);
    fType = typeof(fieldM);
    xM = nothing;

    fieldV = [getproperty(oVecV[iObj], pn)  for iObj = 1 : n];
    if isa(fieldV[1], AbstractArray)
        xM = reduce_array_vector(fieldV, reduceFct);
    else
        xM = reduce_scalar_vector(fieldV, reduceFct);
    end
    return xM

    # if (fType <: AbstractFloat)  ||  (fType ∈ fieldTypes)
    #     # Scalar field
    #     if hasmethod(reduceFct, (Vector{fType},))
    #         xM = reduceFct([getproperty(oVecV[iObj], pn) for iObj = 1 : n]);
    #     end

    # elseif isa(fieldM, AbstractArray)  &&  (eltype(fieldM) <: AbstractFloat)
    #     # Array field
    #     if hasmethod(reduceFct, (Vector{eltype(fieldM)},))
    #         xM = similar(fieldM);
    #         # Loop over elements in the field
    #         for j = 1 : length(xM)
    #             xM[j] = reduceFct([getproperty(oVecV[iObj], pn)[j] for iObj = 1 : n]);
    #         end
    #         @assert size(xM) == size(fieldM)  "Size mismatch: $(size(xM)) vs $(size(fieldM))"
    #     end
    # end

    # if !isnothing(xM)
    #     # Result can only be used if the resulting object matches `fieldM`
    #     @argcheck (typeof(xM) == fType)
    #     if hasmethod(size, (fType,))
    #         @argcheck (size(xM) == size(fieldM))
    #     end
    # end
    # return xM
end

"""
	$(SIGNATURES)

Reduce a `Vector{Array}` element by element. `reduceFct` is applied to the vector of `x[j]` for each index `j` where all the `x` from `fieldV` are stacked into a vector.

If `reduceFct` returns a different type, nothing is assigned. For example, the mean of `Int` values is not assigned.
"""
function reduce_array_vector(fieldV :: AbstractVector{Array{F1, N}},
    reduceFct) where {F1, N}

    n = length(fieldV);
    fType = typeof(fieldV[1]);

    # Array field
    if hasmethod(reduceFct, (Vector{eltype(fType)},))
        for iObj = 1 : n
            @argcheck size(fieldV[iObj]) == size(fieldV[1])
        end

        # Loop over elements in the field
        xM = similar(fieldV[1]);
        for j = eachindex(xM)
            # Reduce element by element
            valueV = [fieldV[iObj][j] for iObj = 1 : n];
            xr = reduceFct(valueV);
            if xr isa eltype(xM)
                xM[j] = xr;
            else
                xM = nothing;
                break;
            end
        end
    else
        xM = nothing;
    end
    if !isnothing(xM)
        @argcheck size(xM) == size(fieldV[1])
        @argcheck (typeof(xM) == fType)
    end
    return xM
end

"""
	$(SIGNATURES)

Reduce a vector of "scalars" so that `reduceFct(fieldV)` is called directly.

If `reduceFct` returns a different type, nothing is assigned.
"""
function reduce_scalar_vector(fieldV :: AbstractVector{F1}, reduceFct) where F1
    fType = typeof(fieldV[1]);
    if hasmethod(reduceFct, (Vector{fType},))
        xr = reduceFct(fieldV);
        if xr isa fType
            xM = xr;
        else
            xM = nothing;
        end
    else
        xM = nothing;
    end
    return xM
end

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


text:
function test_model_ad(model, logp_manual)
    vi = VarInfo(model)
    model(vi, SampleFromPrior())
    x = DynamicPPL.getall(vi)

    # Log probabilities using the model.
    function logp_model(x)
        new_vi = VarInfo(vi, SampleFromPrior(), x)
        model(new_vi, SampleFromPrior())
        return getlogp(new_vi)
    end

    # Check that both functions return the same values.
    lp = logp_manual(x)
    @test logp_model(x) ≈ lp

    # Gradients based on the manual implementation.
    grad = ForwardDiff.gradient(logp_manual, x)

    y, back = Tracker.forward(logp_manual, x)
    @test Tracker.data(y) ≈ lp
    @test Tracker.data(back(1)[1]) ≈ grad

    y, back = Zygote.pullback(logp_manual, x)
    @test y ≈ lp
    @test back(1)[1] ≈ grad

    # Gradients based on the model.
    @test ForwardDiff.gradient(logp_model, x) ≈ grad

    y, back = Tracker.forward(logp_model, x)
    @test Tracker.data(y) ≈ lp
    @test Tracker.data(back(1)[1]) ≈ grad

    y, back = Zygote.pullback(logp_model, x)
    @test y ≈ lp
    @test back(1)[1] ≈ grad
end


text:
mutable struct ClearDistances
    cardinal::SVector{4, Int}
    diagonal::SVector{4, Int}
end

function find_distances(f::Floor, obstacles::Set{Coord}, s::LTState)
    card = MVector{4, Int}(undef)
    diag = MVector{4, Int}(undef)
    for i in 1:4
        d = 1
        while !opaque(f, obstacles, s, s.robot+d*CARDINALS[i])
            d += 1
        end
        card[i] = d-1
    end
    for i in 1:4
        d = 1
        while !opaque(f, obstacles, s, s.robot+d*DIAGONALS[i])
            d += 1
        end
        diag[i] = d-1
    end
    return ClearDistances(card, diag)
end

mutable struct LTDistanceCache
    floor::Floor
    distances::Vector{ClearDistances}
end

function LTDistanceCache(f::Floor, obstacles::Set{Coord})
    dists = Array{ClearDistances}(undef,n_pos(f)^2)
    visited = falses(n_pos(f)^2)
    for i in 1:f.n_cols, j in 1:f.n_rows, k in 1:f.n_cols, l in 1:f.n_rows
        s = LTState(Coord(i,j), Coord(k,l), false)
        ii = stateindex(f, s)
        visited[ii] = true
        dists[ii] = find_distances(f, obstacles, s)
    end
    @assert all(visited)
    push!(dists, ClearDistances(zeros(4), zeros(4)))
    return LTDistanceCache(f, dists)
end

Base.getindex(c::LTDistanceCache, s::LTState) = c.distances[stateindex(c.floor, s)]

function n_clear_cells(d::ClearDistances, dir::Int)
    if dir <= 4
        return d.cardinal[dir]
    else
        return d.diagonal[dir-4]
    end
end

struct ReadingCDF
    cardcdf::Matrix{Float64}
    diagcdf::Matrix{Float64}
end

# reading CDF
function ReadingCDF(f::Floor,
                    std::Float64,
                    shortonly::Bool=false,
                    maxread::Int=ceil(Int, max_diag(f)+4*std))
    maxclear = max(f.n_rows, f.n_cols) - 1
    cardcdf = Array{Float64}(undef, maxclear + 1, maxread + 1)
    diagcdf = Array{Float64}(undef, maxclear + 1, maxread + 1)

    for c in 0:maxclear
        for r in 0:maxread
            cardcdf[c+1, r+1] = (1+erf((r+1.0-c)/(std*sqrt(2))))/2
        end
    end

    for c in 0:maxclear
        for r in 0:maxread
            diagcdf[c+1, r+1] = (1+erf((r+1.0-sqrt(2)*c)/(std*sqrt(2))))/2
        end
    end

    for c in 0:maxclear
        @assert abs(sum(diff(cardcdf[c, :])) - 1.0) < 1e-5
        @assert abs(sum(diff(diagcdf[c, :])) - 1.0) < 1e-5
    end

    return ReadingCDF(cardcdf, diagcdf)
end

function cdf(c::ReadingCDF, dir::Int, clear::Int, reading::Int)
    if dir <= 4 # cardinal
        if reading > size(c.cardcdf, 2)
            return 1.0
        else
            return c.cardcdf[clear + 1, reading + 1]
        end
    else # diagonal
        if reading > size(c.diagcdf, 2)
            return 1.0
        else
            return c.diagcdf[clear + 1, reading + 1]
        end
    end
end


text:
function gsd_bounds(fsim, h=[0.00025], ngrid=10)
    [get_bounds(fsim; type1=0.05, type2=0.2; ngrid=ngrid, h=hh) for hh in h]
end

function xyz_gsd_sim(iter, df, M; gen_data=ph_biweibull_1, model_fun=getBiPHWeibullModel)
    x = gsd_sim(iter, df, M=M, delta=(0.0, 0.0), gen_data=gen_data, model_fun=model_fun)
    y = gsd_sim(iter, df, M=M, delta=(0.0, log(0.6)), gen_data=gen_data, model_fun=model_fun)
    z = gsd_sim(iter, df, M=M, delta=(log(0.8), log(0.6)), gen_data=gen_data, model_fun=model_fun)
    x, y, z
end

function full_aed_sim(iter, bounds_df; h=[(0.00025, true), (0.000250001, false)], objs=[AED.public_utility, AED.sponsor_utility])
    bnds = [[r.a1 r.b1; r.a2 r.b2] for r in eachrow(bounds_df)]
    [xyz_aed_sim(iter, ab[1], h=h, objective=ab[2]) for ab in zip(bnds, objs)]
end

function xzy_aed_sim(iter, bounds; h=[(0.00025, true), (0.00025, false)], objective=AED.public_utility)
    x = AED.weibull_aed_sim(iter, bounds, delta=(0.0, 0.0), h=h, objective=objective)
    y = AED.weibull_aed_sim(iter, bounds, delta=(0.0, log(0.6)), h=h, objective=objective)
    z = AED.weibull_aed_sim(iter, bounds, delta=(log(0.8), log(0.6)), h=h, objective=objective)
    x,y,z
end

## forward simulations for Weibull proportional hazards model
weibull_fsim = AED.get_fsim(20000, AED.getBiPHWeibullModel, [[100, 50], [200, 100]]; max_n=400, M=10000,
                        delta0=[log(1.0), log(0.6)], delta=(log(0.8), log(0.6)), tau=0.0,
                        gen_data=AED.ph_weibull_1, warmup=5000, type1=0.05, type2=0.2)

## forward simulations with Spline proportional hazards model
spline_fsim = AED.get_fsim(20000, AED.getBiPHSplineModel, [[100, 50], [200, 100]]; max_n=400, M=10000,
                        delta0=[log(1.0), log(0.6)], delta=(log(0.8), log(0.6)), tau=0.0,
                        gen_data=AED.ph_billog_1, warmup=5000, type1=0.05, type2=0.2)

weibull_bounds = gsd_bounds(weibull_fsim)
spline_bounds = gsd_bounds(spline_fsim)

weibull_gsd_sim = xyz_gsd_sim(10000, weibull_bounds, 10000; gen_data=AED.ph_biweibull_1, model_fun=AED.getBiPHWeibullModel)
spline_gsd_sim = xyz_gsd_sim(10000, spline_bounds, 10000; gen_data=AED.ph_billog_1, model_fun=AED.getBiPHSplineModel)

weibull_aed_sim = full_aed_sim(1000, weibull_bounds)


text:
include("CalculusFunctional.jl")
include("Evaluation.jl")


text:
function get_GenericDER_states(Qref_Flag::Int)
    if Qref_Flag == 1
        return [:x1, :x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9], 9
    elseif Qref_Flag == 2 || Qref_Flag == 3
        return [:x2, :x3, :x4, :x5, :x6, :x7, :x8, :x9], 8
    else
        error("Unsupported value of Qref_Flag")
    end
end

function get_AggregateDistributedGenerationA_states(Freq_Flag::Int)
    if Freq_Flag == 0
        return [:Vmeas, :Pmeas, :Q_V, :Iq, :Mult, :Fmeas, :dPord, :Pord, :Ip], 9
    elseif Freq_Flag == 1
        return [:Vmeas, :Pmeas, :Q_V, :Iq, :Mult, :Fmeas, :Power_PI, :dPord, :Pord, :Ip], 10
    else
        error("Unsupported value of Freq_Flag")
    end
end

function calculate_IM_torque_params(A::Float64, B::Float64)
    C = 1.0 - A - B
    if A < 0.0 || B < 0.0 || C < 0.0
        error(
            "Unsupported values of A = $(A), B = $(B) or C = $(C). A, B and C must be positive and add up to 1.0",
        )
    end
    return C
end


text:
@parameters t
D = Differential(t)

function InPin(;name)
  @parameters x=13f0
  ODESystem(Equation[],t,Num[],[x]; name)
end

function OutPin(;name)
  @variables x(t) [output=true]
  ODESystem(Equation[],t,[x],Num[]; name)
end

function create_pins(in::Integer, out::Integer)
  inpins = [InPin(;name=Symbol("x$(i)_InPin")) for i in 1:in]
  outpins = [OutPin(;name=Symbol("x$(i)_OutPin")) for i in 1:out]
  inpins, outpins
end

# mutable struct SimType{T} <: DEDataVector{T}
#   x::Array{T,1}
#   observed::Array{T,1}
# end


text:
using Rotations
using GLMakie
using CubedSphere

"""
    visualize_conformal_mapping(w; x_min, y_min, x_max, y_max, n_lines, n_samples, filepath="mapping.png")

Visualize a complex mapping w(z) from the rectangle `x ∈ [x_min, x_max]`, `y ∈ [y_min, y_max]` using `n_lines`
equally spaced lines in x and y each evaluated at `n_samples` sample points.
"""
function visualize_conformal_mapping(w; x_min, y_min, x_max, y_max, n_lines, n_samples, filepath="mapping.png")
    z₁ = [[x + im * y for x in range(x_min, x_max, length=n_samples)] for y in range(y_min, y_max, length=n_lines)]
    z₂ = [[x + im * y for y in range(y_min, y_max, length=n_samples)] for x in range(x_min, x_max, length=n_lines)]
    zs = cat(z₁, z₂, dims=1)
    ws = [w.(z) for z in zs]

    fig = Figure(resolution=(3840, 2160))
    ax = fig[1, 1] = Axis(fig, xlabel="Re{w}", ylabel="Im{w}")
    [lines!(ax, real(z), imag(z), color=RGBAf0(0, 0, 0, 0.25)) for z in zs]
    [lines!(ax, real(w), imag(w), color="dodgerblue2") for w in ws]
    ax.aspect = AxisAspect(1)
    save(filepath, fig)
end

x = range(-1, 1, length=20)
y = range(-1, 1, length=20)

X = zeros(length(x), length(y))
Y = zeros(length(x), length(y))
Z = zeros(length(x), length(y))

for (i, x′) in enumerate(x), (j, y′) in enumerate(y)
    X[i, j], Y[i, j], Z[i, j] = conformal_cubed_sphere_mapping(x′, y′)
end

fig = wireframe(X, Y, Z, show_axis=false)

rotations = (RotX(π/2), RotX(-π/2), RotY(π/2), RotY(-π/2), RotX(π))

for R in rotations
    X′ = similar(X)
    Y′ = similar(Y)
    Z′ = similar(Z)

    for I in CartesianIndices(X)
        X′[I], Y′[I], Z′[I] = R * [X[I], Y[I], Z[I]]
    end

    wireframe!(X′, Y′, Z′)
end

display(fig)


text:
export fastmks

import ..FastMKSModel

using mlpack._Internal.io

import mlpack_jll
const fastmksLibrary = mlpack_jll.libmlpack_julia_fastmks

# Call the C binding of the mlpack fastmks binding.
function fastmks_mlpackMain()
  success = ccall((:fastmks, fastmksLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module fastmks_internal
  import ..fastmksLibrary

import ...FastMKSModel

# Get the value of a model pointer parameter of type FastMKSModel.
function IOGetParamFastMKSModel(paramName::String)::FastMKSModel
  FastMKSModel(ccall((:IO_GetParamFastMKSModelPtr, fastmksLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type FastMKSModel.
function IOSetParamFastMKSModel(paramName::String, model::FastMKSModel)
  ccall((:IO_SetParamFastMKSModelPtr, fastmksLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeFastMKSModel(stream::IO, model::FastMKSModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeFastMKSModelPtr, fastmksLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeFastMKSModel(stream::IO)::FastMKSModel
  buffer = read(stream)
  FastMKSModel(ccall((:DeserializeFastMKSModelPtr, fastmksLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    fastmks(; [bandwidth, base, degree, input_model, k, kernel, naive, offset, query, reference, scale, single, verbose])

This program will find the k maximum kernels of a set of points, using a query
set and a reference set (which can optionally be the same set). More
specifically, for each point in the query set, the k points in the reference set
with maximum kernel evaluations are found.  The kernel function used is
specified with the `kernel` parameter.

For example, the following command will calculate, for each point in the query
set `query`, the five points in the reference set `reference` with maximum
kernel evaluation using the linear kernel.  The kernel evaluations may be saved
with the  `kernels` output parameter and the indices may be saved with the
`indices` output parameter.

```julia
julia> using CSV
julia> reference = CSV.read("reference.csv")
julia> query = CSV.read("query.csv")
julia> indices, kernels, _ = fastmks(k=5, kernel="linear",
            query=query, reference=reference)
```

The output matrices are organized such that row i and column j in the indices
matrix corresponds to the index of the point in the reference set that has j'th
largest kernel evaluation with the point in the query set with index i.  Row i
and column j in the kernels matrix corresponds to the kernel evaluation between
those two points.

This program performs FastMKS using a cover tree.  The base used to build the
cover tree can be specified with the `base` parameter.

# Arguments

 - `bandwidth::Float64`: Bandwidth (for Gaussian, Epanechnikov, and
      triangular kernels).  Default value `1`.

 - `base::Float64`: Base to use during cover tree construction.  Default
      value `2`.

 - `degree::Float64`: Degree of polynomial kernel.  Default value `2`.

 - `input_model::FastMKSModel`: Input FastMKS model to use.
 - `k::Int`: Number of maximum kernels to find.  Default value `0`.

 - `kernel::String`: Kernel type to use: 'linear', 'polynomial', 'cosine',
      'gaussian', 'epanechnikov', 'triangular', 'hyptan'.  Default value
      `linear`.

 - `naive::Bool`: If true, O(n^2) naive mode is used for computation.
      Default value `false`.

 - `offset::Float64`: Offset of kernel (for polynomial and hyptan
      kernels).  Default value `0`.

 - `query::Array{Float64, 2}`: The query dataset.
 - `reference::Array{Float64, 2}`: The reference dataset.
 - `scale::Float64`: Scale of kernel (for hyptan kernel).  Default value
      `1`.

 - `single::Bool`: If true, single-tree search is used (as opposed to
      dual-tree search.  Default value `false`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `indices::Array{Int, 2}`: Output matrix of indices.
 - `kernels::Array{Float64, 2}`: Output matrix of kernels.
 - `output_model::FastMKSModel`: Output for FastMKS model.

"""
function fastmks(;
                 bandwidth::Union{Float64, Missing} = missing,
                 base::Union{Float64, Missing} = missing,
                 degree::Union{Float64, Missing} = missing,
                 input_model::Union{FastMKSModel, Missing} = missing,
                 k::Union{Int, Missing} = missing,
                 kernel::Union{String, Missing} = missing,
                 naive::Union{Bool, Missing} = missing,
                 offset::Union{Float64, Missing} = missing,
                 query = missing,
                 reference = missing,
                 scale::Union{Float64, Missing} = missing,
                 single::Union{Bool, Missing} = missing,
                 verbose::Union{Bool, Missing} = missing,
                 points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, fastmksLibrary), Nothing, ());

  IORestoreSettings("FastMKS (Fast Max-Kernel Search)")

  # Process each input argument before calling mlpackMain().
  if !ismissing(bandwidth)
    IOSetParam("bandwidth", convert(Float64, bandwidth))
  end
  if !ismissing(base)
    IOSetParam("base", convert(Float64, base))
  end
  if !ismissing(degree)
    IOSetParam("degree", convert(Float64, degree))
  end
  if !ismissing(input_model)
    fastmks_internal.IOSetParamFastMKSModel("input_model", convert(FastMKSModel, input_model))
  end
  if !ismissing(k)
    IOSetParam("k", convert(Int, k))
  end
  if !ismissing(kernel)
    IOSetParam("kernel", convert(String, kernel))
  end
  if !ismissing(naive)
    IOSetParam("naive", convert(Bool, naive))
  end
  if !ismissing(offset)
    IOSetParam("offset", convert(Float64, offset))
  end
  if !ismissing(query)
    IOSetParamMat("query", query, points_are_rows)
  end
  if !ismissing(reference)
    IOSetParamMat("reference", reference, points_are_rows)
  end
  if !ismissing(scale)
    IOSetParam("scale", convert(Float64, scale))
  end
  if !ismissing(single)
    IOSetParam("single", convert(Bool, single))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("indices")
  IOSetPassed("kernels")
  IOSetPassed("output_model")
  # Call the program.
  fastmks_mlpackMain()

  return IOGetParamUMat("indices", points_are_rows),
         IOGetParamMat("kernels", points_are_rows),
         fastmks_internal.IOGetParamFastMKSModel("output_model")
end


text:
#=
Implements an online version of Layer Normalization from https://arxiv.org/abs/1607.06450:
    Ba et al (2016): Layer Normalization

The input a is normalized ((a .- μ) ./ σ), and then re-scaled and shifted
with learnable params gain (g) and bias (b).
=#

#=
Input is a, output is y (same size as a):
    y = g ∘ z + b
    z = (a - μₜ) / σₜ
=#
type LayerNorm{T,P<:Params,W<:Weight} <: Learnable
    n::Int  # maps n --> n
    m::T  # this period's mean(aₜ)
    μ::T  # running average of m = mean(aₜ)
    σ::T  # running average of s = std(aₜ)
    σ̂::T  # the clamped version or σ
    αₜ::T
    wgt::W
    z::Vector{T}
    input::SumNode{T,1}    # a
    output::OutputNode{T,1}  # y
    params::P                  # holds g and b
end

function LayerNorm{T}(::Type{T}, n::Int,
                    θ::AbstractVector = zeros(T, 2n),
                    ∇::AbstractVector = zeros(T, 2n);
                    lookback::Int = 100,
                    α::Float64 = NaN,
                    wgt::Weight =  BoundedEqualWeight(isnan(α) ? lookback : α)
                    )
    input = InputNode(T, n)
    output = OutputNode(T, n)
    params = Params(θ, ∇, ((n,), (n,)))
    g, b = params.views
    fill!(g, one(T))
    fill!(b, zero(T))
    LayerNorm(n, zero(T), zero(T),
        one(T), one(T), zero(T),
        wgt, zeros(T,n),
        input, output, params
    )
end
LayerNorm(n::Int, args...; kw...) = LayerNorm(Float64, n, args...; kw...)

function Base.show(io::IO, t::LayerNorm)
    print(io, "LayerNorm{n=$(t.n), mu=$(t.μ), sigma=$(t.σ)}")
end

# compute output = wx + b
function transform!{T}(layer::LayerNorm{T})
    g, b = layer.params.views
    a = layer.input.val
    y = layer.output.val

    # update the αₜ weighting
    OnlineStats.updatecounter!(layer.wgt)
    layer.αₜ = weight(layer.wgt)

    # normalize the layer
    # TODO: μ/σ should be weighted avgs using BoundedEqualWeight
    layer.m = mean(a)
    layer.μ = layer.αₜ * layer.m + (one(T) - layer.αₜ) * layer.μ
    if isnan(layer.μ)
        warn("layer.μ is NaN:")
        @show layer a g b map(extrema, (a, g, b)) layer.wgt.nobs
    end
    layer.σ = layer.αₜ * std(a) + (one(T) - layer.αₜ) * layer.μ
    layer.σ̂ = if layer.σ == 0
        # @show layer.μ, layer.σ, layer.σ̂
        1e-10
    else
        layer.σ
    end

    # mult by g and add b
    @inbounds for o=1:layer.n
        layer.z[o] = (a[o] - layer.μ) / layer.σ̂
        y[o] = g[o] * layer.z[o] + b[o]
    end
    y
end

function grad!{T}(layer::LayerNorm{T})
    g, b = layer.params.views
    ∇g, ∇b = layer.params.∇_views
    a = layer.input.val
    y = layer.output.val
    ∇a = layer.input.∇
    ∇y = layer.output.∇

    D = layer.n
    Dinv = one(T) / D
    Dm1 = D - one(T)
    # @show layer.μ, layer.σ, layer.σ̂

    # TODO: ∇a, ∇g
    fill!(∇a, zero(T))
    for o=1:layer.n
        # ∇g
        ∇g[o] = ∇y[o] * layer.z[o]

        # ∇a
        base_dzₒdaₚ = Dinv - layer.αₜ * layer.z[o] * (layer.μ - layer.m) / layer.σ̂
        for p=1:layer.n
            dzₒdaₚ = base_dzₒdaₚ + layer.z[o] * layer.z[p] / Dm1
            dzₒdaₚ = (o==p ? one(T) : zero(T)) - layer.αₜ * dzₒdaₚ
            dzₒdaₚ /= layer.σ̂
            ∇a[p] += ∇y[o] * g[o] * dzₒdaₚ
        end
    end

    # ∇b
    copy!(∇b, ∇y)
    return
end


text:
"""

# coeftab_plot

Display plot with model coefficients

$(SIGNATURES)

# Arguments
- `dfs::DataFrame`: list of dataframes with sampled model parameters
- `pars`: list of column names (as symbols) to analyze. If not given, take all columns
- `pars_names`: list of strings to be used instead of column names
- `names`: list of dataframe names to be used in plot. If not given, "1", "2", "3"... will be used
- `perc_prob`: probability for confidence interval calculation
- `kwargs`: rest of arguments will be passed to the scatter function

# Examples
```julia
coeftab_plot(m5_1_df, m5_2_df, m5_3_df; pars=(:bA, :bM), names=["m5.1", "m5.2", "m5.3"])
```
"""
function coeftab_plot(dfs::DataFrame...; pars=missing, pars_names=missing,
    names=missing, perc_prob=0.89, kwargs...)::Plots.Plot
    if ismissing(pars)
        pars = unique(Iterators.flatten(propertynames.(dfs)))
        sort!(pars)
    end
    if ismissing(names)
        names = string.(Base.OneTo(length(dfs)))
    end
    if ismissing(pars_names)
        pars_names = string.(pars)
    end

    x = Vector{Float64}()
    xerr = Vector{Tuple{Float64,Float64}}()
    y = Vector{String}()
    for (p_name, p) ∈ zip(pars_names, pars)
        for (name, df) ∈ zip(names, dfs)
            p ∈ propertynames(df) || continue
            μ = mean(df[!,p])
            err = abs.(PI(df[!,p]; perc_prob) .- μ)
            pushfirst!(x, μ)
            label = length(dfs) == 1 ? "$p_name" : "$p_name: $name"
            pushfirst!(y, label)
            pushfirst!(xerr, Tuple(err))
        end
    end
    scatter(x, y; xerr=xerr, kwargs...)
end

export
    coeftab_plot


text:
# This is a very preliminary test function;
# Just small scale examples, e.g., P6 with 10 random signals. More coming!
using Test, MultiscaleGraphSignalTransforms, LinearAlgebra, SparseArrays, JLD2, Plots
@testset "MultiscalGraphSignalTransforms.jl" begin
##############################################################
# 1. Testing basic GHWT functions on 10 random signals on P6 #
##############################################################
    @testset "1. Testing basic GHWT functions" begin
        println("1. Testing basic GHWT functions")
        JLD2.@load "runtests_data/path6randn10.jld2" G
        G = gpath(6, G["tmp"])
        GP = partition_tree_fiedler(G)
        dc2f = ghwt_analysis!(G, GP=GP)

        @testset "Check the Haar basis" begin
            BH=bs_haar(GP)
            dvec=dmatrix2dvec(dc2f, GP, BH)
            frecon=ghwt_synthesis(dvec,GP,BH)
            println("Relative L2 error of the Haar transform: ", norm(G.f-frecon)/norm(G.f))
            @test norm(G.f-frecon)/norm(G.f) < 10 * eps()
        end

        @testset "Check the best basis and its expansion coefficients" begin
            bbc2f = ghwt_c2f_bestbasis(dc2f, GP)
            levlist = collect(enumerate([4; 4; 3; 3; 3; 3]))
            println("The true BB levlist: ", levlist)
            println("The comp BB levlist: ", bbc2f[2].levlist)
            println("They are equal: ", levlist == bbc2f[2].levlist)
            @test levlist == bbc2f[2].levlist
            JLD2.@load "runtests_data/bbcoef.jld2" tmp2
            println("The relative L2 error of the BB coefs: ", norm(tmp2["bbcoef"]-bbc2f[1])/norm(tmp2["bbcoef"]))
            @test norm(tmp2["bbcoef"]-bbc2f[1])/norm(tmp2["bbcoef"]) < 10 * eps()
        end
        println("\n")
    end

################################################################
# 2. Testing GHWT/eGHWT functions on a simple synthetic signal on P6 #
################################################################
    @testset "2. GHWT/eGHWT tests on P6" begin
        f = Array{Float64}([2. -2. 1. 3. -1. -2.]')
        G = gpath(6, f)
        GP = partition_tree_fiedler(G; method = :Lrw)
        dmatrix = ghwt_analysis!(G, GP=GP)
        println("2. Testing the GHWT functions on a path signal: ", f)
        println("The original signal has L1 norm: ", norm(f,1))

        @testset "The full level GHWT best basis" begin
            dvec,BS = ghwt_bestbasis(dmatrix, GP, cfspec=1.)
            (f, GS) = ghwt_synthesis(dvec, GP, BS, G)
            println("The coefficient vectors of the GHWT best basis has L1 norm: ", norm(dvec,1))
            println("Relative L2 error of the synthesized signal: ", norm(G.f-f)/norm(G.f))
            @test norm(G.f-f)/norm(G.f) < 10 * eps()
        end

        @testset "The GHWT best basis restricted levels" begin
            j_start = 3
            j_end = 4
            dvec_restricted,BS_restricted = ghwt_bestbasis(dmatrix, GP, cfspec=1., j_start = j_start, j_end = j_end)
            (f_restricted, GS_restricted) = ghwt_synthesis(dvec_restricted, GP, BS_restricted, G)
            println("The coefficient vectors of the GHWT best basis restricted from level $(j_start) to level $(j_end) in c2f dictionary has L1 norm: ", norm(dvec_restricted,1))
            println("Relative L2 error of the synthesized signal: ", norm(G.f-f_restricted)/norm(G.f))
            @test norm(G.f-f_restricted)/norm(G.f) < 10 * eps()
        end

        @testset "The eGHWT best basis" begin
            dvec_tf, BS_tf = ghwt_tf_bestbasis(dmatrix, GP, cfspec=1.)
            (f_tf, GS_tf) = ghwt_synthesis(dvec_tf, GP, BS_tf, G)
            println("The coefficient vectors of the time-frequency adapted GHWT best basis has L1 norm: ", norm(dvec_tf,1))
            println("Relative L2 error of the synthesized signal: ", norm(G.f-f_tf)/norm(G.f))
            @test norm(G.f-f_tf)/norm(G.f) < 10 * eps()
        end
        println("\n")
    end

#################################################################################
# 3. Testing time-frequency adapted 2D GHWT functions using a simple 3x3 matrix #
#################################################################################
    @testset "3. 2D GHWT/eGHWT on 3x3 matrix" begin
        println("3. Testing 2D GHWT/eGHWT functions")
        matrix = [ 1.0 2.0 3.0; 4.0 5.0 6.0]
        #matrix = matread("termdoc.mat")["matrix"]

        # expand matrix in 2 directions (rows and cols)
        GProws, GPcols = partition_tree_matrixDhillon(matrix)
        dmatrix = ghwt_analysis_2d(matrix, GProws, GPcols)

        # find the best basis using the time-frequency analysis
        # infovec indicate the location of each coefficient in dmatrix
        dvec_tf, infovec_tf = ghwt_tf_bestbasis_2d(matrix, GProws, GPcols)
        matrix_tf = ghwt_tf_synthesis_2d(dvec_tf, infovec_tf, GProws, GPcols)

        # find the best basis using the GHWT
        dvec_ghwt, BSrows, BScols = ghwt_bestbasis_2d(matrix, GProws, GPcols)
        matrix_ghwt = ghwt_synthesis_2d(dvec_ghwt, GProws, GPcols, BSrows, BScols)

        println("The toy matrix [1 2 3; 4 5 6] has 1-vecnorm as ", norm(matrix,1))
        println("The eGHWT best basis of toy matrix [1 2 3; 4 5 6] has 1-vecnorm as ", norm(dvec_tf,1))
        println("The GHWT best basis of toy matrix [1 2 3; 4 5 6] has 1-vecnorm as ", norm(dvec_ghwt,1))
        println("Relative L2 error of the eGHWT synthesized matrix: ", norm(matrix - matrix_tf, 2)/norm(matrix,2))
        @test norm(matrix - matrix_tf, 2)/norm(matrix,2) < 10 * eps()
        println("Relative L2 error of the GHWT synthesized matrix: ", norm(matrix - matrix_ghwt, 2)/norm(matrix,2))
        @test norm(matrix - matrix_ghwt, 2)/norm(matrix,2) < 10 * eps()
        println("\n")
    end

###################################################################
# 4. Testing HGLET functions and hybrid methods on a real dataset #
###################################################################

    @testset "4. HGLET on Toronto Street Network" begin
        println("4. HGLET on Toronto Street Network")
        #JLD2.@load "runtests_data/Dendrite.jld2" G
        #G = G["G"]
        # convert to SparseMatrixCSC{Float64,Int} where Int is machine dependent.
        #A = Array(G["W"])
        #B = sparse(A)
        #G = GraphSig(B, xy = G["xy"], f = G["f"], name = G["name"])

        JLD2.@load "runtests_data/Toronto_new.jld2"
        tmp1 = toronto["G"]
        G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],plotspecs = tmp1["plotspecs"])

        GP = partition_tree_fiedler(G)

        dmatrixH, dmatrixHrw, dmatrixHsym = HGLET_Analysis_All(G, GP) # expansion coefficients of 3-way HGLET bases
        dmatrixG = ghwt_analysis!(G,GP = GP) # expansion coefficients of GHWT

        dvec5,BS5,trans5,p5 = HGLET_GHWT_BestBasis_minrelerror(GP,G,dmatrixH = dmatrixH, dmatrixG = dmatrixG,
                                                                                dmatrixHrw = dmatrixHrw, dmatrixHsym = dmatrixHsym) # best-basis among all combinations of bases

        fS5, GS5 = HGLET_GHWT_Synthesis(reshape(dvec5,(size(dvec5)[1],1)),GP,BS5,trans5,G)

        println("The original signal has L1 norm: ", norm(G.f,1))
        println("The coefficients of best-basis selected from hybrid method has L1 norm: ", norm(dvec5,1))
        println("Relative L2 error of the synthesized signal: ", norm(G.f-fS5)/norm(G.f))
        @test norm(G.f-fS5)/norm(G.f) < 10 * eps()
        println("\n")
    end

###########################################################################
# 5. Testing GraphSig_Plot on mutilated Gaussian signal on the MN roadmap #
###########################################################################
    @testset "5. GraphSig_Plot on the MN roadmap" begin
        println("5. GraphSig_Plot on the MN roadmap")
        println("... this may not display the plot if you run it via the package test.")
        JLD2.@load "runtests_data/MN_MutGauss.jld2" G
        tmp1 = G["G"]
        G=GraphSig(tmp1["W"],xy=tmp1["xy"],f=tmp1["f"],name =tmp1["name"],plotspecs = tmp1["plotspecs"])
        G = Adj2InvEuc(G)
        plt = GraphSig_Plot(G)
        display(plt)
        @test typeof(plt) == Plots.Plot{Plots.GRBackend}
        println("\n")
    end

###################################################
# 6. Testing LP-HGLET functions on a real dataset #
###################################################

    @testset "6. LP-HGLET on RGC100" begin
        println("6. LP-HGLET on RGC100")
        JLD2.@load "runtests_data/Dendrite.jld2" G
        G = G["G_3D"]
        f = G["f"]
        G_Sig = GraphSig(1.0 * G["W"]; xy = G["xy"], f = f)
        G_Sig = Adj2InvEuc(G_Sig)
        GP = partition_tree_fiedler(G_Sig; swapRegion = false)

        dmatrixlH, _ = LPHGLET_Analysis_All(G_Sig, GP; ϵ = 0.3)
        dvec_lphglet, BS_lphglet, _ = HGLET_GHWT_BestBasis(GP, dmatrixH = dmatrixlH)

        fS, GS = LPHGLET_Synthesis(reshape(dvec_lphglet, (length(dvec_lphglet), 1)), GP, BS_lphglet, G_Sig; ϵ = 0.3)
        println("The original signal has L1 norm: ", norm(f, 1))
        println("The coefficients of LP-HGLET best-basis has L1 norm: ", norm(dvec_lphglet, 1))
        println("Relative L2 error of the synthesized signal: ", norm(f - fS) / norm(f))
        @test norm(f - fS) / norm(f) < 20 * eps()
        println("\n")
    end

#########################################################
# 7. Testing VM-NGWP, PC-NGWP, LP-NGWP functions on P64 #
#########################################################

    @testset "7. testing VM-NGWP/PC-NGWP/LP-NGWP on P64" begin
        println("7. testing VM-NGWP/PC-NGWP/LP-NGWP on P64")
        N = 64
        G = gpath(N)
        # use Chebyshev polynomial T₅(x) (x ∈ [0, 1]) as the path signal
        G.f = reshape([16 * x^5 - 20 * x^3 + 5 * x for x in LinRange(0, 1, N)], (N, 1))
        # compute graph Laplacian eigenvectors
        W = G.W
        L = diagm(sum(W; dims = 1)[:]) - W
        𝛌, 𝚽 = eigen(L); 𝚽 = 𝚽 .* sign.(𝚽[1,:])'
        # build Gstar
        Gstar = GraphSig(W)
        GP = partition_tree_fiedler(G; swapRegion = false)
        GP_dual = partition_tree_fiedler(Gstar; swapRegion = false)
        GP_primal = pairclustering(𝚽, GP_dual) # for PC-NGWP
        # construct NGWP dictionaries
        VM_NGWP = vm_ngwp(𝚽, GP_dual)
        PC_NGWP = pc_ngwp(𝚽, GP_dual, GP_primal)
        LP_NGWP = lp_ngwp(𝚽, W, GP_dual; ϵ = 0.3)
        # NGWP analysis
        dmatrix_VM = ngwp_analysis(G, VM_NGWP)
        dvec_vm_ngwp, BS_vm_ngwp = ngwp_bestbasis(dmatrix_VM, GP_dual)
        dmatrix_PC = ngwp_analysis(G, PC_NGWP)
        dvec_pc_ngwp, BS_pc_ngwp = ngwp_bestbasis(dmatrix_PC, GP_dual)
        dmatrix_LP = ngwp_analysis(G, LP_NGWP)
        dvec_lp_ngwp, BS_lp_ngwp = ngwp_bestbasis(dmatrix_LP, GP_dual)

        println("The original signal has L2 norm: ", norm(G.f))
        println("The coefficients of VM-NGWP best-basis has L2 norm: ", norm(dvec_vm_ngwp))
        @test abs(norm(G.f) - norm(dvec_vm_ngwp)) / norm(G.f) < 10 * eps()
        println("The coefficients of PC-NGWP best-basis has L2 norm: ", norm(dvec_pc_ngwp))
        @test abs(norm(G.f) - norm(dvec_pc_ngwp)) / norm(G.f) < 10 * eps()
        println("The coefficients of LP-NGWP best-basis has L2 norm: ", norm(dvec_lp_ngwp))
        @test abs(norm(G.f) - norm(dvec_lp_ngwp)) / norm(G.f) < 10 * eps()
        println("\n")
    end

end
# End of runtests.jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Desription
# ==============================================================================
#
#   Tests related to conversion from quaternions to Euler angles.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# File: ./src/conversions/quaternion_to_angle.jl
# ==============================================

# Functions: quat_to_angle
# ------------------------

@testset "Quaternion => Euler angles (Float64)" begin
    T = Float64
    # We do not need comprehensive test here because `quat_to_angle` first
    # converts a quaternion to DCM and then to Euler angles. Those two
    # operations are already heavily tested.

    q = Quaternion(cosd(45 / 2), sind(45 / 2), 0, 0)
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 45 * pi / 180
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(45 / 2), 0, sind(45 / 2), 0)
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 45 * pi / 180
    @test ea.a3 ≈ 0
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(45 / 2), 0, 0, sind(45 / 2))
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 45 * pi / 180
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 0
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(45 / 2), 0, 0, sind(45 / 2))
    ea = quat_to_angle(q, :YXZ)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 45 * pi / 180
    @test ea.rot_seq === :YXZ
end

@testset "Quaternion => Euler angles (Float64)" begin
    T = Float32
    # We do not need comprehensive test here because `quat_to_angle` first
    # converts a quaternion to DCM and then to Euler angles. Those two
    # operations are already heavily tested.

    q = Quaternion(cosd(T(45 / 2)), sind(T(45 / 2)), 0, 0)
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 45 * pi / 180
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(T(45 / 2)), 0, sind(T(45 / 2)), 0)
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 45 * pi / 180
    @test ea.a3 ≈ 0
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(T(45 / 2)), 0, 0, sind(T(45 / 2)))
    ea = quat_to_angle(q, :ZYX)
    @test eltype(ea) === T
    @test ea.a1 ≈ 45 * pi / 180
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 0
    @test ea.rot_seq === :ZYX

    q = Quaternion(cosd(T(45 / 2)), 0, 0, sind(T(45 / 2)))
    ea = quat_to_angle(q, :YXZ)
    @test eltype(ea) === T
    @test ea.a1 ≈ 0
    @test ea.a2 ≈ 0
    @test ea.a3 ≈ 45 * pi / 180
    @test ea.rot_seq === :YXZ
end


text:
@doc raw"""
    createMosaic(indices, colr, constrnt)

Renders domino mosaic from optimization output

# Resources
Watch Bob Bosch of Oberlin College talk at Google about domino mosaics, and his interpretation of this art form as an integer program.
https://www.youtube.com/watch?v=g3CiVrN-BnY

IP modeling details can be found in this article.
http://www.optimization-online.org/DB_FILE/2003/09/722.pdf

# Examples
```julia-repl
julia> include("ModuleDominoes.jl")
julia> using .Dominoes
julia> createMosaic("Frankenstein", "black", "constrained")
```

# Author
Kalyani Nagaraj
March 2021
"""
function createMosaic(imgname, colr, constraint)

   # CALL OPTIMIZER:
   # ---------------
   dominoes(imgname, colr, constraint)
   imgname * "_" * colr * "_" * constraint

   # 1. LOAD 55 DOMINO IMAGES FROM FILE:
   # ------------------------------------
   file = matopen("../../Data/Domino_Set_Images/dominoes_with_" * colr * "_dots_V7.mat") # -V7 is compatible with Julia when used with package MAT
   f = read(file, "c") # m is read as a array of 55 arrays, each of type UInt8
   d = Int64.(ones(332,162,55)) #Make each matrix of size (332, 162)
   for i = 1:55
      d[:,:,i] = Int64.(f[i][12:end,179:340])
   end
   r = size(d[:,:,1],1)
   c = size(d[:,:,1],2)
   offset = r-2*c

   # 2. READ SOLUTION INDICES FROM FILE:
   # -----------------------------------
   layout = readdlm("../../Data/Optimization_Output/Indices/" * imgname * "_" * colr * "_" * constraint * ".txt", ' ', Int, skipstart = 1)
   R = maximum(layout[:, 4])
   C = maximum(layout[:, 5])

   # 3. FOR EACH ENTRY IN LAYOUT, WRITE DOMINO IMAGE TO MATRIX I:
   # ------------------------------------------------------------
   if colr == "white"
      I = Int64.(255*ones(R*(offset + c), C*(offset + c)))
   elseif colr == "black"
      I = Int64.(zeros(R*(offset + c), C*(offset + c)))
   end

   for domino = 1:size(layout,1)

      m = layout[domino, 1]
      n = layout[domino, 2]
      o = layout[domino, 3]
      i = layout[domino, 4]
      j = layout[domino, 5]

      # 4. Get domino cell number:
      # --------------------------
      k = 10*m + n - sum(1:m) + 1

      if o == 1
         # place domino as is at location (i,j)
         dd = d[:,:,k];

      elseif o == 2
         # swap m and n
         dd = [d[Int64(r/2 + 1):end, :,k]; d[1:Int64(r/2), :, k]]

      elseif o == 3
         # rotate counterclockwise by 90 degrees
         # dd = [d{k}'(end:-1:c/2+1, :); d{k}'(c/2:-1:1, :)];
         dt = d[:,:,k]';  # transpose(d(k)) to get size(dt) = (c, r)
         dd = [dt[end:-1:Int64(c/2+1), :]; dt[Int64(c/2):-1:1, :]] # swap top and bottom half as mirror images

      elseif o == 4
         # first swap m and n, then rotate by 90 degrees
         # dd = [d{k}(r/2 + 1:end, :); d{k}(1:r/2, :)];
         # dd = [dd'(end:-1:c/2+1, :); dd'(c/2:-1:1, :)];
         dt = [d[Int64(r/2 + 1):end, :,k]; d[1:Int64(r/2), :, k]]'    # swap m and n, and perform a matrix transpose
         dd = [dt[end:-1:Int64(c/2+1), :]; dt[Int64(c/2):-1:1, :]]     # swap top and bottom half as mirror images

      end

      I[(i-1)*(c+offset)+1 : (i-1)*(c+offset)+size(dd,1) , (j-1)*(c+offset)+1 : (j-1)*(c+offset)+size(dd,2)] = dd;

   end

   # 5. Save image to file:
   # ----------------------
   I = colorview(Gray, I/256)
   save("../../Data/Mosaic_Images/" * imgname * "_" * colr * "_" * constraint * ".png", I)

end

max_stars_repo_name: kalyaninagaraj/Domino-Mosaics

max_issues_repo_name: kalyaninagaraj/Domino-Mosaics

max_forks_repo_name: kalyaninagaraj/Domino-Mosaics


text:
## Discretized Asset types ond methods ##
type DiscretizedAssetCommon{L <: Lattice}
  time::Float64
  values::Vector{Float64}
  latestPreAdjustment::Float64
  latestPostAdjustment::Float64
  method::L

  # call(::Type{DiscretizedAssetCommon}, t::Float64, v::Vector{Float64}, lpa::Float64, lpoa::Float64) =
  #     new{Lattice}(t, v, lpa, lpoa)
end

DiscretizedAssetCommon() = DiscretizedAssetCommon{NullLattice}(0.0, zeros(0), eps(), eps(), NullLattice())

clone(das::DiscretizedAssetCommon, method::Lattice = das.method) = DiscretizedAssetCommon{typeof(method)}(das.time, das.values, das.latestPreAdjustment, das.latestPostAdjustment, method)

set_time!(a::DiscretizedAsset, t::Float64) = a.common.time = t

get_values(a::DiscretizedAsset) = a.common.values

function set_method!(a::DiscretizedAsset, method::Lattice)
  new_common = clone(a.common, method)
  a.common = new_common

  return a
end

function adjust_values!(dAsset::DiscretizedAsset)
  pre_adjust_values!(dAsset)
  post_adjust_values!(dAsset)

  return dAsset
end

function pre_adjust_values!(dAsset::DiscretizedAsset)
  if ~QuantLib.Math.close_enough(dAsset.common.time, dAsset.common.latestPreAdjustment)
    pre_adjust_values_impl!(dAsset)
    dAsset.common.latestPreAdjustment = dAsset.common.time
  end

  return dAsset
end

function post_adjust_values!(dAsset::DiscretizedAsset)
  if ~QuantLib.Math.close_enough(dAsset.common.time, dAsset.common.latestPostAdjustment)
    post_adjust_values_impl!(dAsset)
    dAsset.common.latestPostAdjustment = dAsset.common.time
  end

  return dAsset
end

function partial_rollback!(dAsset::DiscretizedAsset, t::Float64)
  partial_rollback!(dAsset.common.method, dAsset, t)
  return dAsset
end

function is_on_time(dAsset::DiscretizedAsset, t::Float64)
  grid = dAsset.common.method.tg
  return QuantLib.Math.close_enough(grid.times[return_index(grid, t)], dAsset.common.time)
end

function rollback!(dAsset::DiscretizedAsset, t::Float64)
  rollback!(dAsset.common.method, dAsset, t)

  return dAsset
end

pre_adjust_values_impl!(dAsset::DiscretizedAsset) = dAsset # do nothing
post_adjust_values_impl!(dAsset::DiscretizedAsset) = dAsset # do nothing

function general_reset!(dOption::DiscretizedOption, sz::Int)
  # check methods in main and underlying
  dOption.common.values = zeros(sz)
  adjust_values!(dOption)
  return dOption
end

present_value(dAsset::DiscretizedAsset) = present_value(dAsset.common.method, dAsset)

function post_adjust_values_impl!(dOption::DiscretizedOption)
  partial_rollback!(dOption.underlying, dOption.common.time)
  pre_adjust_values!(dOption.underlying)

  if isa(dOption.exercise, AmericanExercise)
    if dOption.common.time >= dOption.exerciseTimes[1] && dOption.common.time <= dOption.exerciseTimes[1]
      apply_exercise_condition!(dOption)
    end
  else
    @simd for i = 1:length(dOption.exerciseTimes)
      @inbounds t = dOption.exerciseTimes[i]
      if t >= 0.0 && is_on_time(dOption, t)
        apply_exercise_condition!(dOption)
      end
    end
  end

  post_adjust_values!(dOption.underlying)

  return dOption
end

function apply_exercise_condition!(dOption::DiscretizedOption)
  @simd for i = 1:length(dOption.common.values)
    @inbounds dOption.common.values[i] = max(dOption.underlying.common.values[i], dOption.common.values[i])
  end

  return dOption
end

function initialize!(asset::DiscretizedAsset, lattice::Lattice, t::Float64)
  set_method!(asset, lattice)
  initialize!(lattice, asset, t)
  return asset
end