text:
using Random
using NamedArrays

#=
#   Write the team names of the teams that finished first in each group in `seeded`.
#   Do this in order of groups from A to H.
#
#   Then write the names of the teams who finished second in `unseeded` similarly.
#
#   Finally, make sure the teams you listed are included in the dictionaries containing
#   the teams from each country (i.e. `Spain`, `France`, etc.)
#
#   `simdraw(n)` generates an array where the j-th row and k-th column shows the number
#   of times in `n` simulations that the j-th seeded team meets the k-th unseeded team.
#   The i-th seeded (unseeded) team is the winner (runner-up) of the i-th group.
#
#   Run these simulations in Julia using
#
#       include("CLDraw.jl"); n = 10000; simdraw(n)
#
#   `simdraw_multhr(n)` is a multithreaded version of `simdraw(n)`. To simulate a single
#   draw, you can use `getdraw(seeded, unseeded, group, assoc)`.
=#
##
seeded = ["Bayern", "Real Madrid", "Man City", "Liverpool",
          "Chelsea", "Dortmund", "Juventus", "PSG"]
shortseed = ["BAY", "RMA", "MCI", "LIV", "CHE", "DOR", "JUV", "PSG"] # short names for seeded
unseeded = ["Atletico", "Monchengladbach", "Porto", "Atalanta",
            "Sevilla", "Lazio", "Barcelona", "RB Leipzig"]
shortuns = ["ATM", "MGB", "POR", "ATL", "SEV", "LAZ", "BAR", "RBL"] # short names for unseeded
groups = ['A', 'B', 'C', 'D', 'E', 'F', 'G', 'H']

Spain = Dict(team => "ESP" for team in ["Real Madrid", "Barcelona", "Atletico", "Valencia", "Sevilla"])
France = Dict(team => "FRA" for team in ["PSG", "Lyon"])
Germany = Dict(team => "GER" for team in ["Bayern", "RB Leipzig", "Dortmund", "Monchengladbach"])
England = Dict(team => "ENG" for team in  ["Man City", "Liverpool", "Tottenham", "Chelsea"])
Italy = Dict(team => "ITA" for team in ["Juventus", "Napoli", "Atalanta", "Lazio"])
Portugal = Dict(team => "POR" for team in ["Porto"])

group = merge(Dict(seeded[i] => groups[i] for i in eachindex(seeded)),
              Dict(unseeded[i] => groups[i] for i in eachindex(unseeded)))
assoc = merge(Spain, France, Germany, England, Italy, Portugal)

##
function matchable(team, opponent, group, assoc)
    # Checks if `team` and `opponent` are from different groups and associations.
    diffgroup = group[team] != group[opponent]
    diffassoc = assoc[team] != assoc[opponent]
    return diffgroup && diffassoc
end

##
function getmatches(team, seeded, unseeded, group, assoc)
    # Returns the possible opponents of `team`
    opponents = team in seeded ? unseeded : seeded
    # Check if `team` is `seeded` or `unseeded`
    matches = [opp for opp in opponents if matchable(team, opp, group, assoc)]
    return matches
end

##
function draw(team, seeded, unseeded, group, assoc)
    # Chooses an opponent randomly from the potential opponents
    possible = getmatches(team, seeded, unseeded, group, assoc)
    return rand(possible)
end

##
function getdraw(seeded, unseeded, group, assoc)
    # Takes the unmatched teams in `seeded` and `unseeded` and generates a random draw
    # according to CL rules.
    #
    # When the unmatched teams impose no restrictions on the remaining matches,
    # an unseeded team is choosen (uniformly) at random. Then an opponent is chosen
    # (uniformly) at random from the eligible seeded teams.
    #
    # If there are two unseeded teams who can only face the same two seeded teams,
    # these two unseeded teams are matched first.
    if length(seeded) == 1 || length(unseeded) == 1
        # Chooses a random match if there is only one seeded or unseeded team remaining.
        # We have no guarantee that the match is legal (different groups / associations).
        # Uncomment the line before `return Dict(home => away)` to generate a warning
        # for illegal matches.
        home = rand(unseeded)
        away = rand(seeded)
        #matchable(home, away, group, assoc) || @warn("Illegal match.")
        return Dict(home => away)
    end
    seedopps = map(x -> getmatches(x, seeded, unseeded, group, assoc), seeded)
    # Create a vector of the possible opponents of unmatched seeded teams.
    unopps = map(x -> getmatches(x, seeded, unseeded, group, assoc), unseeded)
    # Create a vector of possible opponents of unmatched unseeded teams.
    nseedopps = size.(seedopps, 1)
    # Count the possible opponents of each seeded team.
    nunopps = size.(unopps, 1)
    # Count the possible opponents of each unseeded team.
    oneseedopp = findfirst(isequal(1), nseedopps)
    # Checks if a seeded team has only one possible opponent.
    # Returns `nothing` if a seeded team has multiple possible opponents.
    oneunopp = findfirst(isequal(1), nunopps)
    # Checks if an unseeded team has only one possible opponent.
    # Returns `nothing` if an unseeded team has multiple possible opponents.
    if isnothing(oneseedopp) && isnothing(oneunopp)
        # Checks whether each team has multiple possible opponents.
        twounopps = findall(isequal(2), nunopps)
        if !allunique(unopps[twounopps])
            # If two teams have the same two possible opponents,
            # then these two pairs must face each other.
            #twounopp = findfirst(isequal(2), nunopps)
            #home = unseeded[twounopp]
            # Choose a team at random from the pair of unseeded teams.
            home = rand(unseeded[twounopps])
        else
            # Choose one team randomly from `unseeded`.
            home = rand(unseeded)
        end
        away = draw(home, seeded, unseeded, group, assoc)
    elseif isnothing(oneseedopp)
        # An unseeded team has only one possible opponent.
        home = unseeded[oneunopp]
        away = draw(home, seeded, unseeded, group, assoc)
    else
        # A seeded team has only one possible opponent.
        away = seeded[oneseedopp]
        home = draw(away, seeded, unseeded, group, assoc)
    end
    newunseed = filter(!isequal(home), unseeded)
    newseed = filter(!isequal(away), seeded)
    # Removes the matched teams from `seeded` and `unseeded`.
    tie = Dict(home => away)
    return merge(tie, getdraw(newseed, newunseed, group, assoc))
    # Recursively calls itself on the remaining unmatched teams.
end

##
function simdraw(n, seeded=seeded, unseeded=unseeded, group=group, assoc=assoc)
    count = zeros(Int, length(seeded), length(unseeded))
    for ii in 1:n
        #tempseeded = copy(seeded)
        #tempun = copy(unseeded)
        draw = getdraw(seeded, unseeded, group, assoc)
        for jj in eachindex(unseeded)
            count[:,jj] += draw[unseeded[jj]] .== seeded
        end
    end
    return count
end

##
function simdraw_multhr(n, seeded=seeded, unseeded=unseeded, group=group, assoc=assoc)
    nseeds = length(seeded)
    nunseeds = length(unseeded)
    count = zeros(Int, nseeds, nunseeds, Threads.nthreads())
    Threads.@threads for ii in 1:n
        draw = getdraw(seeded, unseeded, group, assoc)
        for jj in eachindex(unseeded)
            count[:,jj,Threads.threadid()] += draw[unseeded[jj]] .== seeded
        end
    end
    return dropdims(sum(count, dims=3); dims=3)
end

##
function showteams(data, seeded=shortseed, unseeded=shortuns)
    return NamedArray(data, (seeded, unseeded), ("Away", "Home"))
end

max_stars_repo_name: carlocab/CLDraw

max_issues_repo_name: carlocab/CLDraw

max_forks_repo_name: carlocab/CLDraw


text:
# Main application entry point.
# Launch server before launching client.
#
# To run:
# > julia main.jl 127.0.0.1 2001

include("simulation/simulation.jl")

include("app/socket.jl")

using JSON
using .Server

soc_server = Server.SockServer()

@async while true
    input = readavailable(stdin)

    msg = String(input[1:length(input) - 1])
    # println("input: ", msg)
    if msg == "q"  #0x71 # "q" == quit
        println("Quiting...")

        trace = try
            if soc_server.running
                # println("Server is running. Toggling running flag...")
                soc_server.running = false
            end

            nothing
        catch ex
            println("### Exception ###:\n", ex)
            stacktrace(catch_backtrace())
        end

        if trace ≠ nothing
            println(trace)
            println("############################################")
            println("Goodbye.")
            exit(1)
        end

    elseif msg == "help"
        println("--------------------------------------------")
        println("Commands:")
        println("  run")
        println("  run debug")
        println("--------------------------------------------")
    end
end

Server.start(soc_server)

ext: jl


text:
## dcf.jl  decision cost function support for the ROC package.
## (c) 2015 David A. van Leeuwen
##
## Licensed under the MIT software license, see LICENSE.md
##

## effective prior odds
"""
`oeff()`, the effective prior odds given cost function parameters.  Arguments are either:

 - `ptar, cfa, cmiss`: Scalars or Vectors of the prior of a target, the cost of a false positive and the cost of a false negative

 - `dcf::DCF`: a `DCF` object containing `ptar`, `cfa` and `cmiss`.
"""
oeff(ptar=0.5; cfa=1, cmiss=1) = (cmiss ./ cfa) .* (ptar ./ (1-ptar))
oeff(dcf::DCF) = oeff(dcf.ptar; cfa=dcf.cfa, cmiss=dcf.cmiss)
## effective prior
"""
`peff()`, the effective prior given cost function parameters. Arguments are either:

 - `ptar, cfa, cmiss`: Scalars or Vectors of the prior of a target, the cost of a false positive and the cost of a false negative

 - `dcf::DCF`: a `DCF` object containing `ptar`, `cfa` and `cmiss`.
"""
peff(ptar=0.5; cfa=1, cmiss=1) = 1 ./ (1 + (cfa ./ cmiss) * ((1-ptar) ./ ptar))
peff(dcf::DCF) = peff(dcf.ptar, cfa=dcf.cfa, cmiss=dcf.cmiss)
## prior log odds
"""
`plo()`, the prior log-odds given cost function parameters. Arguments are either:

 - `ptar, cfa, cmiss`: Scalars or Vectors of the prior of a target, the cost of a false positive and the cost of a false negative

 - `dcf::DCF`: a `DCF` object containing `ptar`, `cfa` and `cmiss`.
"""
plo(ptar=0.5; cfa=1, cmiss=1) = logit(ptar) + log(cmiss ./ cfa)
plo(dcf::DCF) = plo(dcf.ptar; cfa=dcf.cfa, cmiss=dcf.cmiss)

## Decision Cost Function: dcf = p_tar * C_miss * p_miss + (1-p_tar) * C_fa * p_fa
## Compute this through the Bayes error rate, ber

## Basic definition of what Niko calls Bayes Error Rate for a single prior log-odds
## count errors _at_ the decision boundary both in fa and miss
## For a Bayes Error Rate, scores schould be log-likelihood ratios, and the Bayes decisions
## are made at threshold -plo.
## However, because this routines is also used for dcf(), we allow a separate setting of
## the threshold.
"""
`ber()`, the Bayes Error Rate.  This computes the expected error rate given a set of supervised trials,
the log-likelihood-ratio scores `tar` and `non`, and a threshold based on a given a cost function
specified by its prior log-odds (see `plo()`).

If the scores `tar` and `non` are well-calibrated log-likelihood-ratios, the optimum decision threshold
that minimizes the decision cost function `dcf` is surprisingly simple,

    θ = -plo(dcf).

This function computes the actual cost for such decisions, given the test scores.

Arguments:

 - `tar::Vector`, `non::Vector`, `plo`: target and non-target llr
scores, and prior log odds.  Here, `plo` may be a vector, resulting in
multiple Bayes error rates.

 - `r::Roc`, `plo`: a `Roc` structure computed with `roc(;collapse=false)`.
It is mandatory that the Roc object is not collapsed, because the
actual cost may occur for a threshold between two "corner" points of
the ROC curve.  A `Roc` structure is useful for computing the minimum
Bayes Error Rate for multiple cost functions---in this case the `Roc`
may be `collapse`d.

 - `; thres`: an optional threshold, overriding the theoretical
optimum described above.  This can be used if the scores are not
correctly calibrated.
"""
function ber{T1<:Real}(tar::Vector{T1}, non::Vector{T1}, plo::Real, thres::Real=-plo)
    nmiss = nfa = 0
    for x in tar
        nmiss += x ≤ thres
    end
    for x in non
        nfa += x ≥ thres
    end
    ## the effective prior is sigmoid(plo), use sigmoid(-plo) for p_non for numerical stability
    return sigmoid(plo) * nmiss / length(tar) + sigmoid(-plo) * nfa / length(non)
end
## Both the prior log odds and the threshold may be an array.
## If the plo's are an array, the thresholds---if given---must be an array of the same size,
## because it doesn't really make sense to evaluate different cost functions with the same threshold.
## For larger arrays of thres or plo, it is faster to use the Roc version below
ber{T1<:Real,T2<:Real}(tar::Vector{T1}, non::Vector{T1}, plo::Real, thres::Vector{T2}) = [ber(tar, non, plo, x) for x in thres]
function ber{T1<:Real,T2<:Real}(tar::Vector{T1}, non::Vector{T1}, plo::Vector{T2}, thres::Vector{T2}=-plo)
    size(plo) == size(thres) || error("Inconsistent vector lengths")
    [ber(tar, non, x...) for x in zip(plo,thres)]
end

## Use an _unclollapsed_ ROC for determining the Bayes error rate.
## r = roc(tar, non, collapse=false)

## We can't accurately compute this from a collapsed ROC, because
## the actual threshold might be somewhere among the set of collapsed scores.

