text:
using Distributions, DataFrames, Plots, Turing, StatsPlots, MCMCChains
import Turexas.Environment as e
import Turexas.Agent as a

function step(agent::a.AbstractAgent, env::e.Bandit)
    p = a.policy(agent)
    action = rand(Categorical(p))
    reward = e.step(env, action)
    a.update(agent, action, reward)
    d = DataFrame(action = action, reward = reward)
    return hcat(d, DataFrame(agent.q', :auto))
end

function run(agent::a.AbstractAgent, env::e.Bandit, trial::Int64)
    return reduce(vcat, map(_ -> step(agent, env), 1:trial))
end

# Run MCMC with a small number of samples
function warmup_sampler(result::AbstractDataFrame, model, k::Int64, nchains::Int64 = 4)
    actions = Array{Real, 1}(result.action)
    rewards = Array{Real, 1}(result.reward)
    sample(model(actions, rewards, k), NUTS(), MCMCThreads(), 2, nchains)
end

function fit(result::AbstractDataFrame, model, k::Int64, iter::Int64, nchains::Int64 = 4)
    actions = Array{Real, 1}(result.action)
    rewards = Array{Real, 1}(result.reward)
    sample(model(actions, rewards, k), NUTS(), MCMCThreads(), iter, nchains)
end

# generate an environment which is shared with all simulations
bandit = e.Bandit([0.2, 0.8])
k = length(bandit.p)

# Run simulation and model fitting with Q-learning
@model QLearningModel(actions::Array{Real, 1}, rewards::Array{Real, 1}, k::Int64) = begin
    T = length(actions)
    α ~ Beta(1, 1)
    β ~ Gamma(1, 100)
    agent = a.QLearningAgent(α, β, k)
    for t in 1:T
        action = actions[t]
        reward = rewards[t]
        p = a.policy(agent)
        actions[t] ~ Categorical(p)
        a.update(agent, action, reward)
    end
end

q_agent = a.QLearningAgent(0.05, 2., k)
q_result = run(q_agent, bandit, 500)
# Run MCMC with a small number of samples in advance because it takes a long time to run MCMC the first time.
_ = warmup_sampler(q_result, QLearningModel, k, 4)
chains = fit(q_result, QLearningModel, k, 500, 4)
plot(chains)

# Run simulation and model fitting with particle filter
@model ParticleFilter(actions::Array{Real, 1}, rewards::Array{Real, 1}, k::Int64) = begin
    T = length(actions)
    σ2 ~ truncated(Normal(), 0.001, 5.)
    β ~ Gamma(1, 100)
    agent = a.ParticleFilterAgent(100, σ2, β, k)
    for t in 1:T
        action = actions[t]
        reward = rewards[t]
        p = a.policy(agent)
        actions[t] ~ Categorical(p)
        a.update(agent, action, reward)
    end
end

pf_agent = a.ParticleFilterAgent(100, 0.5, 1., k)
pf_result = run(pf_agent, bandit, 500)
_ = warmup_sampler(pf_result, ParticleFilter, k, 4)
chains = fit(pf_result, ParticleFilter, k, 500, 4)
plot(chains)


text:
module PersistentArrays

export PersistentArray

using Libdl

const SRCDIR = @__DIR__
const PKGDIR = dirname(SRCDIR)
const DEPSDIR = joinpath(PKGDIR, "deps")
const LIBDIR = joinpath(DEPSDIR, "usr", "lib")

# Path to `libpmem`
const libpmem = joinpath(LIBDIR, "libpmem.so")

function __init__()
    # dlopen libpmem so we can `ccall` into it.
    global libpmem
    Libdl.dlopen(libpmem, Libdl.RTLD_GLOBAL)
end

include("lib.jl")
include("array.jl")

end # module

ext: jl


text:
using LinearAlgebra
using StaticArrays
using PyPlot
const Vec2f = SVector{2, Float64}
const Vec4f = SVector{4, Float64}

include("src/geometry/geometry.jl")
include("src/world.jl")
include("src/double_integrator.jl")
include("src/fmt.jl")

obstacle_set = [Rectangle([0.2, 0.3], 0.2, 0.2),
        Rectangle([0.5, 0.5], 0.2, 0.3),
        Rectangle([0.8, 0.3], 0.2, 0.1),
        Rectangle([0.8, 0.6], 0.15, 0.2),
        Rectangle([0.2, 0.7], 0.1, 0.4),
        Rectangle([0.2, 0.7], 0.1, 0.4)]
x_min = [0, 0]
x_max = [1.0, 1.0]
v_min = [-0.5, -0.5]
v_max = [0.5, 0.5]
W = World(x_min, x_max, v_min, v_max, obstacle_set)

s_init = Vec4f([0.1, 0.1, 0.0, 0.0])
s_goal = Vec4f([0.9, 0.9, 0.0, 0.0])

Nsample = 6000
fmt = @time FMTree(s_init, s_goal, Nsample, W)

with_savefig = false
idx_solution = @time solve(fmt, with_savefig)

show(fmt)
for idx in idx_solution
    s0 = fmt.Pset[fmt.parent[idx]]
    s1 = fmt.Pset[idx]
    tau = fmt.time[idx]
    show_trajectory(s0, s1, tau, 20, :blue, 1.5)
end
savefig("./fig/final.png")


text:
module EnvironmentRestore

    export archiveEnvironment, restoreEnvironment, requireGitCommit, requireCleanWorkingTree
    using Pkg

    """
        archiveEnvironment(packages :: Pair{String, String}...)

    Creates a Julia environment that can use your package at a specific commit. The environment
    is defined by the files Project.toml and Manifest.toml.

    Arguments:
    - `packages`: specified in the form `path1 => commitHash1, path2 => commitHash2, ...` where the paths are the file paths to your repositories.

    Keyword arguments are
    - `dir :: String`: the directory in which to store the environment
    - `force :: Bool`: whether or not any existing environment should be overwritten
    """
    function archiveEnvironment(packages :: Pair{String, String}... ;dir=pwd(), force=false)
        projectFiles = ["Project.toml", "Manifest.toml"]
        if !force && !isempty(intersect(readdir(dir), projectFiles))
           print("Directory already contains Julia project files.\nRun with keyword argument ")
           printstyled("force = true", color=:cyan)
           println(" to overwrite.")
           return
        end

        # Start from a copy of the default environment (normally ~/.julia/environments/vX.Y)
        env = dirname(Base.active_project())
        for projFile ∈ projectFiles
            cp(joinpath(env, projFile), joinpath(dir, projFile); force=true)
        end

        # Add packages to the new environment
        try
            Pkg.activate(dir)
            for (path, commit) ∈ packages
                Pkg.add(path=path, rev=commit)
            end
        catch e
            printstyled(stderr, "\nERROR: ", bold=true, color=:red)
            println("Modifying the new environment " * dir * " caused the following error: ")
            throw(e)
        finally
            Pkg.activate(env)
        end

        printstyled(" Successfully archived ", bold=true, color=:green)
        println("environment in ", joinpath(dir, "Manifest.toml"), ".")
    end

    """
        restoreEnvironment(dir=pwd())

    Restore the environment to the specified directory.
    """
    function restoreEnvironment(dir=pwd())
        Pkg.activate(dir)
        Pkg.instantiate(verbose=true)
    end

    """
        requireGitCommit(commit; strict=true)

    Check if the current commit matches the given (full) commit hash.
    If there is a match, this does nothing.
    Otherwise, if strict = true, this causes an error that interrupts the program.
    If strict = false, this asks the user for permission to continue instead of showing an error.
    """
    function requireGitCommit(commit; strict=true)
        currentCommit = chomp(read(`git rev-parse HEAD`, String))
        currentCommit ==  commit && return

        printstyled("┌ Required commit: ", bold=true, color=:yellow)
        println(commit)
        printstyled("└ Active commit:   ", bold=true, color=:yellow)
        println(currentCommit)

        if !strict
            print("Proceed anyway? (y/N) ")
            readline() == "y" && return
        end
        println("")
        error("Commit requirement not satisfied.")
    end

    """
        requireCleanWorkingTree(;strict=true)

    Check if the Git working tree is clean.
    If it is clean, this does nothing.
    Otherwise, if strict = true, this causes an error that interrupts the program.
    If strict = false, this asks the user for permission to continue instead of showing an error.
    """
    function requireCleanWorkingTree(;strict=true)
        cmd = `git status --short`
        modifiedFiles = read(cmd, String)
        isempty(modifiedFiles) && return

        printstyled("Working tree is not clean. The following files have been changed:\n", color=:yellow)
        run(cmd)
        if !strict
            print("Proceed anyway? (y/N) ")
            readline() == "y" && return
        end
        println("")
        error("Working tree requirement not satisfied.")
    end

end


text:
include("AbstractMolecularContainer.jl")
include("features.jl")
include("obtain.jl")
include("clone.jl")
include("comparison.jl")


text:
module printing

export set_printing

import Base: show

using Compat
using ..bases, ..states
using ..operators, ..operators_dense, ..operators_sparse
using ..operators_lazytensor, ..operators_lazysum, ..operators_lazyproduct
using ..spin, ..fock, ..nlevel, ..particle, ..subspace, ..manybody, ..sparsematrix

"""
    QuantumOptics.set_printing(; standard_order, rounding_tol)

Set options for REPL output.

# Arguments
* `standard_order=false`: For performance reasons, the order of the tensor
    product is inverted, i.e. `tensor(a, b)=kron(b, a)`. When changing this
    to `true`, the output shown in the REPL will exhibit the correct order.
* `rounding_tol=1e-17`: Tolerance for floating point errors shown in the output.
"""
function set_printing(; standard_order::Bool=_std_order, rounding_tol::Real=_round_tol)
    global _std_order = standard_order
    global _round_tol = rounding_tol
    global machineprecorder = Int32(round(-log10(_round_tol), 0))
    nothing
end
set_printing(standard_order=false, rounding_tol=1e-17)

function show(stream::IO, x::GenericBasis)
    if length(x.shape) == 1
        write(stream, "Basis(dim=$(x.shape[1]))")
    else
        s = replace(string(x.shape), " ", "")
        write(stream, "Basis(shape=$s)")
    end
end

function show(stream::IO, x::CompositeBasis)
    write(stream, "[")
    for i in 1:length(x.bases)
        show(stream, x.bases[i])
        if i != length(x.bases)
            write(stream, " ⊗ ")
        end
    end
    write(stream, "]")
end

function show(stream::IO, x::SpinBasis)
    d = denominator(x.spinnumber)
    n = numerator(x.spinnumber)
    if d == 1
        write(stream, "Spin($n)")
    else
        write(stream, "Spin($n/$d)")
    end
end

function show(stream::IO, x::FockBasis)
    write(stream, "Fock(cutoff=$(x.N))")
end

function show(stream::IO, x::NLevelBasis)
    write(stream, "NLevel(N=$(x.N))")
end

function show(stream::IO, x::PositionBasis)
    write(stream, "Position(xmin=$(x.xmin), xmax=$(x.xmax), N=$(x.N))")
end

function show(stream::IO, x::MomentumBasis)
    write(stream, "Momentum(pmin=$(x.pmin), pmax=$(x.pmax), N=$(x.N))")
end

function show(stream::IO, x::SubspaceBasis)
    write(stream, "Subspace(superbasis=$(x.superbasis), states:$(length(x.basisstates)))")
end

function show(stream::IO, x::ManyBodyBasis)
    write(stream, "ManyBody(onebodybasis=$(x.onebodybasis), states:$(length(x.occupations)))")
end

function show(stream::IO, x::Ket)
    write(stream, "Ket(dim=$(length(x.basis)))\n  basis: $(x.basis)\n")
    if !_std_order
        Base.showarray(stream, x.data, false; header=false)
    else
        showquantumstatebody(stream, x)
    end
end

function show(stream::IO, x::Bra)
    write(stream, "Bra(dim=$(length(x.basis)))\n  basis: $(x.basis)\n")
    if !_std_order
        Base.showarray(stream, x.data, false; header=false)
    else
        showquantumstatebody(stream, x)
    end
end

function showoperatorheader(stream::IO, x::Operator)
    write(stream, "$(typeof(x).name.name)(dim=$(length(x.basis_l))x$(length(x.basis_r)))\n")
    if bases.samebases(x)
        write(stream, "  basis: ")
        show(stream, basis(x))
    else
        write(stream, "  basis left:  ")
        show(stream, x.basis_l)
        write(stream, "\n  basis right: ")
        show(stream, x.basis_r)
    end
end

show(stream::IO, x::Operator) = showoperatorheader(stream, x)

function showquantumstatebody(stream::IO, x::Union{Ket,Bra})
    #the permutation is used to invert the order A x B = B.data x A.data to A.data x B.data
    perm = collect(length(basis(x).shape):-1:1)
    if length(perm) == 1
        Base.showarray(stream, round.(x.data, machineprecorder), false; header=false)
    else
        Base.showarray(stream,
        round.(permutesystems(x,perm).data, machineprecorder), false; header=false)
    end
end

function permuted_densedata(x::DenseOperator)
    lbn = length(x.basis_l.shape)
    rbn = length(x.basis_r.shape)
    perm = collect(max(lbn,rbn):-1:1)
    #padd the shape with additional x1 subsystems s.t. x has symmetric number of subsystems
    decomp = lbn > rbn ? [x.basis_l.shape; x.basis_r.shape; fill(1,lbn-rbn)] :
                         [x.basis_l.shape; fill(1,rbn-lbn); x.basis_r.shape]

    data = reshape(x.data, decomp...)
    data = permutedims(data, [perm; perm + length(perm)])
    data = reshape(data, length(x.basis_l), length(x.basis_r))

    return round.(data, machineprecorder)
end

function permuted_sparsedata(x::SparseOperator)
    lbn = length(x.basis_l.shape)
    rbn = length(x.basis_r.shape)
    perm = collect(max(lbn,rbn):-1:1)
    #padd the shape with additional x1 subsystems s.t. x has symmetric number of subsystems
    decomp = lbn > rbn ? [x.basis_l.shape; x.basis_r.shape; fill(1,lbn-rbn)] :
                         [x.basis_l.shape; fill(1,rbn-lbn); x.basis_r.shape]

    data = sparsematrix.permutedims(x.data, decomp, [perm; perm + length(perm)])

    return round.(data, machineprecorder)
end

function show(stream::IO, x::DenseOperator)
    showoperatorheader(stream, x)
    write(stream, "\n")
    if !_std_order
        Base.showarray(stream, x.data, false; header=false)
    else
        Base.showarray(stream, permuted_densedata(x), false; header=false)
    end
end

function show(stream::IO, x::SparseOperator)
    showoperatorheader(stream, x)
    if nnz(x.data) == 0
        write(stream, "\n    []")
    else
        if !_std_order
            show(stream, x.data)
        else
            show(stream, permuted_sparsedata(x))
        end
    end
end

function show(stream::IO, x::LazyTensor)
    showoperatorheader(stream, x)
    write(stream, "\n  operators: $(length(x.operators))")
    s = replace(string(x.indices), " ", "")
    write(stream, "\n  indices: $s")
end

function show(stream::IO, x::Union{LazySum, LazyProduct})
    showoperatorheader(stream, x)
    write(stream, "\n  operators: $(length(x.operators))")
end

end # module


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Private functions and variables of the Jacchia-Bowman 2008 Atmospheric
#   Model, a product of the Space Environment Technologies.
#
#   The source code is available at the following git:
#
#       http://sol.spacenvironment.net/jb2008/code.html
#
#   For more information about the model, see:
#
#       http://sol.spacenvironment.net/~JB2006/
#       http://sol.spacenvironment.net/~JB2008/
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Bowman, B. R., Tobiska, W. K., Marcos, F. A., Huang, C. Y., Lin, C. S.,
#       Burke, W. J (2008). A new empirical thermospheric density model JB2008
#       using new solar and geomagnetic indices. AIAA/AAS Astrodynamics
#       Specialist Conference, Honolulu, Hawaii.
#
#   [2] Bowman, B. R., Tobiska, W. K., Marcos, F. A., Valladares, C (2007). The
#       JB2006 empirical thermospheric density model. Journal of Atmospheric and
#       Solar-Terrestrial Physics, v. 70, p. 774-793.

#   [3] Jacchia, L. G (1970). New static models of the thermosphere and
#       exosphere with empirical temperature profiles. SAO Special Report #313.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

################################################################################
#                              Private Constants
################################################################################

# Constants to compute the `ΔTc` correction.
const _jb2008_B =
    (-0.457512297e+01, -0.512114909e+01, -0.693003609e+02,  0.203716701e+03,
      0.703316291e+03, -0.194349234e+04,  0.110651308e+04, -0.174378996e+03,
      0.188594601e+04, -0.709371517e+04,  0.922454523e+04, -0.384508073e+04,
     -0.645841789e+01,  0.409703319e+02, -0.482006560e+03,  0.181870931e+04,
     -0.237389204e+04,  0.996703815e+03,  0.361416936e+02,)

const _jb2008_C =
    (-0.155986211e+02, -0.512114909e+01, -0.693003609e+02,  0.203716701e+03,
      0.703316291e+03, -0.194349234e+04,  0.110651308e+04, -0.220835117e+03,
      0.143256989e+04, -0.318481844e+04,  0.328981513e+04, -0.135332119e+04,
      0.199956489e+02, -0.127093998e+02,  0.212825156e+02, -0.275555432e+01,
      0.110234982e+02,  0.148881951e+03, -0.751640284e+03,  0.637876542e+03,
      0.127093998e+02, -0.212825156e+02,  0.275555432e+01)

# F(z) global model values, 1997 - 2006 fit.
const _jb2008_fzm = (0.2689e+00,-0.1176e-01, 0.2782e-01,-0.2782e-01,0.3470e-03)

# G(t) global model values, 1997 - 2006 fit.
const _jb2008_gtm = (-0.3633e+00, 0.8506e-01, 0.2401e+00,-0.1897e+00,
                     -0.2554e+00,-0.1790e-01, 0.5650e-03,-0.6407e-03,
                     -0.3418e-02,-0.1252e-02)

# Coefficients for high altitude density correction.
const _jb2008_cht = (0.22e0, -0.20e-02, 0.115e-02, -0.211e-05)

################################################################################
#                              Private Functions
################################################################################

"""
    function _jb2008_grav(z::R) where R

Compute the gravity [m/s] at altitude `z` [km] according to the model Jacchia
1971 [3].

"""
@inline function _jb2008_grav(z::R) where R
    # Mean Earth radius [km].
    Re = R(6356.776)

    # Gravity at Earth surface [m/s²].
    g₀ = R(9.80665)

    # Gravity at desired altitude [m/s²].
    g₀*(1 + z/Re)^(-2)
end

"""
    function _jb2008_highaltitude(h::Number, F10ₐ::Number)

Compute the high altitude exospheric density correction factor in altitude `h`
[km] and the averaged 10.7-cm solar flux (81-day centered on input time)
\\[10⁻²² W/(M² Hz)].

This function uses the model in Section 6.2 of [2].

"""
function _jb2008_highaltitude(h::Number, F10ₐ::Number)
    # Auxiliary variables.
    C = _jb2008_cht

    # Compute the high-altitude density correction.
    FρH = 1.0
    @inbounds if 1000 <= h <= 1500
        z = (h - 1000)/500

        # In this case, the `F₁₅₀₀` is the density factor at 1500 km.
        F₁₅₀₀ = C[1] + C[2]*F10ₐ + C[3]*1500 + C[4]*1500*F10ₐ

        ∂F₁₅₀₀∂z = 500( C[3] + C[4]*F10ₐ )

        # In [2], there is an error when computing this coefficients (eq. 11).
        # The partial derivative has the `500` factor and the coefficients have
        # other `500` factor that **is not** present in the source-code.
        c0  = +1
        c1  = +0
        c2  = +3F₁₅₀₀ - ∂F₁₅₀₀∂z - 3
        c3  = -2F₁₅₀₀ + ∂F₁₅₀₀∂z + 2
        FρH = @evalpoly(z, c0, c1, c2, c3)
    elseif h > 1500
        FρH = C[1] + C[2]*F10ₐ + C[3]*h + C[4]*h*F10ₐ
    end

    FρH
end

"""
    function _jb2008_int(z₀::Number, z₁::Number, R::Number, Tx::Number, T∞::Number, δf::Function)

Compute the integral of the function `δf` between `z₀` and `z₁` using the
Newton-Cotes 4th degree method. `R` is a number that defines the step size, `Tx`
is the temperature at the inflection point, and `T∞` is the exospheric
temperature.

The signature of the function `δf` is:

    δf(z, Tx, T∞)

and it must be `_jb2008_δf1` or `_jb2008_δf2`.

This function returns a tuple containing the integral and last value of `z` used
in the numerical algorithm.

"""
function _jb2008_int(z₀::Number, z₁::Number, R::Number, Tx::Number, T∞::Number,
                     δf::Function)
    # Weights for the Newton-Cotes Five-Point Quad. Formula.
    WT = (14/45, 64/45, 24/45, 64/45, 14/45)

    # Compute the number of integration steps.
    #
    # This is computed so that `z₂ = z₁*(zr)^n`. Hence, `zr` is the factor that
    # defines the size of each integration interval.
    al = log(z₁/z₀)
    n  = floor(Int64, al/R) + 1
    zr = exp(al/n)

    # Initialize the integration auxiliary variables.
    zi₁ = z₀
    zj  = 0.0
    Mbj = 0.0
    Tlj = 0.0
    int = 0.0 # This variable stores the integral from `z₀` to `z₁`.

    # For each integration step, use the Newton-Cotes 4th degree formula to
    # integrate (Boole's rule).
    @inbounds for i = 1:n
        zi₀   = zi₁                   # The beginning of the i-th integration step.
        zi₁   = zr*zi₁                # The end of the i-th integration step.
        Δz    = (zi₁ - zi₀)/4         # Step for the i-th integration step.
        int_i = WT[1]*δf(zi₀, Tx, T∞) # First term of the Newton-Cotes 4th degree sum.

        # Compute the Newton-Cotes 4th degree sum.
        for j = 2:5
            zj  = zi₀ + (j-1)*Δz

            # Value of the integrand at `zj`.
            int_i += WT[j] * δf(zj, Tx, T∞)
        end
        int += int_i * Δz
    end

    int, zj
end

"""
    function _jb2008_M(z::R) where R

Compute the mean molecular mass at altitude `z` [km] using the empirical profile
in eq. 1 [3].

"""
@inline function _jb2008_M(z::R) where R
    !(90 <= z < 105.1) && @warn "The empirical model for the mean molecular mass is valid only for 90 <= z <= 105 km."

    @evalpoly(z - 100, +28.15204,
                        -0.085586,
                        +1.2840e-4,
                        -1.0056e-5,
                        -1.0210e-5,
                        +1.5044e-6,
                        +9.9826e-8)
end

"""
    function _jb2008_semiannual(doy::Number, h::Number, F10ₐ::Number, S10ₐ::Number, M10ₐ::Number)

Compute the semiannual variation of the density considering the JB2008 model
[1].

# Args

* `doy`: Day of the year + fraction of the day.
* `h`: Height [km].
* `F10ₐ`: Averaged 10.7-cm flux (81-day centered on input-time)
          \\[10⁻²² W/(M² Hz)].
* `S10ₐ`: EUV 81-day averaged centered index.
* `M10ₐ`: MG2 81-day averaged centered index.

# Returns

* Semiannual F(z) heigh function.
* Semiannual G(t) yearly periodic function.
* Semiannual variation of the density `Δsalog₁₀ρ`.

"""
function _jb2008_semiannual(doy::Number, h::Number, F10ₐ::Number, S10ₐ::Number,
                            M10ₐ::Number)
    # Auxiliary variables.
    B = _jb2008_fzm
    C = _jb2008_gtm
    z = h/1000
    ω = 2π*(doy-1)/365  # See eq. 4 [2].

    # Compute the new 81-day centered solar index for F(z) according to eq. 4
    # [1].
    Fsmjₐ = 1F10ₐ - 0.7S10ₐ - 0.04M10ₐ

    # Compute the semiannual F(z) height function according to eq. 5 [1].
    @inbounds Fz = B[1] + B[2]*Fsmjₐ + B[3]*z*Fsmjₐ + B[4]*z^2*Fsmjₐ + B[5]*z*Fsmjₐ^2

    # Compute the new 81-day centered solar index for G(t) according to eq. 6
    # [1].
    Fsmₐ = 1F10ₐ - 0.75S10ₐ - 0.37M10ₐ

    # Compute the semiannual G(t) yearly periodic function according to eq. 7
    # [1].
    sω,  cω  = sincos(1ω)
    s2ω, c2ω = sincos(2ω)

    @inbounds Gt =        C[1] + C[2]*sω + C[3]*cω + C[4]*s2ω + C[5]*c2ω +
                   Fsmₐ*( C[6] + C[7]*sω + C[8]*cω + C[9]*s2ω + C[10]*c2ω )

    (Fz < 1e-6) && (Fz = 1e-6)

    Δsalog₁₀ρ = Fz*Gt

    (Fz, Gt, Δsalog₁₀ρ)
end

"""
    function _jb2008_T(z::R, Tx::R, T∞::R) where R<:Number

Compute the temperature [K] at height `z` [km] given the temperature `Tx` [K] at
the inflection point, and the exospheric temperature `T∞` [K] according to the
theory of the model Jacchia 1971 [3].

The inflection point is considered to by `z = 125 km`.

"""
function _jb2008_T(z::R, Tx::R, T∞::R) where R<:Number
    # Constants
    # =========

    T₁  = 183    # Temperature at the lower bound [K].
    z₁  = 90     # Altitude of the lower bound [km].
    zx  = 125    # Altitude of the inflection point [km].
    Δz₁ = z₁-zx

    # Check the parameters
    # ====================
    (z  < z₁) && error("The altitude must not be lower than $(z₁) km.")
    (T∞ < 0 ) && error("The exospheric temperature must be positive.")

    # Compute the temperature gradient at the inflection point.
    Gx = R(1.9)*( Tx - T₁ )/( zx - z₁ )

    # Compute the temperature at desire altitude
    # ==========================================

    Δz = z - zx

    if z <= zx
        c₁ = Gx
        c₂ = R(0)
        c₃ = R(-5.1)/(R(5.7)*Δz₁^2)*Gx
        c₄ = R( 0.8)/(R(1.9)*Δz₁^3)*Gx
        T  = @evalpoly(Δz, Tx, c₁, c₂, c₃, c₄)
    else
        A  = R(2/π) * ( T∞ - Tx )
        T  = Tx + A*atan( Gx/A * Δz * ( 1 + R(4.5e-6) * Δz^2.5 ) )
    end

    T
end

"""
    function _jb2008_ΔTc(F10::Number, lst::Number, glat::Number, h::Number)

Compute the correction in the `Tc` for Jacchia-Bowman model.

This correction is mention in [2]. However, the equations do not seem to match
those in the source-code. The ones implemented here are exactly the same as in
the source-code.

# Args

* `F10`: F10.7 flux.
* `lst`: Local solar time (0 - 24) \\[hr].
* `glat`: Geocentric latitude [rad].
* `h`: Altitude [km].

# Returns

The correction `ΔTc` [K].

"""
function _jb2008_ΔTc(F10::Number, lst::Number, glat::Number, h::Number)
    # Auxiliary variables according to [2, p.  784].
    B  = _jb2008_B
    C  = _jb2008_C
    F  = (F10 - 100)/100
    θ  = lst/24
    θ² = θ^2
    θ³ = θ^3
    θ⁴ = θ^4
    θ⁵ = θ^5
    ϕ  = cos(glat)

    ΔTc = 0.0

    # Compute the temperature variation given the altitude.
    @inbounds if 120 <= h <= 200
        ΔTc200 = C[17]     + C[18]*θ*ϕ   + C[19]*θ²*ϕ   + C[20]*θ³*ϕ +
                 C[21]*F*ϕ + C[22]*θ*F*ϕ + C[23]*θ²*F*ϕ

        ΔTc200Δz = C[1]      + B[2]*F     + C[3]*θ*F    + C[4]*θ²*F    +
                   C[5]*θ³*F + C[6]*θ⁴*F  + C[7]*θ⁵*F   + C[8]*θ*ϕ     +
                   C[9]*θ²*ϕ + C[10]*θ³*ϕ + C[11]*θ⁴*ϕ  + C[12]*θ⁵*ϕ   +
                   C[13]*ϕ   + C[14]*F*ϕ  + C[15]*θ*F*ϕ + C[16]*θ²*F*ϕ

        zp  = (h - 120)/80
        ΔTc = (3ΔTc200 - ΔTc200Δz)*zp^2 + (ΔTc200Δz - 2ΔTc200)*zp^3

    elseif 200 < h <= 240
        H = (h - 200)/50

        ΔTc = C[ 1]*H      + B[2]*F*H     + C[ 3]*θ*F*H   + C[ 4]*θ²*F*H   +
              C[ 5]*θ³*F*H + C[6]*θ⁴*F*H  + C[ 7]*θ⁵*F*H  + C[ 8]*θ*ϕ*H    +
              C[9]*θ²*ϕ*H  + C[10]*θ³*ϕ*H + C[11]*θ⁴*ϕ*H  + C[12]*θ⁵*ϕ*H   +
              C[13]*ϕ*H    + C[14]*F*ϕ*H  + C[15]*θ*F*ϕ*H + C[16]*θ²*F*ϕ*H +
              C[17]        + C[18]*θ*ϕ    + C[19]*θ²*ϕ    + C[20]*θ³*ϕ     +
              C[21]*F*ϕ    + C[22]*θ*F*ϕ  + C[23]*θ²*F*ϕ

    elseif 240 < h <= 300
        H = 40/50

        aux1 = C[ 1]*H      + B[ 2]*F*H    + C[ 3]*θ*F*H   + C[ 4]*θ²*F*H   +
               C[ 5]*θ³*F*H + C[ 6]*θ⁴*F*H + C[ 7]*θ⁵*F*H  + C[ 8]*θ*ϕ*H    +
               C[ 9]*θ²*ϕ*H + C[10]*θ³*ϕ*H + C[11]*θ⁴*ϕ*H  + C[12]*θ⁵*ϕ*H   +
               C[13]*ϕ*H    + C[14]*F*ϕ*H  + C[15]*θ*F*ϕ*H + C[16]*θ²*F*ϕ*H +
               C[17]        + C[18]*θ*ϕ    + C[19]*θ²*ϕ    + C[20]*θ³*ϕ     +
               C[21]*F*ϕ    + C[22]*θ*F*ϕ  + C[23]*θ²*F*ϕ

        aux2 = C[ 1]        + B[ 2]*F    + C[ 3]*θ*F + C[ 4]*θ²*F + C[ 5]*θ³*F  +
               C[ 6]*θ⁴*F   + C[ 7]*θ⁵*F + C[ 8]*θ*ϕ + C[ 9]*θ²*ϕ + C[10]*θ³*ϕ  +
               C[11]*θ⁴*ϕ   + C[12]*θ⁵*ϕ + C[13]*ϕ   + C[14]*F*ϕ  + C[15]*θ*F*ϕ +
               C[16]*θ²*F*ϕ

        H = 300/100

        ΔTc300 = B[ 1]        + B[ 2]*F      + B[ 3]*θ*F    + B[ 4]*θ²*F   +
                 B[ 5]*θ³*F   + B[ 6]*θ⁴*F   + B[ 7]*θ⁵*F   + B[ 8]*θ*ϕ    +
                 B[ 9]*θ²*ϕ   + B[10]*θ³*ϕ   + B[11]*θ⁴*ϕ   + B[12]*θ⁵*ϕ   +
                 B[13]*H*ϕ    + B[14]*θ*H*ϕ  + B[15]*θ²*H*ϕ + B[16]*θ³*H*ϕ +
                 B[17]*θ⁴*H*ϕ + B[18]*θ⁵*H*ϕ + B[19]*ϕ

        ΔTc300Δz = B[13]*ϕ + B[14]*θ*ϕ + B[15]*θ²*ϕ + B[16]*θ³*ϕ + B[17]*θ⁴*ϕ +
                   B[18]*θ⁵*ϕ

        aux3 = 3ΔTc300 - ΔTc300Δz - 3aux1 - 2aux2
        aux4 = ΔTc300 - aux1 - aux2 - aux3
        zp   = (h - 240)/60
        ΔTc  = aux1 + aux2*zp + aux3*zp^2 + aux4*zp^3

    elseif 300 < h <= 600
         H  = h/100
        ΔTc = B[ 1]        + B[ 2]*F      + B[ 3]*θ*F    + B[ 4]*θ²*F   +
              B[ 5]*θ³*F   + B[ 6]*θ⁴*F   + B[ 7]*θ⁵*F   + B[ 8]*θ*ϕ    +
              B[ 9]*θ²*ϕ   + B[10]*θ³*ϕ   + B[11]*θ⁴*ϕ   + B[12]*θ⁵*ϕ   +
              B[13]*H*ϕ    + B[14]*θ*H*ϕ  + B[15]*θ²*H*ϕ + B[16]*θ³*H*ϕ +
              B[17]*θ⁴*H*ϕ + B[18]*θ⁵*H*ϕ + B[19]*ϕ

    elseif 600 < h <= 800
        zp   = (h - 600)/100
        hp   = 600/100
        aux1 = B[ 1]         + B[ 2]*F       + B[ 3]*θ*F     + B[ 4]*θ²*F    +
               B[ 5]*θ³*F    + B[ 6]*θ⁴*F    + B[ 7]*θ⁵*F    + B[ 8]*θ*ϕ     +
               B[ 9]*θ²*ϕ    + B[10]*θ³*ϕ    + B[11]*θ⁴*ϕ    + B[12]*θ⁵*ϕ    +
               B[13]*hp*ϕ    + B[14]*θ*hp*ϕ  + B[15]*θ²*hp*ϕ + B[16]*θ³*hp*ϕ +
               B[17]*θ⁴*hp*ϕ + B[18]*θ⁵*hp*ϕ + B[19]*ϕ

        aux2 = B[13]*ϕ    + B[14]*θ*ϕ + B[15]*θ²*ϕ + B[16]*θ³*ϕ + B[17]*θ⁴*ϕ +
               B[18]*θ⁵*ϕ

        aux3 = -(3aux1 + 4aux2)/4
        aux4 =  ( aux1 +  aux2)/4
        ΔTc  = aux1 + aux2*zp + aux3*zp^2 + aux4*zp^3
    end

    ΔTc
end

"""
    function _jb2008_δf1(z, Tx, T∞)

Auxiliary function to compute the integrand in `_jb2008_int`.

"""
@inline function _jb2008_δf1(z, Tx, T∞)
    Mb = _jb2008_M(z)
    Tl = _jb2008_T(z, Tx, T∞)
    g  = _jb2008_grav(z)
    Mb * g / Tl
end

"""
    function _jb2008_δf2(z, Tx, T∞)

Auxiliary function to compute the integrand in `_jb2008_int`.

"""
@inline function _jb2008_δf2(z, Tx, T∞)
    Tl = _jb2008_T(z, Tx, T∞)
    g  = _jb2008_grav(z)
    g / Tl
end


text:
#__precompile__()
module CustomREPLCompletions

# Function that looks for custom completions for specific function. Look at examples.jl for examples.
function custom_completions end

should_evaluate(ex)=true #Numbers, symbols, QuoteNodes etc
function should_evaluate(ex::Expr)
  if ex.head==:(.)
    should_evaluate(ex.args[2])
  end
  return false
end

function custom_complete_methods(ex_org::Expr)
  should_evaluate(ex_org.args[1]) || return false, UTF8String[]
  func, found = Base.REPLCompletions.get_value(ex_org.args[1], Main)
    (!found || (found && !isa(func,Function))) && return false, UTF8String[]
    args_ex = DataType[Base.REPLCompletions.get_type(ex_org.args[1],Main)[1]]
    for ex in ex_org.args[2:end]
      if should_evaluate(ex)
        val, found = Base.REPLCompletions.get_type(ex, Main)
        push!(args_ex, val)
      else
        push!(args_ex,Void)
      end
    end
    t_in = Tuple{args_ex...} # Input types
    if method_exists(custom_completions,t_in)
        args_val=[should_evaluate(ex_org.args[i]) ? Base.REPLCompletions.get_value(ex_org.args[i],Main)[1] : nothing for i=2:length(args_ex)]
        return true, custom_completions(func,args_val...)
    end
    return false, UTF8String[]
end

#This is a modified version of Base that also finds square brackets
function find_start_brace_or_bracket(s::AbstractString)
    braces = 0
    brackets = 0
    r = RevString(s)
    i = start(r)
    in_single_quotes = false
    in_double_quotes = false
    in_back_ticks = false
    in_brackets  = false
    while !done(r, i)
        c, i = next(r, i)
        if !in_single_quotes && !in_double_quotes && !in_back_ticks
            if c == '('
                braces += 1
            elseif c == ')'
                braces -= 1
            elseif c == '['
                brackets += 1
            elseif c == ']'
                brackets -= 1
            elseif c == '\''
                in_single_quotes = true
            elseif c == '"'
                in_double_quotes = true
            elseif c == '`'
                in_back_ticks = true
            end
        else
            if !in_back_ticks && !in_double_quotes && c == '\''
                in_single_quotes = !in_single_quotes
            elseif !in_back_ticks && !in_single_quotes && c == '"'
                in_double_quotes = !in_double_quotes
            elseif !in_single_quotes && !in_double_quotes && c == '`'
                in_back_ticks = !in_back_ticks
            end
        end
        braces == 1 && break
        brackets == 1 && break
    end
    braces != 1 && brackets!= 1 && return 0:-1, -1, true
    method_name_end = reverseind(r, i)
    startind = nextind(s, rsearch(s, Base.REPLCompletions.non_identifier_chars, method_name_end))
    return startind:endof(s), method_name_end, braces==1 ? true : false
end

#This is nearly the same function as in Base.
function Base.REPLCompletions.completions(string::AbstractString, pos)
    # First parse everything up to the current position
    partial = string[1:pos]
    inc_tag = Base.syntax_deprecation_warnings(false) do
        Base.incomplete_tag(parse(partial, raise=false))
    end

    #Here starts the extra code for this package, we are looking especially for cases where we are either in a function or in a non-closed string.
    #This is the extra check if the string ends with " inside function call
    if inc_tag in [:other,:string]
        if inc_tag==:other
            newpartial=endswith(partial,"\"") ? partial*"," : partial
            newinc_tag=inc_tag
        else
            newpartial=partial*"\","
            newinc_tag = Base.syntax_deprecation_warnings(false) do
                Base.incomplete_tag(parse(newpartial, raise=false))
            end
        end
        frange, method_name_end, is_brace = find_start_brace_or_bracket(newpartial)
        if frange!=0:-1
        ex = Base.syntax_deprecation_warnings(false) do
            is_brace ? parse(newpartial[frange] * ")", raise=false) : parse(newpartial[frange] * "]", raise=false)
        end
        #If ex is a :ref change to getindex call
        if isa(ex,Expr) && ex.head==:ref
          ex.head=:call
          unshift!(ex.args,:getindex)
        end
        if isa(ex, Expr) && ex.head==:call
            success,suggestions=custom_complete_methods(ex)
            if success
                m = match(r"[\t\n\r\"'/`@\$><=;|&\{]| (?!\\)", reverse(partial))
                startpos = nextind(partial, reverseind(partial, m.offset))
                r = startpos:pos
                pshort=partial[r]
                filter!(x->startswith(x,pshort),suggestions)
                isempty(suggestions) || return suggestions, r, true
            end
        end
      end
    end

    if inc_tag in [:cmd, :string]
        m = match(r"[\t\n\r\"'`@\$><=;|&\{]| (?!\\)", reverse(partial))
        startpos = nextind(partial, reverseind(partial, m.offset))
        r = startpos:pos
        paths, r, success = Base.REPLCompletions.complete_path(replace(string[r], r"\\ ", " "), pos)
        if inc_tag == :string &&
           length(paths) == 1 &&                              # Only close if there's a single choice,
           !isdir(expanduser(replace(string[startpos:start(r)-1] * paths[1], r"\\ ", " "))) &&  # except if it's a directory
           (length(string) <= pos || string[pos+1] != '"')    # or there's already a " at the cursor.
            paths[1] *= "\""
        end
        #Latex symbols can be completed for strings
        (success || inc_tag==:cmd) && return sort(paths), r, success
    end

    ok, ret = Base.REPLCompletions.bslash_completions(string, pos)
    ok && return ret

    # Make sure that only bslash_completions is working on strings
    inc_tag==:string && return UTF8String[], 0:-1, false

    if inc_tag == :other && Base.REPLCompletions.should_method_complete(partial)
        frange, method_name_end = Base.REPLCompletions.find_start_brace(partial)
        ex = Base.syntax_deprecation_warnings(false) do
            parse(partial[frange] * ")", raise=false)
        end
        if isa(ex, Expr) && ex.head==:call
            return Base.REPLCompletions.complete_methods(ex), start(frange):method_name_end, false
        end
    elseif inc_tag == :comment
        return UTF8String[], 0:-1, false
    end

    dotpos = rsearch(string, '.', pos)
    startpos = nextind(string, rsearch(string, Base.REPLCompletions.non_identifier_chars, pos))

    ffunc = (mod,x)->true
    suggestions = UTF8String[]
    comp_keywords = true
    if Base.REPLCompletions.afterusing(string, startpos)
        # We're right after using or import. Let's look only for packages
        # and modules we can reach from here

        # If there's no dot, we're in toplevel, so we should
        # also search for packages
        s = string[startpos:pos]
        if dotpos <= startpos
            for dir in [Pkg.dir(); LOAD_PATH; pwd()]
                isdir(dir) || continue
                for pname in readdir(dir)
                    if pname[1] != '.' && pname != "METADATA" &&
                        pname != "REQUIRE" && startswith(pname, s)
                        # Valid file paths are
                        #   <Mod>.jl
                        #   <Mod>/src/<Mod>.jl
                        #   <Mod>.jl/src/<Mod>.jl
                        if isfile(joinpath(dir, pname))
                            endswith(pname, ".jl") && push!(suggestions, pname[1:end-3])
                        else
                            mod_name = if endswith(pname, ".jl")
                                pname[1:end - 3]
                            else
                                pname
                            end
                            if isfile(joinpath(dir, pname, "src",
                                               "$mod_name.jl"))
                                push!(suggestions, mod_name)
                            end
                        end
                    end
                end
            end
        end
        ffunc = (mod,x)->(isdefined(mod, x) && isa(mod.(x), Module))
        comp_keywords = false
    end
    startpos == 0 && (pos = -1)
    dotpos <= startpos && (dotpos = startpos - 1)
    s = string[startpos:pos]
    comp_keywords && append!(suggestions, Base.REPLCompletions.complete_keyword(s))
    append!(suggestions, Base.REPLCompletions.complete_symbol(s, ffunc))
    return sort(unique(suggestions)), (dotpos+1):pos, true
end

include("examples.jl")
@load_examples

end # module


text:
"""
    structure(anatomy)
"""
function structure(a::AbstractAnatomy) end

"""
"""
function structuretype(::AbstractAnatomy) end

"""
    coordinates(anatomy)
"""
function positions(a::AbstractAnatomy) end

"""
    weights
"""
function weights(a::AbstractAnatomy) end

"""
    weighttype(anatomy)
"""
function weighttype(::AbstractAnatomy) end

"""
    in(a::AbstractAnatomy, b::AbstractAnatomy) -> Bool

"""
function Base.in(a::AbstractAnatomy, b::AbstractAnatomy) end

"""
connect() -> Synapse
"""
function GeometryBasics.connect() end


text:
#
# Framing the input context.
#

FrameQuery(q::Query) =
    Query(
        FrameSig(),
        [q],
        input(q) |> setrelative(false),
        output(q))

immutable FrameSig <: AbstractSignature
end

ev(sig::FrameSig, args::Vector{Query}, ity::Input, oty::Output, iflow::InputFlow) =
    ev(sig, args..., ity, oty, iflow)

ev(::FrameSig, arg::Query, ::Input, oty::Output, iflow::InputFlow) =
    let iflow =
        InputFlow(
            iflow.ctx,
            domain(iflow),
            iflow.vals,
            InputFrame(OneTo(iflow.len+1)),
            iflow.slotflows)
        OutputFlow(
            oty,
            column(ev(arg, iflow)))
    end

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
module ReactiveMPTest

using Test, Documenter, ReactiveMP
using TestSetExtensions

using Aqua
Aqua.test_all(ReactiveMP; ambiguities=false)

include("test_helpers.jl")

using .ReactiveMPTestingHelpers

# doctest(ReactiveMP)

@testset ExtendedTestSet "ReactiveMP" begin

    enabled_tests = lowercase.(ARGS)

    function key_to_filename(key)
        splitted = split(key, ":")
        return length(splitted) === 1 ? string("test_", first(splitted), ".jl") : string(join(splitted[1:end - 1], "/"), "/test_", splitted[end], ".jl")
    end

    function filename_to_key(filename)
        splitted   = split(filename, "/")
        if length(splitted) === 1
            return replace(replace(first(splitted), ".jl" => ""), "test_" => "")
        else
            path, name = splitted[1:end - 1], splitted[end]
            return string(join(path, ":"), ":", replace(replace(name, ".jl" => ""), "test_" => ""))
        end
    end

    function addtests(filename)
        key = filename_to_key(filename)
        if isempty(enabled_tests) || key in enabled_tests
            include(filename)
        end
    end

    @testset "Testset helpers" begin
        @test key_to_filename(filename_to_key("distributions/test_normal_mean_variance.jl")) == "distributions/test_normal_mean_variance.jl"
        @test filename_to_key(key_to_filename("distributions:normal_mean_variance")) == "distributions:normal_mean_variance"
        @test key_to_filename(filename_to_key("test_message.jl")) == "test_message.jl"
        @test filename_to_key(key_to_filename("message")) == "message"
    end

    addtests("algebra/test_correction.jl")

    addtests("test_math.jl")

    addtests("constraints/prod/test_prod_final.jl")

    addtests("test_distributions.jl")
    addtests("distributions/test_common.jl")
    addtests("distributions/test_bernoulli.jl")
    addtests("distributions/test_beta.jl")
    addtests("distributions/test_categorical.jl")
    addtests("distributions/test_contingency.jl")
    addtests("distributions/test_dirichlet_matrix.jl")
    addtests("distributions/test_dirichlet.jl")
    addtests("distributions/test_gamma.jl")
    addtests("distributions/test_mv_normal_mean_covariance.jl")
    addtests("distributions/test_mv_normal_mean_precision.jl")
    addtests("distributions/test_mv_normal_weighted_mean_precision.jl")
    addtests("distributions/test_normal_mean_variance.jl")
    addtests("distributions/test_normal_mean_precision.jl")
    addtests("distributions/test_normal_weighted_mean_precision.jl")
    addtests("distributions/test_normal.jl")
    addtests("distributions/test_pointmass.jl")
    addtests("distributions/test_wishart.jl")
    addtests("distributions/test_sample_list.jl")

    addtests("test_message.jl")

    addtests("test_variable.jl")

    addtests("test_node.jl")
    addtests("nodes/test_addition.jl")
    addtests("nodes/test_bifm.jl")
    addtests("nodes/test_bifm_helper.jl")
    addtests("nodes/test_subtraction.jl")
    addtests("nodes/test_probit.jl")

    addtests("rules/addition/test_marginals.jl")
    addtests("rules/addition/test_in1.jl")
    addtests("rules/addition/test_in2.jl")
    addtests("rules/addition/test_out.jl")

    addtests("rules/bifm/test_marginals.jl")
    addtests("rules/bifm/test_in.jl")
    addtests("rules/bifm/test_out.jl")
    addtests("rules/bifm/test_zprev.jl")
    addtests("rules/bifm/test_znext.jl")

    addtests("rules/bifm_helper/test_in.jl")
    addtests("rules/bifm_helper/test_out.jl")

    addtests("rules/subtraction/test_marginals.jl")
    addtests("rules/subtraction/test_in1.jl")
    addtests("rules/subtraction/test_in2.jl")
    addtests("rules/subtraction/test_out.jl")

    addtests("rules/bernoulli/test_out.jl")
    addtests("rules/bernoulli/test_p.jl")
    addtests("rules/bernoulli/test_marginals.jl")

    addtests("rules/beta/test_out.jl")
    addtests("rules/beta/test_marginals.jl")

    addtests("rules/normal_mean_variance/test_out.jl")
    addtests("rules/normal_mean_variance/test_mean.jl")
    addtests("rules/normal_mean_variance/test_var.jl")

    addtests("rules/normal_mean_precision/test_out.jl")
    addtests("rules/normal_mean_precision/test_mean.jl")
    addtests("rules/normal_mean_precision/test_precision.jl")

    addtests("rules/mv_normal_mean_covariance/test_out.jl")
    addtests("rules/mv_normal_mean_covariance/test_mean.jl")

    addtests("rules/mv_normal_mean_precision/test_out.jl")
    addtests("rules/mv_normal_mean_precision/test_mean.jl")
    addtests("rules/mv_normal_mean_precision/test_precision.jl")

    addtests("rules/probit/test_out.jl")
    addtests("rules/probit/test_in.jl")

    addtests("rules/wishart/test_marginals.jl")
    addtests("rules/wishart/test_out.jl")

end

end


text:
mutable struct ROCK2ConstantCache{T,T2,zType} <: OrdinaryDiffEqConstantCache
  ms::SVector{46, Int}
  fp1::SVector{46, T}
  fp2::SVector{46, T}
  recf::Vector{T2}
  zprev::zType
  mdegprev::Int
  mdeg::Int
  recind::Int
end
@cache struct ROCK2Cache{uType,rateType,uNoUnitsType} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  gprev::uType
  gprev2::uType
  tmp::uType
  atmp::uNoUnitsType
  fsalfirst::rateType
  k::rateType
  k2::rateType
  constantcache::ROCK2ConstantCache
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  constantcache = ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u) # WIP: not sure about what type to use in here
  gprev = similar(u)
  gprev2 = similar(u)
  tmp = similar(u)
  atmp = similar(u,uEltypeNoUnits)
  fsalfirst = zero(rate_prototype)
  k = zero(rate_prototype)
  k2 = zero(rate_prototype)
  ROCK2Cache(u, uprev, gprev, gprev2, tmp, atmp, fsalfirst, k, k2, constantcache)
end

function alg_cache(alg::ROCK2,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  ROCK2ConstantCache(uEltypeNoUnits, uEltypeNoUnits, u) # WIP: not sure about what type to use in here
end


text:
# cash_flows.jl
# module CF

using QuantLib.Time, QuantLib.Math

# BlackIborCouponPricer() = BlackIborCouponPricer(0.0, 0.0, false)

type CouponMixin{DC <: DayCount}
  accrualStartDate::Date
  accrualEndDate::Date
  refPeriodStart::Date
  refPeriodEnd::Date
  dc::DC
  accrualPeriod::Float64
end

accrual_start_date(coup::Coupon) = coup.couponMixin.accrualStartDate
accrual_end_date(coup::Coupon) = coup.couponMixin.accrualEndDate
ref_period_start(coup::Coupon) = coup.couponMixin.refPeriodStart
ref_period_end(coup::Coupon) = coup.couponMixin.refPeriodEnd
get_dc(coup::Coupon) = coup.couponMixin.dc

# accural_period{C <: Coupon}(coup::C) = coup.couponMixin.accrualPeriod
accrual_period!(coup::Coupon, val::Float64) = coup.couponMixin.accrualPeriod = val
function accrual_period(coup::Coupon)
  if coup.couponMixin.accrualPeriod == -1.0
    p = year_fraction(get_dc(coup), accrual_start_date(coup), accrual_end_date(coup), ref_period_start(coup), ref_period_end(coup))
    accrual_period!(coup, p)
  end

  return coup.couponMixin.accrualPeriod
end

## types of cash flows ##
type SimpleCashFlow <: CashFlow
  amount::Float64
  date::Date
end

amount(cf::SimpleCashFlow) = cf.amount
date(cf::SimpleCashFlow) = cf.date
date_accrual_end(cf::SimpleCashFlow) = cf.date

date(coup::Coupon) = coup.paymentDate
date_accrual_end(coup::Coupon) = accrual_end_date(coup)

type Dividend <: CashFlow
  amount::Float64
  date::Date
end

amount(div::Dividend) = div.amount
date(div::Dividend) = div.date
date_accrual_end(div::Dividend) = div.date

# legs to build cash flows
abstract Leg <: CashFlows

type ZeroCouponLeg <: Leg
  redemption::SimpleCashFlow
end

## Function wrapper for solvers ##
type IRRFinder{L <: Leg, DC <: DayCount, C <: CompoundingType, F <: Frequency}
  leg::L
  npv::Float64
  dc::DC
  comp::C
  freq::F
  includeSettlementDateFlows::Bool
  settlementDate::Date
  npvDate::Date
end

## this function can pass itself or its derivative ##
function operator(finder::IRRFinder)
  function _inner(y::Float64)
    yld = InterestRate(y, finder.dc, finder.comp, finder.freq)
    NPV = npv(finder.leg, yld, finder.includeSettlementDateFlows, finder.settlementDate, finder.npvDate)
    return finder.npv - NPV
  end

  # derivative
  function _inner(y::Float64, ::Derivative)
    yld = InterestRate(y, finder.dc, finder.comp, finder.freq)
    return duration(ModifiedDuration(), finder.leg, yld, finder.dc, finder.includeSettlementDateFlows, finder.settlementDate, finder.npvDate)
  end

  return _inner
end

get_latest_coupon{L <: Leg}(leg::L) = get_latest_coupon(leg, leg.coupons[end])
get_latest_coupon{L <: Leg}(leg::L, simp::SimpleCashFlow) = leg.coupons[end - 1]
get_latest_coupon{L <: Leg, C <: Coupon}(leg::L, coup::C) = coup

check_coupon{C <: CashFlow}(x::C) = isa(x, Coupon)

get_pay_dates{C <: Coupon}(coups::Vector{C}) = Date[date(coup) for coup in coups]

get_reset_dates{C <: Coupon}(coups::Vector{C}) = Date[accrual_start_date(coup) for coup in coups]

## Pricer Methods ##
# function initialize!(pricer::IborCouponPricer, coupon::IborCoupon)
#   # stuff
#   if !pricer.initialized
#     payment_date = date(coupon)
#     if payment_date > coupon.iborIndex.yts.referenceDate
#       pricer.discount = discount(yts, payment_date)
#     else
#       pricer.discount = 1.0
#     end
#     pricer.spreadLegValue = coupon.spread * accrual_period(coupon) * pricer.discount
#
#     pricer.initialized = true
#   end
#
#   return pricer
# end
#
# function adjusted_fixing(pricer::BlackIborCouponPricer, coupon::IborCoupon, yts::YieldTermStructure, cap_vol::OptionletVolatilityStructure,
#                         fixing::Float64 = index_fixing(coupon, yts))
#   # stuff
#   if !coupon.isInArrears
#     return fixing
#   end
#
#   d1 = coupon.fixingDate
#   ref_date = cap_vol.referenceDate
#
#   if d1 <= ref_date
#     return fixing
#   end
#
#   d2 = value_date(coupon.iborIndex, d1)
#   d3 = maturity_date(coupon.iborIndex, d2)
#   tau = year_fraction(coupon.iborIndex.dc, d2, d3)
#   variance = black_variance(cap_vol, d1, fixing)
#   adjustment = fixing * fixing * variance * tau / (1.0 + fixing * tau)
#   return fixing + adjustment
# end
#
# function swaplet_price(pricer::BlackIborCouponPricer, coupon::IborCoupon, yts::YieldTermStructure, cap_vol::OptionletVolatilityStructure)
#   swapl_price = adjusted_fixing(pricer, coupon, yts, cap_vol) * accrual_period(coupon) * pricer.discount
#   return coupon.gearing * swapl_price + pricer.spreadLegValue
# end
#
# swaplet_rate(pricer::BlackIborCouponPricer, coupon::IborCoupon, yts::YieldTermStructure, cap_vol::OptionletVolatilityStructure) =
#   swaplet_price(pricer, coupon, yts) / (accrual_period(coupon) * pricer.discount)

## Floating Rate Coupon Methods ##
# function calc_rate(coupon::IborCoupon, yts::YieldTermStructure, cap_vol::OptionletVolatilityStructure)
#   # first initialize pricer
#   initialize!(coupon.pricer, coupon, yts)
#
#   # then get swaplet rate
#   return swaplet_rate(coupon.pricer, coupon, yts, cap_vol)
# end

## NPV METHODS ##
function npv(leg::Leg, yts::YieldTermStructure, settlement_date::Date, npv_date::Date)
  # stuff
  totalNPV = 0.0
  if length(leg.coupons) == 0
    return totalNPV
  end

  for i in leg
    if has_occurred(i, settlement_date)
      continue
    end
    @inbounds totalNPV += amount(i) * discount(yts, date(i))
  end

  # redemption - not needed with new iterator
  # totalNPV += amount(leg.redemption) * discount(yts, leg.redemption.date)
  return totalNPV / discount(yts, npv_date)
end

function npv(leg::Leg, y::InterestRate, include_settlement_cf::Bool, settlement_date::Date, npv_date::Date)
  if length(leg.coupons) == 0
    return 0.0
  end

  totalNPV = 0.0
  discount_ = 1.0
  last_date = npv_date

  for cp in leg
    if has_occurred(cp, settlement_date)
      continue
    end
    coupon_date = date(cp)
    amount_ = amount(cp)

    if isa(cp, Coupon)
      ref_start_date = ref_period_start(cp)
      ref_end_date = ref_period_end(cp)
    else
      if last_date == npv_date
        ref_start_date = coupon_date - Dates.Year(1)
      else
        ref_start_date = last_date
      end
      ref_end_date = coupon_date
    end

    b = discount_factor(y, last_date, coupon_date, ref_start_date, ref_end_date)

    discount_ *= b
    last_date = coupon_date

    totalNPV += amount_ * discount_
  end

  # redemption - not needed with iterator
  #redempt_date = leg.redemption.date
  # amount = amount(leg.redemption)

  # b = discount_factor(leg.redemption, last_date, redempt_date, last_date, redempt_date)
  # discount *= b
  # totalNPV += amount * discount

  # now we return total npv
  return totalNPV
end

function npv(leg::ZeroCouponLeg, yts::YieldTermStructure, settlement_date::Date, npv_date::Date)
  if amount(leg.redemption) == 0.0
    return 0.0
  end

  totalNPV = amount(leg.redemption) * discount(yts, date(leg.redemption))
  return totalNPV / discount(yts, npv_date)
end

function npv(leg::ZeroCouponLeg, y::InterestRate, include_settlement_cf::Bool, settlement_date::Date, npv_date::Date)
  if amount(leg.redemption) == 0.0
    return 0.0
  end

  redempt_date = date(leg.redemption)

  ref_start_date = redempt_date - Dates.Year(1)
  ref_end_date = redempt_date

  discount_ = discount_factor(y, npv_date, redempt_date, ref_start_date, ref_end_date)

  totalNPV = amount(leg.redemption) * discount_

  return totalNPV
end

function npvbps(leg::Leg, yts::YieldTermStructure, settlement_date::Date, npv_date::Date, includeSettlementDateFlows::Bool = true)
  npv = 0.0
  bps = 0.0

  for cp in leg
    if has_occurred(cp, settlement_date, includeSettlementDateFlows)
      continue
    end

    df = discount(yts, date(cp))
    npv += amount(cp) * df
    if isa(cp, Coupon)
      bps += cp.nominal * accrual_period(cp) * df
    end
    # if isa(cp, FixedRateCoupon)
    #   println("###########################")
    #   println("df ", df)
    #   println("bps ", bps)
    #   println("Date ", date(cp))
    #   println("nom ", cp.nominal)
    #   println("acru p ", accrual_period(cp))
    #   println("DayCount ", cp.dc)
    # end
  end
  d = discount(yts, npv_date)
  npv /= d
  bps = BASIS_POINT * bps / d

  return npv, bps
end

## Duration Calculations ##
modified_duration_calc(::SimpleCompounding, c::Float64, B::Float64, t::Float64, ::Float64, ::Frequency) = c * B * B * t
modified_duration_calc(::CompoundedCompounding, c::Float64, B::Float64, t::Float64, r::Float64, N::Frequency) = c * t * B / (1 + r / QuantLib.Time.value(N))
modified_duration_calc(::ContinuousCompounding, c::Float64, B::Float64, t::Float64, ::Float64, ::Frequency) = c * B * t
modified_duration_calc(::SimpleThenCompounded, c::Float64, B::Float64, t::Float64, r::Float64, N::Frequency) =
  t <= 1.0 / N ? modified_duration_calc(Simple(), c, B, t, r, N) : modified_duration_calc(CompoundedCompounding(), c, B, t, r, N)

function duration(::ModifiedDuration, leg::Leg, y::InterestRate, dc::DayCount, include_settlement_cf::Bool, settlement_date::Date, npv_date::Date = Date())
  if length(leg.coupons) == 0 # TODO make this applicable to redemption too
    return 0.0
  end

  if npv_date == Date()
    npv_date = settlement_date
  end

  P = 0.0
  t = 0.0
  dPdy = 0.0
  r = y.rate
  N = y.freq
  last_date = npv_date

  for cp in leg
    if has_occurred(cp, settlement_date)
      continue
    end

    c = amount(cp)
    coupon_date = date(cp)
    if isa(cp, FixedRateCoupon)
      ref_start_date = ref_period_start(cp)
      ref_end_date = ref_period_end(cp)
    else
      if last_date == npv_date
        ref_start_date = coupon_date - Dates.Year(1)
      else
        ref_start_date = last_date
      end
      ref_end_date = coupon_date
    end

    t += year_fraction(dc, last_date, coupon_date, ref_start_date, ref_end_date)

    B = discount_factor(y, t)
    P += c * B

    dPdy -= modified_duration_calc(y.comp, c, B, t, r, N)

    last_date = coupon_date
  end

  if P == 0.0
    return 0.0 # no cashflows
  end

  return -dPdy / P # reverse derivative sign
end

function duration(::ModifiedDuration, leg::ZeroCouponLeg, y::InterestRate, dc::DayCount, include_settlement_cf::Bool, settlement_date::Date, npv_date::Date = Date())
  if amount(leg.redemption) == 0.0 || has_occurred(leg.redemption, settlement_date)
    return 0.0
  end

  redempt_date = date(leg.redemption)

  if npv_date == Date()
    npv_date = settlement_date
  end

  r = y.rate
  N = y.freq
  c = amount(leg.redemption)

  ref_start_date = redempt_date - Dates.Year(1)
  ref_end_date = redempt_date

  t = year_fraction(dc, npv_date, redempt_date, ref_start_date, ref_end_date)

  B = discount_factor(y, t)
  P = c * B

  dPdy = modified_duration_calc(y.comp, c, B, t, r, N)

  if P == 0.0
    return 0.0
  end

  return -dPdy / P # reverse derivative sign
end

function aggregate_rate(cf::Leg, _first::Int, _last::Int)
  if cf.coupons[_first] == cf.coupons[_last]
    return 0.0
  end

  _next_cf = cf.coupons[_first]
  paymentDate = date(_next_cf)
  result = 0.0
  i = _first
  while i < length(cf.coupons) && date(_next_cf) == paymentDate
    if isa(_next_cf, Coupon)
      result += calc_rate(_next_cf)
    end

    i += 1
    _next_cf = cf.coupons[i]
  end

  return result
end

# functions for sorting and finding
sort_cashflow(cf::Coupon) = accrual_end_date(cf) #cf.couponMixin.accrualEndDate
sort_cashflow(simp::SimpleCashFlow) = simp.date

prev_cf(cf::Coupon, d::Date) = d > cf.paymentDate
prev_cf(simp::SimpleCashFlow, d::Date) = d > simp.date

next_cf(cf::Coupon, d::Date) = d < cf.paymentDate
next_cf(simp::SimpleCashFlow, d::Date) = d < simp.date

function previous_cashflow_date(cf::Leg, settlement_date::Date)
  # right now we can assume cashflows are sorted by date because of schedule
  prev_cashflow_idx = findprev(prev_cf, cf.coupons, length(cf.coupons), settlement_date)

  return prev_cashflow_idx == 0 ? 0 : cf.coupons[prev_cashflow_idx].paymentDate
end

function accrual_days(cf::CashFlows, dc::DayCount, settlement_date::Date)
  last_payment = previous_cashflow_date(cf, settlement_date)

  if last_payment == 0
    return 0.0
  else
    return day_count(dc, last_payment, settlement_date)
  end
end

function next_cashflow(cf::Leg, settlement_date::Date)
  if settlement_date > date(cf.coupons[end])
    return length(cf.coupons)
  end

  return findnext(next_cf, cf.coupons, 1, settlement_date)
end

function next_coupon_rate(cf::Leg, settlement_date::Date)
  next_cf_idx = next_cashflow(cf, settlement_date)
  return aggregate_rate(cf, next_cf_idx, length(cf.coupons))
end

function accrued_amount(cf::Leg, settlement_date::Date, include_settlement_cf::Bool = false)
  next_cf_idx = next_cashflow(cf, settlement_date)
  if cf.coupons[next_cf_idx] == length(cf.coupons)
    return 0.0
  end

  _next_cf = cf.coupons[next_cf_idx]
  paymentDate = date(_next_cf)
  result = 0.0
  i = next_cf_idx
  while i < length(cf.coupons) && date(_next_cf) == paymentDate
    result += accrued_amount(_next_cf, settlement_date)
    i += 1
    _next_cf = cf.coupons[i]
  end

  return result
end

accrued_amount(cf::ZeroCouponLeg, ::Date, ::Bool= false) = 0.0
accrued_amount(simp::SimpleCashFlow, ::Date, ::Bool = false) = 0.0

# accrual_period!{C <: Coupon}(coup::C, p::Float64) = coup.accrualPeriod = p
# function accrual_period{C <: Coupon}(coup::C)
#   if accrual_period(coup) == -1.0
#     p = year_fraction(get_dc(coup), accrual_start_date(coup), accrual_end_date(coup), ref_period_start(coup), ref_period_end(coup))
#     accrual_period!(coup, p)
#   end
#
#   return p
# end

function has_occurred(cf::CashFlow, ref_date::Date, include_settlement_cf::Bool = true)
  # will need to expand this
  if ref_date < date(cf) || (ref_date == date(cf) && include_settlement_cf)
    return false
  else
    return true
  end
end

function maturity_date(leg::Leg)
  # sort cashflows
  sorted_cf = sort(leg.coupons, by=sort_cashflow)
  return date_accrual_end(sorted_cf[end])
end

# Yield Calculations ##
function yield(leg::Leg, npv::Float64, dc::DayCount, compounding::CompoundingType, freq::Frequency, include_settlement_cf::Bool, settlement_date::Date,
              npv_date::Date, accuracy::Float64, max_iter::Int, guess::Float64)
  solver = NewtonSolver(max_iter)
  obj_fun = IRRFinder(leg, npv, dc, compounding, freq, include_settlement_cf, settlement_date, npv_date)

  return solve(solver, QuantLib.operator(obj_fun), accuracy, guess, guess / 10.0)
end

## ITERATORS ##
# this is to iterate through cash flows and redemption
# Base.start(f::FixedRateLeg) = 1
# function Base.next(f::FixedRateLeg, state)
#   if state > length(f.coupons)
#     f.redemption, state + 1
#   else
#     f.coupons[state], state + 1
#   end
# end
#
# Base.done(f::FixedRateLeg, state) = length(f.coupons) + 1 < state

Base.start(f::Leg) = 1
Base.next(f::Leg, state) = f.coupons[state], state + 1
Base.done(f::Leg, state) = length(f.coupons) == state - 1

Base.getindex(f::Leg, i::Int) = f.coupons[i]
Base.endof(f::Leg) = endof(f.coupons)
# end


text:
__precompile__()

module SkyCoords
using StaticArrays

export AbstractSkyCoords, ICRSCoords, GalCoords, FK5Coords, separation

include("types.jl")

"""
    SkyCoords.str2rad(::AbstractString, force_ha=false)

Parses strings that specify common astronomical coordinates into radians. A simplified version of the regex allowed for both

* hour angle - `"xx [h] xx ['′m:] xx [\\\"″s]?"`
* degree - `"xx [°d:] xx ['′m:] xx [\\\"″s]?"`

So, you can see that this will parse an hour angle if the leading unit contains an `h` otherwise it will parse a degree if the leading unit contains a `°`, `d` or `:`. Also, you can forego the final unit, since it will be parsed automatically. The values (`xx`) will be parsed as `Float64`.

If `force_ha` is `true` then the input will be parsed as an hour angle using a looser regex that allows the leading unit to be a `:`

* loose hour angle - `"xx [h:] xx ['′m:] xx [\\\"″s]?"`

# Examples
```jldoctest
julia> import SkyCoords: str2rad

julia> str2rad("12h52m64.300s")
3.3731614843033575

julia> str2rad("0°12'5\\\"") # Note that you have to escape the "
0.003514899188044136

julia> str2rad("-  10° 02 ′  10.885 ″") # Whitespace does not matter
-0.17389837681291273

julia> str2rad("12:0:0", true)
3.141592653589793
```
"""
function str2rad(input::AbstractString, force_ha = false)
    input = replace(strip(input), r"\s" => "")
    ha_r = r"([+-]?\d+\.?\d*)h(\d+\.?\d*)['′m:](\d+\.?\d*)[\"″s]?"
    ha_r_loose = r"([+-]?\d+\.?\d*)[h:](\d+\.?\d*)['′m:](\d+\.?\d*)[\"″s]?"
    deg_r = r"([+-]?\d+\.?\d*)[°d:](\d+\.?\d*)['′m:](\d+\.?\d*)[\"″s]?"
    if force_ha
        m = match(ha_r_loose, input)
        m === nothing && error("Could not parse $input as hour angle")
        h, m, s = parse.(Float64, m.captures)
        rad = hr2rad(h)
        rad += hr2rad(m / 60)
        rad += hr2rad(s / 3600)
        return rad
    elseif occursin(ha_r, input)
        m = match(ha_r, input)
        h, m, s = parse.(Float64, m.captures)
        rad = hr2rad(h)
        rad += hr2rad(m / 60)
        rad += hr2rad(s / 3600)
        return rad
    elseif occursin(deg_r, input)
        m = match(deg_r, input)
        d, m, s = parse.(Float64, m.captures)
        rad = deg2rad(d)
        rad += deg2rad(m / 60)
        rad += deg2rad(s / 3600)
        return rad
    else
        error("Could not parse $input to sky coordinates")
    end
end

@inline hr2rad(h::Number) = π * h / 12

# -----------------------------------------------------------------------------
# Helper functions: Create rotation matrix about a given axis (x, y, z)

function xrotmat(angle)
    s, c = sincos(angle)
    SMatrix{3,3}(1, 0, 0, 0, c, -s, 0, s, c)
end

function yrotmat(angle)
    s, c = sincos(angle)
    SMatrix{3,3}(c, 0, s, 0, 1, 0, -s, 0, c)
end

function zrotmat(angle)
    s, c = sincos(angle)
    SMatrix{3,3}(c, -s, 0, s, c, 0, 0, 0, 1)
end

# (lon, lat) -> [x, y, z] unit vector
function coords2cart(lon, lat)
    sinlon, coslon = sincos(lon)
    sinlat, coslat = sincos(lat)
    SVector{3}(coslat * coslon, coslat * sinlon, sinlat)
end

# [x, y, z] unit vector -> (lon, lat)
cart2coords(v) = atan(v[2], v[1]), atan(v[3], sqrt(v[1] * v[1] + v[2] * v[2]))

# -----------------------------------------------------------------------------
# Constant rotation matricies and precession matrix function

# ICRS --> FK5 at J2000 (See USNO Circular 179, section 3.5)
eta0 = deg2rad(-19.9 / 3600000.0)
xi0 = deg2rad(9.1 / 3600000.0)
da0 = deg2rad(-22.9 / 3600000.0)
const ICRS_TO_FK5J2000 = xrotmat(-eta0) * yrotmat(xi0) * zrotmat(da0)
const FK5J2000_TO_ICRS = ICRS_TO_FK5J2000'

# FK5J2000 --> Gal
# Note that galactic pole and zeropoint of l are somewhat arbitrary
# and not officially defined (as far as I know). The values below are
# from astropy.coordinates, which includes the following comment:
# | "This gives better consistency with other codes than using the values
# |  from Reid & Brunthaler 2004 and the best self-consistency between FK5
# |  -> Galactic and FK5 -> FK4 -> Galactic. The lon0 angle was found by
# |  optimizing the self-consistency."
ngp_fk5j2000_ra = deg2rad(192.8594812065348)
ngp_fk5j2000_dec = deg2rad(27.12825118085622)
lon0_fk5j2000 = deg2rad(122.9319185680026)
const FK5J2000_TO_GAL = (zrotmat(pi - lon0_fk5j2000) *
                         yrotmat(pi / 2.0 - ngp_fk5j2000_dec) * zrotmat(ngp_fk5j2000_ra))
const GAL_TO_FK5J2000 = FK5J2000_TO_GAL'

# Gal --> ICRS: simply chain through FK5J2000
const GAL_TO_ICRS = FK5J2000_TO_ICRS * GAL_TO_FK5J2000
const ICRS_TO_GAL = GAL_TO_ICRS'

# FK5J2000 --> FK5{epoch}
# Computes the precession matrix from J2000 to the given Julian equinox.
# Expression from from Capitaine et al. 2003 as expressed in the USNO
# Circular 179.  This should match the IAU 2006 standard from SOFA.
const pzeta = [2.650545, 2306.083227, 0.2988499, 0.01801828, -0.000005971, -0.0000003173]
const pz = [-2.650545, 2306.077181, 1.0927348, 0.01826837, -0.000028596, -0.0000002904]
const ptheta = [0.0, 2004.191903, -0.4294934, -0.04182264, -0.000007089, -0.0000001274]
function precess_from_j2000(equinox)
    t = (equinox - 2000.0) / 100.0
    tn = 1.0
    zeta = pzeta[1]
    z = pz[1]
    theta = ptheta[1]
    for i = 2:6
        tn *= t
        zeta += pzeta[i] * tn
        z += pz[i] * tn
        theta += ptheta[i] * tn
    end
    zeta = deg2rad(zeta / 3600.0)
    z = deg2rad(z / 3600.0)
    theta = deg2rad(theta / 3600.0)
    zrotmat(-z) * yrotmat(theta) * zrotmat(-zeta)
end

# -----------------------------------------------------------------------------
# Type-dependent methods

lon(c::GalCoords) = c.l
lat(c::GalCoords) = c.b
lon(c::AbstractSkyCoords) = c.ra
lat(c::AbstractSkyCoords) = c.dec

# Abstract away specific field names (ra, dec vs l, b)
coords2cart(c::AbstractSkyCoords) = coords2cart(lon(c), lat(c))

# Rotation matrix between coordinate systems: `rotmat(to, from)`
# Note that all of these return SMatrix{3,3}{Float64}, regardless of
# element type of input coordinates.
rotmat(::Type{T1}, ::Type{T2}) where {T1<:GalCoords,T2<:ICRSCoords} = ICRS_TO_GAL
rotmat(::Type{T1}, ::Type{T2}) where {T1<:ICRSCoords,T2<:GalCoords} = GAL_TO_ICRS

# Define both these so that `convert(FK5Coords{e}, ...)` and
# `convert(FK5Coords{e,T}, ...)` both work. Similar with other
# FK5Coords rotmat methods below.
@generated rotmat(::Type{FK5Coords{e1}}, ::Type{T2}) where {e1,T2<:ICRSCoords} =
    precess_from_j2000(e1) * ICRS_TO_FK5J2000
@generated rotmat(::Type{FK5Coords{e1,T1}}, ::Type{T2}) where {e1,T1,T2<:ICRSCoords} =
    precess_from_j2000(e1) * ICRS_TO_FK5J2000

@generated rotmat(::Type{FK5Coords{e1}}, ::Type{T2}) where {e1,T2<:GalCoords} =
    precess_from_j2000(e1) * GAL_TO_FK5J2000
@generated rotmat(::Type{FK5Coords{e1,T1}}, ::Type{T2}) where {e1,T1,T2<:GalCoords} =
    precess_from_j2000(e1) * GAL_TO_FK5J2000

@generated rotmat(::Type{T1}, ::Type{FK5Coords{e2,T2}}) where {T1<:ICRSCoords,e2,T2} =
    FK5J2000_TO_ICRS * precess_from_j2000(e2)'

@generated rotmat(::Type{T1}, ::Type{FK5Coords{e2,T2}}) where {T1<:GalCoords,e2,T2} =
    FK5J2000_TO_GAL * precess_from_j2000(e2)'

@generated rotmat(::Type{FK5Coords{e1}}, ::Type{FK5Coords{e2,T2}}) where {e1,e2,T2} =
    precess_from_j2000(e1) * precess_from_j2000(e2)'
@generated rotmat(::Type{FK5Coords{e1,T1}}, ::Type{FK5Coords{e2,T2}}) where {e1,T1,e2,T2} =
    precess_from_j2000(e1) * precess_from_j2000(e2)'

# ------------------------------------------------------------------------------
# Distance between coordinates

function _separation(λ_1, ϕ_1, λ_2, ϕ_2)
    Δλ = λ_2 - λ_1
    sin_Δλ, cos_Δλ = sincos(Δλ)
    sin_ϕ1, cos_ϕ1 = sincos(ϕ_1)
    sin_ϕ2, cos_ϕ2 = sincos(ϕ_2)
    return atan(
        hypot(cos_ϕ2 * sin_Δλ, cos_ϕ1 * sin_ϕ2 - sin_ϕ1 * cos_ϕ2 * cos_Δλ),
        sin_ϕ1 * sin_ϕ2 + cos_ϕ1 * cos_ϕ2 * cos_Δλ,
    )
end

"""
    separation(c1::AbstractSkyCoords, c2::AbstractSkyCoords) -> distance

Return angular separation between two sky coordinates, in radians.

The angular separation is calculated using the [Vincenty formula](http://en.wikipedia.org/wiki/Great-circle_distance), which is slightly more
complex and computationally expensive than some alternatives, but is stable at
at all distances, including the poles and antipodes.
"""
separation(c1::T, c2::T) where {T<:AbstractSkyCoords} =
    _separation(lon(c1), lat(c1), lon(c2), lat(c2))

separation(c1::T1, c2::T2) where {T1<:AbstractSkyCoords,T2<:AbstractSkyCoords} =
    separation(c1, convert(T1, c2))

end # module


text:
include(joinpath(@__DIR__, "staticrbd.jl"))
using .StaticRBD
using RigidBodyDynamics
using LinearAlgebra

BLAS.set_num_threads(1)
GC.enable(true)

# repeated info from BenchmarkInputDataGen module, but don't really want to depend on that here.

struct URDFConfig
    basename::String
    floating::Bool
end

const URDF_DIR = joinpath(@__DIR__, "..", "..", "..", "description", "urdf")
const URDF_CONFIGS = [URDFConfig("iiwa.urdf", false), URDFConfig("atlas.urdf", true), URDFConfig("hyq.urdf", true)]
const CSV_DIR = joinpath(@__DIR__, "..", "..", "csv", "rbdjl")

for config in URDF_CONFIGS
    urdf = joinpath(URDF_DIR, config.basename)
    robotname, _ = splitext(config.basename)
    csv = joinpath(CSV_DIR, "rbdjl", robotname, ".csv")
    floating = config.floating

    c_floating = Cint(floating)
    for c_scalar_type in [Cint(1), Cint(2)]
        mechanism = GC.@preserve urdf csv begin
            c_urdf = Base.unsafe_convert(Cstring, urdf)
            c_csv = Base.unsafe_convert(Cstring, csv)
            StaticRBD.create_mechanism(c_urdf, c_floating, c_scalar_type)
        end
        state = StaticRBD.create_state(mechanism)
        result = StaticRBD.create_dynamics_result(mechanism)

        jointwrenches = result.jointwrenches
        accelerations = result.accelerations
        M = result.massmatrix
        vd = result.v̇

        q = configuration(state)
        v = velocity(state)
        vd_desired = similar(v)
        tau = similar(v)

        StaticRBD.inverse_dynamics(tau, jointwrenches, accelerations, state, vd_desired)
        StaticRBD.mass_matrix(M, state)
        StaticRBD.dynamics(result, state, tau)
    end
end

max_stars_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_issues_repo_name: CobbledSteel/rbd-benchmarks-riscv

max_forks_repo_name: CobbledSteel/rbd-benchmarks-riscv


text:
@testset "dict_utils" begin
    setindex_no_overwrite! = DiagnosisClassification.setindex_no_overwrite!
    setindex_same_value!   = DiagnosisClassification.setindex_same_value!

    @testset "setindex_no_overwrite!" begin
        dict = Dict{String, String}()
        @test setindex_no_overwrite!(dict, "hello", "world") isa Nothing
        @test_throws ErrorException setindex_no_overwrite!(dict, "hello", "world2")
    end

    @testset "setindex_same_value!" begin
        dict = Dict{String, String}()
        @test setindex_same_value!(dict, "hello", "world") isa Nothing
        @test setindex_same_value!(dict, "hello", "world") isa Nothing
        @test setindex_same_value!(dict, "hello", "world") isa Nothing
        @test_throws ErrorException setindex_same_value!(dict, "hello", "world2")
    end
end


text:
using Test, LinearMaps
import LinearMaps: FiveArg, ThreeArg

const matrixstyle = VERSION ≥ v"1.3.0-alpha.115" ? FiveArg() : ThreeArg()

const testallocs = VERSION ≥ v"1.4-"

include("linearmaps.jl")

include("transpose.jl")

include("functionmap.jl")

include("scaledmap.jl")

include("composition.jl")

include("linearcombination.jl")

include("wrappedmap.jl")

include("uniformscalingmap.jl")

include("numbertypes.jl")

include("blockmap.jl")

include("kronecker.jl")

include("conversion.jl")

include("left.jl")


text:
module RungeKutta

    using DelimitedFiles
    using Markdown
    using PrettyTables
    using Reexport

    import GenericLinearAlgebra
    import LinearAlgebra: istril
    import Polynomials
    import Polynomials: Polynomial

    include("utils.jl")

    include("tableau.jl")
    include("tableau_partitioned.jl")

    export Tableau, PartitionedTableau

    export isexplicit,
           isimplicit,
           isdiagnonallyimplicit,
           isfullyimplicit

    include("order_conditions.jl")

    export check_order_conditions_b,
           check_order_conditions_c,
           check_order_conditions_d,
           satisfies_simplifying_assumption_b,
           satisfies_simplifying_assumption_c,
           satisfies_simplifying_assumption_d,
           solve_simplifying_assumption_b,
           solve_simplifying_assumption_c,
           solve_simplifying_assumption_d

    include("symmetry.jl")

    export check_symmetry,
           issymmetric

    include("symplecticity.jl")

    export check_symplecticity,
           compute_symplecticity_error,
           issymplectic,
           get_symplectic_conjugate_coefficients,
           get_symplectic_conjugate_coefficients!,
           symplecticize

    include("Tableaus.jl")
    @reexport using .Tableaus

    include("PartitionedTableaus.jl")
    @reexport using .PartitionedTableaus

end


text:
#=
grid_linear_to_cart:
- Julia version: 1.1.0
- Author: Jirik
- Date: 2019-03-17
=#

grid_size= [2, 3, 4]

cube_carthesian_position = [2,2,2]
i, j, k = cube_carthesian_position
# i, j, k = [2,2,2]

function grid_voxel_cart_to_id(grid_size, carthesian)
    sz1,sz2,sz3 = grid_size
    trf = 0
    ind = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k
    return ind
end

function grid_voxel_id_to_cart(grid_size, ind)
#     ind = (sz2 * sz3) * (i - 1)  + (j - 1) * sz3 + k
    sz1,sz2,sz3 = grid_size
    layer = sz2*sz3
    ir = div(ind, layer) + 1
    rest1 = mod(ind, layer)
#     println("rest1 ", rest1)
    row = sz3
    jr = div(rest1, row) + 1
    rest2 = mod(rest1, row)
    kr = rest2
    return [ir, jr, kr]
end

function grid_x_face_to_carthesian(grid_size, ind)
    return grid_voxel_id_to_cart(grid_size, ind)
end

function grid_y_face_to_carthesian(grid_size, ind)
#     f20 = nax1 +
#         (sz2 + 1) * sz3 * (i - 1 + trf)  + (j - 1 + trf) * sz3 + k
    sz1,sz2,sz3 = grid_size
    nax1 = (1 + sz1) * sz2 * sz3
    ind = ind - nax1

    layer = (sz2 + 1) * sz3
    ir = div(ind, layer) + 1
    rest1 = mod(ind, layer)

    jr = div(rest1, sz3) + 1
    rest2 = mod(rest1, sz3)
    kr = rest2
    return [ir, jr, kr]
end

function grid_z_face_to_carthesian(grid_size, ind)
#     nax3_layer = (sz3 + 1) * sz2 * (i - 1)
#     nax3_row = (j - 1) * (sz3 + 1)
#     f30 = nax1 +  nax2 +
#         nax3_layer + nax3_row  + k + trf
    sz1,sz2,sz3 = grid_size
    nax1 = (1 + sz1) * sz2 * sz3
    nax2 = sz1 * (1 + sz2) * sz3
    ind = ind - nax1 - nax2

    layer = (sz3 + 1) * sz2
    ir = div(ind, layer) + 1
    rest1 = mod(ind, layer)

    row = (sz3 + 1)
    jr = div(rest1, row) + 1
    rest2 = mod(rest1, row)
    kr = rest2
    return [ir, jr, kr]
end

# grid_size = sz
# sz1,sz2,sz3 = grid_size
# trf = 0
# ind = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k

# nax1 = (1 + sz1) * sz2 * sz3
ind = grid_voxel_cart_to_id(sz, [2,2,2])
println("linear index: ", ind)

ir, jr, kr = grid_voxel_id_to_cart(sz, ind)

println("carthesian index: ", ir, " ", jr, " ", kr)

#####################################

function grid_face_id_to_cartesian(grid_size, fid)
    sz1,sz2,sz3 = grid_size
    nax1 = (1 + sz1) * sz2 * sz3
    nax2 = sz1 * (1 + sz2) * sz3

    if fid <= nax1
        # x-face it is the same index as voxel index
        axis = 1
        voxel_cart = grid_x_face_to_carthesian(grid_size, fid)

    elseif fid <= (nax1 + nax2)
        # y-face
        axis = 2
        voxel_cart = grid_y_face_to_carthesian(grid_size, fid)
    else
        # z-face
        axis = 3
        voxel_cart = grid_z_face_to_carthesian(grid_size, fid)
    end
    return voxel_cart, axis
end

# this is for grid_size = [2,3,4] description of voxel_cart = [2,2,2]
grid_size= [2, 3, 4]
fid = 18
fid = 90
fid = 58

voxel_cart, axis = grid_face_id_to_cartesian(grid_size, fid)
println("voxel cart: ", voxel_cart, ", axis: ", axis)

# i, j, k = cube_carthesian_position
#
# f10 = (sz2 * sz3) * (i - 1 + trf)  + (j - 1) * sz3 + k
# @debug ("number of 1st axis faces: ", nax1, ", ")
# f20 = nax1 +
#     (sz2 + 1) * sz3 * (i - 1)  + (j - 1 + trf) * sz3 + k
#
# @debug ("2st axis faces: ", nax2, ", ")
# nax3_layer = (sz3 + 1) * sz2 * (i - 1)
# nax3_row = (j - 1) * (sz3 + 1)
# @debug ("3st axis faces in one layer and in one row: ",
#     nax3_layer, " ", nax3_row,  "\n")
# f30 = nax1 +  nax2 +
#     nax3_layer + nax3_row  + k + trf


text:
function compute()
	"""
	seen=Set{BigInt}()
	for a=2:100
		for b=2:100
			d=BigInt(a);
			push!(seen,d^b)
		end
	end
	return length(seen)
	"""

	seen=Set(a^b for a::BigInt in 2:100 for b::BigInt in 2:100)

	return length(seen)

end

println(compute())


text:
using Trixi, OrdinaryDiffEq

surface_flux = FluxLaxFriedrichs()
volume_flux = flux_ranocha
dg = DGMulti(polydeg = 3, element_type = Line(), approximation_type = Polynomial(),
             surface_integral = SurfaceIntegralWeakForm(surface_flux),
             volume_integral = VolumeIntegralFluxDifferencing(volume_flux))

equations = CompressibleEulerEquations1D(1.4)

initial_condition = initial_condition_convergence_test
source_terms = source_terms_convergence_test

mesh = DGMultiMesh(dg, cells_per_dimension=(8,),
                   coordinates_min=(-1.0,), coordinates_max=(1.0,), periodicity=true)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, dg;
                                    source_terms=source_terms)

tspan = (0.0, 1.1)
ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()
alive_callback = AliveCallback(alive_interval=10)
analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, uEltype=real(dg))
callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        alive_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt = 0.5 * estimate_dt(mesh, dg), save_everystep=false, callback=callbacks);

