text:
# This file is inspired from Polyhedra.jl

const ABSZTOL(N::Type{<:AbstractFloat}) = N(10) * sqrt(eps(N))
const ABSZTOL(N::Type{Rational{INNER}}) where {INNER} = zero(N)

"""
    _leq(x::N, y::N; [kwargs...]) where {N<:Real}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`      -- number
- `y`      -- another number (of the same numeric type as `x`)
- `kwargs` -- not used

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

This is a fallback implementation for numbers of type `Real`. If the arguments
are floating point numbers, see `_leq(x::AbstractFloat, y::AbstractFloat)`.
"""
_leq(x::N, y::N; kwargs...) where {N<:Real} = x <= y

"""
    _leq(x::N, y::M; [kwargs...]) where {N<:Real, M<:Real}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`      -- number
- `y`      -- another number (of possibly different numeric type than `x`)
- `kwargs` -- optional arguments; see `?_leq` for the available options

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

This implementation calls Julia's `promote(x, y)` function, which converts all
arguments to a common numeric type, returning them as a tuple. The conversion
is such that the common type to which the values are converted can represent
them as faithfully as possible.
"""
_leq(x::N, y::M; kwargs...) where {N<:Real, M<:Real} =
    _leq(promote(x, y)...; kwargs...)

"""
    _geq(x::Real, y::Real; [kwargs...])

Determine if `x` is greater than or equal to `y`.

### Input

- `x` -- number
- `y` -- another number (of possibly different numeric type than `x`)

### Output

A boolean that is `true` iff `x >= y`.

### Algorithm

This function falls back to `_leq(y, x)`, with type promotion if needed. See the
documentation of `_leq` for further details.
"""
_geq(x::Real, y::Real; kwargs...) = _leq(y, x; kwargs...)

"""
    isapproxzero(x::N; ztol::Real=ABSZTOL(N)) where {N<:Real}

Determine if `x` is approximately zero.

### Input

- `x`    -- number
- `ztol` -- (optional, default: `ABSZTOL`) tolerance against zero

### Output

A boolean that is `true` iff `x ≈ 0`.

### Algorithm

It is considered that `x ≈ 0` whenever `x` (in absolute value) is smaller than
the tolerance for zero, `ztol`.
"""
function isapproxzero(x::N; ztol::Real=ABSZTOL(N)) where {N<:Real}
    return abs(x) <= ztol
end

"""
    _isapprox(x::N, y::N;
              rtol::Real=Base.rtoldefault(N),
              ztol::Real=ABSZTOL(N),
              atol::Real=zero(N)) where {N<:Real}

Determine if `x` is approximately equal to `y`.

### Input

- `x`    -- number
- `y`    -- another number (of the same numeric type as `x`)
- `rtol` -- (optional, default: `Base.rtoldefault(N)`) relative tolerance
- `ztol` -- (optional, default: `ABSZTOL(N)`) absolute tolerance for comparison
            against zero
- `atol` -- (optional, default: `zero(N)`) absolute tolerance

### Output

A boolean that is `true` iff `x ≈ y`.

### Algorithm

We first check if `x` and `y` are both approximately zero, using
`isapproxzero(x, y)`.
If that fails, we check if `x ≈ y`, using Julia's `isapprox(x, y)`.
In the latter check we use `atol` absolute tolerance and `rtol` relative
tolerance.

Comparing to zero with default tolerances is a special case in Julia's
`isapprox`, see the last paragraph in `?isapprox`. This function tries to
combine `isapprox` with its default values and a branch for `x ≈ y ≈ 0` which
includes `x == y == 0` but also admits a tolerance `ztol`.

Note that if `x = ztol` and `y = -ztol`, then `|x-y| = 2*ztol` and still
`_isapprox` returns `true`.
"""
function _isapprox(x::N, y::N;
                   rtol::Real=Base.rtoldefault(N),
                   ztol::Real=ABSZTOL(N),
                   atol::Real=zero(N)) where {N<:Real}
    if isapproxzero(x, ztol=ztol) && isapproxzero(y, ztol=ztol)
        return true
    else
        return isapprox(x, y, rtol=rtol, atol=atol)
    end
end

# generic "dense"
function _isapprox(x::AbstractVector{N}, y::AbstractVector{N};
                   rtol::Real=Base.rtoldefault(N),
                   ztol::Real=ABSZTOL(N),
                   atol::Real=zero(N)) where {N<:Real}
    n = length(x)
    if length(x) != length(y)
        return false
    end
    @inbounds for i in 1:n
        if !_isapprox(x[i], y[i], rtol=rtol, ztol=ztol, atol=atol)
            return false
        end
    end
    return true
end

# sparse
function _isapprox(x::SparseVector{N}, y::SparseVector{N};
                   rtol::Real=Base.rtoldefault(N),
                   ztol::Real=ABSZTOL(N),
                   atol::Real=zero(N)) where {N<:Real}
    @assert length(x) == length(y)
    return x.nzind == y.nzind && _isapprox(x.nzval, y.nzval, rtol=rtol, ztol=ztol, atol=atol)
end

"""
    ispermutation(u::AbstractVector{T}, v::AbstractVector)::Bool where {T}

Check that two vectors contain the same elements up to reordering.

### Input

- `u` -- first vector
- `v` -- second vector

### Output

`true` iff the vectors are identical up to reordering.

### Examples

```jldoctest
julia> LazySets.ispermutation([1, 2, 2], [2, 2, 1])
true

julia> LazySets.ispermutation([1, 2, 2], [1, 1, 2])
false
```

### Notes

Containment check is performed using `LazySets._in(e, v)`, so in the case of
floating point numbers, the precision to which the check is made is determined
by the type of elements in `v`. See `_in` and `_isapprox` for more information.
"""
function ispermutation(u::AbstractVector{T}, v::AbstractVector)::Bool where {T}
    if length(u) != length(v)
        return false
    end
    occurrence_map = Dict{T, Int}()
    has_duplicates = false
    for e in u
        if !_in(e, v)
            return false
        end
        if haskey(occurrence_map, e)
            occurrence_map[e] += 1
            has_duplicates = true
        else
            occurrence_map[e] = 1
        end
    end
    if has_duplicates
        for e in v
            if !haskey(occurrence_map, e) || occurrence_map[e] == 0
                return false
            end
            occurrence_map[e] -= 1
        end
    end
    return true
end

# alias for Julia's containment check
function _in(x, itr) where {T}
    return x ∈ itr
end

# approximate containment check for numbers in floating point
function _in(x::AbstractVector{T}, itr) where {T<:AbstractFloat}
    return any(y -> _isapprox(x, y), itr)
end

"""
    _leq(x::N, y::N;
         rtol::Real=Base.rtoldefault(N),
         ztol::Real=ABSZTOL(N),
         atol::Real=zero(N)) where {N<:AbstractFloat}

Determine if `x` is smaller than or equal to `y`.

### Input

- `x`    -- number
- `y`    -- another number (of the same numeric type as `x`)
- `rtol` -- (optional, default: `Base.rtoldefault(N)`) relative tolerance
- `ztol` -- (optional, default: `ABSZTOL(N)`) absolute tolerance for comparison
            against zero
- `atol` -- absolute tolerance

### Output

A boolean that is `true` iff `x <= y`.

### Algorithm

The `x <= y` comparison is split into `x < y` or `x ≈ y`; the latter is
implemented by extending Juila's built-in `isapprox(x, y)` with an absolute
tolerance that is used to compare against zero.
"""
function _leq(x::N, y::N;
              rtol::Real=Base.rtoldefault(N),
              ztol::Real=ABSZTOL(N),
              atol::Real=zero(N)) where {N<:AbstractFloat}
    return x <= y || _isapprox(x, y, rtol=rtol, ztol=ztol, atol=atol)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e1afd90c-7594-11eb-2738-2551799bdcea
using VegaLite, DataFrames, Query, VegaDatasets

# ╔═╡ 476b4810-7595-11eb-3c93-618fe1364786
cars = dataset("cars")

# ╔═╡ 54b2b76c-7595-11eb-36b1-cdb5851dfb07
cars |> @select(:Acceleration, :Name) |> collect

# ╔═╡ 82514288-7595-11eb-0268-e3e5cfb72ddc
function foo(data, origin)
    df = data |> @filter(_.Origin==origin) |> DataFrame

    return df |> @vlplot(:point, :Acceleration, :Miles_per_Gallon)
end

# ╔═╡ 8707c7d4-7595-11eb-27b2-dd31906a13b9
p = foo(cars, "USA")

# ╔═╡ 992e8ff6-7595-11eb-0a6e-596def39ca9d
p |> save("foo.png")

# ╔═╡ a93d61b0-7595-11eb-0982-3d302e8e7e7d

# ╔═╡ Cell order:
# ╠═e1afd90c-7594-11eb-2738-2551799bdcea
# ╠═476b4810-7595-11eb-3c93-618fe1364786
# ╠═54b2b76c-7595-11eb-36b1-cdb5851dfb07
# ╠═82514288-7595-11eb-0268-e3e5cfb72ddc
# ╠═8707c7d4-7595-11eb-27b2-dd31906a13b9
# ╠═992e8ff6-7595-11eb-0a6e-596def39ca9d
# ╠═a93d61b0-7595-11eb-0982-3d302e8e7e7d

ext: jl


text:
using DelimitedFiles
bingodata = readdlm("bingodata.txt",',',Int32)
@show bingodata

bingoboards = readdlm("bingoboards.txt",' ',Int32,'\n')
@show bingoboards
@show typeof(bingoboards)

function get_bingo_boards(bingoboarddata)
    index = 1

    global boardcount = Int(length(bingoboarddata)/25)

    boards = Vector{Matrix{Int32}}(undef,boardcount)
    println( "# of boards: $boardcount")

    for b = 1:boardcount
        #println("Get board $b")
        boards[b] = bingoboarddata[index:index+4, :]

        #println("boards[$b] = $(boards[b])")
        index += 5
    end

    return boards
end

boards = get_bingo_boards(bingoboards)
@show boards

function draw_ball(ballnumber, boards)
    for bi in eachindex(boards)
        if boards[bi][1] == 10_000
            continue
        end

        #println("boards[$bi] = $(boards[bi])")
        for bb in eachindex(boards[bi])
            if boards[bi][bb] == ballnumber
                #println("Replacing $ballnumber")
                boards[bi][bb] = -1
            end
        end
    end
end

function check_for_bingo(boards)
    bingo = fill(-1,5)
    for bi in eachindex(boards)
        if boards[bi][1] == 10_000
            continue
        end

        for a in 1:5
            if boards[bi][a,:] == bingo
                #println("BINGO! $(boards[bi])")
                return bi # boards[bi]
            end
        end
        for a in 1:5
            if boards[bi][:,a] == bingo
                #println("BINGO! $(boards[bi])")
                return bi #boards[bi]
            end
        end
    end

    return nothing
end

#function find_winner(boardcount, boards)
    losers = falses(boardcount)
    winners = trues(boardcount)

    for ball in bingodata
        draw_ball(ball, boards)
        #println("ball= $ball, boards = $boards")
        while   (bx = check_for_bingo(boards)) != nothing
            println("We found it! board # $bx")
            losers[bx] = true
            #if winners == losers
                @show "found last one: $bx"
                bingo = boards[bx]
                bingo[bingo.==-1] .= 0
                global lastbingo = deepcopy(bingo)
            @show lastbingo, bx
                global lastball = ball
            #return
            #end
            boards[bx][1] = 10_000
        end
    end
#end

#@show boards[2]
@show lastbingo
println("Sum of bingo card: $(sum(lastbingo))")
println("Times ball number: $(sum(lastbingo)*lastball)")


text:
#unitarylis.jl
#Algorithm 2.5 of Random Eigenvalues by Alan Edelman

#Experiment:    Generate random orthogonal/unitary matrices
#Theory:        Counts longest increasing subsequence statistics

## Parameters
t = 200000      # Number of trials
n = 4           # permutation size
k = 2           # length of longest increasing subsequence

include("patiencesort.jl")

function unitarylis(t, n, k)
    v = zeros(t) # samples
    ## Experiment
    for i = 1:t
        X = QR(complex(randn(k,k), randn(k,k)))[:Q]
        X = X*diagm(sign(complex(randn(k),randn(k))))
        v[i] = abs(trace(X)) ^ (2n)
    end
    z = mean(v)
    c = 0
    for i=1:factorial(n)
        c = c + int(patiencesort(nthperm([1:n],i))<=k)
    end
    return (z, c)
end

println(unitarylis(t, n, k))


text:
"""
  daxpy(alpha, x, y)

BLAS daxpy: y = alpha*x + y
"""
function daxpy(a::Float64, x::Vector{Float64}, y::Vector{Float64})
    @inbounds for i in 1:length(x)
        y[i] = a*x[i] + y[i]
    end
end

"""
    axpy(alpha, x, y)

Template version of BLAS axpy: y = alpha*x + y
"""
function axpy(a::T, x::Vector{T}, y::Vector{T}) where
    {T <: Real}
    # - This would do but slower - #
    # y = a*x + y

    # - This is equivalent to the for-loop below
    @. y = a*x + y # This is eaiser to read
    # or
    #y .= a.*x .+ y

    # - Keep this for doc reason - #
    # @inbounds for i in 1:length(x)
    #     y[i] = a*x[i] + y[i]
    # end