## Bayes false alarm and miss rates, we need this for the NBE plot.
## First define a helper function, scalar and array versions
function ber_famiss(r::Roc, field::Symbol, plo::Real, thres::Real=-plo)
    i = binsearch(thres, getfield(r, field)) + 1
    return sigmoid(-plo) * r.pfa[i], sigmoid(plo) * r.pmiss[i]
end
function ber_famiss{T<:Real}(r::Roc, f::Symbol, plo::Real, thres::Array{T})
    tuples = [ber_famiss(r, f, x) for x in thres]
    return [ [ x[y] for x in tuples] for y in 1:2]
end
function ber_famiss{T<:Real}(r::Roc, f::Symbol, plo::Array{T}, thres::Array{T}=-plo)
    size(plo) == size(thres) || error("Inconsistent vector lengths")
    tuples = [ber_famiss(r, f, x) for x in plo]
    return [ [ x[y] for x in tuples] for y in 1:2]
end

## actual errors
ber_famiss(r::Roc, plo::ArrayOrReal, thres=-plo) = ber_famiss(r, :θ, plo, thres)

ber(r::Roc, plo::Real, thres::Real=-plo) = +(ber_famiss(r, plo, thres)...)
ber{T<:Real}(r::Roc, plo::Real, thres::Array{T}) = [ber(r, plo, x) for x in thres]
ber{T<:Real}(r::Roc, plo::Array{T}, thres::Array{T}=-plo) = [ber(r, x...) for x in zip(plo,thres)]

## minimum ber is best computed using a Roc structure
minber_famiss(r::Roc, plo::ArrayOrReal) = ber_famiss(r, :llr, plo)
"""
`minber()`: The minimum Byes Error Rate for a set of socres, found
after an optimal score-to-llr (log-likelihood-ratio) transformation.
The optimal transformation corresponds to the convex hull of the ROC,
where the optimal llr values correspont to the negative log of the
slope of the appropriate line segment spanning the ROC.  This is
equivalent to running the pool-adjacent-violators algorithm on the
test data.

In order to compute `minber()` for multiple cost functions, as in a Normalized Bayes Error Rate plot,
it is advantageous to compute a `Roc` object once.

Arguments: see `ber()`.
"""
minber(r::Roc, plo::Real) = +(minber_famiss(r, plo)...)
minber{T<:Real}(r::Roc, plo::Array{T}) = [minber(r, x) for x in plo]

## compute ROC if scores are given
function minber{T<:Real}(tar::Vector{T}, non::Vector{T}, plo::ArrayOrReal)
    r = roc(tar, non)
    return minber(r, plo)
end

## factor to normalize the Bayes error rate, for normalized bayes error rate (or normalized cdet)
normfactor(plo) = 1 + exp(abs(plo))

## factor to convert Bayes error rates to actual / minimum costs.
costfactor(d::DCF) = d.ptar .* d.cmiss + (1-d.ptar) .* d.cfa

## to norm or not to norm
function applyfactor(d::DCF, er, norm::Bool)
    lo = plo(d)
    if norm
        return normfactor(lo) .* er
    else
        return costfactor(d) .* er
    end
end

## Decision cost functions

## Allow a default DCF to be set
Base.copy(d::DCF) = DCF(d.ptar, d.cfa, d.cmiss)
"""
`setdcf()` allows a global decision cost function to be set.  This DCF can be used in `dcf()` and `mindcf()`.   The decision cost function is defined as

    dcf = ptar Cmiss pmiss + (1-ptar) Cfa pfa.

and is defined by the parameters `ptar` (the target prior), `Cfa` (the cost of a false positive) and `Cmiss` (the cost of a false negative.  Arguments:

 - `ptar`. The target prior, default `0.5`.

 - `cfa`.  The cost of a false positive (false alarm), default `1`

 - `cmiss`.  The cost of a false negative (miss), default `1`.

Multiple simultaneous cost functions can be set by specifying any, or a combination, of these parameters as Vectors.

The current values of the parameters of the DCF can be found using `getdcf()`.
"""
function setdcf(;ptar=0.5, cfa=1, cmiss=1, d=DCF(ptar, cfa, cmiss))
    global global_dcf = copy(d)
end
## set a default DCF
setdcf()
"""
`getdcf()` retrieves the current valua of the DCF parameters.
"""
getdcf() = global_dcf

"""
`dcf()`.  Compute the classical decision cost function:

    dcf = ptar Cmiss pmiss + (1-ptar) Cfa pfa.

The _relative_ value of the cost function is only determined by the ratio

    oeff = (ptar/(1-ptar) (Cmiss/Cfa).

Therefore, this function is computed using `ber()` and scaling accordingly.  Arguments:

 - `tar`, `non`: target and non-target scores

 - `tnt::TNT`: a `TNT` object containing target and non-target scores.

 - `r::Roc`: a `Roc` object, the result of `roc()`

 - `; d::DCF`: a decision cost function, default `getdcf()`.  This can be a vector or DCFs.

 - `; thres`: the threshold used to make decisions, default `-plo(d)`

 - `; norm`: are the costs normalized to a trivial system deciding using the prior only, default: `false`

"""
dcf(tar::Vector, non::Vector; d::DCF=getdcf(), thres=-plo(d), norm=false) = applyfactor(d, ber(tar, non, plo(d), thres), norm)
dcf(tnt::TNT; d::DCF=getdcf(), thres=-plo(d), norm=false) = dcf(tnt.tar, tnt.non, d=d, thres=thres, norm=norm)
dcf(r::Roc; d::DCF=getdcf(), thres=-plo(d), norm=false) = applyfactor(d, ber(r, plo(d), thres), norm)

"""

`mindcf()` computes the minimum costs of the classical decision cost
function (see `dcf()`) that can be obtained by varying the threshold.
Arguments are the same as for `dcf()`, except that there is no
threshold.

This function uses `minber()` and scales accordingly.
"""
mindcf(r::Roc; d::DCF=getdcf(), norm=false) = applyfactor(d, minber(r, plo(d)), norm)
mindcf(tar::Vector, non::Vector; d::DCF=getdcf(), norm=false) = applyfactor(d, minber(roc(tar, non), plo(d)), norm)
mindcf(tnt::TNT; d::DCF=getdcf(), norm=false) = mindcf(tnt.tar, tnt.non, d=d, norm=norm)

Base.show(io::IO, ::MIME"text/plain", dcf::DCF) = show(io, a)
if VERSION < v"0.5.0-dev+4356"
    Base.writemime(io::IO, ::MIME"text/plain", dcf::DCF) = show(io, dcf)
end

function Base.show(io::IO, dcf::DCF)
    println(io, "Ptar = ", dcf.ptar, ", Cfa = ", dcf.cfa, ", Cmiss = ", dcf.cmiss)
    println(io, " prior log-odds = ", plo(dcf))
    println(io, " effective prior odds = ", oeff(dcf))
    println(io, " effective prior = ", peff(dcf))
end


text:
# Replace or uncomment a line below to provide the path to a directory of spectra on your system
#path_to_spectra = "C:\\Users\\eford\\Box Sync\\SOAP simulations\\Aug2016_workshop\\SOAP_Spectra\\planet_10ms_150k"
path_to_spectra = joinpath(homedir(),"SOAP_Spectra/planet_10ms_150k")

max_stars_repo_name: eford/RVSpectraKitLearn

max_issues_repo_name: eford/RVSpectraKitLearn

max_forks_repo_name: eford/RVSpectraKitLearn


text:
export  ProHermite, degree,
        FamilyProHermite, FamilyScaledProHermite,
        DProPolyHermite,
        FamilyDProPolyHermite, FamilyDScaledProPolyHermite,
        FamilyD2ProPolyHermite, FamilyD2ScaledProPolyHermite
        derivative!, derivative,
        evaluate!, evaluate,
        vander!, vander

"""
ProHermite <: Hermite

An immutable structure for probabilistic Hermite functions defined as ψem(x) = Hem(x)*exp(-x^2/4).

## Fields
-  `m` : order of the function
-  `Poly` : probabilistic Hermite polynomial of order m
-  `scaled` : with rescaling to have unitary norm

## Constructors
- `ProHermite(m, Poly, scaled)`
- `ProHermite(m; scaled = false)`
"""
struct ProHermite <: Hermite
    m::Int64
    Poly::ProPolyHermite
    scaled::Bool
end
#
function Base.show(io::IO, P::ProHermite)
    println(io, string(P.m)*"-th order probabilistic Hermite function, scaled = "*string(P.scaled))
end

ProHermite(m::Int64; scaled::Bool = false) = ProHermite(m, ProPolyHermite(m; scaled = scaled), scaled)

"""
$(TYPEDSIGNATURES)

Return the degree of the physicist Hermite function `P`.
"""
degree(P::ProHermite) = P.m

"""
$(TYPEDSIGNATURES)

Evaluates the physicist Hermite function `P`.
"""
(P::ProHermite)(x) = P.Poly.P(x)*exp(-x^2/4)

const FamilyProHermite = ntuple(i->ProHermite(i-1),CstMaxDegree+1)
const FamilyScaledProHermite = ntuple(i->ProHermite(i-1; scaled = true),CstMaxDegree+1)

# Store Pe′_n - Pe_n * X/2 with Pe_n the n-th Probabilistic Hermite Polynomial

"""
$(TYPEDSIGNATURES)

Store Pe'ₙ - Peₙ * X/2 with Peₙ the n-th Probabilistic Hermite Polynomial
"""
function DProPolyHermite(m::Int64; scaled::Bool)
    k = 1
    factor = exp(loggamma(m+1) - loggamma(m+1-k))
    P = ProPolyHermite(m; scaled = false)
    Pprime = derivative(P, 1)
    if scaled == false
        return Pprime - ImmutablePolynomial((0.0, 0.5))*P.P
    else
        C = 1/Cpro(m)
        return C*(Pprime - ImmutablePolynomial((0.0, 0.5))*P.P)
    end
end

const FamilyDProPolyHermite = ntuple(i->AbstractProHermite(DProPolyHermite(i-1; scaled = false), false), CstMaxDegree+1)
const FamilyDScaledProPolyHermite = ntuple(i->AbstractProHermite(DProPolyHermite(i-1; scaled = true), true), CstMaxDegree+1)

# Store Pe″n -  Pe′n * X  + Pen * (X^2/2 - 1/2) with Pe_n the n-th Probabilistic Hermite Polynomial
"""
$(TYPEDSIGNATURES)

Store Pe″ₙ  - Pe′ₙ * X  + Peₙ * (X^2/2 - 1/2) with Peₙ the n-th Probabilistic Hermite Polynomial
"""
function D2ProPolyHermite(m::Int64; scaled::Bool)
    k = 2
    factor = exp(loggamma(m+1) - loggamma(m+1-k))
    P = ProPolyHermite(m; scaled = false)
    Pprime  = derivative(P, 1)
    Ppprime = derivative(P, 2)
    if scaled == false
        return Ppprime - ImmutablePolynomial((0.0, 1.0))*Pprime +
               P.P*ImmutablePolynomial((-0.5, 0.0, 0.25))
    else
        C = 1/Cpro(m)
        return C*(Ppprime - ImmutablePolynomial((0.0, 1.0))*Pprime +
                  P.P*ImmutablePolynomial((-0.5, 0.0, 0.25)))
    end
end

const FamilyD2ProPolyHermite = ntuple(i->AbstractProHermite(D2ProPolyHermite(i-1; scaled = false), false), CstMaxDegree+1)
const FamilyD2ScaledProPolyHermite = ntuple(i->AbstractProHermite(D2ProPolyHermite(i-1; scaled = true), true), CstMaxDegree+1)

"""
$(TYPEDSIGNATURES)

Computes in-place the k-th derivative of the probabilistic Hermite function `F`.
"""
function derivative!(dF, F::ProHermite, k::Int64, x)
    m = F.m
    @assert k>-2 "anti-derivative is not implemented for k<-1"
    @assert size(dF,1) == size(x,1) "Size of dF and x don't match"
    N = size(x,1)
    if k==0
        # map!(F, dF, x)
        @avx @. dF = F.Poly.P(x)*exp(-x^2/4)
        return  dF
    elseif k==1
        if F.scaled ==false
            Pprime = FamilyDProPolyHermite[m+1]
            @avx @. dF = Pprime.P(x)*exp(-x^2/4)
            return dF
        else
            Pprime = FamilyDScaledProPolyHermite[m+1]
            @avx @. dF = Pprime.P(x)*exp(-x^2/4)
            return dF
        end
    elseif k==2
        if F.scaled ==false
            Ppprime = FamilyD2ProPolyHermite[m+1]
            @avx @. dF = Ppprime.P(x)*exp(-x^2/4)
            return dF
        else
            Ppprime = FamilyD2ScaledProPolyHermite[m+1]
            @avx @. dF = Ppprime.P(x)*exp(-x^2/4)
            return dF
        end
        # map!(xi->ForwardDiff.derivative(y->ForwardDiff.derivative(z->F(z), y), xi), dF, x)
    elseif k==-1
        # Call derivative function for PhyHermite{m}
        dF .= derivative(FamilyPhyHermite[m+1], k, (1/sqrt(2))*x)
        rmul!(dF, 1/sqrt(2^m))
        rmul!(dF, 1/sqrt(2)^k)
        if F.scaled == true
            rmul!(dF, 1/Cpro(m))
        end
        return dF
    end
end

"""
$(TYPEDSIGNATURES)

Computes in-place the k-th derivative of the probabilistic Hermite function `F`.
"""
derivative(F::ProHermite, k::Int64, x::Array{Float64,1}) = derivative!(zero(x), F, k, x)