summary_callback() # print the timer summary


text:
using FastGroupBy, Compat
import Base.Forward
import SortingAlgorithms: RADIX_MASK, RADIX_SIZE, uint_mapping

"""
sumby that uses much less memory but is slower

It saves memory by not creating a separate copy of by and val vectors but instead
follows a chain of assignment for the swap operations.

"""
function sumby_lessmem_chain{T, S<:Number}(by::AbstractVector{T},  val::AbstractVector{S})
  o = Forward
  lo = 1
  hi = length(by)

  if lo >= hi;  return by;  end

  if !isbits(T)
      error("Radix sort only sorts bits types (got $T)")
  end

  # Init
  iters = ceil(Integer, sizeof(T)*8/RADIX_SIZE)
  bin = zeros(UInt32, 2^RADIX_SIZE, iters)
  if lo > 1;  bin[1,:] = lo-1;  end

  # Histogram for each element, radix
  for i = lo:hi
    v = uint_mapping(o, by[i])
    for j = 1:iters
        idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
        @inbounds bin[idx,j] += 1
    end
  end

  # Sort!
  len = hi-lo+1
  for j = 1:iters
      # Unroll first data iteration, check for degenerate case
      v = uint_mapping(o, by[hi])
      idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

      # are all values the same at this radix?
      if bin[idx,j] == len;  continue;  end

      cbin = cumsum(bin[:,j])

      # chain assignment to save memory
      touched = BitArray(hi) # to keep track of which cell has been touched
      i = hi
      while i > 0
          begin_index = i # this is where the chain starts
          tmp = (by[i], val[i])
          v = uint_mapping(o, tmp[1])
          idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1

          while true
              ci = cbin[idx]
              new_tmp = (by[ci], val[ci])
              by[ci], val[ci] = tmp[1], tmp[2]
              tmp = new_tmp
              touched[ci] = true # keep track of which number is done
              cbin[idx] -= 1

              if ci == begin_index # the chain has ended
                  break;
              end
              v = uint_mapping(o, tmp[1])
              idx = @compat(Int((v >> (j-1)*RADIX_SIZE) & RADIX_MASK)) + 1
          end

          while i > 0 && touched[i]
              i -= 1
          end
      end
  end

  sumby_contiguous(by, val)