end


text:
# Copied from VectorLH

using StatsBase  # for countmap

"""
count_elem()

Count the number of times each element occurs in a vector
NaN is kept as a separate value unless `omitNan` is set

OUT
    outM
        matrix with values and counts as columns
"""
function count_elem(xV :: Array{T,1}, omitNan = false) where T <: Any
    if omitNan
        dictV = countmap(filter(y -> ~isnan(y), xV));
    else
        dictV = countmap(xV);
    end

    # Make dict into vector of counts
    countM = hcat([[key, val] for (key, val) in dictV]...)';
    # countV = zeros(T, size(xV));
    # for (key, val) in dictV
    #     countV[key] = val;
    # end
    return sortslices(countM, dims=1)
end

"""
count_indices()

Count the number of times each element occurs in a vector
Return vector of counts
Assumes that values are integer indices
"""
function count_indices(xV :: Array{T,1}) where T <: Integer
    countM = count_elem(xV, false);
    countV = zeros(Int64, maximum(countM[:,1]));
    countV[countM[:,1]] = countM[:, 2];
    return countV;
end

# ----------


text:
export
    ERFA_DPI,
    ERFA_D2PI,
    ERFA_DR2D,
    ERFA_DD2R,
    ERFA_DR2AS,
    ERFA_DAS2R,
    ERFA_DS2R,
    ERFA_TURNAS,
    ERFA_DMAS2R,
    ERFA_DTY,
    ERFA_DAYSEC,
    ERFA_DJ,
    ERFA_DJC,
    ERFA_DJM,
    ERFA_DJ00,
    ERFA_DJM0,
    ERFA_DJM00,
    ERFA_DJM77,
    ERFA_TTMTAI,
    ERFA_DAU,
    ERFA_CMP,
    ERFA_AULT,
    ERFA_DC,
    ERFA_ELG,
    ERFA_ELB,
    ERFA_TDB0,
    ERFA_SRS,
    ERFA_WGS84,
    ERFA_GRS80,
    ERFA_WGS72

export
    eraASTROM,
    eraLDBODY,
    eraA2af,
    eraA2tf,
    eraAb,
    eraAf2a,
    eraAnp,
    eraAnpm,
    eraApcg,
    eraApcg13,
    eraApci,
    eraApci13,
    eraApco,
    eraApco13,
    eraApcs,
    eraApcs13,
    eraAper,
    eraAper13,
    eraApio,
    eraApio13,
    eraAtci13,
    eraAtciq,
    eraAtciqn,
    eraAtciqz,
    eraAtco13,
    eraAtic13,
    eraAticq,
    eraAticqn,
    eraAtio13,
    eraAtioq,
    eraAtoc13,
    eraAtoi13,
    eraAtoiq,
    eraBi00,
    eraBp00,
    eraBp06,
    eraBpn2xy,
    eraCal2jd,
    eraC2s,
    eraC2i00a,
    eraC2i00b,
    eraC2i06a,
    eraC2ibpn,
    eraC2ixy,
    eraC2ixys,
    eraC2t00a,
    eraC2t00b,
    eraC2t06a,
    eraC2tcio,
    eraC2teqx,
    eraC2tpe,
    eraC2txy,
    eraCp,
    eraCpv,
    eraCr,
    eraDat,
    eraD2dtf,
    eraD2tf,
    eraDtdb,
    eraDtf2d,
    eraEceq06,
    eraEcm06,
    eraEe00,
    eraEe00a,
    eraEe00b,
    eraEe06a,
    eraEect00,
    eraEform,
    eraEo06a,
    eraEors,
    eraEpb,
    eraEpb2jd,
    eraEpj,
    eraEpj2jd,
    eraEpv00,
    eraEqec06,
    eraEqeq94,
    eraEra00,
    eraFad03,
    eraFae03,
    eraFaf03,
    eraFaju03,
    eraFal03,
    eraFalp03,
    eraFama03,
    eraFame03,
    eraFane03,
    eraFaom03,
    eraFapa03,
    eraFasa03,
    eraFaur03,
    eraFave03,
    eraFk52h,
    eraFk5hip,
    eraFk5hz,
    eraFw2m,
    eraFw2xy,
    eraG2icrs,
    eraGc2gd,
    eraGc2gde,
    eraGd2gc,
    eraGd2gce,
    eraGmst00,
    eraGmst06,
    eraGmst82,
    eraGst00a,
    eraGst00b,
    eraGst06,
    eraGst06a,
    eraGst94,
    eraH2fk5,
    eraHfk5z,
    eraIcrs2g,
    eraIr,
    eraJd2cal,
    eraJdcalf,
    eraLd,
    eraLdn,
    eraLdsun,
    eraLteceq,
    eraLtecm,
    eraLteqec,
    eraLtp,
    eraLtpb,
    eraLtpecl,
    eraLtpequ,
    eraNum00a,
    eraNum00b,
    eraNum06a,
    eraNumat,
    eraNut00a,
    eraNut00b,
    eraNut06a,
    eraNut80,
    eraNutm80,
    eraObl06,
    eraObl80,
    eraP06e,
    eraP2pv,
    eraP2s,
    eraPap,
    eraPas,
    eraPb06,
    eraPdp,
    eraPfw06,
    eraPlan94,
    eraPm,
    eraPmp,
    eraPmat00,
    eraPmat06,
    eraPmat76,
    eraPmpx,
    eraPmsafe,
    eraPn,
    eraPn00,
    eraPn00a,
    eraPn00b,
    eraPn06,
    eraPn06a,
    eraPnm00a,
    eraPnm00b,
    eraPnm06a,
    eraPnm80,
    eraPom00,
    eraPpp,
    eraPpsp,
    eraPr00,
    eraPrec76,
    eraPv2p,
    eraPv2s,
    eraPvdpv,
    eraPvm,
    eraPvmpv,
    eraPvppv,
    eraPvstar,
    eraPvtob,
    eraPvu,
    eraPvup,
    eraPvxpv,
    eraPxp,
    eraRefco,
    eraRm2v,
    eraRv2m,
    eraRx,
    eraRy,
    eraRz,
    eraRxp,
    eraRxpv,
    eraRxr,
    eraS00,
    eraS00a,
    eraS00b,
    eraS06,
    eraS06a,
    eraS2c,
    eraS2p,
    eraS2pv,
    eraS2xpv,
    eraSp00,
    eraSepp,
    eraSeps,
    eraStarpm,
    eraStarpv,
    eraSxp,
    eraSxpv,
    eraTaitt,
    eraTaiut1,
    eraTaiutc,
    eraTcbtdb,
    eraTcgtt,
    eraTdbtcb,
    eraTdbtt,
    eraTf2a,
    eraTf2d,
    eraTr,
    eraTrxp,
    eraTrxpv,
    eraTttai,
    eraTttcg,
    eraTttdb,
    eraTtut1,
    eraUt1tai,
    eraUt1tt,
    eraUt1utc,
    eraUtctai,
    eraUtcut1,
    eraXy06,
    eraXys00a,
    eraXys00b,
    eraXys06a,
    eraZp,
    eraZpv,
    eraZr

@deprecate ERFA_DPI ERFA.DPI
@deprecate ERFA_D2PI ERFA.D2PI
@deprecate ERFA_DR2D ERFA.DR2D
@deprecate ERFA_DD2R ERFA.DD2R
@deprecate ERFA_DR2AS ERFA.DR2AS
@deprecate ERFA_DAS2R ERFA.DAS2R
@deprecate ERFA_DS2R ERFA.DS2R
@deprecate ERFA_TURNAS ERFA.TURNAS
@deprecate ERFA_DMAS2R ERFA.DMAS2R
@deprecate ERFA_DTY ERFA.DTY
@deprecate ERFA_DAYSEC ERFA.DAYSEC
@deprecate ERFA_DJ ERFA.DJ
@deprecate ERFA_DJC ERFA.DJC
@deprecate ERFA_DJM ERFA.DJM
@deprecate ERFA_DJ00 ERFA.DJ00
@deprecate ERFA_DJM0 ERFA.DJM0
@deprecate ERFA_DJM00 ERFA.DJM00
@deprecate ERFA_DJM77 ERFA.DJM77
@deprecate ERFA_TTMTAI ERFA.TTMTAI
@deprecate ERFA_DAU ERFA.DAU
@deprecate ERFA_CMP ERFA.CMP
@deprecate ERFA_AULT ERFA.AULT
@deprecate ERFA_DC ERFA.DC
@deprecate ERFA_ELG ERFA.ELG
@deprecate ERFA_ELB ERFA.ELB
@deprecate ERFA_TDB0 ERFA.TDB0
@deprecate ERFA_SRS ERFA.SRS
@deprecate ERFA_WGS84 ERFA.WGS84
@deprecate ERFA_GRS80 ERFA.GRS80
@deprecate ERFA_WGS72 ERFA.WGS72