"""
$(TYPEDSIGNATURES)

Evaluates in-place the probabilistic Hermite function `P`.
"""
function evaluate!(dV, P::ProHermite, x)
    m = P.m
    N = size(x,1)
    @assert size(dV) == (N, m+1) "Wrong dimension of the Vander matrix"

    # H_0(x) = 1, H_1(x) = x
    col0 = view(dV,:,1)
    @avx @. col0 = exp(-0.25*x^2)

    if P.scaled
        rmul!(col0, 1/Cpro(0))
    end
    if m == 0
        return dV
    end

    if P.scaled
        rmul!(col0, Cpro(0))
    end

    col = view(dV,:,2)
    @avx @. col = x*col0

    if P.scaled
        rmul!(col0, 1/Cpro(0))
        rmul!(col, 1/Cpro(1))
    end

    if m == 1
        return dV
    end

    if P.scaled
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = (Cpro(i-1)* x * col - Cpro(i-2)*(i-1)*colm1)
            # dV[:,i+1] = Cpro(i-1)*x .* dV[:,i] - Cpro(i-2)*i*dV[:,i-1]
            rmul!(colp1, 1/Cpro(i))
        end
    else
        @inbounds for i=2:m
            colp1 = view(dV,:,i+1)
            col   = view(dV,:,i)
            colm1 = view(dV,:,i-1)
            @avx @. colp1 = (x * col - (i-1)*colm1)
        end
    end
    return dV
end

"""
$(TYPEDSIGNATURES)

Evaluates the probabilistic Hermite function `P`.
"""
evaluate(P::ProHermite, x::Array{Float64,1}) = evaluate!(zeros(size(x,1), P.m+1), P, x)

# Use ψe^{(k)}_n(x) = 1/√(2^n)*ψ^{k}_n(x/√2)
"""
$(TYPEDSIGNATURES)

Computes in-place the Vandermonde matrix of the k-th derivative of the probabilistic Hermite function `P`.
"""
function vander!(dV, P::ProHermite, k::Int64, x)
    m = P.m
    if k==0
        evaluate!(dV, P, x)
        return dV
    else
        # x is an Array{Float64,1} or a view of it
        C = 1/√(2.0)
        vander!(dV, FamilyPhyHermite[m+1], k, C*x)

        if P.scaled==true
            @inbounds for i=0:m
                col = view(dV,:,i+1)
                @avx @. col *= 1.0/(Cpro(i)*√(2.0^i)*√(2.0)^k)
            end
        else
            @inbounds for i=0:m
                col = view(dV,:,i+1)
                @avx @. col *= 1.0/(√(2.0^i)*√(2.0)^k)
            end
        end
        return dV
    end
end

"""
$(TYPEDSIGNATURES)

Computes the Vandermonde matrix of the k-th derivative of the probabilistic Hermite function `P`.
"""
vander(P::ProHermite, k::Int64, x::Array{Float64,1}) = vander!(zeros(size(x,1), P.m+1), P, k, x)


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 6e3ef066-d115-11eb-2338-013a707dfe8a
begin
	using PlutoUI
	using DataFrames
	using Statistics
	using RollingFunctions
	using HTTP
	using JSON
	using Chain
	using Tables
	using PerceptualColourMaps
	using Gadfly
	using MLStyle
	Gadfly.set_default_plot_size(16cm, 14cm)
end

# ╔═╡ 0da5a4ba-6b23-4f23-8165-1bd6183abbc1
html"<button onclick='present()'>Present</button>"

# ╔═╡ ba990146-058e-440e-b172-4d35b2b63985
md"""
## Looking Back

### Updates from last time
  - Compared performance on current source code to previous source code
  - Showed graphs colored by `num_embeddings` parameter but results were inconclusive because parameters were randomly sampled.

### Goals set last time
- Compare pretrained GPT to GPT with randomized parameters
- Choose best hyperparameters (one set across all games for each `num_embeddings`).
- Compare `num_embeddings` runs using single best parameter-set per `num_embeddings` value.
- Lingustic Analysis: We need to re-run to record the necessary data. Any particular runs we should prioritize?
- Generalization-Focused Tasks ([BabyAI](https://github.com/mila-iqia/babyai), [ProcGen](https://openai.com/blog/procgen-benchmark/))
  - Finish implementation of few-shot learning algorithm
  - Implement learning-from-tokens algorithm
"""

# ╔═╡ 087ee9a9-e5e6-4607-a547-824d2fea470e
md"""
However, if we extend both lines, we see that the actual performance curves are nearly identcal, measured against time.

"""

# ╔═╡ 07b1a576-f7a1-404d-9043-2827da1f4d07
md"""
Furthermore, PPO learns much more quickly than A2C when performance is measured against time steps.
"""

# ╔═╡ 0be25163-f9c1-41bc-82cf-9f63ff8e609a
md"## Compare GPT with randomly initialized GPT"

# ╔═╡ 654ad80a-a802-49d5-8373-b0b4056fd8f1
md"""
The following graph compares the performance of our GPT architecture with and without randomized (but still frozen) GPT parameters. The randomized parameter architecture is still running, so results are currently inconclusive.
"""

# ╔═╡ 942d8c47-bb2f-410c-aeea-cb2e0190cfcf
md"""
| ` num_embeddings` | pretrained | random |
|-------------------|------------|--------|
| Breakout          | 8          | 2      |
| Beamrider         | 4          | 4      |
"""

# ╔═╡ ab745b36-164c-4534-be1f-a703a2010f3e
md"""
## Choose best hyperparameters (one set across all games for each num_embeddings)
Currently the parameter sweep is still in progress. I am performing grid-search on the following:
-  `env`:
    - Seaquest-v0
    - BeamRider-v0
    - BreakoutNoFrameskip-v0
-  `gpt_size`: medium
-  `ppo_epoch`  (number of PPO updates per rollout): 3
-  `clip_param` (PPO clip parameter):
    - 0.1
    - 0.2
    - 0.3
-  `num_env_steps` (cut-off for training in steps): 10000000
-  `num_processes` (number of parallel processes):
    - 8
    - 16
-  `lr` (learning-rate):
    - 0.00025
    - 0.0003
-  `num_embeddings` (number of input embeddings to GPT):
    - 1
    - 2
    - 4
    - 8
Grid search is not complete, but there does seem to already be some correlation between performance and `num_embeddings`.
"""

# ╔═╡ 43695fbb-f8f2-4c05-a149-7ae545b2c738
md"##"

# ╔═╡ a0b13a5c-0337-43ba-89b1-e277954d4d94
md"##"

# ╔═╡ 7d96dcb9-3d63-4d14-b092-411109a7fac0
md"##"

# ╔═╡ ce675911-0344-43b6-8a8d-6750ddbb1924
max_returns = Dict(
	"BeamRider-v0" => 1590,
	"PongNoFrameskip-v0" => 20.7,
	"BreakoutNoFrameskip-v0" => 274.8,
	"Seaquest-v0" => 1204.5,
	"Qbert-v0" => 14293.3
) ;# from PPO paper

# ╔═╡ 42b94309-0cf6-4220-9c91-1fa9ba1b37af
EPISODE_RETURN = "episode return";

# ╔═╡ ab24c6bd-2721-4162-a8c2-fd7aabab1155
md"##"

# ╔═╡ 3838787a-f138-45ec-9a58-7879aeab9b99
md"""
## Linguistic Analysis
We performed analysis on the best BeamRider run using GPT-2 `medium`, which can be viewed [here](https://colab.research.google.com/drive/1heCWKd8oyOSMLtaVxpezIaSQZ9ODnyeP#scrollTo=nSrYCATp_oKa).

A high level summary of that we found:
 - a
"""

# ╔═╡ 8798b46c-1b0b-4d34-93ef-1d5062f8a632
md"""
## Architectures
This report compares the performance of a standard PPO baseline against a "GPT-in-the-middle" variant. The architectures are as follows:

#### PPO baseline

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- Linear layer with output-size $512$
- Parallel value head and action head

#### GPT
(differences in **bold**)

- Convolution with
  - output-size $32$
  - kernel-shape $8\times 8$
  - stride $4\times 4$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $4\times 4$
  - stride $2\times 2$
- ReLU
- Convolution with
  - output-size $64$
  - kernel-shape $3\times 3$
  - stride $1\times 1$
- ReLU
- **Linear layer with output-size = $n \times e$ where $n$ is the number of embeddings (a hyperparameter) and $e$ is the GPT embedding size (see table below).**
- **GPT torso**
- **one `hidden_size` layer**
- one layer projecting to action and value

| Architecture | parameters | embedding size ($e$) |
|--------------|------------|----------------|
| small        | 117M       | 768            |
| medium       | 345M       | 1024           |
| large        | 774M       | 1280           |
| xl           | 1558M      | 1600           |
"""

# ╔═╡ 33909132-16c2-4eb3-97d3-607010a82e52
md"""
## Next Steps
- How to choose best hyperparameters (one set across all games for each `num_embeddings`).
- Graph comparing the speed of different GPT sizes (e.g. run 1 million samples for each size and see how long it takes).
- "Time Steps" -> "Steps"
- Display interesting hyperparameters for the runs below each graph (e.g. how many embs for the first few graphs in this notebook).
"""

# ╔═╡ 623ed8e1-7f9e-4092-bf39-b9a3d90a7d23
# Colormap Gallery: https://colorcet.com/gallery.html
colors(n) = cmap("I1"; N=n);

# ╔═╡ 16ce1ff4-3a0e-48e5-ae84-30afb618804c
begin
	num_embeddings = [1,2,4,8]
	color_scale = Scale.color_discrete_hue(colors, levels=reverse(num_embeddings))
end;

# ╔═╡ 5d3bfc57-6a5a-41e5-8f63-6902e4958936
begin
	struct HTMLDocument
		embedded
	end
	function Base.show(io::IO, mime::MIME"text/html", doc::HTMLDocument)
		println(io, "<html>")
		show(io, mime, doc.embedded)
		println(io, "</html>")
	end
end;

# ╔═╡ 9af550be-4ef4-4938-a710-41f96392b05a

# ╔═╡ 41bd06a2-e7bd-46c6-9249-6e69223b0e11
begin
	HASRUA_ENDPOINT_URL = "http://rldl12.eecs.umich.edu:1200/v1/graphql"
	HASURA_ADMIN_SECRET = "hNuxFSTJMk28GMiFDpZDbviKtelFQamcb20UJiYZSIZ4LYsyLCCwQFYA6Y8HRSUo"
end;

# ╔═╡ 94185c13-b6d2-4337-b5ce-336f5e128032
function gql_query(query:: String; variables:: Dict = nothing)
	r = HTTP.request(
		"POST",
		HASRUA_ENDPOINT_URL;
		verbose=3,
		body= JSON.json(Dict("query" => query, "variables" => variables)) ,
		headers=Dict("x-hasura-admin-secret" => HASURA_ADMIN_SECRET),
	)
	@chain r begin
		_.body
		String
		JSON.parse(_, null=missing)
		_["data"]
	end
end;

# ╔═╡ ce12c840-ece3-48bd-8696-40c1f9802695
for (x,y) in [1=>2]
	x
end

# ╔═╡ c5b77c36-436f-4682-8d52-945ae25fb47c
gadfly_theme = :default;

# ╔═╡ 0f54629c-e045-47bf-90a5-36e44c05b8f0
function sweep_runs(sweep_ids::AbstractVector{Int}, max_step::Int)
	set_default_plot_size(18cm, 15cm)

	query = """
		query getSweepRuns(\$ids: [Int!], \$max_step: Int!) {
		  logs_less_than_step(args: {max_step: \$max_step}, where: {run: {sweep_id: {_in: \$ids}}}) {
			log
			run_id
			run {
			  metadata
			  sweep_id
			}
		  }
		}
  	"""
	rows = @chain gql_query(query; variables=Dict("ids" => sweep_ids, "max_step" => max_step)) begin
		_["logs_less_than_step"]
		map(d -> Dict(
				"run_id" => d["run_id"],
				"sweep_id" => d["run"]["sweep_id"],
				d["log"]...,
				d["run"]["metadata"]["parameters"]...
				), _)
		map(d -> Dict(
				d...,
				[k => v for (k1, v1, k2, v2) in [
							(
								"hours", get(d, "time-delta", 0) / 3600,
								"time-delta", get(d, "hours", 0) * 3600,
							),
							(
								"env", get(d, "env_name", nothing),
								"env_name", get(d, "env", nothing),
							)
						]
						for (k, v) in [
								(k1, get(d, k1, v1)),
								(k2, get(d, k2, v2)),
								]]...,
				[name => get(d, name, false) for name in [
							"randomize_parameters"
						]]...,
				[name => get(d, name, nothing) for name in [
							"action_hidden_size",
							"gpt",
							"time",
							"gae",
							"gradient_clip",
							"nonlinearity",
							"normalize_observation",
							"normalize_torso_output",
							"optimizer",
							"num_embeddings",
							"save_interval",
							"save_path",
							"config",
							# This is where the compatibility with old gpt stuff starts:
							"action loss",
							"alpha",
							"clip_param",
							"cuda",
							"linguistic_analysis_path",
							"entropy",
							"eps",
							"eval_interval",
							"gae_lambda",
							"gpt_size",
							"linear_lr_decay",
							"lr",
							"max_grad_norm",
							"num_env_steps",
							"num_mini_batch",
							"num_processes",
							"ppo_epoch",
							"recurrent_policy",
							"use_proper_time_limits",
							"value loss",
							"value_coef",
							"async_envs",
							"charts_path",
							"episode length",
							"epsilon",
							"jit",
							"lambda_",
							"learning_rate",
							"logger",
							"loss",
							"num_envs",
							"resnet",
							"rnn",
							"training_steps",
							"graphql_endpoint",
							"host_machine",
						]]...
				), _)
		collect
	end
	vcat(DataFrame.(rows)...)

end;