end


text:
using Test
import MultivariateBases
using DynamicPolynomials

function quadratic_test(
    optimizer, config::MOIT.Config,
    cone::SumOfSquares.PolyJuMP.PolynomialSet, basis, bivariate::Bool)
    atol = config.atol
    rtol = config.rtol

    model = _model(optimizer)

    @variable(model, α)

    @polyvar x
    if bivariate
        @polyvar y
        poly = x^2 + α*x*y + y^2
        cert_monos = [x, y]
        monos = [x^2, x*y, y^2]
    else
        poly = x^2 + α*x + 1
        cert_monos = [x, 1]
        monos = [x^2, x, 1]
    end
    cref = @constraint(model, poly in cone, basis = basis)

    @objective(model, Max, α)
    optimize!(model)

    if basis == ChebyshevBasis
        err = ErrorException("`certificate_monomials` is not supported with `$(basis{typeof(x + 1.0)})`, use `certificate_basis` instead.")
        @test_throws err certificate_monomials(cref)
    else
        @test certificate_monomials(cref) == cert_monos
    end

    @test termination_status(model) == MOI.OPTIMAL
    @test objective_value(model) ≈ 2.0 atol=atol rtol=rtol

    @test primal_status(model) == MOI.FEASIBLE_POINT
    @test value(α) ≈ 2.0 atol=atol rtol=rtol

    test_constraint_primal(cref, value(poly))

    p = gram_matrix(cref)
    @test getmat(p) ≈ ones(2, 2) atol=atol rtol=rtol
    if basis == ChebyshevBasis
        @test p.basis.polynomials == cert_monos
    else
        @test p.basis.monomials == cert_monos
    end

    a = moment_value.(moments(dual(cref)))
    @test a[2] ≈ -1.0 atol=atol rtol=rtol
    @test a[1] + a[3] ≈ 2.0 atol=atol rtol=rtol

    @test dual_status(model) == MOI.FEASIBLE_POINT
    for μ in [dual(cref), moments(cref)]
        @test μ isa AbstractMeasure{Float64}
        @test length(moments(μ)) == 3
        @test a ≈ moment_value.(moments(μ)) atol=atol rtol=rtol
        @test monomial.(moments(μ)) == monos
    end

    ν = moment_matrix(cref)
    @test getmat(ν) ≈ [a[1] a[2]
                       a[2] a[3]] atol=atol rtol=rtol
    if basis == ChebyshevBasis
        @test p.basis.polynomials == cert_monos
    else
        @test ν.basis.monomials == cert_monos
    end

    S = SumOfSquares.SOSPolynomialSet{
        SumOfSquares.FullSpace, Monomial{true}, MonomialVector{true}, SumOfSquares.Certificate.Newton{typeof(cone), basis, Tuple{}}
    }
    @test list_of_constraint_types(model) == [(Vector{AffExpr}, S)]
    test_delete_bridge(
        model, cref, 1,
        ((MOI.VectorOfVariables, MOI.Nonnegatives, 0),
         (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 0),
         (MOI.VectorAffineFunction{Float64},
          SumOfSquares.PolyJuMP.ZeroPolynomialSet{
              SumOfSquares.FullSpace, basis,
              Monomial{true}, MonomialVector{true}},
          0)))