@deprecate ASTROM ERFA.ASTROM
@deprecate LDBODY ERFA.LDBODY
@deprecate eraA2af ERFA.a2af
@deprecate eraA2tf ERFA.a2tf
@deprecate eraAb ERFA.ab
@deprecate eraAf2a ERFA.af2a
@deprecate eraAnp ERFA.anp
@deprecate eraAnpm ERFA.anpm
@deprecate eraApcg ERFA.apcg
@deprecate eraApcg13 ERFA.apcg13
@deprecate eraApci ERFA.apci
@deprecate eraApci13 ERFA.apci13
@deprecate eraApco ERFA.apco
@deprecate eraApco13 ERFA.apco13
@deprecate eraApcs ERFA.apcs
@deprecate eraApcs13 ERFA.apcs13
@deprecate eraAper ERFA.aper
@deprecate eraAper13 ERFA.aper13
@deprecate eraApio ERFA.apio
@deprecate eraApio13 ERFA.apio13
@deprecate eraAtci13 ERFA.atci13
@deprecate eraAtciq ERFA.atciq
@deprecate eraAtciqn ERFA.atciqn
@deprecate eraAtciqz ERFA.atciqz
@deprecate eraAtco13 ERFA.atco13
@deprecate eraAtic13 ERFA.atic13
@deprecate eraAticq ERFA.aticq
@deprecate eraAticqn ERFA.aticqn
@deprecate eraAtio13 ERFA.atio13
@deprecate eraAtioq ERFA.atioq
@deprecate eraAtoc13 ERFA.atoc13
@deprecate eraAtoi13 ERFA.atoi13
@deprecate eraAtoiq ERFA.atoiq
@deprecate eraBi00 ERFA.bi00
@deprecate eraBp00 ERFA.bp00
@deprecate eraBp06 ERFA.bp06
@deprecate eraBpn2xy ERFA.bpn2xy
@deprecate eraCal2jd ERFA.cal2jd
@deprecate eraC2s ERFA.c2s
@deprecate eraC2i00a ERFA.c2i00a
@deprecate eraC2i00b ERFA.c2i00b
@deprecate eraC2i06a ERFA.c2i06a
@deprecate eraC2ibpn ERFA.c2ibpn
@deprecate eraC2ixy ERFA.c2ixy
@deprecate eraC2ixys ERFA.c2ixys
@deprecate eraC2t00a ERFA.c2t00a
@deprecate eraC2t00b ERFA.c2t00b
@deprecate eraC2t06a ERFA.c2t06a
@deprecate eraC2tcio ERFA.c2tcio
@deprecate eraC2teqx ERFA.c2teqx
@deprecate eraC2tpe ERFA.c2tpe
@deprecate eraC2txy ERFA.c2txy
@deprecate eraCp Base.copy
@deprecate eraCpv Base.copy
@deprecate eraCr Base.copy
@deprecate eraDat ERFA.dat
@deprecate eraD2dtf ERFA.d2dtf
@deprecate eraD2tf ERFA.d2tf
@deprecate eraDtdb ERFA.dtdb
@deprecate eraDtf2d ERFA.dtf2d
@deprecate eraEceq06 ERFA.eceq06
@deprecate eraEcm06 ERFA.ecm06
@deprecate eraEe00 ERFA.ee00
@deprecate eraEe00a ERFA.ee00a
@deprecate eraEe00b ERFA.ee00b
@deprecate eraEe06a ERFA.ee06a
@deprecate eraEect00 ERFA.eect00
@deprecate eraEform ERFA.eform
@deprecate eraEo06a ERFA.eo06a
@deprecate eraEors ERFA.eors
@deprecate eraEpb ERFA.epb
@deprecate eraEpb2jd ERFA.epb2jd
@deprecate eraEpj ERFA.epj
@deprecate eraEpj2jd ERFA.epj2jd
@deprecate eraEpv00 ERFA.epv00
@deprecate eraEqec06 ERFA.eqec06
@deprecate eraEqeq94 ERFA.eqeq94
@deprecate eraEra00 ERFA.era00
@deprecate eraFad03 ERFA.fad03
@deprecate eraFae03 ERFA.fae03
@deprecate eraFaf03 ERFA.faf03
@deprecate eraFaju03 ERFA.faju03
@deprecate eraFal03 ERFA.fal03
@deprecate eraFalp03 ERFA.falp03
@deprecate eraFama03 ERFA.fama03
@deprecate eraFame03 ERFA.fame03
@deprecate eraFane03 ERFA.fane03
@deprecate eraFaom03 ERFA.faom03
@deprecate eraFapa03 ERFA.fapa03
@deprecate eraFasa03 ERFA.fasa03
@deprecate eraFaur03 ERFA.faur03
@deprecate eraFave03 ERFA.fave03
@deprecate eraFk52h ERFA.fk52h
@deprecate eraFk5hip ERFA.fk5hip
@deprecate eraFk5hz ERFA.fk5hz
@deprecate eraFw2m ERFA.fw2m
@deprecate eraFw2xy ERFA.fw2xy
@deprecate eraG2icrs ERFA.g2icrs
@deprecate eraGc2gd ERFA.gc2gd
@deprecate eraGc2gde ERFA.gc2gde
@deprecate eraGd2gc ERFA.gd2gc
@deprecate eraGd2gce ERFA.gd2gce
@deprecate eraGmst00 ERFA.gmst00
@deprecate eraGmst06 ERFA.gmst06
@deprecate eraGmst82 ERFA.gmst82
@deprecate eraGst00a ERFA.gst00a
@deprecate eraGst00b ERFA.gst00b
@deprecate eraGst06 ERFA.gst06
@deprecate eraGst06a ERFA.gst06a
@deprecate eraGst94 ERFA.gst94
@deprecate eraH2fk5 ERFA.h2fk5
@deprecate eraHfk5z ERFA.hfk5z
@deprecate eraIcrs2g ERFA.icrs2g
@deprecate eraIr() Base.eye(3)
@deprecate eraJd2cal ERFA.jd2cal
@deprecate eraJdcalf ERFA.jdcalf
@deprecate eraLd ERFA.ld
@deprecate eraLdn ERFA.ldn
@deprecate eraLdsun ERFA.ldsun
@deprecate eraLteceq ERFA.lteceq
@deprecate eraLtecm ERFA.ltecm
@deprecate eraLteqec ERFA.lteqec
@deprecate eraLtp ERFA.ltp
@deprecate eraLtpb ERFA.ltpb
@deprecate eraLtpecl ERFA.ltpecl
@deprecate eraLtpequ ERFA.ltpequ
@deprecate eraNum00a ERFA.num00a
@deprecate eraNum00b ERFA.num00b
@deprecate eraNum06a ERFA.num06a
@deprecate eraNumat ERFA.numat
@deprecate eraNut00a ERFA.nut00a
@deprecate eraNut00b ERFA.nut00b
@deprecate eraNut06a ERFA.nut06a
@deprecate eraNut80 ERFA.nut80
@deprecate eraNutm80 ERFA.nutm80
@deprecate eraObl06 ERFA.obl06
@deprecate eraObl80 ERFA.obl80
@deprecate eraP06e ERFA.p06e
@deprecate eraP2pv ERFA.p2pv
@deprecate eraP2s ERFA.p2s
@deprecate eraPap ERFA.pap
@deprecate eraPas ERFA.pas
@deprecate eraPb06 ERFA.pb06
@deprecate eraPdp ERFA.pdp
@deprecate eraPfw06 ERFA.pfw06
@deprecate eraPlan94 ERFA.plan94
@deprecate eraPm ERFA.pm
@deprecate eraPmp ERFA.pmp
@deprecate eraPmat00 ERFA.pmat00
@deprecate eraPmat06 ERFA.pmat06
@deprecate eraPmat76 ERFA.pmat76
@deprecate eraPmpx ERFA.pmpx
@deprecate eraPmsafe ERFA.pmsafe
@deprecate eraPn ERFA.pn
@deprecate eraPn00 ERFA.pn00
@deprecate eraPn00a ERFA.pn00a
@deprecate eraPn00b ERFA.pn00b
@deprecate eraPn06 ERFA.pn06
@deprecate eraPn06a ERFA.pn06a
@deprecate eraPnm00a ERFA.pnm00a
@deprecate eraPnm00b ERFA.pnm00b
@deprecate eraPnm06a ERFA.pnm06a
@deprecate eraPnm80 ERFA.pnm80
@deprecate eraPom00 ERFA.pom00
@deprecate eraPpp ERFA.ppp
@deprecate eraPpsp ERFA.ppsp
@deprecate eraPr00 ERFA.pr00
@deprecate eraPrec76 ERFA.prec76
@deprecate eraPv2p ERFA.pv2p
@deprecate eraPv2s ERFA.pv2s
@deprecate eraPvdpv ERFA.pvdpv
@deprecate eraPvm ERFA.pvm
@deprecate eraPvmpv ERFA.pvmpv
@deprecate eraPvppv ERFA.pvppv
@deprecate eraPvstar ERFA.pvstar
@deprecate eraPvtob ERFA.pvtob
@deprecate eraPvu ERFA.pvu
@deprecate eraPvup ERFA.pvup
@deprecate eraPvxpv ERFA.pvxpv
@deprecate eraPxp ERFA.pxp
@deprecate eraRefco ERFA.refco
@deprecate eraRm2v ERFA.rm2v
@deprecate eraRv2m ERFA.rv2m
@deprecate eraRx ERFA.rx
@deprecate eraRy ERFA.ry
@deprecate eraRz ERFA.rz
@deprecate eraRxp ERFA.rxp
@deprecate eraRxpv ERFA.rxpv
@deprecate eraRxr ERFA.rxr
@deprecate eraS00 ERFA.s00
@deprecate eraS00a ERFA.s00a
@deprecate eraS00b ERFA.s00b
@deprecate eraS06 ERFA.s06
@deprecate eraS06a ERFA.s06a
@deprecate eraS2c ERFA.s2c
@deprecate eraS2p ERFA.s2p
@deprecate eraS2pv ERFA.s2pv
@deprecate eraS2xpv ERFA.s2xpv
@deprecate eraSp00 ERFA.sp00
@deprecate eraSepp ERFA.sepp
@deprecate eraSeps ERFA.seps
@deprecate eraStarpm ERFA.starpm
@deprecate eraStarpv ERFA.starpv
@deprecate eraSxp ERFA.sxp
@deprecate eraSxpv ERFA.sxpv
@deprecate eraTaitt ERFA.taitt
@deprecate eraTaiut1 ERFA.taiut1
@deprecate eraTaiutc ERFA.taiutc
@deprecate eraTcbtdb ERFA.tcbtdb
@deprecate eraTcgtt ERFA.tcgtt
@deprecate eraTdbtcb ERFA.tdbtcb
@deprecate eraTdbtt ERFA.tdbtt
@deprecate eraTf2a ERFA.tf2a
@deprecate eraTf2d ERFA.tf2d
@deprecate eraTr ERFA.tr
@deprecate eraTrxp ERFA.trxp
@deprecate eraTrxpv ERFA.trxpv
@deprecate eraTttai ERFA.tttai
@deprecate eraTttcg ERFA.tttcg
@deprecate eraTttdb ERFA.tttdb
@deprecate eraTtut1 ERFA.ttut1
@deprecate eraUt1tai ERFA.ut1tai
@deprecate eraUt1tt ERFA.ut1tt
@deprecate eraUt1utc ERFA.ut1utc
@deprecate eraUtctai ERFA.utctai
@deprecate eraUtcut1 ERFA.utcut1
@deprecate eraXy06 ERFA.xy06
@deprecate eraXys00a ERFA.xys00a
@deprecate eraXys00b ERFA.xys00b
@deprecate eraXys06a ERFA.xys06a
@deprecate eraZp() Base.zeros(3)
@deprecate eraZpv() Base.zeros(6)
@deprecate eraZr() Base.zeros(3, 3)


text:
using Dice
using Dice: num_flips, num_nodes

for n = 7:7
    code = @dice begin

        function uniform(bitrange::Int, bitwidth::Int)
            ans = fill(DistBool(dicecontext(), false), bitwidth)
            for i = 1:bitrange
                ans[i] = flip(0.5)
            end
            DistInt(ans)
        end

        function gcd(a::DistInt, b::DistInt)
            for _ = 1 : 1 + max_bits(b) ÷ log2(MathConstants.golden)
                t = b
                converged = prob_equals(b, 0)
                amb, _ = (a % b)
                b = if converged
                        b
                    else
                        amb
                    end
                a = if converged
                        a
                    else
                        t
                    end
            end
            return a
        end

        a, _ = uniform(n, n+1) + 1
        b, _ = uniform(n, n+1) + 1
        g = gcd(a, b)
        prob_equals(g, 1)
    end

    @show n
    @show 2^n

    bdd = compile(code)
    @show num_flips(bdd)
    @show num_nodes(bdd)

    @time @show infer(code, :bdd)
end


text:
# The Dockerfile installs Pluto in a system-wide depot.
PLUTO_DEPOT = Base.DEPOT_PATH[2]

# For each package environment in our custom depot, we force-copy the
# environment to the first-listed depot, which lives on /var because it must be
# writable by the user.
import Pkg
Pkg.activate("pluto-env", shared=true)
origdir = Pkg.envdir(PLUTO_DEPOT)
envdir = Pkg.envdir()
mkpath(envdir)
for dir in readdir(origdir)
    rm("$envdir/$dir"; force=true, recursive=true)
    run(`cp -a $origdir/$dir $envdir/$dir`)
end

# Start Pluto with appropriate options
import Pluto
Pluto.run(
    port=8000,
    launch_browser=false,
    require_secret_for_access=false,
    dismiss_update_notification=true,
    show_file_system=false,
)

max_stars_repo_name: garrison/sandstorm-pluto-jl

max_issues_repo_name: garrison/sandstorm-pluto-jl

max_forks_repo_name: garrison/sandstorm-pluto-jl


text:
using Test

using iLQGames: SystemTrajectory, time_cont2disc, timepoints, zero!

h, ΔT, nx, nu = 100, 0.1, 4, 2
traj = Base.zero(SystemTrajectory{h, ΔT, nx, nu})

@testset "time" begin
    @test time_cont2disc(traj, 5.0) ==  51
    @test collect(timepoints(traj)) ≈ collect(0:0.1:9.9)
end

@testset "zero" begin
    traj_pertubed = deepcopy(traj)
    for i in eachindex(traj_pertubed.x)
        traj_pertubed.x[i] = randn(typeof(traj_pertubed.x[i]))
        traj_pertubed.u[i] = randn(typeof(traj_pertubed.u[i]))
    end
    @assert any(!iszero(x) for x in traj_pertubed.x)
    @assert any(!iszero(u) for u in traj_pertubed.u)

    function test_all_zero(traj)
        @test all(iszero(x) for x in zero(traj.x))
        @test all(iszero(u) for u in zero(traj.u))
    end

    test_all_zero(traj)
    test_all_zero(zero(traj_pertubed))
    zero!(traj_pertubed)
    test_all_zero(traj_pertubed)
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from quaternion to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export quat_to_dcm

"""
    quat_to_dcm(q::Quaternion)

Convert the quaternion `q` to a Direction Cosine Matrix (DCM).

# Examples

```jldoctest
julia> q = Quaternion(cosd(45/2), sind(45/2), 0, 0);

julia> quat_to_dcm(q)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
 1.0   0.0       0.0
 0.0   0.707107  0.707107
 0.0  -0.707107  0.707107
```
"""
function quat_to_dcm(q::Quaternion)
    # Auxiliary variables.
    q0 = q.q0
    q1 = q.q1
    q2 = q.q2
    q3 = q.q3

    return DCM(
        q0^2+q1^2-q2^2-q3^2,   2(q1*q2+q0*q3)   ,   2(q1*q3-q0*q2),
          2(q1*q2-q0*q3)   , q0^2-q1^2+q2^2-q3^2,   2(q2*q3+q0*q1),
          2(q1*q3+q0*q2)   ,   2(q2*q3-q0*q1)   , q0^2-q1^2-q2^2+q3^2
    )'
end


text:
cfg = NeuroEvolution.get_config("../cfg/test.yaml")
cfg["n_out"] = 1
cfg["n_in"] = 2
cfg["log_gen"] = 0
cfg["save_gen"] = 0

@testset "Evolution creation" begin
    # test_fitness::Function = x::NEATIndiv -> fitness_xor(x, cfg["n_in"])
    e::Evolution = NEAT(cfg, fitness_xor, cfg["n_in"])
    # @test length(e.population) == cfg["n_population"]
    @test all(typeof.(e.population) .== NEATIndividual)
end

@testset "Evaluate" begin
    # test_fitness::Function = x::NEATIndiv -> fitness_xor(x, cfg["n_in"])
    e::Evolution = NEAT(cfg, fitness_xor, cfg["n_in"])
    e.evaluate(e)
    @test all(getfield.(e.population, :fitness) .!= [-Inf])
end