# ╔═╡ 03f1dfc7-970e-4920-9df9-79dd9f048e65
function a2c_vs_ppo(a2c_steps, ppo_steps)
	a2c = sweep_runs([85], a2c_steps * 4)
	time_min = minimum(a2c[!, :time])
	seconds =  (a2c[!, :time] .- time_min) ./ 1000000
	hours =  seconds ./ 3600
	a2c[!, :hours] = hours
	a2c[!, :step] /= 4
	insertcols!(a2c, :algorithm=>"A2C")
	ppo = sweep_runs([672], ppo_steps)
	insertcols!(ppo, :algorithm=>"PPO")
	vcat(a2c, ppo)
end;

# ╔═╡ 67944dbb-0ebb-44d9-b6f3-79e8d5610f61
begin
	Gadfly.with_theme(gadfly_theme) do
		plot(
			sweep_runs([672, 687], 50000000),
			x=:step, y="episode return",
			yintercept=[274.8],
			group=:run_id,
			color=:randomize_parameters,
			Guide.xlabel("Step"),
			Guide.ylabel("Episode Return"),
			Geom.line,
			Geom.hline(style=:dot,color="#b88fff"),
			Scale.color_discrete(colors),
			Guide.colorkey(title="Random Parameters"),
			Guide.title("Breakout"),
			alpha=[0.5]
		) |> HTMLDocument
	end
end

# ╔═╡ de9db738-822c-4bd6-8516-a35bd7369929
begin
	Gadfly.with_theme(gadfly_theme) do
		plot(
			sweep_runs([690, 787], 50000000),
			x=:step, y="episode return",
			yintercept=[274.8],
			group=:run_id,
			color=:randomize_parameters,
			Guide.xlabel("Step"),
			Guide.ylabel("Episode Return"),
			Geom.line,
			Geom.hline(style=:dot,color="#b88fff"),
			Scale.color_discrete(colors),
			Guide.colorkey(title="Random Parameters"),
			Guide.title("BeamRider"),
			alpha=[0.5]
		) |> HTMLDocument
	end
end

# ╔═╡ c85805a8-725d-4a84-8e8a-e1f3f1ddeae5
begin
	Gadfly.with_theme(gadfly_theme) do
		plot(
			filter(sweep_runs([784], 10000000)) do row
				row.env == "Seaquest-v0"
			end,
			x=:step, y="episode return",
			yintercept=[1204.5],
			group=:run_id,
			color="num_embeddings",
			Guide.xlabel("Step"),
			Guide.ylabel("Episode Return"),
			Geom.line,
			Geom.hline(style=:dot,color="purple"),
			color_scale,
			Guide.colorkey(title="#embeddings"),
			Guide.title("Seaquest"),
			alpha=[0.5]
		) |> HTMLDocument
	end
end

# ╔═╡ e8572b0c-6556-4c22-aa0f-cddbbda8b954
begin
	Gadfly.with_theme(gadfly_theme) do
		plot(
			filter(sweep_runs([784], 10000000)) do row
				row.env == "BeamRider-v0"
			end,
			x=:step, y="episode return",
			yintercept=[1590],
			group=:run_id,
			color="num_embeddings",
			Guide.xlabel("Step"),
			Guide.ylabel("Episode Return"),
			Geom.line,
			Geom.hline(style=:dot,color="purple"),
			color_scale,
			Guide.colorkey(title="#embeddings"),
			Guide.title("BeamRider"),
			alpha=[0.5]
		) |> HTMLDocument
	end
end

# ╔═╡ acdc2e20-8bb2-41f7-95db-973fd969490b
begin
	Gadfly.with_theme(gadfly_theme) do
		plot(
			filter(sweep_runs([784], 10000000)) do row
				row.env == "BreakoutNoFrameskip-v0"
			end,
			x=:step, y="episode return",
			yintercept=[274.8],
			group=:run_id,
			color="num_embeddings",
			Guide.xlabel("Step"),
			Guide.ylabel("Episode Return"),
			Geom.line,
			Geom.hline(style=:dot,color="purple"),
			color_scale,
			Guide.colorkey(title="#embeddings"),
			Guide.title("Breakout"),
			alpha=[0.5]
		) |> HTMLDocument
	end
end

# ╔═╡ 4e37c894-cdd4-48b2-a25a-cf7cda2d6d7b
sweeps = sweep_runs([784], 10000000) ;

# ╔═╡ df070706-0482-4804-84cd-03f9655dac88
min_returns = @chain sweeps begin
	dropmissing(_, EPISODE_RETURN)
	groupby(_, [:env])
	combine(_, EPISODE_RETURN => minimum)
	Dict(k=>v for (k,v) in eachrow(_))
end;

# ╔═╡ 132deac0-130a-4bab-8ae9-460f5a66a776
dframe = @chain sweeps begin
	filter(:step => >=(8000000), _)
	groupby(_, [:env])
	transform(_, ["env", EPISODE_RETURN] =>
		function (envs, ret)
			@match [Set(envs)...] (
				[env] => (ret .- min_returns[env]) ./ max_returns[env]
			)
		end => [EPISODE_RETURN])
	_[!, filter(names(_)) do name
			!(name in [
				"action loss",
				"config",
				"cuda",
				"entropy",
				"env",
				"eval_interval",
				"fps",
				"gradient norm",
				"hours",
				"log_interval",
				"log_level",
				"num_env_steps",
				"recurrent_policy",
				"run_id",
				"save_interval",
				"save_path",
				"step",
				"subcommand",
				"sweep_id",
				"time",
				"time-delta",
				"value loss",
			])
		end]
end;

# ╔═╡ cadd8b9d-791d-41d6-9f6c-b549d3bbd45f
function filter_by_type(ty)
	filter(name -> eltype(dframe[:, name]) == ty, names(dframe))
end;

# ╔═╡ b0d1c3b7-3cd6-4cd6-81b5-ad09f3d5cd10
df = @chain dframe begin
	groupby(_, "run ID")
	combine(_,
		filter_by_type(Float64) .=> first,
		filter_by_type(Int64) .=> first,
		filter_by_type(Bool) .=> first,
		EPISODE_RETURN .=> mean .=> :episode_return_mean,
		"env_name" .=> first .=> :env,
		)
	_[!, filter(n -> !(n in  ["run ID", "episode return_first"]), names(_))]
	sort!(_, [:episode_return_mean], rev=true)
	rename(name -> replace(name, "_first" => ""), _)
end;

# ╔═╡ 2a1fcf4d-19bb-40d9-90d1-62155843b853
begin
	bool_df = DataFrame()
	for name in names(df)
		if !occursin("episode_return", name)
			for value in df[:, name]
				new_name = string(replace(name, "_first"=>""), " = ", value)
				bool_df[:, :($new_name)] = df[:, name] .== value
			end
		end
	end
    bool_df.episode_return_mean = df.episode_return_mean
	bool_df
end;

# ╔═╡ 98c0fea9-f348-4e30-81c5-4e2f1e5e239e
function get_cor_df(env)

	filtered = @chain df begin
		filter(row -> row.env == env, _)
		_[!, filter(name -> name != "env", names(_))]
	end

	bool_df = DataFrame()
	for name in names(df)
		if !occursin("episode_return", name)
			for value in df[:, name]
				new_name = string(replace(name, "_first"=>""), " = ", value)
				bool_df[:, :($new_name)] = df[:, name] .== value
			end
		end
	end
    bool_df.episode_return_mean = df.episode_return_mean

	cor_mat = cor(Matrix(bool_df))
	return_index = findfirst(name -> name == "episode_return_mean", names(bool_df))

	correlation = cor_mat[:, return_index]

	name = names(bool_df)
	mapping = Dict(zip(name, correlation))
	[
		mapping["num_embeddings = 1"],
		mapping["num_embeddings = 2"],
		mapping["num_embeddings = 4"],
		mapping["num_embeddings = 8"],
	]
end;

# ╔═╡ 8bf7862f-9be4-44de-9d14-b8056d9d27c2
cor_df = vcat([
		DataFrame(
			num_embeddings=[1, 2, 4, 8],
			env=repeat([env], 4),
			correlation=get_cor_df(env)
			) for env in [
				"Seaquest-v0",
				"BreakoutNoFrameskip-v4",
				"BeamRider-v0",
			]
		]...)

# ╔═╡ 5a968482-35cc-4afa-b658-0ab130a47508
plot(cor_df, x=:num_embeddings, y=:correlation, color=:env, Geom.point, Geom.line)

# ╔═╡ bbed8889-d401-4167-9370-9927aabbc83b
# @bind window_size Select(string.(collapse_runs ? (100:100:500) : (1:10)), default=collapse_runs ? "200" : "5")
window_size = "11";

# ╔═╡ Cell order:
# ╟─0da5a4ba-6b23-4f23-8165-1bd6183abbc1
# ╟─ba990146-058e-440e-b172-4d35b2b63985
# ╟─087ee9a9-e5e6-4607-a547-824d2fea470e
# ╟─07b1a576-f7a1-404d-9043-2827da1f4d07
# ╟─03f1dfc7-970e-4920-9df9-79dd9f048e65
# ╟─0be25163-f9c1-41bc-82cf-9f63ff8e609a
# ╟─654ad80a-a802-49d5-8373-b0b4056fd8f1
# ╟─67944dbb-0ebb-44d9-b6f3-79e8d5610f61
# ╠═de9db738-822c-4bd6-8516-a35bd7369929
# ╟─942d8c47-bb2f-410c-aeea-cb2e0190cfcf
# ╟─ab745b36-164c-4534-be1f-a703a2010f3e
# ╟─43695fbb-f8f2-4c05-a149-7ae545b2c738
# ╠═c85805a8-725d-4a84-8e8a-e1f3f1ddeae5
# ╟─a0b13a5c-0337-43ba-89b1-e277954d4d94
# ╟─e8572b0c-6556-4c22-aa0f-cddbbda8b954
# ╟─7d96dcb9-3d63-4d14-b092-411109a7fac0
# ╟─acdc2e20-8bb2-41f7-95db-973fd969490b
# ╟─4e37c894-cdd4-48b2-a25a-cf7cda2d6d7b
# ╟─df070706-0482-4804-84cd-03f9655dac88
# ╟─ce675911-0344-43b6-8a8d-6750ddbb1924
# ╟─cadd8b9d-791d-41d6-9f6c-b549d3bbd45f
# ╟─42b94309-0cf6-4220-9c91-1fa9ba1b37af
# ╟─132deac0-130a-4bab-8ae9-460f5a66a776
# ╟─b0d1c3b7-3cd6-4cd6-81b5-ad09f3d5cd10
# ╟─2a1fcf4d-19bb-40d9-90d1-62155843b853
# ╟─98c0fea9-f348-4e30-81c5-4e2f1e5e239e
# ╟─ab24c6bd-2721-4162-a8c2-fd7aabab1155
# ╠═8bf7862f-9be4-44de-9d14-b8056d9d27c2
# ╠═5a968482-35cc-4afa-b658-0ab130a47508
# ╠═3838787a-f138-45ec-9a58-7879aeab9b99
# ╟─8798b46c-1b0b-4d34-93ef-1d5062f8a632
# ╠═33909132-16c2-4eb3-97d3-607010a82e52
# ╟─16ce1ff4-3a0e-48e5-ae84-30afb618804c
# ╟─623ed8e1-7f9e-4092-bf39-b9a3d90a7d23
# ╟─5d3bfc57-6a5a-41e5-8f63-6902e4958936
# ╟─6e3ef066-d115-11eb-2338-013a707dfe8a
# ╟─9af550be-4ef4-4938-a710-41f96392b05a
# ╟─41bd06a2-e7bd-46c6-9249-6e69223b0e11
# ╟─94185c13-b6d2-4337-b5ce-336f5e128032
# ╟─ce12c840-ece3-48bd-8696-40c1f9802695
# ╟─c5b77c36-436f-4682-8d52-945ae25fb47c
# ╟─0f54629c-e045-47bf-90a5-36e44c05b8f0
# ╟─bbed8889-d401-4167-9370-9927aabbc83b


text:
"""
    FixedRankTaxonomy

    What makes this a `FixedRankTaxonomy` is that within the
    ClassificationTree `tree`, the depth of each node `d` corresponds
    to the index `d` in `ranks`.

    To use a taxonomy that enables classifications that do not correspond
    to a rank in `ranks`, e.g. super-class, sub-genera, sub-sub-family, etc,
    see `FlexibleTaxonomy`

"""
struct FixedRankTaxonomy <: AbstractTaxonomy
    ranks::Vector{TaxonomicRank}
    root::Union{TaxonNode, Nothing}
    FixedRankTaxonomy(root::TaxonNode; ranks = defaultranks()) = new(ranks, root)
end

"""
    FlexibleRankTaxonomy

    A `FlexibleTaxonomy` contains a set of `TaxonomicRanks`,
    but does not enfore that the depth of each node `d` correspond
    to a rank; instead "between"-ranks separating
    named ranks can exist, and are described as being arbitrary numbers of either `super` or `sub`
    from the nearest named rank.

"""
struct FlexibleRankTaxonomy <: AbstractTaxonomy
    ranks::Vector{TaxonomicRank}
    root::Union{TaxonNode, Nothing}
    FlexibleRankTaxonomy(root::TaxonNode; ranks = defaultranks()) = new(ranks, root)
end

tree(t::AbstractTaxonomy) = t.tree
ranks(t::AbstractTaxonomy) = t.ranks
root(t::AbstractTaxonomy) = t.root

Base.show(io::IO,t::FlexibleRankTaxonomy) = Base.show(io, AbstractTrees.print_tree(root(t)))
Base.show(io::IO, t::FixedRankTaxonomy) = Base.show(io, AbstractTrees.print_tree(root(t)))

defaultranks() = [TaxonomicRank(Symbol(s),i-1) for (i,s) in enumerate(["domain", "kingdom", "phylum", "class", "order", "family", "genus", "species"])]