end
sos_univariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SOSCone(), MonomialBasis, false)
sd_tests["sos_univariate_quadratic"] = sos_univariate_quadratic_test
sos_bivariate_quadratic_test(optimizer, config)    = quadratic_test(optimizer, config, SOSCone(), MonomialBasis, true)
sd_tests["sos_bivariate_quadratic"] = sos_bivariate_quadratic_test
sos_scaled_univariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SOSCone(), ScaledMonomialBasis, false)
sd_tests["sos_scaled_univariate_quadratic"] = sos_scaled_univariate_quadratic_test
sos_scaled_bivariate_quadratic_test(optimizer, config)    = quadratic_test(optimizer, config, SOSCone(), ScaledMonomialBasis, true)
sd_tests["sos_scaled_bivariate_quadratic"] = sos_scaled_bivariate_quadratic_test
sos_cheby_univariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SOSCone(), ChebyshevBasis, false)
sd_tests["sos_cheby_univariate_quadratic"] = sos_cheby_univariate_quadratic_test
sos_cheby_bivariate_quadratic_test(optimizer, config)    = quadratic_test(optimizer, config, SOSCone(), ChebyshevBasis, true)
sd_tests["sos_scaled_bivariate_quadratic"] = sos_scaled_bivariate_quadratic_test