@testset "Populate" begin
    cfg = NeuroEvolution.get_config("../cfg/test.yaml")
    # test_fitness::Function = x::NEATIndiv -> fitness_xor(x, cfg["n_in"])
    e::Evolution = NEAT(cfg, fitness_xor, cfg["n_in"])
    e.evaluate(e)
    e.populate(e)
    @test length(e.population) > 0
end

@testset "NEAT on XOR" begin
    cfg = NeuroEvolution.get_config("../cfg/test.yaml")
    cfg["verbose"]=true
    # Run evolution
    # test_fitness::Function = x::NEATIndiv -> fitness_xor(x, cfg["n_in"])
    e = NEAT(cfg, fitness_xor, cfg["n_in"])
    Cambrian.run!(e)

    @test length(e.population) <= 1.1 * cfg["n_population"]

    X, y = xor_dataset(cfg["n_in"], 100)
    max_f = sum(log_fitness.(y, y)) / length(X)
    y_false = []
    for i in y
        push!(y_false, 1 .- i)
    end
    min_f = sum(log_fitness.(y, y_false)) / length(X)
    if cfg["verbose"]
        println("\n-- NEAT on XOR --")
        println("Worst possible fitness: ", min_f)
        println("Best  possible fitness: ", max_f)
    end
    # Analyse results
    best = sort(e.population, rev=true)
    @test best[1].fitness[1] <= max_f
    @test best[1].fitness[1] > min_f

    y = process(best[1], [1., 1.])
    @test typeof(y[1])==Float64

    if cfg["verbose"]
        println("Final fitness: ", best[1].fitness[1])
        println("Population size: ", length(e.population))
        println("Genome size: ", length(best[1].genes))
        println("Max innovation: ", e.cfg["innovation_max"])
        println("Species: ", length(e.cfg["Species"]))
        println("\nBest individuals: ", getfield.(best[1:10], :fitness), "\n")
        println("\n\n", best[1])
    end
end

@testset "Cambrian.GA x NEAT" begin
    cfg = NeuroEvolution.get_config("../cfg/ga.yaml")

    e = GA_NEAT(NEATIndividual, cfg, fitness_xor; id="test")
    step!(e)
    @test length(e.population) == cfg["n_population"]
    run!(e)
    best = sort(e.population, rev=true)

    X, y = xor_dataset(cfg["n_in"], 100)
    max_f = sum(log_fitness.(y, y)) / length(X)
    y_false = []
    for i in y
        push!(y_false, 1 .- i)
    end
    min_f = sum(log_fitness.(y, y_false)) / length(X)
    @test best[1].fitness[1] <= max_f
    @test best[1].fitness[1] > min_f

    y = process(best[1], [1., 1.])
    @test typeof(y[1])==Float64
end


text:
using Test
using ConstraintSolver
const CS = ConstraintSolver

import MathOptInterface
const MOI = MathOptInterface

import ForwardDiff
using JuMP
using LinearAlgebra: dot