text:
module Turf

import Base: contains

using Reexport
@reexport using GeoInterface
@reexport using GeoJSON

export radians_to_length, length_to_radians, length_to_degrees, convert_length, convert_area,
    angle_adjacent, distance, distance_to_segment, rhumb_distance, rhumb_bearing, rhumb_destination,
    destination, bearing, bearing_to_azimuth, bbox, bbox_polygon, center, centroid, transform_rotate, transform_scale,
    bezier, concave, clockwise, nearestpoint, earth_radius, area_factors, units_factors, parallel,
    Spline, linesegment, point_on_line, point_in_polygon, mediancenter, masscenter, meancenter, circle,
    pnorm_distance, distance_weight, ellipse, explode, midpoint, square, flip, lineclip, polygonclip,
    linearc, sector, planepoint, to_WGS84, to_mercator, point_grid, rectangle_grid, hexgrid, polygon_tangents,
    area, contains, within, disjoint, line_intersects, polygon_to_line, transform_translate, valid, crosses, convert_to,
    combine, scale, square_grid, triangle_grid, tag, point_to_line_distance, simplify, simplify!, clean, clean!, overlap

include("Constants.jl")
include("Utils.jl")
include("lib/Angle.jl")
include("lib/Distance.jl")
include("lib/Destination.jl")
include("lib/Bearing.jl")
include("lib/BBox.jl")
include("lib/Lines.jl")
include("lib/Centering.jl")
include("lib/Transformations.jl")
include("lib/Splines.jl")
include("lib/Circle.jl")
include("lib/Ellipse.jl")
include("lib/Planes.jl")
include("lib/Square.jl")
include("lib/Grids.jl")
include("lib/BezierSpline.jl")
include("lib/Booleans.jl")

end


text:
using FillArrays, LinearAlgebra, SparseArrays, Random, Test
import FillArrays: AbstractFill