sdsos_univariate_quadratic_test(optimizer, config)         = quadratic_test(optimizer, config, SDSOSCone(), MonomialBasis, false)
soc_tests["sdsos_univariate_quadratic"] = sdsos_univariate_quadratic_test
sdsos_bivariate_quadratic_test(optimizer, config)          = quadratic_test(optimizer, config, SDSOSCone(), MonomialBasis, true)
soc_tests["sdsos_bivariate_quadratic"] = sdsos_bivariate_quadratic_test
sdsos_scaled_univariate_quadratic_test(optimizer, config)  = quadratic_test(optimizer, config, SDSOSCone(), ScaledMonomialBasis, false)
soc_tests["sdsos_scaled_univariate_quadratic"] = sdsos_scaled_univariate_quadratic_test
sdsos_scaled_bivariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SDSOSCone(), ScaledMonomialBasis, true)
soc_tests["sdsos_scaled_bivariate_quadratic"] = sdsos_scaled_bivariate_quadratic_test
sdsos_cheby_univariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, SDSOSCone(), ChebyshevBasis, false)
soc_tests["sdsos_cheby_univariate_quadratic"] = sdsos_cheby_univariate_quadratic_test
sdsos_cheby_bivariate_quadratic_test(optimizer, config)    = quadratic_test(optimizer, config, SDSOSCone(), ChebyshevBasis, true)
soc_tests["sdsos_scaled_bivariate_quadratic"] = sdsos_scaled_bivariate_quadratic_test