@testset "StableSet" begin
    @testset "Weighted stable set" begin
        matrix = [
            0 1 1 0
            1 0 1 0
            1 1 0 1
            0 0 1 0
        ]
        model = CS.Optimizer(logging = [], keep_logs = true)
        x = [MOI.add_constrained_variable(model, MOI.ZeroOne()) for _ in 1:4]
        for i in 1:4, j in 1:4
            if matrix[i, j] == 1 && i < j
                (z, _) = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
                MOI.add_constraint(model, z, MOI.Integer())
                MOI.add_constraint(model, z, MOI.LessThan(1.0))
                f = MOI.ScalarAffineFunction(
                    [
                        MOI.ScalarAffineTerm(1.0, x[i][1]),
                        MOI.ScalarAffineTerm(1.0, x[j][1]),
                        MOI.ScalarAffineTerm(1.0, z),
                    ],
                    0.0,
                )
                MOI.add_constraint(model, f, MOI.EqualTo(1.0))
            end
        end
        weights = [0.2, 0.1, 0.2, 0.1]
        terms = [MOI.ScalarAffineTerm(weights[i], x[i][1]) for i in eachindex(x)]
        objective = MOI.ScalarAffineFunction(terms, 0.0)
        MOI.set(model, MOI.ObjectiveFunction{typeof(objective)}(), objective)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        MOI.optimize!(model)
        var_x = [x[i][1] for i in 1:4]
        CS.save_logs(CS.get_inner_model(model), "stable_set.json", :x => var_x)
        rm("stable_set.json")

        @test MOI.get(model, MOI.VariablePrimal(), x[4][1]) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x[1][1]) == 1
        @test MOI.get(model, MOI.ObjectiveValue()) ≈ 0.3
    end

    @testset "Weighted stable set JuMP" begin
        matrix = [
            0 1 1 0
            1 0 1 0
            1 1 0 1
            0 0 1 0
        ]
        m = Model(CSJuMPTestOptimizer())
        x = @variable(m, x[1:4], Bin)
        for i in 1:4, j in (i + 1):4
            if matrix[i, j] == 1
                zcomp = @variable(m)
                JuMP.set_binary(zcomp)
                @constraint(m, x[i] + x[j] + zcomp == 1)
            end
        end
        w = [0.2, 0.1, 0.2, 0.1]
        @objective(m, Max, dot(w, x))
        optimize!(m)
        @test JuMP.value(x[4]) ≈ 1
        @test JuMP.value(x[1]) ≈ 1
        @test JuMP.value(x[2]) ≈ 0
        @test JuMP.value(x[3]) ≈ 0
        @test JuMP.objective_value(m) ≈ 0.3
    end

    @testset "Bigger stable set using <= " begin
        matrix = [
            0 1 0 1 0 0 1 0 0 1
            1 0 1 0 1 0 0 0 0 0
            0 1 0 1 0 0 0 0 1 0
            1 0 1 0 1 1 0 0 0 0
            0 1 0 1 0 1 1 1 0 0
            0 0 0 1 1 0 1 0 0 0
            1 0 0 0 1 1 0 1 0 0
            0 0 0 0 1 0 1 0 1 0
            0 0 1 0 0 0 0 1 0 1
            1 0 0 0 0 0 0 0 1 0
        ]
        n = 10
        m = Model(CSJuMPTestOptimizer())
        @variable(m, x[1:n], Bin)
        nconstraints = 0
        for i in 1:n, j in (i + 1):n
            if matrix[i, j] == 1
                @constraint(m, x[i] + x[j] <= 1)
                nconstraints += 1
            end
        end
        w = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
        @objective(m, Max, dot(w, x))
        optimize!(m)
        com = CS.get_inner_model(m)
        @test is_solved(com)
        @test com.info.n_constraint_types.inequality == length(com.constraints)
        @test com.info.n_constraint_types.inequality == nconstraints

        @test MOI.get(m, MOI.ObjectiveValue()) ≈ 2.7
    end

    @testset "Bigger stable set using with Cbc" begin
        matrix = [
            0 1 0 1 0 0 1 0 0 1
            1 0 1 0 1 0 0 0 0 0
            0 1 0 1 0 0 0 0 1 0
            1 0 1 0 1 1 0 0 0 0
            0 1 0 1 0 1 1 1 0 0
            0 0 0 1 1 0 1 0 0 0
            1 0 0 0 1 1 0 1 0 0
            0 0 0 0 1 0 1 0 1 0
            0 0 1 0 0 0 0 1 0 1
            1 0 0 0 0 0 0 0 1 0
        ]
        n = 10
        cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
        m = Model(optimizer_with_attributes(
            CS.Optimizer,
            "lp_optimizer" => cbc_optimizer,
            "logging" => [],
        ))
        @variable(m, x[1:n], Bin)
        nconstraints = 0
        for i in 1:n, j in (i + 1):n
            if matrix[i, j] == 1
                @constraint(m, x[i] + x[j] <= 1)
                nconstraints += 1
            end
        end
        w = [0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
        @objective(m, Max, dot(w, x))
        optimize!(m)
        com = CS.get_inner_model(m)
        @test com.info.n_constraint_types.inequality == length(com.constraints)
        @test com.info.n_constraint_types.inequality == nconstraints

        @test MOI.get(m, MOI.ObjectiveValue()) ≈ 2.7
    end

    function weighted_stable_set(w)
        matrix = [
            0 1 1 0
            1 0 1 0
            1 1 0 1
            0 0 1 0
        ]
        model = CS.Optimizer(solution_type = Real, logging = [])
        x = [MOI.add_constrained_variable(model, MOI.ZeroOne()) for _ in 1:4]
        for i in 1:4, j in 1:4
            if matrix[i, j] == 1 && i < j
                (z, _) = MOI.add_constrained_variable(model, MOI.GreaterThan(0.0))
                MOI.add_constraint(model, z, MOI.Integer())
                MOI.add_constraint(model, z, MOI.LessThan(1.0))
                f = MOI.ScalarAffineFunction(
                    [
                        MOI.ScalarAffineTerm(1.0, x[i][1]),
                        MOI.ScalarAffineTerm(1.0, x[j][1]),
                        MOI.ScalarAffineTerm(1.0, z),
                    ],
                    0.0,
                )
                MOI.add_constraint(model, f, MOI.EqualTo(1.0))
            end
        end
        # sum(w[i] * x[i] for i in V) - stable_set == 0
        terms = [MOI.ScalarAffineTerm(w[i], x[i][1]) for i in eachindex(x)]
        objective = MOI.ScalarAffineFunction(terms, zero(eltype(w)))
        MOI.set(model, MOI.ObjectiveFunction{typeof(objective)}(), objective)
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
        MOI.optimize!(model)
        return MOI.get(model, MOI.ObjectiveValue())
    end

    # type piracy for the greater good
    # see https://github.com/JuliaDiff/ForwardDiff.jl/issues/318
    function Base.typemin(::Type{ForwardDiff.Dual{T,V,N}}) where {T,V,N}
        ForwardDiff.Dual{T,V,N}(typemin(V))
    end

    @testset "Differentiating stable set" begin
        weights = [0.2, 0.1, 0.2, 0.1]
        # ∇w = Zygote.gradient(weighted_stable_set, weights)
        ∇w = ForwardDiff.gradient(weighted_stable_set, weights)
        @test ∇w[1] ≈ 1
        @test ∇w[4] ≈ 1
        @test ∇w[2] ≈ ∇w[3] ≈ 0
    end
end


text:
# template log-on and fetch graph

using Caesar
using IncrementalInference

include(joinpath(dirname(@__FILE__),"blandauthremote.jl"))
addrdict["session"] = "SESSSHARK_16_11_14"  # "SESSROX"
addrdict["robot"] = "robot"
addrdict["user"] = "user"
cloudGraph, addrdict = standardcloudgraphsetup(addrdict=addrdict)

fg = initfg(cloudgraph=cloudGraph, sessionname=addrdict["session"])
fullLocalGraphCopy!(fg)
savejld(fg, file="sandshark_testing_seg417_iNA.jld")

plotKDE(fg, :x1710, dims=[1;2], api=localapi)

plotKDE(fg, :x1723, dims=[1;2], api=localapi)
plotKDE(fg, :x1724, dims=[3], api=localapi)
plotKDE(fg, :x1725, dims=[3], api=localapi)

plotKDE(fg, :x1662, dims=[1;2], api=localapi)

plotKDE(fg, :l200050, dims=[1;2], api=localapi)
plotLocalProduct(fg, :l200050, dims=[1;2], api=localapi)

# pp, = localProduct(fg, :x1710, api=localapi)

plotLocalProduct(fg, :x1710, dims=[1;2], api=localapi)
plotLocalProduct(fg, :x1723, dims=[1;2], api=localapi, levels=1)
pl = plotLocalProduct(fg, :x1724, dims=[1;2], api=localapi, levels=1)

using KernelDensityEstimate
pl = plotKDEresiduals(fg, :x1710l200050, fill=true)

fncvar = getfnctype(getVert(fg,fg.fIDs[:x1710l200050]))

using Gadfly, Colors

Gadfly.draw(PDF("residual_x1710l200050_NA.pdf",15cm,10cm),pl)

h1 = plotKDE(fncvar.range, title=string("Likelihood density estimate, $(:x1710l200050)\n", fncvar), xlbl="Range [m]", fill=true )
h2 = plotKDE(fncvar.bearing, xlbl="Bearing [rad]", fill=true )

Gadfly.draw(PDF("factor_x1710l200050_stack.pdf",15cm,20cm),vstack(h1,h2))

# TODO: Please fix with correct session for example.
ret = whosNear2D(cloudGraph, "SESSSHARK_16_11_14", "robot", "user" x=21.7, y=-54.0)

sort(collect(keys(ret)))

#


text:
module thread_test

using Base.Threads
using Logging
@Logging.configure(level=DEBUG)

function threaded_call(A)
    tid = threadid()
    A[tid] = tid*100
end

function test_loop()
    A = zeros(Float64, 2,3)

    @threads all for i=1:nthreads()
        tid = threadid()
        A[i] = tid*10
    end
    A
end

function test_block()
    A = zeros(Float64,2,3)
    @threads all begin
        tid = threadid()
        A[tid] = tid
    end
    @show A
end

function test_atomic_add()
    A = zeros(Int64,2,3)
    a = Atomic()
    @threads all for i=1:6
        atomic_add!(a,1)
    end
    @show a
end

function test1()
    a = rand(Int16, 10)
    @show a
    @threads all threaded_call(a)

end

function threaded_add_using_spinlock()
    s = SpinLock()
    x = 0
    n = 4
    @threads all for i=1:n
        lock!(s)
        x = x + 1
        unlock!(s)
    end
    @show x
end

function threaded_add_using_mutex()
    m = Mutex()
    x = 0
    n = 4
    @threads all for i=1:n
        lock!(m)
        x += 1
        unlock!(m)
    end

    @show x
end

function test_spinlocks()

end

end


text:
export evalObj, evalObjAndGrad

function myMap(f::Function,Θ::AbstractArray)
    return f(Θ)
end

function myMap(f::Function,Θ::Tuple)
    fΘ = Array{Any}(undef,length(Θ))
    for k=1:length(Θ)
        fΘ[k] = myMap(f,Θ[k])
    end
    return tuple(fΘ...)
end

function evalObjAndGrad(J,Θ::Vector,parms,ps)
    # put theta into vectors
    parms = vec2param!(Θ,parms)
    Jc,back = Zygote.Tracker.forward(()->J(parms), ps)
    gc = back(1)

    dJ = Θ .* 0.0
    cnt = 0;
    for p in ps
        gp = vec(gc[p].data)
        dJ[cnt+1:cnt+length(gp)] = gp
        cnt +=length(gp)
    end
    return Jc.data,dJ
end

function evalObj(J,Θ::Vector,parms,ps)
    # put theta into vectors
    parms = vec2param!(Θ,parms)
    Jc = J(parms)
    return Jc.data
end


text:
@testset "Testing random_mc.jl" begin
    @testset "Test random_markov_chain" begin
        n, k = 5, 3
        mc_dicts = (Dict("P" => random_markov_chain(n).p, "k" => n),
                    Dict("P" => random_markov_chain(n, k).p, "k" => k))
        for d in mc_dicts
            P = d["P"]
            @test size(P) == (n, n)
            @test all(x->(countnz(x)==d["k"]),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end
    end

    @testset "Test random_stochastic_matrix" begin
        n, k = 5, 3
        Ps = (random_stochastic_matrix(n), random_stochastic_matrix(n, k))
        for P in Ps
            @test all(P .>= 0) == true
            @test all(x->isapprox(sum(x), 1),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end
    end

    @testset "Test random_stochastic_matrix with k=1" begin
        n, k = 3, 1
        P = random_stochastic_matrix(n, k)
        @test all((P .== 0) .| (P .== 1))
        @test all(x->isequal(sum(x), 1),
                  [P[i, :] for i in 1:size(P)[1]]) == true
    end

    @testset "Test errors properly thrown" begin
        # n <= 0
        @test_throws ArgumentError random_markov_chain(0)

        # k <= 0
        @test_throws ArgumentError random_markov_chain(2, 0)

        # k > n
        @test_throws ArgumentError random_markov_chain(2, 3)
    end

    @testset "Test random_discrete_dp" begin
        num_states, num_actions = 5, 4
        num_sa = num_states * num_actions
        k = 3
        ddp = random_discrete_dp(num_states, num_actions; k=k)

        # Check shapes
        @test size(ddp.R) == (num_states, num_actions)
        @test size(ddp.Q) == (num_states, num_actions, num_states)

        # Check ddp.Q[:, a, :] is a stochastic matrix for all actions `a`
        @test all(ddp.Q .>= 0) == true
        for a in 1:num_actions
            P = reshape(ddp.Q[:, a, :], (num_states, num_states))
            @test all(x->isapprox(sum(x), 1),
                      [P[i, :] for i in 1:size(P)[1]]) == true
        end

        # Check number of nonzero entries for each state-action pair
        @test sum(ddp.Q .> 0, 3) ==
            ones(Int, (num_states, num_actions, 1)) * k
    end
end  # @testset


text:
using Distributed, PortAudio

# Modified from Jiahao Chen's example in the obsolete AudioIO module.
# Will use first output device found in system's listing or DEFAULTDEVICE if set below
const DEFAULTDEVICE = -1

function paudio()
    devs = PortAudio.devices()
    if DEFAULTDEVICE < 0
        devnum = findfirst(x -> x.maxoutchans > 0, devs)
        (devnum == nothing) && error("No output device for audio found")
    else
        devnum = DEFAULTDEVICE + 1
    end
    return ostream = PortAudioStream(devs[devnum].name, 0, 2)
end

play(ostream, sample::Array{Float64, 1}) = write(ostream, sample)
play(ostr, sample::Array{Int64, 1}) = play(ostr, Float64.(sample))

struct Note{S <: Real, T <: Real}
    pitch::S
    duration::T
    sustained::Bool
end

function play(
    ostream,
    A::Note,
    samplingfreq::Real = 44100,
    shape::Function = t -> 0.6sin(t) + 0.2sin(2t) + 0.05 * sin(8t),
)
    timesamples = 0:(1 / samplingfreq):(A.duration * (A.sustained ? 0.98 : 0.9))
    v = Float64[shape(2π * A.pitch * t) for t in timesamples]
    if !A.sustained
        decay_length = div(length(timesamples), 5)
        v[(end - decay_length):(end - 1)] =
            v[(end - decay_length):(end - 1)] .* LinRange(1, 0, decay_length)
    end
    play(ostream, v)
    sleep(A.duration)
end

function parsevoice(melody::String; tempo = 132, beatunit = 4, lyrics = nothing)
    ostream = paudio() # initialize audio for output
    lyrics_syllables = lyrics == nothing ? nothing : split(lyrics)
    lyrics_syllables != nothing && (lyrics_syllables[end] *= "\n")
    note_idx = 1
    oldduration = 4
    for line in split(melody, '\n')
        percent_idx = findfirst('%', line) # Trim comment
        percent_idx == nothing || (line = line[1:(percent_idx - 1)])
        for token in split(line)
            pitch, duration, dotted, sustained = parsetoken(token)
            duration == nothing && (duration = oldduration)
            oldduration = duration
            dotted && (duration *= 1.5)
            if lyrics_syllables != nothing && 1 <= note_idx <= length(lyrics_syllables)
                # Print the lyrics, omitting hyphens
                if lyrics_syllables[note_idx][end] == '-'
                    print(join(split(lyrics_syllables[note_idx][:], "")[1:(end - 1)]), "")
                else
                    print(lyrics_syllables[note_idx], ' ')
                end
            end
            play(ostream, Note(pitch, (beatunit / duration) * (60 / tempo), sustained))
            note_idx += 1
        end
    end
end

function parsetoken(token, Atuning::Real = 220)
    state = :findpitch
    pitch = 0.0
    sustain = dotted = false
    lengthbuf = Char[]
    for char in token
        if state == :findpitch
            scale_idx =
                something(findfirst(char, String(collect('a':'g'))), 0) +
                something(findfirst(char, String(collect('A':'G'))), 0)
            if scale_idx != 0
                halfsteps = [12, 14, 3, 5, 7, 8, 10]
                pitch = Atuning * 2^(halfsteps[scale_idx] / 12)
                state = :findlength
            elseif char == 'r'
                pitch, state = 0, :findlength
            else
                error("unknown pitch: $char")
            end
        elseif state == :findlength
            if char == '#'
                pitch *= 2^(1 / 12) # sharp
            elseif char == 'b'
                pitch /= 2^(1 / 12) # flat
            elseif char == '\''
                pitch *= 2          # higher octave
            elseif char == ','
                pitch /= 2          # lower octave
            elseif char == '.'
                dotted = true       # dotted note
            elseif char == '~'
                sustain = true      # tied note
            else
                push!(lengthbuf, char)
                # Check for "is" and "es" suffixes for sharps and flats
                if length(lengthbuf) >= 2
                    if lengthbuf[(end - 1):end] == "is"
                        pitch *= 2^(1 / 12)
                        lengthbuf = lengthbuf[1:(end - 2)]
                    elseif lengthbuf[(end - 1):end] == "es"
                        pitch /= 2^(1 / 12)
                        lengthbuf = lengthbuf[1:(end - 2)]
                    end
                end
            end
        end
    end
    #finalize length
    lengthstr = String(lengthbuf)
    duration = isempty(lengthstr) ? nothing : tryparse(Int, lengthstr)
    return (pitch, duration, sustain, dotted)
end

parsevoice(
    """
f# f# g a a g f# e d d e f# f#~ f#8 e e2
f#4 f# g a a g f# e d d e f# e~ e8 d d2
e4 e f# d e f#8~ g8 f#4 d e f#8~ g f#4 e d e a,
f#2 f#4 g a a g f# e d d e f# e~ e8 d8 d2""",
    lyrics = """
      Freu- de, schö- ner Göt- ter- fun- ken, Toch- ter aus E- li- - si- um!
      Wir be- tre- ten feu- er- trun- ken, Himm- li- sche, dein Hei- - lig- thum!
      Dei- ne Zau- ber bin den - wie- der, was die - Mo- de streng ge- theilt,
      al- le mensch- en wer- den Brü- der wo dein sanf- ter Flü- - gel weilt.
      """,
)

# And now with harmony!

soprano = @spawn parsevoice(
    """
f'#. f'#. g'. a'. a'. g'. f'#. e'~ e'8 d.'4 d.' e.' f#'. f#'.~ f#' e'8 e'4~ e'2
""",
    lyrics = "Freu- de, schö- ner Göt- ter- fun- ken, Toch- ter aus E- li- - si- um!",
)
alto = @spawn parsevoice("""
a. a. a. a.  a.  a. a. a~ g8 f#.4 a.  a.  a. a.~ a a8 a4~ a2
""")
tenor = @spawn parsevoice("""
d. d. e. f#. f#. e. d. d~ e8 f#.4 f#. a,. d. d.~ d c#8 c#4 c#2
""")
bass = @spawn parsevoice("""
d. d. d. d. a,. a,. a,. b,~ c8 d. a., a., a., a., a, a8, a,4 a,2
""")
wait(soprano)
wait(alto)
wait(tenor)
wait(bass)

soprano = @spawn parsevoice(
    """
f'#.4 f'#. g'. a'. a'. g'. f'#. e'. d'. d'. e'. f'#. e'.~ e' d'8 d'4~ d'2
""",
    lyrics = "Wir be- tre- ten feu- er- trun- ken, Himm- li- sche, dein Hei- - lig- thum!",
)
alto = @spawn parsevoice("""
a.4 a. b. c'. c'. b. a. g. f#. f#. g. f#. g.~ g4 f#8 f#~ f#2
""")
tenor = @spawn parsevoice("""
d.4 d. d. d. d. d. d. d. d. d. c#. d. c#.~ c# d8 d d2
""")
bass = @spawn parsevoice("""
d.4 d. d. d. a,. a,. a,. a., a., a., a., a., a.,~ a, a,8 d, d,2
""")
wait(soprano)
wait(alto)
wait(tenor)
wait(bass)


text:
function cutensornetGetCudartVersion()
    ccall((:cutensornetGetCudartVersion, libcuquantum), Csize_t, ())
end

@checked function cutensornetContraction(handle, plan, rawDataIn, rawDataOut, workspace, workspaceSize, sliceId, stream)
    initialize_context()
    ccall((:cutensornetContraction, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionPlan_t, Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{Cvoid}, UInt64, Int64, CUstream), handle, plan, rawDataIn, rawDataOut, workspace, workspaceSize, sliceId, stream)
end

@checked function cutensornetCreateContractionOptimizerConfig(handle, optimizerConfig)
    initialize_context()
    ccall((:cutensornetCreateContractionOptimizerConfig, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, Ptr{cutensornetContractionOptimizerConfig_t}), handle, optimizerConfig)
end

@checked function cutensornetCreate(handle)
    initialize_context()
    ccall((:cutensornetCreate, libcuquantum), cutensornetStatus_t, (Ptr{cutensornetHandle_t},), handle)
end

@checked function cutensornetContractionGetWorkspaceSize(handle, descNet, optimizerInfo, workspaceSize)
    initialize_context()
    ccall((:cutensornetContractionGetWorkspaceSize, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetNetworkDescriptor_t, cutensornetContractionOptimizerInfo_t, Ptr{UInt64}), handle, descNet, optimizerInfo, workspaceSize)
end

@checked function cutensornetDestroyContractionOptimizerInfo(optimizerInfo)
    initialize_context()
    ccall((:cutensornetDestroyContractionOptimizerInfo, libcuquantum), cutensornetStatus_t, (cutensornetContractionOptimizerInfo_t,), optimizerInfo)
end

@checked function cutensornetLoggerOpenFile(logFile)
    initialize_context()
    ccall((:cutensornetLoggerOpenFile, libcuquantum), cutensornetStatus_t, (Cstring,), logFile)
end

@checked function cutensornetContractionAutotunePreferenceSetAttribute(handle, autotunePreference, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionAutotunePreferenceSetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionAutotunePreference_t, cutensornetContractionAutotunePreferenceAttributes_t, Ptr{Cvoid}, Csize_t), handle, autotunePreference, attr, buf, sizeInBytes)
end

@checked function cutensornetLoggerSetCallback(callback)
    initialize_context()
    ccall((:cutensornetLoggerSetCallback, libcuquantum), cutensornetStatus_t, (cutensornetLoggerCallback_t,), callback)
end

@checked function cutensornetLoggerForceDisable()
    initialize_context()
    ccall((:cutensornetLoggerForceDisable, libcuquantum), cutensornetStatus_t, ())
end

@checked function cutensornetLoggerSetFile(file)
    initialize_context()
    ccall((:cutensornetLoggerSetFile, libcuquantum), cutensornetStatus_t, (Ptr{FILE},), file)
end

@checked function cutensornetContractionOptimize(handle, descNet, optimizerConfig, workspaceSizeConstraint, optimizerInfo)
    initialize_context()
    ccall((:cutensornetContractionOptimize, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetNetworkDescriptor_t, cutensornetContractionOptimizerConfig_t, UInt64, cutensornetContractionOptimizerInfo_t), handle, descNet, optimizerConfig, workspaceSizeConstraint, optimizerInfo)
end

@checked function cutensornetCreateContractionAutotunePreference(handle, autotunePreference)
    initialize_context()
    ccall((:cutensornetCreateContractionAutotunePreference, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, Ptr{cutensornetContractionAutotunePreference_t}), handle, autotunePreference)
end

@checked function cutensornetCreateNetworkDescriptor(handle, numInputs, numModesIn, extentsIn, stridesIn, modesIn, alignmentRequirementsIn, numModesOut, extentsOut, stridesOut, modesOut, alignmentRequirementsOut, dataType, computeType, descNet)
    initialize_context()
    ccall((:cutensornetCreateNetworkDescriptor, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, Int32, Ptr{Int32}, Ptr{Ptr{Int64}}, Ptr{Ptr{Int64}}, Ptr{Ptr{Int32}}, Ptr{UInt32}, Int32, Ptr{Int64}, Ptr{Int64}, Ptr{Int32}, UInt32, cudaDataType_t, cutensornetComputeType_t, Ptr{cutensornetNetworkDescriptor_t}), handle, numInputs, numModesIn, extentsIn, stridesIn, modesIn, alignmentRequirementsIn, numModesOut, extentsOut, stridesOut, modesOut, alignmentRequirementsOut, dataType, computeType, descNet)
end

@checked function cutensornetContractionOptimizerConfigGetAttribute(handle, optimizerConfig, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionOptimizerConfigGetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionOptimizerConfig_t, cutensornetContractionOptimizerConfigAttributes_t, Ptr{Cvoid}, Csize_t), handle, optimizerConfig, attr, buf, sizeInBytes)
end

@checked function cutensornetContractionOptimizerInfoGetAttribute(handle, optimizerInfo, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionOptimizerInfoGetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionOptimizerInfo_t, cutensornetContractionOptimizerInfoAttributes_t, Ptr{Cvoid}, Csize_t), handle, optimizerInfo, attr, buf, sizeInBytes)
end