@testset "fill array constructors and convert" begin
    for (Typ, funcs) in ((:Zeros, :zeros), (:Ones, :ones))
        @eval begin
            @test $Typ((-1,5)) == $Typ((0,5))
            @test $Typ(5) isa AbstractVector{Float64}
            @test $Typ(5,5) isa AbstractMatrix{Float64}
            @test $Typ(5) == $Typ((5,))
            @test $Typ(5,5) == $Typ((5,5))
            @test eltype($Typ(5,5)) == Float64

            for T in (Int, Float64)
                Z = $Typ{T}(5)
                @test eltype(Z) == T
                @test Array(Z) == $funcs(T,5)
                @test Array{T}(Z) == $funcs(T,5)
                @test Array{T,1}(Z) == $funcs(T,5)

                @test convert(AbstractArray,Z) ≡ Z
                @test convert(AbstractArray{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                @test convert(AbstractVector{T},Z) ≡ AbstractVector{T}(Z) ≡ Z

                @test $Typ{T,1}(2ones(T,5)) == Z
                @test $Typ{T}(2ones(T,5)) == Z
                @test $Typ(2ones(T,5)) == Z

                Z = $Typ{T}(5, 5)
                @test eltype(Z) == T
                @test Array(Z) == $funcs(T,5,5)
                @test Array{T}(Z) == $funcs(T,5,5)
                @test Array{T,2}(Z) == $funcs(T,5,5)

                @test convert(AbstractArray,Z) ≡ convert(AbstractFill,Z) ≡ Z
                @test convert(AbstractArray{T},Z) ≡ convert(AbstractFill{T},Z) ≡ AbstractArray{T}(Z) ≡ Z
                @test convert(AbstractMatrix{T},Z) ≡ convert(AbstractFill{T,2},Z) ≡ AbstractMatrix{T}(Z) ≡ Z

                @test_throws Exception convert(Fill{Float64}, [1,1,2])
                @test_throws Exception convert(Fill, [])
                @test convert(Fill{Float64}, [1,1,1]) ≡ Fill(1.0, 3)
                @test convert(Fill, Float64[1,1,1]) ≡ Fill(1.0, 3)
                @test convert(Fill{Float64}, Fill(1.0,2)) ≡ Fill(1.0, 2) # was ambiguous
                @test convert(Fill{Int}, Ones(20)) ≡ Fill(1, 20)

                @test $Typ{T,2}(2ones(T,5,5)) == Z
                @test $Typ{T}(2ones(T,5,5)) == Z
                @test $Typ(2ones(T,5,5)) == Z

                @test AbstractArray{Float32}(Z) == $Typ{Float32}(5,5)
                @test AbstractArray{Float32,2}(Z) == $Typ{Float32}(5,5)
            end
        end
    end

    @test Fill(1,(-1,5)) ≡ Fill(1,(0,5))
    @test Fill(1.0,5) isa AbstractVector{Float64}
    @test Fill(1.0,5,5) isa AbstractMatrix{Float64}
    @test Fill(1,5) == Fill(1,(5,))
    @test Fill(1,5,5) == Fill(1,(5,5))
    @test eltype(Fill(1.0,5,5)) == Float64

    @test Matrix{Float64}(Zeros{ComplexF64}(10,10)) == zeros(10,10)
    @test_throws InexactError Matrix{Float64}(Fill(1.0+1.0im,10,10))

    for T in (Int, Float64)
        F = Fill{T}(one(T), 5)

        @test eltype(F) == T
        @test Array(F) == fill(one(T),5)
        @test Array{T}(F) == fill(one(T),5)
        @test Array{T,1}(F) == fill(one(T),5)

        F = Fill{T}(one(T), 5, 5)
        @test eltype(F) == T
        @test Array(F) == fill(one(T),5,5)
        @test Array{T}(F) == fill(one(T),5,5)
        @test Array{T,2}(F) == fill(one(T),5,5)

        @test convert(AbstractArray,F) ≡ F
        @test convert(AbstractArray{T},F) ≡ AbstractArray{T}(F) ≡ F
        @test convert(AbstractMatrix{T},F) ≡ AbstractMatrix{T}(F) ≡ F

        @test convert(AbstractArray{Float32},F) == AbstractArray{Float32}(F) ==
                Fill{Float32}(one(Float32),5,5)
        @test convert(AbstractMatrix{Float32},F) == AbstractMatrix{Float32}(F) ==
                Fill{Float32}(one(Float32),5,5)
    end

    @test Eye(5) isa Diagonal{Float64}
    @test Eye(5) == Eye{Float64}(5)
    @test eltype(Eye(5)) == Float64

    for T in (Int, Float64)
        E = Eye{T}(5)
        M = Matrix{T}(I, 5, 5)

        @test eltype(E) == T
        @test Array(E) == M
        @test Array{T}(E) == M
        @test Array{T,2}(E) == M

        @test convert(AbstractArray,E) === E
        @test convert(AbstractArray{T},E) === E
        @test convert(AbstractMatrix{T},E) === E

        @test AbstractArray{Float32}(E) == Eye{Float32}(5)
    end

    @testset "Bool should change type" begin
        x = Fill(true,5)
        y = x + x
        @test y isa Fill{Int,1}
        @test y[1] == 2

        x = Ones{Bool}(5)
        y = x + x
        @test y isa Fill{Int,1}
        @test y[1] == 2
        @test x + Zeros{Bool}(5) ≡ x
        @test x - Zeros{Bool}(5) ≡ x
        @test Zeros{Bool}(5) + x ≡ x
        @test -x ≡ Fill(-1,5)
    end
end

# Check that all pair-wise combinations of + / - elements of As and Bs yield the correct
# type, and produce numerically correct results.
function test_addition_and_subtraction(As, Bs, Tout::Type)
    for A in As, B in Bs
        @test A + B isa Tout{promote_type(eltype(A), eltype(B))}
        @test Array(A + B) == Array(A) + Array(B)

        @test A - B isa Tout{promote_type(eltype(A), eltype(B))}
        @test Array(A - B) == Array(A) - Array(B)

        @test B + A isa Tout{promote_type(eltype(B), eltype(A))}
        @test Array(B + A) == Array(B) + Array(A)

        @test B - A isa Tout{promote_type(eltype(B), eltype(A))}
        @test Array(B - A) == Array(B) - Array(A)
    end
end

# Check that all permutations of + / - throw a `DimensionMismatch` exception.
function test_addition_and_subtraction_dim_mismatch(a, b)
    @test_throws DimensionMismatch a + b
    @test_throws DimensionMismatch a - b
    @test_throws DimensionMismatch b + a
    @test_throws DimensionMismatch b - a
end

@testset "FillArray addition and subtraction" begin
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros(6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros{Int}(6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros(6,6))
    test_addition_and_subtraction_dim_mismatch(Zeros(5), Zeros{Int}(6,5))

    # Construct FillArray for repeated use.
    rng = MersenneTwister(123456)
    A_fill, B_fill = Fill(randn(rng, Float64), 5), Fill(4, 5)

    # Unary +/- constructs a new FillArray.
    @test +A_fill === A_fill
    @test -A_fill === Fill(-A_fill.value, 5)

    # FillArray +/- FillArray should construct a new FillArray.
    test_addition_and_subtraction([A_fill, B_fill], [A_fill, B_fill], Fill)
    test_addition_and_subtraction_dim_mismatch(A_fill, Fill(randn(rng), 5, 2))

    # FillArray + Array (etc) should construct a new Array using `getindex`.
    A_dense, B_dense = randn(rng, 5), [5, 4, 3, 2, 1]
    test_addition_and_subtraction([A_fill, B_fill], [A_dense, B_dense], Array)
    test_addition_and_subtraction_dim_mismatch(A_fill, randn(rng, 5, 2))

    # FillArray + StepLenRange / UnitRange (etc) should yield an AbstractRange.
    A_ur, B_ur = 1.0:5.0, 6:10
    test_addition_and_subtraction([A_fill, B_fill], (A_ur, B_ur), AbstractRange)
    test_addition_and_subtraction_dim_mismatch(A_fill, 1.0:6.0)
    test_addition_and_subtraction_dim_mismatch(A_fill, 5:10)
end

@testset "Other matrix types" begin
    @test Diagonal(Zeros(5)) == Diagonal(zeros(5))

    @test Diagonal(Zeros(8,5)) == Diagonal(zeros(5))
    @test convert(Diagonal, Zeros(5,5)) == Diagonal(zeros(5))
    @test_throws BoundsError convert(Diagonal, Zeros(8,5))

    @test convert(Diagonal{Int}, Zeros(5,5)) == Diagonal(zeros(Int,5))
    @test_throws BoundsError convert(Diagonal{Int}, Zeros(8,5))

    @test convert(Diagonal, Eye(5)) == Diagonal(ones(5))
    @test convert(Diagonal{Int}, Eye(5)) == Diagonal(ones(Int,5))
end

@testset "Sparse vectors and matrices" begin
    @test SparseVector(Zeros(5)) ==
            SparseVector{Float64}(Zeros(5)) ==
            SparseVector{Float64,Int}(Zeros(5)) ==
            convert(AbstractSparseArray,Zeros(5)) ==
            convert(AbstractSparseVector,Zeros(5)) ==
            convert(AbstractSparseArray{Float64},Zeros(5)) ==
            convert(AbstractSparseVector{Float64},Zeros(5)) ==
            convert(AbstractSparseVector{Float64,Int},Zeros(5)) ==
            spzeros(5)

    for (Mat, SMat) in ((Zeros(5,5), spzeros(5,5)), (Zeros(6,5), spzeros(6,5)),
                        (Eye(5), sparse(I,5,5)))
        @test SparseMatrixCSC(Mat) ==
                SparseMatrixCSC{Float64}(Mat) ==
                SparseMatrixCSC{Float64,Int}(Mat) ==
                convert(AbstractSparseArray,Mat) ==
                convert(AbstractSparseMatrix,Mat) ==
                convert(AbstractSparseArray{Float64},Mat) ==
                convert(AbstractSparseArray{Float64,Int},Mat) ==
                convert(AbstractSparseMatrix{Float64},Mat) ==
                convert(AbstractSparseMatrix{Float64,Int},Mat) ==
                SMat
    end
end

@testset "Rank" begin
    @test rank(Zeros(5,4)) == 0
    @test rank(Ones(5,4)) == 1
    @test rank(Fill(2,5,4)) == 1
    @test rank(Fill(0,5,4)) == 0
    @test rank(Eye(2)) == 2
end

@testset "BigInt indices" begin
    for A in (Zeros(BigInt(100)), Ones(BigInt(100)), Fill(2, BigInt(100)))
        @test length(A) isa BigInt
        @test axes(A) == tuple(Base.OneTo{BigInt}(BigInt(100)))
        @test size(A) isa Tuple{BigInt}
    end
    let A = Eye(BigInt(100))
        @test length(A) isa BigInt
        @test axes(A) == tuple(Base.OneTo{BigInt}(BigInt(100)),Base.OneTo{BigInt}(BigInt(100)))
        @test size(A) isa Tuple{BigInt,BigInt}
    end
    for A in (Zeros(BigInt(10), 10), Ones(BigInt(10), 10), Fill(2.0, (BigInt(10), 10)))
        @test size(A) isa Tuple{BigInt,Int}
    end

end

@testset "IndexStyle" begin
    @test IndexStyle(Zeros(5,5)) == IndexStyle(typeof(Zeros(5,5))) == IndexLinear()
end

@testset "Identities" begin
    @test Zeros(3,4) * randn(4,5) === randn(3,4) * Zeros(4,5) === Zeros(3, 5)
    @test_throws DimensionMismatch randn(3,4) * Zeros(3, 3)
    @test eltype(Zeros{Int}(3,4) * fill(1, 4, 5)) == Int
    @test eltype(Zeros{Int}(3,4) * fill(3.4, 4, 5)) == Float64
    @test Zeros(3, 4) * randn(4) == Zeros(3, 4) * Zeros(4) == Zeros(3)
    @test Zeros(3, 4) * Zeros(4, 5) === Zeros(3, 5)

    @test [1,2,3]*Zeros(1) ≡ Zeros(3)
    @test [1,2,3]*Zeros(1,3) ≡ Zeros(3,3)
    @test_throws DimensionMismatch [1,2,3]*Zeros(3)

    # Check multiplication by Adjoint vectors works as expected.
    @test randn(4, 3)' * Zeros(4) === Zeros(3)
    @test randn(4)' * Zeros(4) === zero(Float64)
    @test [1, 2, 3]' * Zeros{Int}(3) === zero(Int)
    @test_throws DimensionMismatch randn(4)' * Zeros(3)

    # Check multiplication by Transpose-d vectors works as expected.
    @test transpose(randn(4, 3)) * Zeros(4) === Zeros(3)
    @test transpose(randn(4)) * Zeros(4) === zero(Float64)
    @test transpose([1, 2, 3]) * Zeros{Int}(3) === zero(Int)
    @test_throws DimensionMismatch transpose(randn(4)) * Zeros(3)

    @test +(Zeros{Float64}(3, 5)) === Zeros{Float64}(3, 5)
    @test -(Zeros{Float32}(5, 2)) === Zeros{Float32}(5, 2)

    # `Zeros` are closed under addition and subtraction (both unary and binary).
    z1, z2 = Zeros{Float64}(4), Zeros{Int}(4)
    @test +(z1) === z1
    @test -(z1) === z1

    test_addition_and_subtraction([z1, z2], [z1, z2], Zeros)
    test_addition_and_subtraction_dim_mismatch(z1, Zeros{Float64}(4, 2))

    # `Zeros` +/- `Fill`s should yield `Fills`.
    fill1, fill2 = Fill(5.0, 4), Fill(5, 4)
    test_addition_and_subtraction([z1, z2], [fill1, fill2], Fill)
    test_addition_and_subtraction_dim_mismatch(z1, Fill(5, 5))

    X = randn(3, 5)
    for op in [+, -]

        # Addition / subtraction with same eltypes.
        @test op(Zeros(6, 4), Zeros(6, 4)) === Zeros(6, 4)
        @test_throws DimensionMismatch op(X, Zeros(4, 6))
        @test eltype(op(Zeros(3, 5), X)) == Float64

        # Different eltypes, the other way around.
        @test op(X, Zeros{Float32}(3, 5)) isa Matrix{Float64}
        @test !(op(X, Zeros{Float32}(3, 5)) === X)
        @test op(X, Zeros{Float32}(3, 5)) == X
        @test !(op(X, Zeros{ComplexF64}(3, 5)) === X)
        @test op(X, Zeros{ComplexF64}(3, 5)) == X

        # Addition / subtraction of Zeros.
        @test eltype(op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5))) == Float64
        @test eltype(op(Zeros{Int}(5, 4), Zeros{Float32}(5, 4))) == Float32
        @test op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5)) isa Zeros{Float64}
        @test op(Zeros{Float64}(4, 5), Zeros{Int}(4, 5)) === Zeros{Float64}(4, 5)
    end

    # Zeros +/- dense where + / - have different results.
    @test +(Zeros(3, 5), X) == X && +(X, Zeros(3, 5)) == X
    @test !(Zeros(3, 5) + X === X) && !(X + Zeros(3, 5) === X)
    @test -(Zeros(3, 5), X) == -X

    # Addition with different eltypes.
    @test +(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
    @test !(+(Zeros{Float32}(3, 5), X) === X)
    @test +(Zeros{Float32}(3, 5), X) == X
    @test !(+(Zeros{ComplexF64}(3, 5), X) === X)
    @test +(Zeros{ComplexF64}(3, 5), X) == X

    # Subtraction with different eltypes.
    @test -(Zeros{Float32}(3, 5), X) isa Matrix{Float64}
    @test -(Zeros{Float32}(3, 5), X) == -X
    @test -(Zeros{ComplexF64}(3, 5), X) == -X

    # Tests for ranges.
    X = randn(5)
    @test !(Zeros(5) + X === X)
    @test Zeros{Int}(5) + (1:5) === (1:5) && (1:5) + Zeros{Int}(5) === (1:5)
    @test Zeros(5) + (1:5) === (1.0:1.0:5.0) && (1:5) + Zeros(5) === (1.0:1.0:5.0)
    @test (1:5) - Zeros{Int}(5) === (1:5)
    @test Zeros{Int}(5) - (1:5) === -1:-1:-5
    @test Zeros(5) - (1:5) === -1.0:-1.0:-5.0
end

@testset "maximum/minimum/svd/sort" begin
    @test maximum(Fill(1, 1_000_000_000)) == minimum(Fill(1, 1_000_000_000)) == 1
    @test svdvals(fill(2,5,6)) ≈ svdvals(Fill(2,5,6))
    @test svdvals(Eye(5)) === Fill(1.0,5)
    @test sort(Ones(5)) == sort!(Ones(5))
end

@testset "Cumsum" begin
    @test sum(Fill(3,10)) ≡ 30
    @test cumsum(Fill(3,10)) ≡ 3:3:30

    @test sum(Ones(10)) ≡ 10.0
    @test cumsum(Ones(10)) ≡ 1.0:10.0

    @test sum(Ones{Int}(10)) ≡ 10
    @test cumsum(Ones{Int}(10)) ≡ Base.OneTo(10)

    @test sum(Zeros(10)) ≡ 0.0
    @test cumsum(Zeros(10)) ≡ Zeros(10)

    @test sum(Zeros{Int}(10)) ≡ 0
    @test cumsum(Zeros{Int}(10)) ≡ Zeros{Int}(10)

    @test cumsum(Zeros{Bool}(10)) ≡ Zeros{Bool}(10)
    @test cumsum(Ones{Bool}(10)) ≡ Base.OneTo{Int}(10)
    @test cumsum(Fill(true,10)) ≡ 1:1:10
end

@testset "Broadcast" begin
    x = Fill(5,5)
    @test (.+)(x) ≡ x
    @test (.-)(x) ≡ -x
    @test exp.(x) ≡ Fill(exp(5),5)
    @test x .+ 1 ≡ Fill(6,5)
    @test x .+ x ≡ Fill(10,5)
    @test x .+ Ones(5) ≡ Fill(6.0,5)
    f = (x,y) -> cos(x*y)
    @test f.(x, Ones(5)) ≡ Fill(f(5,1.0),5)

    y = Ones(5,5)
    @test (.+)(y) ≡ Fill(1.0,5,5)
    @test (.-)(y) ≡ Fill(-1.0,5,5)
    @test exp.(y) ≡ Fill(exp(1),5,5)
    @test y .+ 1 ≡ Fill(2.0,5,5)
    @test y .+ y ≡ Fill(2.0,5,5)
    @test y .* y ≡ y ./ y ≡ y .\ y ≡ y

    @test Zeros{Int}(5) .+ Zeros(5) isa Zeros{Float64}
end

@testset "Sub-arrays" begin
    A = Fill(3.0,5)
    @test A[1:3] ≡ Fill(3.0,3)
    @test A[1:3,1:1] ≡ Fill(3.0,3,1)
    @test_broken A[1:3,2] ≡ Zeros{Int}(3)
    @test_throws BoundsError A[1:26]
    @test A[[true, false, true, false, false]] ≡ Fill(3.0, 2)
    A = Fill(3.0, 2, 2)
    @test A[[true true; true false]] ≡ Fill(3.0, 3)
    @test_throws DimensionMismatch A[[true, false]]

    A = Ones{Int}(5,5)
    @test A[1:3] ≡ Ones{Int}(3)
    @test A[1:3,1:2] ≡ Ones{Int}(3,2)
    @test_broken A[1:3,2] ≡ Ones{Int}(3)
    @test_throws BoundsError A[1:26]
    A = Ones{Int}(2,2)
    @test A[[true false; true false]] ≡ Ones{Int}(2)
    @test A[[true, false, true, false]] ≡ Ones{Int}(2)
    @test_throws DimensionMismatch A[[true false false; true false false]]

    A = Zeros{Int}(5,5)
    @test A[1:3] ≡ Zeros{Int}(3)
    @test A[1:3,1:2] ≡ Zeros{Int}(3,2)
    @test_broken A[1:3,2] ≡ Zeros{Int}(3)
    @test_throws BoundsError A[1:26]
    A = Zeros{Int}(2,2)
    @test A[[true false; true false]] ≡ Zeros{Int}(2)
    @test A[[true, false, true, false]] ≡ Zeros{Int}(2)
    @test_throws DimensionMismatch A[[true false false; true false false]]
end

@testset "Offset indexing" begin
    A = Fill(3, (Base.Slice(-1:1),))
    @test axes(A)  == (Base.Slice(-1:1),)
    @test A[0] == 3
    @test_throws BoundsError A[2]
    @test_throws BoundsError A[-2]

    A = Zeros((Base.Slice(-1:1),))
    @test axes(A)  == (Base.Slice(-1:1),)
    @test A[0] == 0
    @test_throws BoundsError A[2]
    @test_throws BoundsError A[-2]
end

@testset "0-dimensional" begin
    A = Fill{Int,0,Tuple{}}(3, ())

    @test A[] ≡ A[1] ≡ 3
    @test A ≡ Fill{Int,0}(3, ()) ≡ Fill(3, ()) ≡ Fill(3)
    @test size(A) == ()
    @test axes(A) == ()

    A = Ones{Int,0,Tuple{}}(())
    @test A[] ≡ A[1] ≡ 1
    @test A ≡ Ones{Int,0}(()) ≡ Ones{Int}(()) ≡ Ones{Int}()

    A = Zeros{Int,0,Tuple{}}(())
    @test A[] ≡ A[1] ≡ 0
    @test A ≡ Zeros{Int,0}(()) ≡ Zeros{Int}(()) ≡ Zeros{Int}()
end

@testset "unique" begin
    @test unique(Fill(12, 20)) == unique(fill(12, 20))
    @test unique(Fill(1, 0)) == []
    @test unique(Zeros(0)) isa Vector{Float64}
    @test !allunique(Fill("a", 2))
    @test allunique(Ones(0))
end

@testset "Zero .*" begin
    @test (1:10) .* Zeros(10) ≡ Zeros(10) .* (1:10) ≡ Zeros(10) .* randn(10) ≡ Zeros(10)
    @test (1:10) .* Zeros{Int}(10) ≡ Zeros{Int}(10)
    @test_throws DimensionMismatch (1:11) .* Zeros(10)
end

@testset "iterate" begin
    for d in (0, 1, 2, 100)
        for T in (Float64, Int)
            m = Eye(d)
            mcp = [x for x in m]
            @test mcp == m
            @test eltype(mcp) == eltype(m)
        end
    end
end

@testset "properties" begin
    for d in (0, 1, 2, 100)
        @test isone(Eye(d))
    end
end

@testset "any all iszero isone" begin
    for T in (Int, Float64, ComplexF64)
        for d in (0, )
            m = Eye{T}(d)
            @test ! any(isone, m)
            @test ! any(iszero, m)
            @test ! all(iszero, m)
            @test ! all(isone, m)
        end
        for d in (1, )
            m = Eye{T}(d)
            @test ! any(iszero, m)
            @test ! all(iszero, m)
            @test any(isone, m)
            @test all(isone, m)

            onem = Ones{T}(d, d)
            @test isone(onem)
            @test ! iszero(onem)

            zerom = Zeros{T}(d, d)
            @test ! isone(zerom)
            @test  iszero(zerom)

            fillm0 = Fill(T(0), d, d)
            @test ! isone(fillm0)
            @test   iszero(fillm0)

            fillm1 = Fill(T(1), d, d)
            @test isone(fillm1)
            @test ! iszero(fillm1)

            fillm2 = Fill(T(2), d, d)
            @test ! isone(fillm2)
            @test ! iszero(fillm2)
        end
        for d in (2, 3)
            m = Eye{T}(d)
            @test any(iszero, m)
            @test ! all(iszero, m)
            @test any(isone, m)
            @test ! all(isone, m)

            m1 = Ones{T}(d, d)
            @test ! isone(m1)
            @test ! iszero(m1)
            @test all(isone, m1)
            @test ! all(iszero, m1)

            m2 = Zeros{T}(d, d)
            @test ! isone(m2)
            @test iszero(m2)
            @test ! all(isone, m2)
            @test  all(iszero, m2)

            m3 = Fill(T(2), d, d)
            @test ! isone(m3)
            @test ! iszero(m3)
            @test ! all(isone, m3)
            @test ! all(iszero, m3)
            @test ! any(iszero, m3)

            m4 = Fill(T(1), d, d)
            @test ! isone(m4)
            @test ! iszero(m4)
        end
    end
end


text:
cd(dirname(@__FILE__))
download("https://github.com/vermaseren/form/releases/download/v4.1-20131025/form-4.1-x86_64-linux.tar.gz", "./form-4.1-x86_64-linux.tar.gz")
run(`tar xzvf form-4.1-x86_64-linux.tar.gz`)
run(`mv form-4.1-x86_64-linux/tform bin/`)
run(`mv form-4.1-x86_64-linux/form bin/`)
run(`rmdir form-4.1-x86_64-linux`)
run(`rm form-4.1-x86_64-linux.tar.gz`)


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

function value(obj::TimeVectorFunction, arg0::FieldAbsoluteDate)
    return jcall(obj, "value", Vector{RealFieldElement}, (FieldAbsoluteDate,), arg0)
end


text:
for i=1:x ## loop able to produce more than 1 gigapixel per hour per core on my laptop if fed correctly
    for j=1:y
        a = los(pincam,i,j)
        #frame[i,j]
        temp = raymarch(a[1],a[2],scene,1000.0) ## upto 3x over head due to excessive copying
        ## Stencils on this ^ would be fun :)

        myPic[i,j] = ShaderArray[temp.Shader](temp) ##TODO:reflections have no place in this model
        ## some improvement possible with the shader
    end
end


text:
using Documenter, SIAMFANLEquations, DocumenterTools
push!(LOAD_PATH,"../src/")
makedocs(sitename="SIAMFANLEquations.jl",
authors="C. T. Kelley",
format = Documenter.HTML(
               prettyurls = get(ENV, "CI", nothing) == "true"
           ),
pages = Any[
     "Home" => "index.md",
     "Solvers" => Any[
       "functions/nsol.md",
       "functions/ptcsol.md",
       "functions/nsoli.md",
       "functions/ptcsoli.md",
       ],
     "Scalar Equations" => Any[
       "functions/nsolsc.md",
       "functions/ptcsolsc.md",
       "functions/secant.md",
       ],
     "Linear Solvers" => Any[
       "functions/kl_gmres.md",
       "functions/kl_bicgstab.md",
       ]
]
)
deploydocs(
     repo="github.com/ctkelley/SIAMFANLEquations.jl.git"
)


text:
"""
Charged Particle in an axisymmetric magnetic field of the form
``B(x,y,z) = (1 + x^2 + y^2) e_z``.
"""
module SymmetricField

    import ElectromagneticFields.SymmetricQuadratic

    export pauli_particle_3d_pode, hamiltonian, toroidal_momentum

    SymmetricQuadratic.@code() # inject magnetic field code

    const qᵢ = [1., 0., 0.]
    const vᵢ = [0., 1., 1.]
    const parameters = (μ = 1E-2,)

    include("pauli_particle_3d.jl")

end


text:
using Distributions
using Evolutionary

"""

Return the expected fraction of the dataset processed per iteration. Here, `θs` is the fraction of
the dataset stored, `ps` is the number of partitions, and `Ls` is the probability of finishing
among the `nwait` fastest workers, for each worker.
"""
function fraction_processed(ps::AbstractVector, θs::AbstractVector, Ls::AbstractVector)
    rv = 0.0
    for i in 1:length(Ls)
        rv += Ls[i] * θs[i] / ps[i]
    end
    rv
end

"""Helper function for running event-drive simulations"""
function simulate(;θs, ps, cms_comp, cvs_comp, ds_comm, nwait, niterations=100, nsamples=10, update_latency=0)
    comp_means = cms_comp .* θs ./ ps
    comp_vars = cvs_comp .* θs ./ ps
    ds_comp = DSAGAnalysis.distribution_from_mean_variance.(Gamma, comp_means, comp_vars)
    nworkers = length(ps)
    latency = 0.0
    Ls = zeros(nworkers)
    for _ in 1:nsamples
        dfi, Lsi = simulate_iterations(;nwait, niterations, ds_comm, ds_comp, update_latency)
        latency += mean(dfi.latency)
        Ls .+= Lsi
    end
    latency /= nsamples
    Ls ./= nsamples
    latency, Ls
end

"""

Minimize the variance of the vector `Ls .* θs ./ ps` using a genetic optimization algorithm.
"""
function balance_contribution(ps0::AbstractVector{<:Real}, min_processed_fraction::Real; θs::AbstractVector{<:Real}, comp_means, comp_vars, comm_means, comm_vars, nwait::Integer, populationSize::Integer=100, tournamentSize::Integer=10, mutationRate::Real=1.0, time_limit::Real=10.0)
    nworkers = length(ps0)
    length(θs) == nworkers || throw(DimensionMismatch("θs has dimension $(length(θs)), but nworkers is $nworkers"))
    length(comp_means) == nworkers || throw(DimensionMismatch("comp_means has dimension $(length(comp_means)), but nworkers is $nworkers"))
    length(comp_vars) == nworkers || throw(DimensionMismatch("comp_vars has dimension $(length(comp_vars)), but nworkers is $nworkers"))
    length(comm_means) == nworkers || throw(DimensionMismatch("comm_means has dimension $(length(comp_means)), but nworkers is $nworkers"))
    length(comm_vars) == nworkers || throw(DimensionMismatch("comm_vars has dimension $(length(comp_vars)), but nworkers is $nworkers"))
    all((x)->0<x<=1, θs) || throw(ArgumentError("The entries of θs must be in (0, 1], but got $θs"))
    all((x)->0<x, ps0) || throw(ArgumentError("The entries of ps must be positive, but got $ps0"))

    # initialization
    @info "comp" comp_means comp_vars
    @info "comm" comm_means comm_vars
    ds_comm = DSAGAnalysis.distribution_from_mean_variance.(Gamma, comm_means, comm_vars)
    cms_comp = comp_means ./ (θs ./ ps0)
    cvs_comp = comp_vars ./ (θs ./ ps0)

    # constraint
    # (the total expected conbtribution must be above some threshold)
    function g(ps, Ls)
        fraction_processed(ps, θs, Ls) - min_processed_fraction
    end

    # objective function
    # (the variance of the contribution between workers)
    fworst = 0.0
    function f(ps)
        _, Ls = simulate(;θs, ps, cms_comp, cvs_comp, ds_comm, nwait)
        c = g(ps, Ls)
        if c < 0
            rv = fworst + abs(c)
        else
            rv = var(Ls .* θs ./ ps)
            fworst = max(fworst, rv)
        end
        rv
    end

    # evolutionary algorithm setup
    nworkers = length(ps0)
    selection = Evolutionary.tournament(tournamentSize)
    crossover = Evolutionary.LX()
    lower = max.(ones(Int, nworkers), ceil.(Int, ps0 ./ 2))
    upper = 2 .* ps0
    mutation = Evolutionary.domainrange((lower .- upper) ./ 10) # as recommended in the BGA paper

    # integer-output mutation that wraps another mutation
    function integer_mutation(x)
        if isone(mutationRate) || rand() < mutationRate
            mutation(x)
        end
        for i in 1:length(x)
            if rand() < 0.5
                x[i] = floor(x[i])
            else
                x[i] = ceil(x[i])
            end
        end
        x
    end

    # optimization algorithm
    opt = Evolutionary.GA(;populationSize, mutationRate=1.0, selection, crossover, mutation=integer_mutation)
    options = Evolutionary.Options(;time_limit, Evolutionary.default_options(opt)...)
    Evolutionary.optimize(f, lower, upper, float.(ps0), opt, options)
end

function load_balancer(chin::Channel; chout::Channel, ps0::AbstractVector, θs, min_processed_fraction::Real, nwait::Integer)
    @info "load_balancer task started"
    nworkers = length(ps0)
    0 < nworkers || throw(ArgumentError("nworkers must be positive, but is $nworkers"))
    length(θs) == nworkers || throw(DimensionMismatch("θs has dimension $(length(θs)), but nworkers is $nworkers"))
    ps = Vector{Int}(ps0)
    comp_means = fill(NaN, nworkers)
    comp_vars = fill(NaN, nworkers)
    comm_means = fill(NaN, nworkers)
    comm_vars = fill(NaN, nworkers)

    function process_sample(v)
        0 < v.worker <= nworkers || throw(ArgumentError("v.worker is $(v.worker), but nworkers is $nworkers"))
        isnan(v.comp_mean) || (comp_means[v.worker] = v.comp_mean)
        isnan(v.comp_var) || (comp_vars[v.worker] = v.comp_var)
        isnan(v.comm_mean) || (comm_means[v.worker] = v.comm_mean)
        isnan(v.comm_var) || (comm_vars[v.worker] = v.comm_var)
        return
    end

    # helper to check if there is any missing latency data
    all_populated = false
    function check_populated()
        if all_populated
            return all_populated
        end
        all_populated = iszero(count(isnan, comp_means))
        if all_populated
            all_populated = all_populated && iszero(count(isnan, comp_vars))
        end
        if all_populated
            all_populated = all_populated && iszero(count(isnan, comm_means))
        end
        if all_populated
            all_populated = all_populated && iszero(count(isnan, comm_vars))
        end
        all_populated
    end

    while isopen(chin)

        # consume all values currently in the channel
        try
            vin = take!(chin)
            process_sample(vin)
        catch e
            if e isa InvalidStateException
                break
            else
                rethrow()
            end
        end
        while isready(chin)
            try
                vin = take!(chin)
                process_sample(vin)
            catch e
                if e isa InvalidStateException
                    break
                else
                    rethrow()
                end
            end
        end

        # verify that we have complete latency information for all workers
        if !check_populated()
            continue
        end

        # new = balance_contribution(ps, min_processed_fraction; θs, ds_comm, cms_comp, cvs_comp, nwait)
        @info "started load-balancing optimization"
        result = balance_contribution(ps, min_processed_fraction; θs, comp_means, comp_vars, comm_means, comm_vars, nwait)
        new = result.minimizer

        # push any changes into the output channel
        for i in 1:nworkers
            if new[i] != ps[i]
                vout = @NamedTuple{worker::Int,p::Int}((i, new[i]))
                @info vout
                push!(chout, vout)
                ps[i] = new[i]
            end
        end
    end
    @info "load_balancer task finished"
    return
end


text:
# functions to call out to OSRM. See osrmjl.cpp for the C++ side.
# Since Julia can only call plain C functions, not C++, we have a thin C++ wrapper with
# extern "C" functions to allocate a new OSRM engine, compute a distance matrix, and destroy
# the OSRM engine when you're done with it.

# To use the OSRM module, libosrmjl.so (or libosrmjl.dylib on macOS) must be in your LD_LIBRARY_PATH,
# which likely means you have to build it first, and may need to build OSRM as well. To build libosrmjl,
# run cmake .. && cmake --build . in the cxx/build directory, and put the resulting shared object file somewhere.
# You may need to build OSRM from source so you have the requisitie libraries available. osrm-backend from Homebrew
# is not compatible with Boost from Homebrew.

module OSRM

struct Coordinate
    latitude::Float64
    longitude::Float64
end

const osrmjl = "libosrmjl"

mutable struct OSRMInstance
    _engine::Ptr{Any}
    file_path::String
    algorithm::String
    running::Bool
end

# Start OSRM, with the file path to an already built OSRM graph, and an algorithm
# specification which is mld for multi-level Dijkstra, and ch for contraction hierarchies.
function start_osrm(file_path::String, algorithm::String)::OSRMInstance
    algorithm = lowercase(algorithm)

    if (algorithm != "mld" && algorithm != "ch")
        error("Algorithm must be 'mld' for Multi-Level Dijkstra, or 'ch' for Contraction Hierarchies.")
    end

    ptr = @ccall osrmjl.init_osrm(file_path::Cstring, algorithm::Cstring)::Ptr{Any}
    return OSRMInstance(ptr, file_path, algorithm, true)
end

function distance_matrix(osrm::OSRMInstance, origins::Vector{Coordinate}, destinations::Vector{Coordinate})
    if !osrm.running
        error("OSRM is not running!")
    end

    n_origins::Csize_t = length(origins)
    n_destinations::Csize_t = length(destinations)
    origin_lats::Vector{Float64} = map(c -> c.latitude, origins)
    origin_lons::Vector{Float64} = map(c -> c.longitude, origins)
    destination_lats::Vector{Float64} = map(c -> c.latitude, destinations)
    destination_lons::Vector{Float64} = map(c -> c.longitude, destinations)

    durations::Array{Float64, 2} = fill(-1.0, (n_origins, n_destinations))::Array{Float64, 2}
    distances::Array{Float64, 2} = fill(-1.0, (n_origins, n_destinations))::Array{Float64, 2}

    stat = @ccall osrmjl.distance_matrix(
        osrm._engine::Ptr{Any},
        n_origins::Csize_t,
        origin_lats::Ptr{Float64},
        origin_lons::Ptr{Float64},
        n_destinations::Csize_t,
        destination_lats::Ptr{Float64},
        destination_lons::Ptr{Float64},
        durations::Ptr{Float64},
        distances::Ptr{Float64}
    )::Cint

    if (stat != 0)
        error("osrm failed, status $stat")
    end

    return (durations=durations, distances=distances)
end

function stop_osrm!(osrm::OSRMInstance)
    if osrm.running
        @ccall osrmjl.stop_osrm(osrm._engine::Ptr{Any})::Cvoid
        osrm.running = false
    else
        @warn "stop_osrm! called on already stopped OSRM instance"
    end
end

export start_osrm, distance_matrix, stop_osrm!, Coordinate, OSRMInstance

end


text:
mutable struct HopOnOffSingleCarSimulator
    timestep::Float64
    time_to_finish::Distributions.Normal{Float64}
    params::Parameters
    rng::AbstractRNG
end

function HopOnOffSingleCarSimulator(params::Parameters, rng::AbstractRNG=Random.GLOBAL_RNG)
    return HopOnOffSingleCarSimulator(params.time_params.MDP_TIMESTEP,
                                      Distributions.Normal(0.0,params.time_params.CAR_TIME_STD), params, rng)
end

function reset_sim(sim::HopOnOffSingleCarSimulator,meanval::Float64,stdval::Float64)
    sim.time_to_finish = Distributions.Normal(meanval,stdval)
end

function reset_sim(sim::HopOnOffSingleCarSimulator)
    # randomly sample the initial time for car to reach waypoint
    val = rand(sim.rng,Uniform(0.5*sim.params.time_params.MDP_TIMESTEP*sim.params.time_params.HORIZON_LIM,
                               1.25*sim.params.time_params.MDP_TIMESTEP*sim.params.time_params.HORIZON_LIM))
    sim.time_to_finish = Distributions.Normal(val,sim.params.time_params.CAR_TIME_STD)
end

function sample_finish_time(sim::HopOnOffSingleCarSimulator)
    return rand(sim.rng,sim.time_to_finish)
end

function step_sim(sim::HopOnOffSingleCarSimulator)

    # These params are just for testing simulator
    DELAY_SPEEDUP_PROB = 0.4
    MAX_DELAY_SPEEDUP = sim.timestep*2.0

    old_mean = mean(sim.time_to_finish)
    car_std = std(sim.time_to_finish)

    # Split 1 into 3 bins - 0.0: delay_speedup_prob/2 : 1.0-dsp/2 : 1.0
    probval = rand(sim.rng)
    if probval < DELAY_SPEEDUP_PROB/2.0
        is_delay = -1
    elseif probval < 1.0 - DELAY_SPEEDUP_PROB/2.0
        is_delay = 0
    else # probval > 1.0 - dsp/2.0
        is_delay = 1
    end

    if is_delay == -1
        # Delay
        delay = rand(sim.rng,Uniform(0.0,MAX_DELAY_SPEEDUP))
        sim.time_to_finish = Distributions.Normal(old_mean+delay, car_std)
    elseif is_delay == 0
        # Normal
        sim.time_to_finish = Distributions.Normal(old_mean - sim.timestep, car_std)
    else
        # Speedup
        # Can't sample more than time left
        max_speedup_val = min(old_mean, MAX_DELAY_SPEEDUP)
        speedup = rand(sim.rng,Uniform(0.0,max_speedup_val))

        if sim.timestep + speedup > old_mean
            sim.time_to_finish = Distributions.Normal(0.0,car_std)
        else
            sim.time_to_finish = Distributions.Normal(old_mean - (sim.timestep+speedup), car_std)
        end
    end

    # If too little time left, car is at end
    is_done =  (mean(sim.time_to_finish) < std(sim.time_to_finish))

    return is_done
end


text:
"""
Draw a ```tile``` using ```Luxor.jl```.
If the ```tile``` has no subtiles, it colors it with a random color from ```palette```.

This function does not set up any ```Luxor.Drawing()``` environment by itself.
Regard it more as a ```Luxor.rect()``` function.

Optional arguments:
    - ```origin```: a ```Tuple``` of two numbers two specify the coordinates of the top-left corner of the tile. Defaults to ```(0,0)```;
    - ```palette```: a ```Vector{Vector{Int}}``` containing rgba values of the colors to use to paint the ```tile```. Defaults to ```palette_StampingItOut```.
"""
function draw(tile::TileNode; origin=(0,0), palette=palette_StampingItOut)
    col = rand(palette/255)
    setcolor(col...)
    rect(origin[1], origin[2], tile.b, tile.h, :fillstroke)
end
function draw(tile::TileLeaf; origin=(0,0), palette=palette_StampingItOut)
    b = base(tile)
    h = height(tile)
    draw(tile.subtiles[1]; origin=(origin[1],origin[2]), palette=palette)
    if tile.orientation == H
        draw(tile.subtiles[2]; origin=(origin[1]+base(tile.subtiles[1]), origin[2]), palette=palette)
    else
        draw(tile.subtiles[2]; origin=(origin[1], origin[2]+height(tile.subtiles[1])), palette=palette)
    end
end


text:
@testset "DouglasRachford" begin
    # Though this seems a strange way, it is a way to compute the mid point
    M = Sphere(2)
    p = [1., 0., 0.]
    q = [0., 1., 0.]
    r = [0., 0., 1.]
    start = [0.,0.,1.]
    result = geodesic(M, p, q, 0.5)
    F(x) = distance(M, x, p)^2 + distance(M, x, q)^2
    prox1 = (η,x) -> prox_distance(M, η, p, x)
    prox2 = (η,x) -> prox_distance(M, η, q, x)
    @test_throws ErrorException DouglasRachford(M,F,Array{Function,1}([prox1,]),start) # we need more than one prox
    xHat = DouglasRachford(M,F,[prox1,prox2],start)
    @test F(start) > F(xHat)
    @test_broken distance(M, xHat, result) ≈ 0
    # but we can also compute the riemannian center of mass (locally) on Sn
    # though also this is not that useful, but easy to test that DR works
    F2(x) = distance(M,x,p)^2 + distance(M,x,q)^2 + distance(M,x,r)^2
    prox3 = (η,x) -> prox_distance(M,η,r,x)
    o = DouglasRachford(M,F2,[prox1,prox2,prox3],start;
    debug = [DebugCost(), DebugIterate(), DebugProximalParameter(),100],
    record = [RecordCost(), RecordProximalParameter()],
    return_options=true
    )
    xHat2 = get_solver_result(o)
    drec2 = get_record(o)
    xCmp = 1/sqrt(3)*ones(3)
    # since the default does not run that long -> rough estimate
    @test_broken distance(M,xHat2,xCmp) ≈ 0 atol = 10^(-5)
end


text:
@test isequal(Float64(Float256(42.0)), 42.0)
@test isequal(Float64(Float256(42.0f0)), 42.0)
@test isequal(Float64(Float256(42)), 42.0)
@test isequal(Float64(Float256(0x2a)), 42.0)

@test isequal(Float128(Float256(42.0)), Float128(42.0))
@test isequal(Float128(Float256(42.0f0)), Float128(42.0))
@test isequal(Float128(Float256(42)), Float128(42.0))
@test isequal(Float128(Float256(0x2a)), Float128(42.0))

@test Float64(eps(Float256)) <= 1.0e-60

@test isequal(big(zero(Float256)), big(0))
@test isequal(big(one(Float256)), big(1))
@test abs(big(Float256(π)) - big(π)) <= eps(Float256)
@test abs(big(Float256(ℯ)) - big(ℯ)) <= eps(Float256)

pi256 = f256"3.141592653589793238462643383279502884197169399375105820974944592307816406286198"
@test isequal(big(pi256), Float256(π))
@test isequal(parse(Float256, string(pi256)), pi256)

for n in 1:1000
    x64 = rand(Float64)
    y64 = rand(Float64)

    @test isequal(Float64(Float256(x64)), x64)

    x128 = Float128(rand(BigFloat))
    y128 = Float128(rand(BigFloat))

    @test isequal(Float128(Float256(x128)), x128)

    xbig = rand(BigFloat)
    @test isequal(big(Float256(x64)), big(x64))
    @test isequal(big(Float256(x128)), big(x128))
    @test abs(big(Float256(xbig)) - xbig) <= eps(Float256)

    x = Float256(rand(BigFloat))
    y = Float256(rand(BigFloat))

    @test isequal(big(+ x), + big(x))
    @test isequal(big(- x), - big(x))
    @test isequal(big(abs(x)), abs(big(x)))
    @test approxeq(Float256, 2, [x], big(inv(x)),  inv(big(x)))

    @test abs(big(x + y) - (big(x) + big(y))) <= eps(Float256)
    @test abs(big(x128 + y) - (big(x128) + big(y))) <= eps(Float256)
    @test abs(big(x64 + y) - (big(x64) + big(y))) <= eps(Float256)
    @test abs(big(x + y128) - (big(x) + big(y128))) <= eps(Float256)
    @test abs(big(x + y64) - (big(x) + big(y64))) <= eps(Float256)

    @test abs(big(x - y) - (big(x) - big(y))) <= eps(Float256)
    @test abs(big(x128 - y) - (big(x128) - big(y))) <= eps(Float256)
    @test abs(big(x64 - y) - (big(x64) - big(y))) <= eps(Float256)
    @test abs(big(x - y128) - (big(x) - big(y128))) <= eps(Float256)
    @test abs(big(x - y64) - (big(x) - big(y64))) <= eps(Float256)

    @test abs(big(x * y) - (big(x) * big(y))) <= eps(Float256)
    @test abs(big(x128 * y) - (big(x128) * big(y))) <= eps(Float256)
    @test abs(big(x64 * y) - (big(x64) * big(y))) <= eps(Float256)
    @test abs(big(x * y128) - (big(x) * big(y128))) <= eps(Float256)
    @test abs(big(x * y64) - (big(x) * big(y64))) <= eps(Float256)

    @test approxeq(Float256, 2, [x, y], big(x / y), big(x) / big(y))
    @test approxeq(Float256, 2, [x, y], big(x128 / y), big(x128) / big(y))
    @test approxeq(Float256, 2, [x, y], big(x64 / y), big(x64) / big(y))
    @test approxeq(Float256, 2, [x, y], big(x / y128), big(x) / big(y128))
    @test approxeq(Float256, 2, [x, y], big(x / y64), big(x) / big(y64))

    @test approxeq(Float256, 2, [x, y], big(x \ y), big(x) \ big(y))
    @test approxeq(Float256, 2, [x, y], big(x128 \ y), big(x128) \ big(y))
    @test approxeq(Float256, 2, [x, y], big(x64 \ y), big(x64) \ big(y))
    @test approxeq(Float256, 2, [x, y], big(x \ y128), big(x) \ big(y128))
    @test approxeq(Float256, 2, [x, y], big(x \ y64), big(x) \ big(y64))

    @test isequal(x == x, big(x) == big(x))
    @test isequal(x != x, big(x) != big(x))
    @test isequal(x < x, big(x) < big(x))
    @test isequal(x <= x, big(x) <= big(x))
    @test isequal(x > x, big(x) > big(x))
    @test isequal(x >= x, big(x) >= big(x))

    @test isequal(x == y, big(x) == big(y))
    @test isequal(x != y, big(x) != big(y))
    @test isequal(x < y, big(x) < big(y))
    @test isequal(x <= y, big(x) <= big(y))
    @test isequal(x > y, big(x) > big(y))
    @test isequal(x >= y, big(x) >= big(y))

    @test isequal(x64 == y, big(x64) == big(y))
    @test isequal(x64 != y, big(x64) != big(y))
    @test isequal(x64 < y, big(x64) < big(y))
    @test isequal(x64 <= y, big(x64) <= big(y))
    @test isequal(x64 > y, big(x64) > big(y))
    @test isequal(x64 >= y, big(x64) >= big(y))

    @test isequal(x == y64, big(x) == big(y64))
    @test isequal(x != y64, big(x) != big(y64))
    @test isequal(x < y64, big(x) < big(y64))
    @test isequal(x <= y64, big(x) <= big(y64))
    @test isequal(x > y64, big(x) > big(y64))
    @test isequal(x >= y64, big(x) >= big(y64))

    @test isequal(big(floor(1.0e10 * x)), floor(big(1.0e10 * x)))
    @test isequal(big(ceil(1.0e10 * x)), ceil(big(1.0e10 * x)))
    @test isequal(big(round(1.0e10 * x)), round(big(1.0e10 * x)))
    @test isequal(big(trunc(1.0e10 * x)), trunc(big(1.0e10 * x)))

    @test abs(big(x ^ 2) - big(x) ^ 2) <= eps(Float256)
    @test approxeq(Float256, 10, [x], big(x ^ 5), big(x) ^ 5)
    @test approxeq(Float256, 1, [], big(sqrt(abs(x))), sqrt(big(abs(x))))
    @test approxeq(Float256, 10, [abs(x)],
                   big(nroot(abs(x), 5)), big(abs(x)) ^ (big(1)/5))

    @test approxeq(Float256, 10, [x], big(exp(x)), exp(big(x)))
    @test approxeq(Float256, 10, [x], big(exp10(x)), exp10(big(x)))
    @test approxeq(Float256, 10, [x], big(exp2(x)), exp2(big(x)))
    @test approxeq(Float256, 10,
                   [abs(x)], big(log(abs(x))), log(big(abs(x))))
    @test approxeq(Float256, 10,
                   [abs(x)], big(log10(abs(x))), log10(big(abs(x))))
    @test approxeq(Float256, 10,
                   [abs(x)], big(log2(abs(x))), log2(big(abs(x))))

    @test approxeq(Float256, 10, [x], big(sin(x)), sin(big(x)))
    @test approxeq(Float256, 10, [x], big(cos(x)), cos(big(x)))
    @test approxeq(Float256, 10, [x], big(tan(x)), tan(big(x)))
    @test approxeq(Float256, 10, [x], big(sincos(x)[1]), sin(big(x)))
    @test approxeq(Float256, 10, [x], big(sincos(x)[2]), cos(big(x)))

    @test approxeq(Float256, 10, [x], big(asin(x)), asin(big(x)))
    @test approxeq(Float256, 100, [x], big(acos(x)), acos(big(x)))
    @test approxeq(Float256, 10, [x], big(atan(x)), atan(big(x)))
    @test approxeq(Float256, 10, [x, y],
                   big(atan(x, y)), atan(big(x), big(y)))
    @test approxeq(Float256, 10, [x, y],
                   big(hypot(x, y)), hypot(big(x), big(y)))

    @test approxeq(Float256, 10, [x], big(sinh(x)), sinh(big(x)))
    @test approxeq(Float256, 10, [x], big(cosh(x)), cosh(big(x)))
    @test approxeq(Float256, 10, [x], big(tanh(x)), tanh(big(x)))
    @test approxeq(Float256, 10, [x], big(sincosh(x)[1]), sinh(big(x)))
    @test approxeq(Float256, 10, [x], big(sincosh(x)[2]), cosh(big(x)))

    @test approxeq(Float256, 10, [1, x], big(asinh(x)), asinh(big(x)))
    @test approxeq(Float256, 10, [x + 2], big(acosh(x + 2)), acosh(big(x + 2)))
    @test approxeq(Float256, 10, [1, x], big(atanh(x)), atanh(big(x)))
end

n = 100
A256 = rand(Float256, n, n)
B256 = inv(A256)
@test norm(B256 * A256 - I) <= n * norm(A256) * norm(B256) * eps(Float256)

Abig = rand(BigFloat, n, n);
Bbig = inv(Abig);
@test norm(Bbig * Abig - I) <= n * norm(Abig) * norm(Bbig) * eps(BigFloat)