dsos_univariate_quadratic_test(optimizer, config)        = quadratic_test(optimizer, config, DSOSCone(), MonomialBasis, false)
linear_tests["dsos_univariate_quadratic"] = dsos_univariate_quadratic_test
dsos_bivariate_quadratic_test(optimizer, config)         = quadratic_test(optimizer, config, DSOSCone(), MonomialBasis, true)
linear_tests["dsos_bivariate_quadratic"] = dsos_bivariate_quadratic_test
dsos_scaled_univariate_quadratic_test(optimizer, config) = quadratic_test(optimizer, config, DSOSCone(), ScaledMonomialBasis, false)
linear_tests["dsos_scaled_univariate_quadratic"] = dsos_scaled_univariate_quadratic_test
dsos_scaled_bivariate_quadratic_test(optimizer, config)  = quadratic_test(optimizer, config, DSOSCone(), ScaledMonomialBasis, true)
linear_tests["dsos_scaled_bivariate_quadratic"] = dsos_scaled_bivariate_quadratic_test
dsos_cheby_univariate_quadratic_test(optimizer, config)  = quadratic_test(optimizer, config, DSOSCone(), ChebyshevBasis, false)
linear_tests["dsos_cheby_univariate_quadratic"] = dsos_cheby_univariate_quadratic_test
dsos_cheby_bivariate_quadratic_test(optimizer, config)   = quadratic_test(optimizer, config, DSOSCone(), ChebyshevBasis, true)
linear_tests["dsos_cheby_bivariate_quadratic"] = dsos_cheby_bivariate_quadratic_test

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Synchronizers
using Synchronizers: stack, unstack, Record, Ledger, serve, isready,  Synchronizer, sync
using Test