@checked function cutensornetContractionOptimizerInfoSetAttribute(handle, optimizerInfo, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionOptimizerInfoSetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionOptimizerInfo_t, cutensornetContractionOptimizerInfoAttributes_t, Ptr{Cvoid}, Csize_t), handle, optimizerInfo, attr, buf, sizeInBytes)
end

@checked function cutensornetDestroyContractionPlan(plan)
    initialize_context()
    ccall((:cutensornetDestroyContractionPlan, libcuquantum), cutensornetStatus_t, (cutensornetContractionPlan_t,), plan)
end

@checked function cutensornetLoggerSetLevel(level)
    initialize_context()
    ccall((:cutensornetLoggerSetLevel, libcuquantum), cutensornetStatus_t, (Int32,), level)
end

@checked function cutensornetDestroy(handle)
    initialize_context()
    ccall((:cutensornetDestroy, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t,), handle)
end

@checked function cutensornetDestroyNetworkDescriptor(desc)
    initialize_context()
    ccall((:cutensornetDestroyNetworkDescriptor, libcuquantum), cutensornetStatus_t, (cutensornetNetworkDescriptor_t,), desc)
end

@checked function cutensornetContractionAutotune(handle, plan, rawDataIn, rawDataOut, workspace, workspaceSize, pref, stream)
    initialize_context()
    ccall((:cutensornetContractionAutotune, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionPlan_t, Ptr{Ptr{Cvoid}}, Ptr{Cvoid}, Ptr{Cvoid}, UInt64, cutensornetContractionAutotunePreference_t, CUstream), handle, plan, rawDataIn, rawDataOut, workspace, workspaceSize, pref, stream)
end

@checked function cutensornetDestroyContractionAutotunePreference(autotunePreference)
    initialize_context()
    ccall((:cutensornetDestroyContractionAutotunePreference, libcuquantum), cutensornetStatus_t, (cutensornetContractionAutotunePreference_t,), autotunePreference)
end

@checked function cutensornetCreateContractionOptimizerInfo(handle, descNet, optimizerInfo)
    initialize_context()
    ccall((:cutensornetCreateContractionOptimizerInfo, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetNetworkDescriptor_t, Ptr{cutensornetContractionOptimizerInfo_t}), handle, descNet, optimizerInfo)
end

function cutensornetGetVersion()
    ccall((:cutensornetGetVersion, libcuquantum), Csize_t, ())
end

function cutensornetGetErrorString(error)
    ccall((:cutensornetGetErrorString, libcuquantum), Cstring, (cutensornetStatus_t,), error)
end

@checked function cutensornetDestroyContractionOptimizerConfig(optimizerConfig)
    initialize_context()
    ccall((:cutensornetDestroyContractionOptimizerConfig, libcuquantum), cutensornetStatus_t, (cutensornetContractionOptimizerConfig_t,), optimizerConfig)
end

@checked function cutensornetContractionOptimizerConfigSetAttribute(handle, optimizerConfig, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionOptimizerConfigSetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionOptimizerConfig_t, cutensornetContractionOptimizerConfigAttributes_t, Ptr{Cvoid}, Csize_t), handle, optimizerConfig, attr, buf, sizeInBytes)
end

@checked function cutensornetLoggerSetMask(mask)
    initialize_context()
    ccall((:cutensornetLoggerSetMask, libcuquantum), cutensornetStatus_t, (Int32,), mask)
end

@checked function cutensornetContractionAutotunePreferenceGetAttribute(handle, autotunePreference, attr, buf, sizeInBytes)
    initialize_context()
    ccall((:cutensornetContractionAutotunePreferenceGetAttribute, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetContractionAutotunePreference_t, cutensornetContractionAutotunePreferenceAttributes_t, Ptr{Cvoid}, Csize_t), handle, autotunePreference, attr, buf, sizeInBytes)
end

@checked function cutensornetCreateContractionPlan(handle, descNet, optimizerInfo, workspaceSize, plan)
    initialize_context()
    ccall((:cutensornetCreateContractionPlan, libcuquantum), cutensornetStatus_t, (cutensornetHandle_t, cutensornetNetworkDescriptor_t, cutensornetContractionOptimizerInfo_t, UInt64, Ptr{cutensornetContractionPlan_t}), handle, descNet, optimizerInfo, workspaceSize, plan)
end

function custatevecApplyGeneralizedPermutationMatrix(handle, sv, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, adjoint, basisBits, nBasisBits, maskBitString, maskOrdering, maskLen, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyGeneralizedPermutationMatrix, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecIndex_t}, Ptr{Cvoid}, cudaDataType_t, Int32, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, adjoint, basisBits, nBasisBits, maskBitString, maskOrdering, maskLen, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecAccessor_create(handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_create, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecAccessorDescriptor_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecExpectationsOnPauliBasis(handle, sv, svDataType, nIndexBits, expectationValues, pauliOperatorsArray, basisBitsArray, nBasisBitsArray, nPauliOperatorArrays)
    initialize_context()
    ccall((:custatevecExpectationsOnPauliBasis, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Ptr{custatevecPauli_t}}, Ptr{Ptr{Int32}}, Ptr{UInt32}, UInt32), handle, sv, svDataType, nIndexBits, expectationValues, pauliOperatorsArray, basisBitsArray, nBasisBitsArray, nPauliOperatorArrays)
end

function custatevecSampler_create(handle, sv, svDataType, nIndexBits, sampler, nMaxShots, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSampler_create, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecSamplerDescriptor_t}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, sampler, nMaxShots, extraWorkspaceSizeInBytes)
end

function custatevecExpectation(handle, sv, svDataType, nIndexBits, expectationValue, expectationDataType, residualNorm, matrix, matrixDataType, layout, basisBits, nBasisBits, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecExpectation, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, Ptr{Cdouble}, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Ptr{Int32}, UInt32, custatevecComputeType_t, Ptr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, expectationValue, expectationDataType, residualNorm, matrix, matrixDataType, layout, basisBits, nBasisBits, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecBatchMeasure(handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, randnum, collapse)
    initialize_context()
    ccall((:custatevecBatchMeasure, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble, custatevecCollapseOp_t), handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, randnum, collapse)
end

function custatevecExpectation_bufferSize(handle, svDataType, nIndexBits, matrix, matrixDataType, layout, nBasisBits, computeType, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecExpectation_bufferSize, libcuquantum), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, UInt32, custatevecComputeType_t, Ptr{Csize_t}), handle, svDataType, nIndexBits, matrix, matrixDataType, layout, nBasisBits, computeType, extraWorkspaceSizeInBytes)
end

function custatevecApplyMatrix_bufferSize(handle, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, nTargets, nControls, computeType, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyMatrix_bufferSize, libcuquantum), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Int32, UInt32, UInt32, custatevecComputeType_t, Ptr{Csize_t}), handle, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, nTargets, nControls, computeType, extraWorkspaceSizeInBytes)
end

function custatevecAccessor_createReadOnly(handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_createReadOnly, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{custatevecAccessorDescriptor_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Ptr{Csize_t}), handle, sv, svDataType, nIndexBits, accessor, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecGetStream(handle, streamId)
    initialize_context()
    ccall((:custatevecGetStream, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{CUstream}), handle, streamId)
end

function custatevecAccessor_get(handle, accessor, externalBuffer, _begin, _end)
    initialize_context()
    ccall((:custatevecAccessor_get, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, custatevecIndex_t, custatevecIndex_t), handle, accessor, externalBuffer, _begin, _end)
end

function custatevecGetVersion()
    ccall((:custatevecGetVersion, libcuquantum), Csize_t, ())
end

function custatevecLoggerForceDisable()
    ccall((:custatevecLoggerForceDisable, libcuquantum), custatevecStatus_t, ())
end

function custatevecSetWorkspace(handle, workspace, workspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSetWorkspace, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, Csize_t), handle, workspace, workspaceSizeInBytes)
end

function custatevecApplyExp(handle, sv, svDataType, nIndexBits, theta, paulis, targets, nTargets, controls, controlBitValues, nControls)
    initialize_context()
    ccall((:custatevecApplyExp, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Cdouble, Ptr{custatevecPauli_t}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, theta, paulis, targets, nTargets, controls, controlBitValues, nControls)
end

function custatevecGetProperty(type, value)
    initialize_context()
    ccall((:custatevecGetProperty, libcuquantum), custatevecStatus_t, (libraryPropertyType, Ptr{Int32}), type, value)
end

function custatevecAbs2SumOnZBasis(handle, sv, svDataType, nIndexBits, abs2sum0, abs2sum1, basisBits, nBasisBits)
    initialize_context()
    ccall((:custatevecAbs2SumOnZBasis, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, abs2sum0, abs2sum1, basisBits, nBasisBits)
end

function custatevecDestroy(handle)
    initialize_context()
    ccall((:custatevecDestroy, libcuquantum), custatevecStatus_t, (custatevecHandle_t,), handle)
end

function custatevecLoggerSetLevel(level)
    initialize_context()
    ccall((:custatevecLoggerSetLevel, libcuquantum), custatevecStatus_t, (Int32,), level)
end

function custatevecApplyGeneralizedPermutationMatrix_bufferSize(handle, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, basisBits, nBasisBits, maskLen, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyGeneralizedPermutationMatrix_bufferSize, libcuquantum), custatevecStatus_t, (custatevecHandle_t, cudaDataType_t, UInt32, Ptr{custatevecIndex_t}, Ptr{Cvoid}, cudaDataType_t, Ptr{Int32}, UInt32, UInt32, Ptr{Csize_t}), handle, svDataType, nIndexBits, permutation, diagonals, diagonalsDataType, basisBits, nBasisBits, maskLen, extraWorkspaceSizeInBytes)
end

function custatevecCollapseByBitString(handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, norm)
    initialize_context()
    ccall((:custatevecCollapseByBitString, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble), handle, sv, svDataType, nIndexBits, bitString, bitOrdering, bitStringLen, norm)
end

function custatevecSampler_sample(handle, sampler, bitStrings, bitOrdering, bitStringLen, randnums, nShots, output)
    initialize_context()
    ccall((:custatevecSampler_sample, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecSamplerDescriptor_t}, Ptr{custatevecIndex_t}, Ptr{Int32}, UInt32, Ptr{Cdouble}, UInt32, custatevecSamplerOutput_t), handle, sampler, bitStrings, bitOrdering, bitStringLen, randnums, nShots, output)