### Testing Int

io = IOBuffer()
stack(io,123)

ionew = IOBuffer(take!(io))
@test unstack(ionew,Int)==123

### Testing Record

record = Record("hello",b"world")

io = IOBuffer()
stack(io,record)

ionew = IOBuffer(take!(io))
@test unstack(ionew,Record)==record

### Testing now the ledger
masterdir = dirname(@__FILE__) * "/master/"

mledger = Ledger(masterdir)
push!(mledger,Record("hellohere",b"world"))

# ### Now let's test synchronization

task = @async serve(2000,mledger)

sleep(1)

slavedir = dirname(@__FILE__) * "/slave/"
sledger = Ledger(slavedir)

s = Synchronizer(2000,sledger)

@test isready(s)==true

sync(s)

@test mledger==sledger

# WIP IO is too fast
#sledger2 = Ledger(sledger)
#@test sledger==sledger2


text:
using GLFW, Test

println(GLFW.GetVersionString())

# GLFWError uses enum for known error code
@test isa(GLFW.GLFWError(GLFW.NOT_INITIALIZED, "").code, GLFW.ErrorCode)
@test GLFW.GLFWError(0x00010001, "").code == GLFW.NOT_INITIALIZED

# GLFWError uses integer for unrecognized error code
@test isa(GLFW.GLFWError(0xDEADBEEF, "").code, Integer)

if !haskey(ENV, "CI")  # AppVeyor and Travis CI don't support OpenGL
	include("windowclose.jl")
end


text:
using LinearAlgebra
using Random
using Statistics: mean
using Test

using Setfield: @lens

using Bifurcations
using Bifurcations.Codim1
using Bifurcations.Codim2
using Bifurcations: BifurcationProblem, special_intervals
using Bifurcations.Codim2: cast_container

macro test_nothrow(ex)
    quote
        @test begin
            $(esc(ex))
            true
        end
    end
end

function nullshow(x)
    show(devnull, x)
end

function nullshow(mime::MIME, x)
    show(devnull, mime, x)
end

function smoke_test_solver_show(solver)
    @test_nothrow nullshow(solver)
    @test_nothrow nullshow(solver.sol)

    for sweep in solver.sol.sweeps
        @test_nothrow nullshow(sweep)
        for point in special_intervals(sweep)
            @test_nothrow nullshow(point)
        end
    end
end

function smoke_test_solver_plot(solver)
    for p in special_intervals(solver)
        @test_nothrow nullshow(plot(p))
    end
    for sweep in solver.sol.sweeps
        @test_nothrow nullshow(plot(sweep))
    end
    @test_nothrow nullshow(plot(solver.sol))
    @test_nothrow nullshow(plot(solver))
end

function _generalize_f(f, A)
    function g(u::U, p, t) where {T, U <: AbstractArray{T}}
        v = cast_container(A, u)
        H = f(v, p, t)
        return cast_container(U, H)
    end

    g(args...) = f(args...)

    return g
end

generalized_f(mod) = _generalize_f(mod.f, typeof(mod.u0))

function uniquify_points(point_list; digits=3)
    uniquified = []
    seen = []
    for point in point_list
        y = round.(point.u, digits=digits)
        for x in seen
            if x == y
                @goto seen
            end
        end
        push!(uniquified, point)
        push!(seen, y)

        @label seen
    end
    return uniquified
end