end

function custatevecGetErrorName(status)
    initialize_context()
    ccall((:custatevecGetErrorName, libcuquantum), Cstring, (custatevecStatus_t,), status)
end

function custatevecCreate(handle)
    initialize_context()
    ccall((:custatevecCreate, libcuquantum), custatevecStatus_t, (Ptr{custatevecHandle_t},), handle)
end

function custatevecAccessor_setExtraWorkspace(handle, accessor, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecAccessor_setExtraWorkspace, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, Csize_t), handle, accessor, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecApplyMatrix(handle, sv, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, targets, nTargets, controls, nControls, controlBitValues, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecApplyMatrix, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cvoid}, cudaDataType_t, custatevecMatrixLayout_t, Int32, Ptr{Int32}, UInt32, Ptr{Int32}, UInt32, Ptr{Int32}, custatevecComputeType_t, Ptr{Cvoid}, Csize_t), handle, sv, svDataType, nIndexBits, matrix, matrixDataType, layout, adjoint, targets, nTargets, controls, nControls, controlBitValues, computeType, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecCollapseOnZBasis(handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, norm)
    initialize_context()
    ccall((:custatevecCollapseOnZBasis, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Int32, Ptr{Int32}, UInt32, Cdouble), handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, norm)
end

function custatevecLoggerSetMask(mask)
    ccall((:custatevecLoggerSetMask, libcuquantum), custatevecStatus_t, (Int32,), mask)
end

function custatevecSetStream(handle, streamId)
    initialize_context()
    ccall((:custatevecSetStream, libcuquantum), custatevecStatus_t, (custatevecHandle_t, CUstream), handle, streamId)
end

function custatevecGetDefaultWorkspaceSize(handle, workspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecGetDefaultWorkspaceSize, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Csize_t}), handle, workspaceSizeInBytes)
end

function custatevecAbs2SumArray(handle, sv, svDataType, nIndexBits, abs2sum, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen)
    initialize_context()
    ccall((:custatevecAbs2SumArray, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Cdouble}, Ptr{Int32}, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32), handle, sv, svDataType, nIndexBits, abs2sum, bitOrdering, bitOrderingLen, maskBitString, maskOrdering, maskLen)
end

function custatevecLoggerOpenFile(logFile)
    ccall((:custatevecLoggerOpenFile, libcuquantum), custatevecStatus_t, (Cstring,), logFile)
end

function custatevecLoggerSetFile(file)
    ccall((:custatevecLoggerSetFile, libcuquantum), custatevecStatus_t, (Ptr{FILE},), file)
end

function custatevecGetErrorString(status)
    ccall((:custatevecGetErrorString, libcuquantum), Cstring, (custatevecStatus_t,), status)
end

function custatevecLoggerSetCallback(callback)
    ccall((:custatevecLoggerSetCallback, libcuquantum), custatevecStatus_t, (custatevecLoggerCallback_t,), callback)
end

function custatevecSampler_preprocess(handle, sampler, extraWorkspace, extraWorkspaceSizeInBytes)
    initialize_context()
    ccall((:custatevecSampler_preprocess, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecSamplerDescriptor_t}, Ptr{Cvoid}, Csize_t), handle, sampler, extraWorkspace, extraWorkspaceSizeInBytes)
end

function custatevecMeasureOnZBasis(handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, randnum, collapse)
    initialize_context()
    ccall((:custatevecMeasureOnZBasis, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{Cvoid}, cudaDataType_t, UInt32, Ptr{Int32}, Ptr{Int32}, UInt32, Cdouble, custatevecCollapseOp_t), handle, sv, svDataType, nIndexBits, parity, basisBits, nBasisBits, randnum, collapse)
end

function custatevecAccessor_set(handle, accessor, externalBuffer, _begin, _end)
    initialize_context()
    ccall((:custatevecAccessor_set, libcuquantum), custatevecStatus_t, (custatevecHandle_t, Ptr{custatevecAccessorDescriptor_t}, Ptr{Cvoid}, custatevecIndex_t, custatevecIndex_t), handle, accessor, externalBuffer, _begin, _end)
end


text:
module TestDoubleDicts

using Test

using VecMathOptInterface
const MOI = VecMathOptInterface

const DoubleDicts = VecMathOptInterface.Utilities.DoubleDicts

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function _test_iterator(dict)
    kk, vv = [], []
    for (k, v) in dict
        push!(kk, k)
        push!(vv, v)
    end
    @test length(kk) == length(keys(dict))
    @test isempty(setdiff(kk, keys(dict)))
    @test length(vv) == length(values(dict))
    @test isempty(setdiff(vv, values(dict)))
    @test dict == Dict(kk .=> vv)
    return
end

function _test_basic_functionality(dict, k_values, v_values)
    @test_throws ErrorException sizehint!(dict, 1)

    @test isempty(dict)
    @test length(dict) == 0
    _test_iterator(dict)

    dict[k_values[1]] = v_values[1]
    delete!(dict, k_values[1])
    _test_iterator(dict)

    dict[k_values[3]] = v_values[3]
    _test_iterator(dict)

    empty!(dict)
    @test isempty(dict)
    @test !haskey(dict, k_values[1])

    dict[k_values[1]] = v_values[1]
    @test haskey(dict, k_values[1])
    @test values(dict) == [v_values[1]]
    @test keys(dict) == [k_values[1]]
    @test dict[k_values[1]] == v_values[1]
    @test get(dict, k_values[1], nothing) == v_values[1]
    @test get(dict, k_values[1], v_values[2]) == v_values[1]
    for i in eachindex(k_values)
        if i != 1
            @test get(dict, k_values[i], nothing) === nothing
            # Test that the implementation does not only work with `nothing`
            @test get(dict, k_values[i], v_values[i]) == v_values[i]
        end
    end

    _test_iterator(dict)

    delete!(dict, k_values[1])
    @test !haskey(dict, k_values[1])
    for i in eachindex(k_values)
        @test get(dict, k_values[i], nothing) === nothing
        # Test that the implementation does not only work with `nothing`
        @test get(dict, k_values[i], v_values[i]) == v_values[i]
    end
    @test isempty(dict)

    dict[k_values[1]] = v_values[1]
    @test length(dict) == 1

    for (k, v) in zip(k_values, v_values)
        dict[k] = v
    end

    _test_iterator(dict)

    empty!(dict)
    @test isempty(dict)
    @test length(dict) == 0

    dict[k_values[1]] = v_values[1]
    idict = dict[MOI.VariableIndex, MOI.Integer]
    bdict = dict[MOI.VariableIndex, MOI.ZeroOne]
    idict_ = dict[MOI.VariableIndex, MOI.Integer]
    @test idict.dict === idict_.dict
    sizehint!(idict, 2)
    @test length(idict) == 1
    @test length(bdict) == 0
    @test values(bdict) == typeof(v_values[3])[]
    @test keys(bdict) == typeof(k_values[3])[]

    @test haskey(idict, k_values[1])
    @test !haskey(idict, k_values[2])
    @test !haskey(bdict, k_values[3])

    @test values(idict) == [v_values[1]]
    @test keys(idict) == [k_values[1]]
    @test idict[k_values[1]] == v_values[1]

    idict[k_values[2]] = v_values[2]
    @test haskey(idict, k_values[2])
    _test_iterator(idict)

    delete!(idict, k_values[2])
    @test !haskey(idict, k_values[2])
    _test_iterator(idict)

    @test !isempty(idict)
    @test isempty(bdict)
    empty!(idict)
    empty!(bdict)
    @test isempty(idict)
    @test isempty(bdict)

    bdict[k_values[3]] = v_values[3]
    length(bdict) == 1

    edict = dict[MOI.VariableIndex, MOI.EqualTo{Bool}]
    ek = MOI.ConstraintIndex{MOI.VariableIndex,MOI.EqualTo{Bool}}(1)
    delete!(edict, ek)
    @test_throws KeyError edict[ek]
    sizehint!(edict, 0)
    @test length(edict) == 0
    @test_throws KeyError edict[ek]
    delete!(edict, ek)
    _test_iterator(edict)
    return
end

function test_DoubleDict()
    dict = DoubleDicts.DoubleDict{Float64}()
    keys = [
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(1),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(2),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(1),
    ]
    vals = [1.0, 2.0, 1.0]
    _test_basic_functionality(dict, keys, vals)
    return
end

function test_IndexDoubleDict()
    dict = DoubleDicts.IndexDoubleDict()
    keys = [
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(1),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.Integer}(2),
        MOI.ConstraintIndex{MOI.VariableIndex,MOI.ZeroOne}(1),
    ]
    vals = keys
    _test_basic_functionality(dict, keys, vals)
    src = DoubleDicts.IndexDoubleDict()
    for (k, v) in zip(keys, vals)
        dict[k] = v
    end
    dest = DoubleDicts.IndexDoubleDict()
    MOI.Utilities._reverse_dict(dest, src)
    for (k, v) in src
        @test dest[v] == k
    end
    return
end

end  # module

TestDoubleDicts.runtests()


text:
using Genie, Genie.Renderer.Html, Stipple, StipplePlotly

pd(name) = PlotData(
  x = ["Jan2019", "Feb2019", "Mar2019", "Apr2019", "May2019",
        "Jun2019", "Jul2019", "Aug2019", "Sep2019", "Oct2019",
        "Nov2019", "Dec2019"],
  y = Int[rand(1:100_000) for x in 1:12],
  plot = StipplePlotly.Charts.PLOT_TYPE_SCATTER,
  name = name
)

Base.@kwdef mutable struct Model <: ReactiveModel
  data::R{Vector{PlotData}} = [pd("Random 1"),pd("Random 2")]
  layout::R{PlotLayout} = PlotLayout(
      plot_bgcolor = "#333",
      title = PlotLayoutTitle(text="Random numbers", font=Font(24))
    )
  config::R{PlotConfig} = PlotConfig()
end

model = Stipple.init(Model())

function ui()
  page(
    vm(model), class="container", [
      plot(:data, layout = :layout, config = :config)
    ]
  ) |> html
end

route("/", ui)

up()


text:
@doc raw"""
    adjoint_differential_bezier_control(M::AbstractManifold, b::BezierSegment, t, η)
    adjoint_differential_bezier_control!(
        M::AbstractManifold,
        Y::BezierSegment,
        b::BezierSegment,
        t,
        η,
    )

evaluate the adjoint of the differential of a Bézier curve on the manifold `M`
with respect to its control points `b` based on a point `t```∈[0,1]`` on the
curve and a tangent vector ``η∈T_{β(t)}\mathcal M``.
This can be computed in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(M::AbstractManifold, b::BezierSegment, t, η)
    n = length(b.pts)
    if n == 2
        return BezierSegment([
            adjoint_differential_geodesic_startpoint(M, b.pts[1], b.pts[2], t, η),
            adjoint_differential_geodesic_endpoint(M, b.pts[1], b.pts[2], t, η),
        ])
    end
    c = [b.pts, [similar.(b.pts[1:l]) for l in (n - 1):-1:2]...]
    for i in 2:(n - 1) # casteljau on the tree -> forward with interims storage
        c[i] .= shortest_geodesic.(Ref(M), c[i - 1][1:(end - 1)], c[i - 1][2:end], Ref(t))
    end
    Y = [η, [similar(η) for i in 1:(n - 1)]...]
    for i in (n - 1):-1:1 # propagate adjoints -> backward without interims storage
        Y[1:(n - i + 1)] .=
            [ # take previous results and add start&end point effects
                adjoint_differential_geodesic_startpoint.(
                    Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y[1:(n - i)]
                )...,
                zero_vector(M, c[i][end]),
            ] .+ [
                zero_vector(M, c[i][1]),
                adjoint_differential_geodesic_endpoint.(
                    Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y[1:(n - i)]
                )...,
            ]
    end
    return BezierSegment(Y)
end
function adjoint_differential_bezier_control!(
    M::AbstractManifold, Y::BezierSegment, b::BezierSegment, t, η
)
    n = length(b.pts)
    if n == 2
        adjoint_differential_geodesic_startpoint!(M, Y.pts[1], b.pts[1], b.pts[2], t, η)
        adjoint_differential_geodesic_endpoint!(M, Y.pts[2], b.pts[1], b.pts[2], t, η)
        return Y
    end
    c = [b.pts, [similar.(b.pts[1:l]) for l in (n - 1):-1:2]...]
    for i in 2:(n - 1) # casteljau on the tree -> forward with interims storage
        c[i] .= shortest_geodesic.(Ref(M), c[i - 1][1:(end - 1)], c[i - 1][2:end], Ref(t))
    end
    Y.pts[1] = η
    for i in (n - 1):-1:1 # propagate adjoints -> backward without interims storage
        Y.pts[1:(n - i + 1)] .=
            [ # take previous results and add start&end point effects
                adjoint_differential_geodesic_startpoint.(
                    Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y.pts[1:(n - i)]
                )...,
                zero_vector(M, c[i][end]),
            ] .+ [
                zero_vector(M, c[i][1]),
                adjoint_differential_geodesic_endpoint.(
                    Ref(M), c[i][1:(end - 1)], c[i][2:end], Ref(t), Y.pts[1:(n - i)]
                )...,
            ]
    end
    return Y
end

@doc raw"""
    adjoint_differential_bezier_control(
        M::AbstractManifold,
        b::BezierSegment,
        t::AbstractVector,
        X::AbstractVector,
    )
    adjoint_differential_bezier_control!(
        M::AbstractManifold,
        Y::BezierSegment,
        b::BezierSegment,
        t::AbstractVector,
        X::AbstractVector,
    )
evaluate the adjoint of the differential of a Bézier curve on the manifold `M`
with respect to its control points `b` based on a points `T```=(t_i)_{i=1}^n`` that
are pointwise in `` t_i∈[0,1]`` on the curve and given corresponding tangential
vectors ``X = (η_i)_{i=1}^n``, ``η_i∈T_{β(t_i)}\mathcal M``
This can be computed in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve and[^BergmannGousenbourger2018].

[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
    > by minimizing the acceleration of a Bézier curve_.
    > Frontiers in Applied Mathematics and Statistics, 2018.
    > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
"""
function adjoint_differential_bezier_control(
    M::AbstractManifold, b::BezierSegment, t::AbstractVector, X::AbstractVector
)
    effects = [bt.pts for bt in adjoint_differential_bezier_control.(Ref(M), Ref(b), t, X)]
    return BezierSegment(sum(effects))
end
function adjoint_differential_bezier_control!(
    M::AbstractManifold,
    Y::BezierSegment,
    b::BezierSegment,
    t::AbstractVector,
    X::AbstractVector,
)
    Z = BezierSegment(similar.(Y.pts))
    fill!.(Y.pts, zero(eltype(first(Y.pts))))
    for i in 1:length(t)
        adjoint_differential_bezier_control!(M, Z, b, t[i], X[i])
        Y.pts .+= Z.pts
    end
    return Y
end

@doc raw"""
    adjoint_differential_bezier_control(
        M::AbstractManifold,
        B::AbstractVector{<:BezierSegment},
        t,
        X
    )
    adjoint_differential_bezier_control!(
        M::AbstractManifold,
        Y::AbstractVector{<:BezierSegment},
        B::AbstractVector{<:BezierSegment},
        t,
        X
    )

evaluate the adjoint of the differential of a composite Bézier curve on the
manifold `M` with respect to its control points `b` based on a points `T```=(t_i)_{i=1}^n``
that are pointwise in ``t_i∈[0,1]`` on the curve and given corresponding tangential
vectors ``X = (η_i)_{i=1}^n``, ``η_i∈T_{β(t_i)}\mathcal M``
This can be computed in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(
    M::AbstractManifold, B::AbstractVector{<:BezierSegment}, t, X
)
    Y = broadcast(b -> BezierSegment(zero_vector.(Ref(M), b.pts)), B) # Double broadcast
    return adjoint_differential_bezier_control!(M, Y, B, t, X)
end
function adjoint_differential_bezier_control!(
    M::AbstractManifold,
    Y::AbstractVector{<:BezierSegment},
    B::AbstractVector{<:BezierSegment},
    t,
    X,
)
    # doubly nested broadbast on the Array(Array) of CPs (note broadcast _and_ .)
    if (0 > t) || (t > length(B))
        error(
            "The parameter ",
            t,
            " to evaluate the composite Bézier curve at is outside the interval [0,",
            length(B),
            "].",
        )
    end
    for y in Y
        fill!.(y.pts, zero(eltype(first(y.pts))))
    end
    seg = max(ceil(Int, t), 1)
    localT = ceil(Int, t) == 0 ? 0.0 : t - seg + 1
    adjoint_differential_bezier_control!(M, Y[seg], B[seg], localT, X)
    return Y
end
@doc raw"""
    adjoint_differential_bezier_control(
        M::AbstractManifold,
        T::AbstractVector,
        X::AbstractVector,
    )
    adjoint_differential_bezier_control!(
        M::AbstractManifold,
        Y::AbstractVector{<:BezierSegment},
        T::AbstractVector,
        X::AbstractVector,
    )

Evaluate the adjoint of the differential with respect to the controlpoints at several times `T`.
This can be computed in place of `Y`.

See [`de_casteljau`](@ref) for more details on the curve.
"""
function adjoint_differential_bezier_control(
    M::AbstractManifold,
    B::AbstractVector{<:BezierSegment},
    T::AbstractVector,
    X::AbstractVector,
)
    Y = broadcast(b -> BezierSegment(zero_vector.(Ref(M), b.pts)), B) # Double broadcast
    return adjoint_differential_bezier_control!(M, Y, B, T, X)
end
function adjoint_differential_bezier_control!(
    M::AbstractManifold,
    Y::AbstractVector{<:BezierSegment},
    B::AbstractVector{<:BezierSegment},
    T::AbstractVector,
    X::AbstractVector,
)
    Z = [BezierSegment(similar.(y.pts)) for y in Y]
    for j in 1:length(T) # for all times
        adjoint_differential_bezier_control!(M, Z, B, T[j], X[j])
        for i in 1:length(Z)
            Y[i].pts .+= Z[i].pts
        end
    end
    return Y
end
@doc raw"""
    adjoint_differential_geodesic_startpoint(M,p, q, t, X)
    adjoint_differential_geodesic_startpoint!(M, Y, p, q, t, X)

Compute the adjoint of ``D_p γ(t; p, q)[X]`` (in place of `Y`).

# See also

[`differential_geodesic_startpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_geodesic_startpoint(M::AbstractManifold, p, q, t, X)
    return adjoint_Jacobi_field(M, p, q, t, X, βdifferential_geodesic_startpoint)
end
function adjoint_differential_geodesic_startpoint!(M::AbstractManifold, Y, p, q, t, X)
    return adjoint_Jacobi_field!(M, Y, p, q, t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    adjoint_differential_geodesic_endpoint(M, p, q, t, X)
    adjoint_differential_geodesic_endpoint!(M, Y, p, q, t, X)

Compute the adjoint of ``D_q γ(t; p, q)[X]`` (in place of `Y`).

# See also

[`differential_geodesic_endpoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_geodesic_endpoint(M::AbstractManifold, p, q, t, X)
    return adjoint_Jacobi_field(M, q, p, 1 - t, X, βdifferential_geodesic_startpoint)
end
function adjoint_differential_geodesic_endpoint!(M::AbstractManifold, Y, p, q, t, X)
    return adjoint_Jacobi_field!(M, Y, q, p, 1 - t, X, βdifferential_geodesic_startpoint)
end

@doc raw"""
    adjoint_differential_exp_basepoint(M, p, X, Y)
    adjoint_differential_exp_basepoint!(M, Z, p, X, Y)

Computes the adjoint of ``D_p \exp_p X[Y]`` (in place of `Z`).

# See also

[`differential_exp_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_basepoint(M::AbstractManifold, p, X, Y)
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_basepoint)
end
function adjoint_differential_exp_basepoint!(M::AbstractManifold, Z, p, X, Y)
    return adjoint_Jacobi_field!(M, Z, p, exp(M, p, X), 1.0, Y, βdifferential_exp_basepoint)
end

@doc raw"""
    adjoint_differential_exp_argument(M, p, X, Y)
    adjoint_differential_exp_argument!(M, Z, p, X, Y)

Compute the adjoint of ``D_X\exp_p X[Y]`` (in place of `Z`).
Note that ``X ∈  T_p(T_p\mathcal M) = T_p\mathcal M`` is still a tangent vector.

# See also

[`differential_exp_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_exp_argument(M::AbstractManifold, p, X, Y)
    return adjoint_Jacobi_field(M, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end
function adjoint_differential_exp_argument!(M::AbstractManifold, Z, p, X, Y)
    return adjoint_Jacobi_field!(M, Z, p, exp(M, p, X), 1.0, Y, βdifferential_exp_argument)
end

@doc raw"""
    adjoint_differential_log_basepoint(M, p, q, X)
    adjoint_differential_log_basepoint!(M, Y, p, q, X)

computes the adjoint of ``D_p log_p q[X]`` (in place of `Y`).

# See also
[`differential_log_basepoint`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_basepoint(M::AbstractManifold, p, q, X)
    return adjoint_Jacobi_field(M, p, q, 0.0, X, βdifferential_log_basepoint)
end
function adjoint_differential_log_basepoint!(M::AbstractManifold, Y, p, q, X)
    return adjoint_Jacobi_field!(M, Y, p, q, 0.0, X, βdifferential_log_basepoint)
end

@doc raw"""
    adjoint_differential_log_argument(M, p, q, X)
    adjoint_differential_log_argument!(M, Y, p, q, X)

Compute the adjoint of ``D_q log_p q[X]`` (in place of `Y`).

# See also
[`differential_log_argument`](@ref), [`adjoint_Jacobi_field`](@ref)
"""
function adjoint_differential_log_argument(M::AbstractManifold, p, q, X)
    # order of p and q has to be reversed in this call, cf. Persch, 2018 Lemma 2.3
    return adjoint_Jacobi_field(M, q, p, 1.0, X, βdifferential_log_argument)
end
function adjoint_differential_log_argument!(M::AbstractManifold, Y, p, q, X)
    return adjoint_Jacobi_field!(M, Y, q, p, 1.0, X, βdifferential_log_argument)
end

@doc raw"""
    Y = adjoint_differential_forward_logs(M, p, X)
    adjoint_differential_forward_logs!(M, Y, p, X)

Compute the adjoint differential of [`forward_logs`](@ref) ``F`` orrucirng,
in the power manifold array `p`, the differential of the function

``F_i(p) = \sum_{j ∈ \mathcal I_i} \log_{p_i} p_j``

where ``i`` runs over all indices of the `PowerManifold` manifold `M` and ``\mathcal I_i``
denotes the forward neighbors of ``i``
Let ``n`` be the number dimensions of the `PowerManifold` manifold (i.e. `length(size(x)`)).
Then the input tangent vector lies on the manifold ``\mathcal M' = \mathcal M^n``.
The adjoint differential can be computed in place of `Y`.

# Input

* `M`     – a `PowerManifold` manifold
* `p`     – an array of points on a manifold
* `X`     – a tangent vector to from the n-fold power of `p`, where n is the `ndims` of `p`

# Ouput

`Y` – resulting tangent vector in ``T_p\mathcal M`` representing the adjoint
  differentials of the logs.
"""
function adjoint_differential_forward_logs(
    M::PowerManifold{𝔽,TM,TSize,TPR}, p, X
) where {𝔽,TM,TSize,TPR}
    Y = zero_vector(M, p)
    return adjoint_differential_forward_logs!(M, Y, p, X)
end
function adjoint_differential_forward_logs!(
    M::PowerManifold{𝔽,TM,TSize,TPR}, Y, p, X
) where {𝔽,TM,TSize,TPR}
    power_size = power_dimensions(M)
    d = length(power_size)
    N = PowerManifold(M.manifold, TPR(), power_size..., d)
    R = CartesianIndices(Tuple(power_size))
    maxInd = last(R).I
    for i in R # iterate over all pixel
        for k in 1:d # for all direction combinations
            I = [i.I...] # array of index
            J = I .+ 1 .* (1:d .== k) #i + e_k is j
            if all(J .<= maxInd) # is this neighbor in range?
                j = CartesianIndex{d}(J...) # neigbbor index as Cartesian Index
                Y[M, I...] =
                    Y[M, I...] + adjoint_differential_log_basepoint(
                        M.manifold, p[M, I...], p[M, J...], X[N, I..., k]
                    )
                Y[M, J...] =
                    Y[M, J...] + adjoint_differential_log_argument(
                        M.manifold, p[M, J...], p[M, I...], X[N, I..., k]
                    )
            end
        end # directions
    end # i in R
    return Y
end


text:
#import Base.print

#=
print : print Assoc in a way that mimics the Sparse Array print.
=#
function print(A::Assoc)
    if !isempty(A)
        r,c,v = find(A)
        [size(A)...; length.(A.col)]

        padR = max(ndigits(size(A)[1]), length.(A.row)...)
        padC = max(ndigits(size(A)[2]), length.(A.col)...)

        println.("  [", rpad.(r, padR), ", ", rpad.(c , padC), "]  =  ", v)
    else
        show(A)
    end
    return nothing
end

#=
printFull : print Assoc in tabular form.
=#
function printFull(A::Assoc)
    if !isempty(A)
        display(full(A))
    else
        show(A)
    end

    return nothing
end

#=
printTriple : return A in triple String form: (r,c) v
=#
function printTriple(A::Assoc)

    if !isempty(A)
        r,c,v = find(A)
        println.("(" .* string.(r) .* "," .* string.(c) .* ")\t".* string.(v))
    else
        show(A)
    end

    return nothing
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
#                    Lauren Milechin (lauren.milechin@mit.edu)
########################################################