text:
# Good candidates: tab10 (glasbey_category10_n256), tableau_10, seaborn_colorblind
const COLORS = ColorSchemes.glasbey_hv_n256.colors

const MARKERS = [:circle, :rect, :dtriangle, :utriangle, :cross, :diamond, :ltriangle, :rtriangle, :pentagon, :xcross, :hexagon]

const CYCLE = Cycle([:color, :marker], covary=true)

# = figure_height / figure_width
const HWRATIO = 0.68


text:
# Code derived from metabolicEP (https://github.com/anna-pa-m/Metabolic-EP)

struct EPFields{T<:AbstractFloat}
    av::Vector{T}
    va::Vector{T}
    a::Vector{T}
    d::Vector{T}
    μ::Vector{T}
    s::Vector{T}
    siteflagave::BitArray{1}
    siteflagvar::BitArray{1}
end

function EPFields(N::Int, expval, T)

    siteflagvar = trues(N)
    siteflagave = trues(N)

    expave, expvar = parseexpval!(expval, siteflagave, siteflagvar)
    av = zeros(T,N)
    var = zeros(T,N)

    for (k,v) in expave
        av[k] = v
    end
    for (k,v) in expvar
        var[k] = v
    end

    EPFields(
        #= av          =# av,
        #= va          =# var,
        #= a           =# zeros(T,N),
        #= d           =# ones(T,N),
        #= μ           =# zeros(T,N),
        #= s           =# ones(T,N),
        #= siteflagave =# siteflagave,
        #= siteflagvar =# siteflagvar
    )
end

EPFields(N::Int, ::Nothing, T) = EPFields(
    #= av          =#  zeros(T,N),
    #= va          =#  zeros(T,N),
    #= a           =#  zeros(T,N),
    #= d           =#  ones(T,N),
    #= μ           =#  zeros(T,N),
    #= s           =#  ones(T,N),
    #= siteflagave =#  trues(N),
    #= siteflagvar =#  trues(N)
)

ext: jl


text:
include("filtering.jl");
##
list = ["PreVehicle",
    "Altanserin",
    "SB242084",
    "Way_100135",
    "Methysergide",
    "Citalopram",
    "SB242084_opt",
    "Saline"]
s = filter(r->r.Treatment in list &&
    r.Trial < 61 &&
    # r.Trial_Travel_to < 10 &&
    r.MouseID in ["pc1","pc2","pc3","pc4","pc5","pc6","pc8","pc9","pc10"],
    streaks)
GroupC_PreVehicle = [Date(2016,07,22),Date(2016,07,27), Date(2016,07,31)]
GroupD_PreVehicle = [Date(2016,07,21),Date(2016,07,25), Date(2016,07,29)]
s[:,:Treatment] = [r.Group == "Group C" && r.Day in GroupC_PreVehicle ? "Control" : r.Treatment for r in eachrow(s)]
s[:,:Treatment] = [r.Group == "Group D" && r.Day in GroupD_PreVehicle ? "Control" : r.Treatment for r in eachrow(s)]
##
optopharma = filter(r ->
    r.Treatment in ["Control", "SB242084_opt"], s)
gd = groupby(optopharma,[:MouseID,:Treatment,:Stim])
df1 = combine([:Num_pokes,:AfterLast,:Trial_Travel_to] => (n,a,t) ->
    (Num_pokes = mean(n), AfterLast = mean(a), Travel = mean(t))
    ,gd)
# select the variable to analyze
df2 = select(df1, Not([:Travel,:Num_pokes]))
df3 = combine(groupby(df2,:Stim)) do dd
    subdf = unstack(dd,:Treatment,:AfterLast)
    # current_drug = Symbol(subdf[1,:Phase])
    # rename!(subdf, current_drug => :Drug)
end
df4 = combine(groupby(df3,[:Stim])) do dd
    wilcoxon(dd,:SB242084_opt, :Control; f = x -> mean(skipmissing(x)))
end
plot_wilcoxon(df4)

open_html_table(df3)
disallowmissing!(df3)
SignedRankTest(df3[df3.Stim,:SB242084_opt],df3[df3.Stim,:Control])
SignedRankTest(df3[.!(df3.Stim),:SB242084_opt],df3[.!(df3.Stim),:Control])


text:
function NBodyProblem(potential,M,v0,u0,tspan,p=nothing; kwargs...)
    # check number of particles
    @assert length(u0) == length(v0)
    ind = i -> div(i-1, dim) + 1
    N   = length(u0)
    dim = div(N, length(M))
    @assert dim*length(M) == N
    @assert dim ∈ (2, 3)

    let fun = FWrapper{typeof(potential),typeof(tspan[1]),typeof(M),typeof(p)}(potential,tspan[1],M,p), cfg = ForwardDiff.GradientConfig(fun, u0)
        function acceleration!(dv,v,x,p,t)
            fun.t = t
            ForwardDiff.gradient!(dv, fun, x, cfg)
            for x in dv.x
                x ./= -M
            end
        end
        return SecondOrderODEProblem{true}(acceleration!, v0, u0, tspan, p; kwargs...)
    end
end

mutable struct FWrapper{F,T,MType,P}
    f::F
    t::T
    M::MType
    p::P
end
(f::FWrapper)(u) = f.f(f.p,f.t,u.x...,f.M)


text:
using NetworkEpidemics
using LightGraphs
using Statistics
using Plots
using ColorSchemes
using LaTeXStrings
using Random

CorrectedMetapopulation(χ::Real, mp::Metapopulation{SI}) = Metapopulation(mp.h, mp.D, SI(χ*mp.dynamics.β))

##

function rewire_network(g::SimpleGraph, p)
    h = SimpleGraph(g)
    for e in edges(g)
        u = src(e)
        v = dst(e)
        if rand() < p
            if rand() < 0.5
                w = u
            else
                w = v
            end
            x = rand(setdiff(vertices(h), [u,v]))
            # If there's already an edge, we just swap them, so the graph is unchanged
            if !has_edge(h,w,x)
                rem_edge!(h,u,v)
                add_edge!(h,w,x)
            end
        end
    end
    return h
end

##
Random.seed!(2020)

N = 1000
M = 10

h = complete_graph(M)
#g = watts_strogatz(N, 100, 0.0)
g = random_regular_graph(N, 50)

k = degree(g)

χ = mean(k)/N # doesn't change so only compute it once

β = 0.1
D = 0.1

mp = Metapopulation(h, D, SI(β))
mpc = CorrectedMetapopulation(χ, mp)

mpx = Metaplex(g, h, D, SI(β))

x0_i = fill(1, N)
x0_i[1:10] .= 2
x0_μ = rand(1:M, N)
x0_mp = zeros(Int, M, 2)
for i in 1:N
    x0_mp[x0_μ[i], x0_i[i]] += 1
end

tmax = 30.0
nmax = 10000
nsims = 1000
nbins = 200

ts_mp, output_mp = gillespie(mp, x0_mp, tmax=tmax, nmax=nmax)
ts_mpx, output_mpx = gillespie(mpx, [x0_i,x0_μ], tmax=tmax, nmax=nmax)

plot(ts_mp, output_mp[2], label="")
plot(ts_mpx, output_mpx[2], label="")

ts = LinRange(0.0, tmax, nbins)
ps = LinRange(0.0, 1.0, 5)

Random.seed!(2020)

ts_av_mpx, u_av_mpx = average(mpx, [x0_i,x0_μ], tmax=tmax, nbins=nbins,nsims=nsims,nmax=nmax)

ts_mf_mp, u_mf_mp = meanfield(mp, x0_mp, tmax=tmax, saveat=ts)
ts_mf_mpc, u_mf_mpc = meanfield(mpc, x0_mp, tmax=tmax, saveat=ts)

u_mf_mpx = zeros(length(ps), length(ts))
d_stds = Array{Float64,2}(undef, length(ps), 100)
for i in 1:length(ps), j in 1:100
    f = rewire_network(g, ps[i])
    d_stds[i,j] = std(degree(f))
    local mpx = Metaplex(f, h, D, SI(β))
    t, u = meanfield(mpx, [x0_i,x0_μ], tmax=tmax, saveat=ts)
    u_mf_mpx[i,:] += sum(u[2], dims=2)
end
u_mf_mpx /= 100

colors = ColorSchemes.GnBu_9;
orange = ColorSchemes.Oranges_3[3];

plot(ts, sum(u_mf_mpc[2], dims=2)/N,
    #label="Meanfield (No correction)",
    label = "Mean degree correction",
    xlabel="Time",
    ylabel="Fraction of infected individuals",
    legend=:bottomright,
    color = orange
    )
    plot!(ts, u_mf_mpx[1,:]/N,
    label = "p=0.0",
    color = colors[5]
    )
    plot!(ts, u_mf_mpx[2,:]/N,
    label = "p=0.25",
    color = colors[6]
    )
    plot!(ts, u_mf_mpx[3,:]/N,
    label = "p=0.5",
    color = colors[7]
    )
    plot!(ts, u_mf_mpx[4,:]/N,
    label = "p=0.75",
    color = colors[8]
    )
    plot!(ts, u_mf_mpx[5,:]/N,
    label = "p=1.0",
    color = colors[9]
    )
    plot!(ts, sum(u_mf_mpc[2], dims=2)/N,
    label="",
    color = orange
    )

# necessary for consistent LaTeX
pyplot()

plot(ps, mean(d_stds, dims=2).^2,
    label="",
    xlabel="p",
    ylabel=L"\sigma^2")

cols = map(x->[x.r, x.g, x.b], colors.colors)
or = [orange.r, orange.g, orange.b]

##
using MATLAB

mat"
figure
xlabel('Time')
ylabel('Fraction of Infected Individuals')
hold on
plot($ts, $(sum(u_mf_mpc[2], dims=2)/N), 'Color', $(or))
plot($ts, $(u_mf_mpx[1,:]/N), 'Color', $(cols[5]))
plot($ts, $(u_mf_mpx[2,:]/N), 'Color', $(cols[6]))
plot($ts, $(u_mf_mpx[3,:]/N), 'Color', $(cols[7]))
plot($ts, $(u_mf_mpx[4,:]/N), 'Color', $(cols[8]))
plot($ts, $(u_mf_mpx[5,:]/N), 'Color', $(cols[9]))
hold off
legend({'Mean degree correction', 'p=0.0', 'p=0.25', 'p=0.5', 'p=0.75', 'p=1.0'}, 'Location','southeast')
"

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
module IndexableBitVectors

export
    # types
    AbstractIndexableBitVector, AbstractBitVector,
    CompactBitVector, SucVector, CSucVector, RRR, RRRNP,
    # query operations
    rank, rank0, rank1, select, select0, select1

import Base: rank, select, show, size, push!, getindex, length, endof, sizeof, convert
using Compat
using Switch

abstract AbstractIndexableBitVector <: DenseArray{Bool,1}
typealias AbstractBitVector Union(BitVector,AbstractIndexableBitVector)

include("common.jl")
include("bitvector.jl")
include("compactbitvector.jl")
include("sucvector.jl")
include("csucvector.jl")
include("rrr.jl")
include("derived.jl")

end # module


text:
@testset "state" begin
    include("state/core.jl")
    include("state/hold.jl")
    include("state/wrap.jl")
    include("state/advance.jl")
    include("state/preserve.jl")
    include("state/tabulate.jl")
    include("state/interpolate.jl")
    include("state/track.jl")
    include("state/remember.jl")
    include("state/provide.jl")
    include("state/drive.jl")
    include("state/call.jl")
    include("state/integrate.jl")
    include("state/accumulate.jl")
    include("state/capture.jl")
    include("state/flag.jl")
    include("state/produce.jl")
    include("state/bisect.jl")
    include("state/solve.jl")
end


text:
# Miscellaneous functions

"""
    logmvgamma(p::Int, a::Real)

Return the logarithm of [multivariate gamma function](https://en.wikipedia.org/wiki/Multivariate_gamma_function) ([DLMF 35.3.1](https://dlmf.nist.gov/35.3.1)).
"""
function logmvgamma(p::Int, a::Real)
    # NOTE: one(a) factors are here to prevent unnecessary promotion of Float32
    res = p * (p - 1) * log(pi * one(a)) / 4
    for ii in 1:p
        res += loggamma(a + (1 - ii) * one(a) / 2)
    end
    return res
end

"""
    logmvbeta(p::Int, a::Real, b::Real)

Return the logarithm of the multivariate beta function ([DLMF 35.3.7](https://dlmf.nist.gov/35.3#E7)).
"""
logmvbeta(p::Int, a::Real, b::Real) = logmvgamma(p, a) + logmvgamma(p, b) - logmvgamma(p, a + b)

"""
    lstirling_asym(x)

The remainder term after
[Stirling's approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation)
to [`loggamma`](@ref):

```math
\\log \\Gamma(x) \\approx x \\log(x) - x + \\log(2\\pi/x)/2 = \\log(x)*(x-1/2) + \\log(2\\pi)/2 - x
```

In Julia syntax, this means:

    lstirling_asym(x) = loggamma(x) + x - (x-0.5)*log(x) - 0.5*log(2π)

For sufficiently large `x`, this can be approximated using the asymptotic
_Stirling's series_ ([DLMF 5.11.1](https://dlmf.nist.gov/5.11.1)):

```math
\\frac{1}{12x} - \\frac{1}{360x^3} + \\frac{1}{1260x^5} - \\frac{1}{1680x^7} + \\ldots
```

The truncation error is bounded by the first omitted term, and is of the same sign.

Relative error of approximation is bounded by
    (174611/125400 x^-19) / (1/12 x^-1 - 1/360 x^-3)
which is < 1/2 ulp for x >= 10.0, and total numeric error appears to be < 2 ulps

# References

* Temme, N. (1996) Special functions: An introduction to the classical functions of
   mathematical physics, Wiley, New York, ISBN: 0-471-11313-1, Chapter 3.6, pp 61-65.
* Weisstein, Eric W. ["Stirling's Series."](http://mathworld.wolfram.com/StirlingsSeries.html).
  MathWorld.
* [OEIS A046968](http://oeis.org/A046968) and [OEIS A046969](http://oeis.org/A046969)
  for the series coefficients
"""
function lstirling_asym end

lstirling_asym(x::BigFloat) = loggamma(x) + x - log(x)*(x - big(0.5)) - log2π/big(2)

lstirling_asym(x::Integer) = lstirling_asym(float(x))

const lstirlingF64 = Float64[lstirling_asym(k) for k in big(1):big(64)]
const lstirlingF32 = Float64[lstirling_asym(k) for k in big(1):big(40)]

function lstirling_asym(x::Float64)
    isinteger(x) && (0 < x ≤ length(lstirlingF64)) && return lstirlingF64[Int(x)]
    t = inv(abs2(x))
    @horner(t,
             8.33333333333333333e-2, #  1/12 x^-1
            -2.77777777777777778e-3, # -1/360 x^-3
             7.93650793650793651e-4, #  1/1260 x^-5
            -5.95238095238095238e-4, # -1/1680 x^-7
             8.41750841750841751e-4, #  1/1188 x^-9
            -1.91752691752691753e-3, # -691/360360 x^-11
             6.41025641025641026e-3, #  1/156 x^-13
            -2.95506535947712418e-2, # -3617/122400 x^-15
             1.79644372368830573e-1)/x #  43867/244188 x^-17
end

function lstirling_asym(x::Float32)
    isinteger(x) && (0 < x ≤ length(lstirlingF32)) && return lstirlingF32[Int(x)]
    t = inv(abs2(x))
    @horner(t,
             8.333333333333f-2, #  1/12 x^-1
            -2.777777777777f-3, # -1/360 x^-3
             7.936507936508f-4, #  1/1260 x^-5
            -5.952380952381f-4, # -1/1680 x^-7
             8.417508417508f-4)/x #  1/1188 x^-9
end

@doc raw"""
owent(h,a) : Returns the value of Owen's T function for (h,a)

Owen's T function:
```math
T(h,a) = \frac{1}{2\pi } \int_{0}^{a} \frac{e^{-\frac{1}{2}h^2(1+x^2)}}{1+x^2}dx\quad(-\infty < h,a < +\infty)
```

For *h* and *a* > 0, *T(h,a)* gives the volume of the uncorrelated bivariate normal distribution with zero mean and unit variance
over the area from *y = ax* and *y = 0* and to the right of *x = h*.

EXAMPLE:
```
julia> owent(0.0625, 0.025)
0.003970281304296922
```

Worst case accuracy is about 2E-16.
"""
function owent(h::T1,a::T1) where {T1 <: Real}

    #*********************
    # shortcut evaluations
    #*********************

    if h < zero(h)
        return owent(abs(h),a)
    end

    if h==zero(h)
        return T1(atan(a)*inv2π)
    end

    if a < zero(a)
        return -owent(h,abs(a))
    end

    if a==zero(a)
        return zero(a)
    end

    if a==one(a)
        return T1(0.125*erfc(-h*invsqrt2)*erfc(h*invsqrt2))
    end

    if a==Inf
        return T1(0.25*erfc(sqrt(h^2)*invsqrt2))
    end

    # below reduces the range from -inf < h,a < +inf to h ≥ 0, 0 ≤ a ≤ 1
    if a > one(a)
        return 0.25*(erfc(-h*invsqrt2) + erfc(-a*h*invsqrt2)) - 0.25*erfc(-h*invsqrt2)*erfc(-a*h*invsqrt2) - owent(a*h,one(a)/a)
    end

    # calculate Owen's T

    if a ≤ 0.999999

        t2(h,a,x) = inv4π*a*exp(-0.5*(h^2)*(one(h)+(a*x)^2))/(one(h)+(a*x)^2)

        x = [-0.9987710072524261, -0.9935301722663508, -0.9841245837228269, -0.9705915925462473, -0.9529877031604309, -0.9313866907065543, -0.9058791367155696
        , -0.8765720202742479, -0.8435882616243935, -0.8070662040294426, -0.7671590325157404, -0.7240341309238146, -0.6778723796326639, -0.6288673967765136
        , -0.5772247260839727, -0.5231609747222331, -0.4669029047509584, -0.4086864819907167, -0.34875588629216075, -0.28736248735545555, -0.22476379039468905
        , -0.1612223560688917, -0.0970046992094627, -0.03238017096286937, 0.03238017096286937, 0.0970046992094627, 0.1612223560688917, 0.22476379039468905
        , 0.28736248735545555, 0.34875588629216075, 0.4086864819907167, 0.4669029047509584, 0.5231609747222331, 0.5772247260839727, 0.6288673967765136
        , 0.6778723796326639, 0.7240341309238146, 0.7671590325157404, 0.8070662040294426, 0.8435882616243935, 0.8765720202742479, 0.9058791367155696
        , 0.9313866907065543, 0.9529877031604309, 0.9705915925462473, 0.9841245837228269, 0.9935301722663508, 0.9987710072524261]

        w = [0.0031533460523059122, 0.0073275539012762885, 0.011477234579234613, 0.015579315722943824, 0.01961616045735561, 0.023570760839324363
        , 0.027426509708356944, 0.031167227832798003, 0.03477722256477052, 0.038241351065830737, 0.04154508294346467, 0.0446745608566943, 0.04761665849249045
        , 0.05035903555385445, 0.05289018948519363, 0.055199503699984116, 0.05727729210040322, 0.05911483969839564, 0.06070443916589387, 0.06203942315989268
        , 0.06311419228625402, 0.06392423858464813, 0.06446616443594998, 0.06473769681268386, 0.06473769681268386, 0.06446616443594998, 0.06392423858464813
        , 0.06311419228625402, 0.06203942315989268, 0.06070443916589387, 0.05911483969839564, 0.05727729210040322, 0.055199503699984116
        , 0.05289018948519363, 0.05035903555385445, 0.04761665849249045, 0.0446745608566943, 0.04154508294346467, 0.038241351065830737, 0.03477722256477052
        , 0.031167227832798003, 0.027426509708356944, 0.023570760839324363, 0.01961616045735561, 0.015579315722943824, 0.011477234579234613, 0.0073275539012762885
        , 0.0031533460523059122]

        towen = dot(w, t2.(h,a,x))

        return T1(towen)
    else
        # a > 0.999999, T6 from paper (quadrature using QuadGK would also work, but be slower)

        j = 0.5*erfc(-h*invsqrt2)
        k = atan((one(a)-a)/(one(a)+a))
        towen = 0.5*j*(one(h)-j)-inv2π*k*exp((-0.5*(one(a)-a)*h^2)/k)

        return T1(towen)
    end
end

owent(h::Real, a::Real) = owent(promote(h,a)...)


text:
"""
    HMC(n_iters::Int, epsilon::Float64, tau::Int)

Hamiltonian Monte Carlo sampler.

Usage:

```julia
HMC(1000, 0.05, 10)
```

Example:

```julia
# Define a simple Normal model with unknown mean and variance.
@model gdemo(x) = begin
  s ~ InverseGamma(2,3)
  m ~ Normal(0, sqrt(s))
  x[1] ~ Normal(m, sqrt(s))
  x[2] ~ Normal(m, sqrt(s))
  return s, m
end

sample(gdemo([1.5, 2]), HMC(1000, 0.05, 10))
```
"""
mutable struct HMC{T} <: Hamiltonian
  n_iters   ::  Int       # number of samples
  epsilon   ::  Float64   # leapfrog step size
  tau       ::  Int       # leapfrog step number
  space     ::  Set{T}    # sampling space, emtpy means all
  gid       ::  Int       # group ID
end
function HMC(epsilon::Float64, tau::Int, space...)
  _space = isa(space, Symbol) ? Set([space]) : Set(space)
  HMC(1, epsilon, tau, _space, 0)
end
function HMC(n_iters::Int, epsilon::Float64, tau::Int)
  HMC(n_iters, epsilon, tau, Set(), 0)
end
function HMC(n_iters::Int, epsilon::Float64, tau::Int, space...)
  _space = isa(space, Symbol) ? Set([space]) : Set(space)
  HMC(n_iters, epsilon, tau, _space, 0)
end
function HMC(alg::HMC, new_gid::Int)
  HMC(alg.n_iters, alg.epsilon, alg.tau, alg.space, new_gid)
end
HMC{T}(alg::HMC, new_gid::Int) where {T} = HMC(alg, new_gid)

# Below is a trick to remove the dependency of Stan by Requires.jl
# Please see https://github.com/TuringLang/Turing.jl/pull/459 for explanations
DEFAULT_ADAPT_CONF_TYPE = Nothing
STAN_DEFAULT_ADAPT_CONF = nothing
#  @init @require Stan="682df890-35be-576f-97d0-3d8c8b33a550" begin
  DEFAULT_ADAPT_CONF_TYPE = Union{DEFAULT_ADAPT_CONF_TYPE,Stan.Adapt}
  STAN_DEFAULT_ADAPT_CONF = Stan.Adapt()
#  end

# NOTE: the implementation of HMC is removed,
#       it now reuses the one of HMCDA
Sampler(alg::HMC) = Sampler(alg::HMC, STAN_DEFAULT_ADAPT_CONF::DEFAULT_ADAPT_CONF_TYPE)
Sampler(alg::HMC, adapt_conf::DEFAULT_ADAPT_CONF_TYPE) = begin
  spl = Sampler(HMCDA(alg.n_iters, 0, 0.0, alg.epsilon * alg.tau, alg.space, alg.gid), adapt_conf)
  spl.info[:pre_set_ϵ] = alg.epsilon
  spl
end

Sampler(alg::Hamiltonian) =  Sampler(alg, STAN_DEFAULT_ADAPT_CONF::DEFAULT_ADAPT_CONF_TYPE)
Sampler(alg::Hamiltonian, adapt_conf::DEFAULT_ADAPT_CONF_TYPE) = begin
  info=Dict{Symbol, Any}()

  # For sampler infomation
  info[:accept_his] = []
  info[:lf_num] = 0
  info[:total_lf_num] = 0
  info[:total_eval_num] = 0

  # For pre-conditioning
  info[:θ_mean] = nothing
  info[:θ_num] = 0
  info[:stds] = nothing
  info[:vars] = nothing
  info[:ad] = Dict()

  # For caching gradient
  info[:grad_cache] = Dict{UInt64,Vector}()
  info[:reverse_diff_cache] = Dict()

  # Adapt configuration
  if adapt_conf != nothing
    info[:adapt_conf] = adapt_conf
  end

  Sampler(alg, info)
end

function sample(model::Function, alg::T;
                                chunk_size=CHUNKSIZE[],               # set temporary chunk size
                                save_state=false,                   # flag for state saving
                                resume_from=nothing,                # chain to continue
                                reuse_spl_n=0,                      # flag for spl re-using
                                adapt_conf=STAN_DEFAULT_ADAPT_CONF  # adapt configuration
                               ) where T<:Hamiltonian

  default_chunk_size = CHUNKSIZE[]  # record global chunk size
  setchunksize(chunk_size)        # set temp chunk size

  spl = reuse_spl_n > 0 ?
        resume_from.info[:spl] :
        Sampler(alg, adapt_conf)

  @assert isa(spl.alg, Hamiltonian) "[Turing] alg type mismatch; please use resume() to re-use spl"

  alg_str = isa(alg, HMC)   ? "HMC"   :
            isa(alg, HMCDA) ? "HMCDA" :
            isa(alg, SGHMC) ? "SGHMC" :
            isa(alg, SGLD)  ? "SGLD"  :
            isa(alg, NUTS)  ? "NUTS"  : "Hamiltonian"

  # Initialization
  time_total = zero(Float64)
  n = reuse_spl_n > 0 ?
      reuse_spl_n :
      alg.n_iters
  samples = Array{Sample}(undef, n)
  weight = 1 / n
  for i = 1:n
    samples[i] = Sample(weight, Dict{Symbol, Any}())
  end

  vi = resume_from == nothing ?
       Base.invokelatest(model, VarInfo(), nothing) :
       deepcopy(resume_from.info[:vi])

  if spl.alg.gid == 0
    link!(vi, spl)
    runmodel(model, vi, spl)
  end

  # HMC steps
  PROGRESS[] && (spl.info[:progress] = ProgressMeter.Progress(n, 1, "[$alg_str] Sampling...", 0))
  for i = 1:n
    @debug "$alg_str stepping..."

    time_elapsed = @elapsed vi = step(model, spl, vi, i == 1)
    time_total += time_elapsed

    if spl.info[:accept_his][end]     # accepted => store the new predcits
      samples[i].value = Sample(vi, spl).value
    else                              # rejected => store the previous predcits
      samples[i] = samples[i - 1]
    end
    samples[i].value[:elapsed] = time_elapsed
    samples[i].value[:lf_eps] = spl.info[:wum][:ϵ][end]

    PROGRESS[] && ProgressMeter.next!(spl.info[:progress])
  end

  println("[$alg_str] Finished with")
  println("  Running time        = $time_total;")
  if ~isa(alg, NUTS)  # accept rate for NUTS is meaningless - so no printing
    accept_rate = sum(spl.info[:accept_his]) / n  # calculate the accept rate
    println("  Accept rate         = $accept_rate;")
  end
  println("  #lf / sample        = $(spl.info[:total_lf_num] / n);")
  println("  #evals / sample     = $(spl.info[:total_eval_num] / n);")
  stds_str = string(spl.info[:wum][:stds])
  stds_str = length(stds_str) >= 32 ? stds_str[1:30]*"..." : stds_str   # only show part of pre-cond
  println("  pre-cond. diag mat  = $(stds_str).")

  setchunksize(default_chunk_size)      # revert global chunk size

  if resume_from != nothing   # concat samples
    pushfirst!(samples, resume_from.value2...)
  end
  c = Chain(0, samples)       # wrap the result by Chain
  if save_state               # save state
    # Convert vi back to X if vi is required to be saved
    if spl.alg.gid == 0 invlink!(vi, spl) end
    spl.info[:grad_cache] = Dict{UInt64,Vector}()
    spl.info[:reverse_diff_cache] = Dict()
    save!(c, spl, model, vi)
  end

  c
end

assume(spl::Sampler{T}, dist::Distribution, vn::VarName, vi::VarInfo) where T<:Hamiltonian = begin
  @debug "assuming..."
  updategid!(vi, vn, spl)
  r = vi[vn]
  # acclogp!(vi, logpdf_with_trans(dist, r, istrans(vi, vn)))
  # r
  @debug "dist = $dist"
  @debug "vn = $vn"
  @debug "r = $r" "typeof(r)=$(typeof(r))"
  r, logpdf_with_trans(dist, r, istrans(vi, vn))
end

assume(spl::Sampler{A}, dists::Vector{D}, vn::VarName, var::Any, vi::VarInfo) where {A<:Hamiltonian,D<:Distribution} = begin
  @assert length(dists) == 1 "[observe] Turing only support vectorizing i.i.d distribution"
  dist = dists[1]
  n = size(var)[end]

  vns = map(i -> copybyindex(vn, "[$i]"), 1:n)

  rs = vi[vns]  # NOTE: inside Turing the Julia conversion should be sticked to

  # acclogp!(vi, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1]))))

  if isa(dist, UnivariateDistribution) || isa(dist, MatrixDistribution)
    @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
    var = rs
  elseif isa(dist, MultivariateDistribution)
    if isa(var, Vector)
      @assert length(var) == size(rs)[2] "[assume] variable and random number dimension unmatched"
      for i = 1:n
        var[i] = rs[:,i]
      end
    elseif isa(var, Matrix)
      @assert size(var) == size(rs) "[assume] variable and random number dimension unmatched"
      var = rs
    else
      error("[Turing] unsupported variable container")
    end
  end

  var, sum(logpdf_with_trans(dist, rs, istrans(vi, vns[1])))
end

observe(spl::Sampler{A}, d::Distribution, value::Any, vi::VarInfo) where A<:Hamiltonian=
  observe(nothing, d, value, vi)

observe(spl::Sampler{A}, ds::Vector{D}, value::Any, vi::VarInfo) where {A<:Hamiltonian,D<:Distribution} =
  observe(nothing, ds, value, vi)


text:
# MadNLP.jl
# Created by Sungho Shin (sungho.shin@wisc.edu)

module MadNLPMa77

import ..MadNLPHSL:
    @kwdef, Logger, @debug, @warn, @error, libhsl,
    SparseMatrixCSC, SparseMatrixCSC, SubVector, StrideOneVector,
    get_tril_to_full, transfer!,
    AbstractOptions, AbstractLinearSolver, set_options!,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia
import ..MadNLPHSL: Mc68

const INPUT_MATRIX_TYPE = :csc

@enum(Ordering::Int,AMD = 1, METIS = 3)

@kwdef mutable struct Options <: AbstractOptions
    ma77_buffer_lpage::Int = 4096
    ma77_buffer_npage::Int = 1600
    ma77_file_size::Int = 2097152
    ma77_maxstore::Int = 0
    ma77_nemin::Int = 8
    ma77_order::Ordering = METIS
    ma77_print_level::Int = -1
    ma77_small::Float64 = 1e-20
    ma77_static::Float64 = 0.
    ma77_u::Float64 = 1e-8
    ma77_umax::Float64 = 1e-4
end

@kwdef mutable struct Control
    f_arrays::Cint = 0
    print_level::Cint = 0
    unit_diagnostics::Cint = 0
    unit_error::Cint = 0
    unit_warning::Cint = 0
    bits::Cint = 0

    buffer_lpage_1::Cint = 0
    buffer_lpage_2::Cint = 0
    buffer_npage_1::Cint = 0
    buffer_npage_2::Cint = 0

    file_size::Clong = 0
    maxstore::Clong = 0

    storage_1::Clong = 0
    storage_2::Clong = 0
    storage_3::Clong = 0

    nemin::Cint = 0
    maxit::Cint = 0
    infnorm::Cint = 0
    thresh::Cdouble = 0.
    nb54::Cint = 0
    action::Cint = 0
    multiplier::Cdouble = 0.
    nb64::Cint = 0
    nbi::Cint = 0
    small::Cdouble = 0.
    static_::Cdouble = 0.
    storage_indef::Clong = 0
    u::Cdouble = 0.
    umin::Cdouble = 0.
    consist_tol::Cdouble = 0.

    ispare_1::Cint = 0
    ispare_2::Cint = 0
    ispare_3::Cint = 0
    ispare_4::Cint = 0
    ispare_5::Cint = 0

    lspare_1::Clong = 0
    lspare_2::Clong = 0
    lspare_3::Clong = 0
    lspare_4::Clong = 0
    lspare_5::Clong = 0

    rspare_1::Cdouble = 0.
    rspare_2::Cdouble = 0.
    rspare_3::Cdouble = 0.
    rspare_4::Cdouble = 0.
    rspare_5::Cdouble = 0.
end

@kwdef mutable struct Info
    detlog::Cdouble = 0.
    detsign::Cint = 0
    flag::Cint = 0
    iostat::Cint = 0
    matrix_dup::Cint = 0
    matrix_rank::Cint = 0
    matrix_outrange::Cint = 0
    maxdepth::Cint = 0
    maxfront::Cint = 0
    minstore::Clong = 0
    ndelay::Cint = 0
    nfactor::Clong = 0
    nflops::Clong = 0
    niter::Cint = 0
    nsup::Cint = 0
    num_neg::Cint = 0
    num_nothresh::Cint = 0
    num_perturbed::Cint = 0
    ntwo::Cint = 0

    stat_1::Cint = 0
    stat_2::Cint = 0
    stat_3::Cint = 0
    stat_4::Cint = 0

    nio_read_1::Clong = 0 # 2
    nio_read_2::Clong = 0 # 2

    nio_write_1::Clong = 0 # 2
    nio_write_2::Clong = 0 # 2

    nwd_read_1::Clong = 0 # 2
    nwd_read_2::Clong = 0 # 2

    nwd_write_1::Clong = 0 # 2
    nwd_write_2::Clong = 0 # 2

    num_file_1::Cint = 0 # 4
    num_file_2::Cint = 0 # 4
    num_file_3::Cint = 0 # 4
    num_file_4::Cint = 0 # 4

    storage_1::Clong = 0 # 4
    storage_2::Clong = 0 # 4
    storage_3::Clong = 0 # 4
    storage_4::Clong = 0 # 4

    tree_nodes::Cint = 0
    unit_restart::Cint = 0
    unused::Cint = 0
    usmall::Cdouble = 0.

    ispare_1::Cint = 0
    ispare_2::Cint = 0
    ispare_3::Cint = 0
    ispare_4::Cint = 0
    ispare_5::Cint = 0

    lspare_1::Clong = 0
    lspare_2::Clong = 0
    lspare_3::Clong = 0
    lspare_4::Clong = 0
    lspare_5::Clong = 0

    rspare_1::Cdouble = 0.
    rspare_2::Cdouble = 0.
    rspare_3::Cdouble = 0.
    rspare_4::Cdouble = 0.
    rspare_5::Cdouble = 0.
end

mutable struct Solver <: AbstractLinearSolver
    tril::SparseMatrixCSC{Float64,Int32}
    full::SparseMatrixCSC{Float64,Int32}
    tril_to_full_view::SubVector{Float64}

    control::Control
    info::Info

    mc68_control::Mc68.Control
    mc68_info::Mc68.Info

    order::Vector{Int32}
    keep::Vector{Ptr{Nothing}}

    opt::Options
    logger::Logger
end

ma77_default_control_d(control::Control) = ccall(
    (:ma77_default_control_d,libhsl),
    Cvoid,
    (Ref{Control},),
    control)
ma77_open_d(n::Cint,fname1::String,fname2::String,fname3::String,fname4::String,
            keep::Vector{Ptr{Cvoid}},control::Control,info::Info) = ccall(
                (:ma77_open_d,libhsl),
                Cvoid,
                (Cint,Ptr{Cchar},Ptr{Cchar},Ptr{Cchar},Ptr{Cchar},
                 Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info}),
                n,fname1,fname2,fname3,fname4,keep,control,info)
ma77_input_vars_d(idx::Cint,nvar::Cint,list::StrideOneVector{Cint},
                  keep::Vector{Ptr{Cvoid}},control::Control,info::Info) = ccall(
                      (:ma77_input_vars_d,libhsl),
                      Cvoid,
                      (Cint,Cint,Ptr{Cint},
                       Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info}),
                      idx,nvar,list,keep,control,info)
ma77_input_reals_d(idx::Cint,length::Cint,reals::StrideOneVector{Cdouble},
                   keep::Vector{Ptr{Cvoid}},control::Control,info::Info) = ccall(
                       (:ma77_input_reals_d,libhsl),
                       Cvoid,
                       (Cint,Cint,Ptr{Cdouble},
                        Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info}),
                       idx,length,reals,keep,control,info)
ma77_analyse_d(order::Vector{Cint},
               keep::Vector{Ptr{Cvoid}},control::Control,info::Info) = ccall(
                   (:ma77_analyse_d,libhsl),
                   Cvoid,
                   (Ptr{Cint},Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info}),
                   order,keep,control,info)
ma77_factor_d(posdef::Cint,keep::Vector{Ptr{Cvoid}},control::Control,info::Info,
              scale::Ptr{Nothing}) = ccall(
                  (:ma77_factor_d,libhsl),
                  Cvoid,
                  (Cint,Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info},Ptr{Nothing}),
                  posdef,keep,control,info,scale)
ma77_factor_solve_d(posdef::Cint,keep::Vector{Ptr{Cvoid}},control::Control,info::Info,
                    scale::Vector{Cdouble},nrhs::Cint,lx::Cint,rhs::Vector{Cdouble}) = ccall(
                        (:ma77_factor_solve_d,libhsl),
                        Cvoid,
                        (Cint,Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info},
                         Ptr{Cdouble},Cint,Cint,Ptr{Cdouble}),
                        posdef,keep,control,info,scale,nrhs,lx,rhs);
ma77_solve_d(job::Cint,nrhs::Cint,lx::Cint,x::Vector{Cdouble},
             keep::Vector{Ptr{Cvoid}},control::Control,info::Info,
             scale::Ptr{Nothing})=ccall(
                 (:ma77_solve_d,libhsl),
                 Cvoid,
                 (Cint,Cint,Cint,Ptr{Float64},
                  Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info},Ptr{Nothing}),
                 job,nrhs,lx,x,keep,control,info,scale);
ma77_finalize_d(keep::Vector{Ptr{Cvoid}},control::Control,info::Info) = ccall(
    (:ma77_finalise_d,libhsl),
    Cvoid,
    (Ptr{Ptr{Cvoid}},Ref{Control},Ref{Info}),
    keep,control,info)

function Solver(csc::SparseMatrixCSC{Float64,Int32};
                option_dict::Dict{Symbol,Any}=Dict{Symbol,Any}(),
                opt=Options(),logger=Logger(),
                kwargs...)

    set_options!(opt,option_dict,kwargs)

    full,tril_to_full_view = get_tril_to_full(csc)
    order = Vector{Int32}(undef,csc.n)

    mc68_info = Mc68.Info()
    mc68_control = Mc68.get_mc68_default_control()

    keep = [C_NULL]

    mc68_control.f_array_in=1
    mc68_control.f_array_out=1
    Mc68.mc68_order_i(Int32(opt.ma77_order),Int32(csc.n),csc.colptr,csc.rowval,order,mc68_control,mc68_info)

    info=Info()
    control=Control()
    ma77_default_control_d(control)
    control.f_arrays = 1
    control.bits = 32
    control.file_size = opt.ma77_file_size
    control.maxstore = opt.ma77_maxstore
    control.print_level = -1

    control.buffer_lpage_1=opt.ma77_buffer_lpage
    control.buffer_lpage_2=opt.ma77_buffer_lpage
    control.buffer_npage_1=opt.ma77_buffer_npage
    control.buffer_npage_2=opt.ma77_buffer_npage

    control.nemin = opt.ma77_nemin
    control.small = opt.ma77_small
    control.static_ = opt.ma77_static
    control.u = opt.ma77_u

    isfile(".ma77_int")   && rm(".ma77_int")
    isfile(".ma77_real")  && rm(".ma77_real")
    isfile(".ma77_work")  && rm(".ma77_work")
    isfile(".ma77_delay") && rm(".ma77_delay")

    ma77_open_d(Int32(full.n),".ma77_int", ".ma77_real", ".ma77_work", ".ma77_delay",
                keep,control,info)

    info.flag < 0 && throw(SymbolicException())

    for i=1:full.n
        ma77_input_vars_d(Int32(i),full.colptr[i+1]-full.colptr[i],
                          view(full.rowval,full.colptr[i]:full.colptr[i+1]-1),
                          keep,control,info);
        info.flag < 0 && throw(LinearSymbolicException())
    end

    ma77_analyse_d(order,keep,control,info)
    info.flag<0 && throw(SymbolicException())

    M = Solver(csc,full,tril_to_full_view,
               control,info,mc68_control,mc68_info,order,keep,opt,logger)
    finalizer(finalize,M)
    return M
end

function factorize!(M::Solver)
    M.full.nzval.=M.tril_to_full_view
    for i=1:M.full.n
        ma77_input_reals_d(Int32(i),M.full.colptr[i+1]-M.full.colptr[i],
                           view(M.full.nzval,M.full.colptr[i]:M.full.colptr[i+1]-1),
                           M.keep,M.control,M.info)
        M.info.flag < 0 && throw(FactorizationException())
    end
    ma77_factor_d(Int32(0),M.keep,M.control,M.info,C_NULL)
    M.info.flag < 0 && throw(FactorizationException())
    return M
end
function solve!(M::Solver,rhs::StrideOneVector{Float64})
    ma77_solve_d(Int32(0),Int32(1),Int32(M.full.n),rhs,M.keep,M.control,M.info,C_NULL);
    M.info.flag < 0 && throw(SolveException())
    return rhs
end

is_inertia(::Solver) = true
function inertia(M::Solver)
    return (M.info.matrix_rank-M.info.num_neg,M.full.n-M.info.matrix_rank,M.info.num_neg)
end

finalize(M::Solver) = ma77_finalize_d(M.keep,M.control,M.info)

function improve!(M::Solver)
    if M.control.u == M.opt.ma77_umax
        @debug(M.logger,"improve quality failed.")
        return false
    end
    M.control.u = min(M.opt.ma77_umax,M.control.u^.75)
    @debug(M.logger,"improved quality: pivtol = $(M.control.u)")
    return true
end

introduce(::Solver)="ma77"

end


text:
@safetestset "PaduaTransforms" begin include("test_PaduaTransforms.jl") end

@safetestset "Differentiation" begin
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation:
        DiffPlan, differentiate!

    @test DiffPlan{Float64}(5).D == [0  1  0  3  0   5;
                                     0  0  4  0  8   0;
                                     0  0  0  6  0  10;
                                     0  0  0  0  8   0;
                                     0  0  0  0  0  10;
                                     0  0  0  0  0   0]

    @test DiffPlan{Float64}(4).D == [0  1  0  3  0;
                                     0  0  4  0  8;
                                     0  0  0  6  0;
                                     0  0  0  0  8;
                                     0  0  0  0  0]

    coeffs = [1  2  3  0;
              4  5  0  0;
              6  0  0  0;
              0  0  0  0]

    P = DiffPlan{Float64}(3)

    ∂x = zeros(4, 4)
    ∂y = zeros(4, 4)

    differentiate!(∂x, ∂y, P, coeffs)

    @test ∂x == [2  12  0  0;
                 5   0  0  0;
                 0   0  0  0;
                 0   0  0  0]

    @test ∂y == [ 4  5  0  0;
                 24  0  0  0;
                  0  0  0  0;
                  0  0  0  0]

    coeffsarr = rand(4, 4, 6)
    ∂xarr = zeros(4, 4, 6)
    ∂yarr = zeros(4, 4, 6)

    differentiate!(∂xarr, ∂yarr, P, coeffsarr)

    for i in 1:6
        differentiate!(∂x, ∂y, P, coeffsarr[:, :, i])
        @test ∂xarr[:, :, i] == ∂x
        @test ∂yarr[:, :, i] == ∂y
    end
end

@safetestset "Integration" begin
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation:
        getintegrator, integrate

    # integrating odd polynomials over symmetric boundary conditions gives 0
    @test all(getintegrator(100)[2:2:end] .== 0)
    @test all(getintegrator(999)[2:2:end] .== 0)

    coeffs = [1 4 7 0;
              2 5 8 0;
              3 6 9 0;
		      0 0 0 0]

    @test integrate(coeffs, getintegrator(3)) ≈ 1 * 4 + 3 * -4/3 + 7 * -4/3 + 9 * 4/9 atol=5eps()
end

@safetestset "getintegrand! with OOPIntPlan" begin
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation.PaduaTransforms
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation:
        DiffPlan, differentiate!, OOPIntPlan, getintegrand!
    using PoincareInvariants.CanonicalSymplecticStructures

    @testset "6 Points in 2 Dimensions" begin
        degree = 2
        D = 2
        phasecoeffs = cat(
            Float64[1 2 3;
                    4 5 0;
                    6 0 0],
            Float64[1 4 6;
                    2 5 0;
                    3 0 0]
        ; dims=3)

        ∂xcoeffs, ∂ycoeffs = differentiate!(similar(phasecoeffs), similar(phasecoeffs),
            DiffPlan{Float64}(degree), phasecoeffs)

        phasepoints = invpaduatransform!(zeros(getpaduanum(degree), D),
            InvPaduaTransformPlan{Float64}(degree), phasecoeffs)

        Ω(v, t, p) = [0 -1; 1  0]

        plan = OOPIntPlan{Float64}(D, degree)

        intcoeffs = zeros(degree+1, degree+1)

        getintegrand!(intcoeffs, plan, Ω, phasepoints, 0, nothing, ∂xcoeffs, ∂ycoeffs)

        @test intcoeffs ≈ [-72   -82  -30;
                           -82  -432    0;
                           -30     0    0] atol=5eps()
    end

    @testset "$N Points in $D Dimensions" for D in [4, 20], N in [53, 11223]
        pointnum = nextpaduanum(N)
        degree = getdegree(pointnum)

        phasecoeffs = zeros(degree+1, degree+1, D)
        phasecoeffs[1, 1, 1] = 0.5  # 0.5
        phasecoeffs[1, 2, 1] = 0.5  # 0.5 x
        phasecoeffs[1, 1, D ÷ 2 + 1] = 0.5  # 0.5
        phasecoeffs[2, 1, D ÷ 2 + 1] = 0.5  # 0.5 y

        ∂xcoeffs, ∂ycoeffs = differentiate!(similar(phasecoeffs), similar(phasecoeffs),
            DiffPlan{Float64}(degree), phasecoeffs)

        phasepoints = invpaduatransform!(zeros(getpaduanum(degree), D),
            InvPaduaTransformPlan{Float64}(degree), phasecoeffs)

        Ω(v, t, p) = CanonicalSymplecticMatrix(D)

        plan = OOPIntPlan{Float64}(D, degree)

        intcoeffs = zeros(degree+1, degree+1)

        getintegrand!(intcoeffs, plan, Ω, phasepoints, 0, nothing, ∂xcoeffs, ∂ycoeffs)

        testcoeffs = zeros(degree+1, degree+1)
        testcoeffs[1, 1] = 0.25

        @test intcoeffs ≈ intcoeffs atol=1eps()
    end
end

@safetestset "compute! with ChebyshevPlan (OOP)" begin
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation.PaduaTransforms
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation:
        ChebyshevPlan, compute!, DiffPlan
    using PoincareInvariants.CanonicalSymplecticStructures

    @testset "$N Points in $D dimensions" for D in [4, 12], N in [20, 1234]
        Ω(v, t, p) = CanonicalSymplecticMatrix(D)
        plan = ChebyshevPlan{Float64}(Ω, D, N, Val(false))

        pointnum = nextpaduanum(N)
        degree = getdegree(pointnum)

        @test plan.degree == degree

        @test plan.paduaplan isa PaduaTransformPlan
        @test plan.paduaplan.degree == degree

        testphasecoeffs = zeros(degree+1, degree+1, D)
        @test plan.phasecoeffs == testphasecoeffs

        @test plan.diffplan isa DiffPlan
        @test size(plan.diffplan.D) == (degree+1, degree+1)

        test∂x = zeros(degree+1, degree+1, D)
        @test size(plan.∂x) == size(test∂x)

        test∂y = zeros(degree+1, degree+1, D)
        @test size(plan.∂y) == size(test∂y)

        testintcoeffs = zeros(degree+1, degree+1)
        @test size(plan.intcoeffs) == size(testintcoeffs)
        @test size(plan.integrator) == (degree+1,)

        phasepoints = getpaduapoints(degree) do x, y
            ntuple(D) do i
                i == 1 && return x
                i == D ÷ 2 + 1 && return y
                return 0
            end
        end

        @test compute!(plan, Ω, phasepoints, 0, nothing) ≈ 4 atol=10eps()

        testphasecoeffs[1, 2, 1] = 1  # 1 x
        testphasecoeffs[2, 1, D ÷ 2 + 1] = 1 # 1 y
        @test maximum(abs, plan.phasecoeffs .- testphasecoeffs) / eps() < 10

        test∂x[1, 1, 1] = 1
        @test maximum(abs, plan.∂x .- test∂x) / eps() < 50

        test∂y[1, 1, D ÷ 2 + 1] = 1
        @test maximum(abs, plan.∂y .- test∂y) / eps() < 50

        testintcoeffs[1, 1] = 1
        @test maximum(abs, plan.intcoeffs .- testintcoeffs) / eps() < 50
    end
end

@safetestset "getpoints and getpointnum" begin
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation.PaduaTransforms
    using PoincareInvariants.SecondPoincareInvariants.ChebyshevImplementation:
        ChebyshevPlan, getpoints, getpointnum

    Ω(v, t, p) = [0 -1; 1 0]
    plan = ChebyshevPlan{Float64}(Ω, 2, 11, Val(false))

    @test getpointnum(plan) == 15
    @test getpoints(plan) == map(getpaduapoints(4)) do v
        (v .+ 1) ./ 2
    end
end


text:
"""
```
get_scenario_input_file(m, scen::Scenario)
```

Get file name of raw scenario targets from `inpath(m, \"scenarios\")`.
"""
function get_scenario_input_file(m::AbstractModel, scen::Scenario)
    basename = string(scen.key) * "_" * scen.vintage * ".jld"
    return inpath(m, "scenarios", basename)
end

"""
```
n_scenario_draws(m, scen::Scenario)
```

Return the number of draws for `scen`, determined using
`get_scenario_input_file(m, scen)`.
"""
function n_scenario_draws(m::AbstractModel, scen::Scenario)
    input_file = get_scenario_input_file(m, scen)
    draws = h5open(input_file, "r") do file
        dataset = HDF5.o_open(file, "arr")
        size(dataset)[1]
    end
    return draws
end

"""
```
load_scenario_targets!(m, scen::Scenario, draw_index)
```

Add the targets from the `draw_index`th draw of the raw scenario targets to
`scen.targets`.
"""
function load_scenario_targets!(m::AbstractModel, scen::Scenario, draw_index::Int)
    path = get_scenario_input_file(m, scen)
    raw_targets = squeeze(h5read(path, "arr", (draw_index, :, :)), 1)
    target_inds = load(path, "target_indices")

    @assert collect(keys(target_inds)) == scen.target_names "Target indices in $path do not match target names in $(scen.key)"

    for (target_name, target_index) in target_inds
        scen.targets[target_name] = raw_targets[target_index, :]
    end

    return scen.targets
end

"""
```
get_scenario_filename(m, scen::AbstractScenario, output_var;
    pathfcn = rawpath, fileformat = :jld, directory = "")
```

Get scenario file name of the form
`pathfcn(m, \"scenarios\", output_var * filestring * string(fileformat))`. If
`directory` is provided (nonempty), then the same file name in that directory
will be returned instead.
"""
function get_scenario_filename(m::AbstractModel, scen::AbstractScenario, output_var::Symbol;
                               pathfcn::Function = rawpath,
                               fileformat::Symbol = :jld,
                               directory::String = "")
    filestring_addl = Vector{String}()
    if isa(scen, SingleScenario)
        push!(filestring_addl, "scen=" * string(scen.key))
    elseif isa(scen, ScenarioAggregate)
        push!(filestring_addl, "sagg=" * string(scen.key))
    end
    push!(filestring_addl, "svin=" * scen.vintage)

    base = string(output_var) * "." * string(fileformat)
    path = pathfcn(m, "scenarios", base, filestring_addl)
    if !isempty(directory)
        path = joinpath(directory, basename(path))
    end
    return path
end

"""
```
get_scenario_output_files(m, scen::SingleScenario, output_vars)
```

Return a `Dict{Symbol, String}` mapping `output_vars` to the raw simulated
scenario outputs for `scen`.
"""
function get_scenario_output_files(m::AbstractModel, scen::SingleScenario,
                                   output_vars::Vector{Symbol})
    output_files = Dict{Symbol, String}()
    for var in output_vars
        output_files[var] = get_scenario_filename(m, scen, var)
    end
    return output_files
end

"""
```
get_scenario_mb_input_file(m, scen::AbstractScenario, output_var)
```

Call `get_scenario_filename` while replacing `forecastut` and `forecast4q` in
`output_var` with `forecast`.
"""
function get_scenario_mb_input_file(m::AbstractModel, scen::AbstractScenario, output_var::Symbol)
    input_file = get_scenario_filename(m, scen, output_var)
    input_file = replace(input_file, "forecastut", "forecast")
    input_file = replace(input_file, "forecast4q", "forecast")
    return input_file
end

"""
```
get_scenario_mb_output_file(m, scen::AbstractScenario, output_var;
    directory = "")
```

Call `get_scenario_filename` while tacking on `\"mb\"` to the front of the base
file name.
"""
function get_scenario_mb_output_file(m::AbstractModel, scen::AbstractScenario, output_var::Symbol;
                                     directory::String = "")
    fullfile = get_scenario_filename(m, scen, output_var, pathfcn = workpath, directory = directory)
    joinpath(dirname(fullfile), "mb" * basename(fullfile))
end

"""
```
get_scenario_mb_metadata(m, scen::SingleScenario, output_var)

get_scenario_mb_metadata(m, agg::ScenarioAggregate, output_var)
```

Return the `MeansBands` metadata dictionary for `scen`.
"""
function get_scenario_mb_metadata(m::AbstractModel, scen::SingleScenario, output_var::Symbol)
    forecast_output_file = get_scenario_mb_input_file(m, scen, output_var)
    metadata = get_mb_metadata(m, :mode, :none, output_var, forecast_output_file)
    metadata[:scenario_key] = scen.key
    metadata[:scenario_vint] = scen.vintage

    return metadata
end

function get_scenario_mb_metadata(m::AbstractModel, agg::ScenarioAggregate, output_var::Symbol)
    forecast_output_file = get_scenario_mb_input_file(m, agg.scenarios[1], output_var)
    metadata = get_mb_metadata(m, :mode, :none, output_var, forecast_output_file)

    # Initialize start and end date
    start_date = date_forecast_start(m)
    end_date   = maximum(keys(metadata[:date_inds]))

    for scen in agg.scenarios
        forecast_output_file = get_scenario_mb_input_file(m, scen, output_var)
        scen_dates = jldopen(forecast_output_file, "r") do file
            read(file, "date_indices")
        end

        # Throw error if start date for this scenario doesn't match
        if map(reverse, scen_dates)[1] != start_date
            error("All scenarios in agg must start from the same date")
        end

        # Update end date if necessary
        end_date = max(end_date, maximum(keys(metadata[:date_inds])))
    end

    dates = quarter_range(start_date, end_date)
    metadata[:date_inds] = OrderedDict{Date, Int}(d => i for (i, d) in enumerate(dates))
    metadata[:scenario_key] = agg.key
    metadata[:scenario_vint] = agg.vintage

    return metadata
end

"""
```
read_scenario_output(m, scen::SingleScenario, class, product, var_name)

read_scenario_output(m, agg::ScenarioAggregate, class, product, var_name)
```

Given either `scen` or `agg`, read in and return all draws of and the
appropriate reverse transform for `var_name`.
"""
function read_scenario_output(m::AbstractModel, scen::SingleScenario, class::Symbol, product::Symbol,
                              var_name::Symbol)
    # Get filename
    filename = get_scenario_mb_input_file(m, scen, Symbol(product, class))

    jldopen(filename, "r") do file
        # Read forecast outputs
        fcast_series = read_forecast_series(file, class, product, var_name)

        # Parse transform
        class_long = get_class_longname(class)
        transforms = read(file, string(class_long) * "_revtransforms")
        transform = parse_transform(transforms[var_name])

        fcast_series, transform
    end
end

function read_scenario_output(m::AbstractModel, agg::ScenarioAggregate, class::Symbol,
                              product::Symbol, var_name::Symbol)
    # Aggregate scenarios
    nscens = length(agg.scenarios)
    agg_draws = Vector{Matrix{Float64}}(nscens)

    # If not sampling, initialize vector to record number of draws in each
    # scenario in order to update `agg.proportions` and `agg.total_draws` at the
    # end
    if !agg.sample
        n_scen_draws = zeros(Int, nscens)
    end

    # Initialize transform so it can be assigned from within the following for
    # loop. Each transform read in from read_scenario_output will be the
    # same. We just want to delegate the transform parsing to the recursive
    # read_scenario_output call.
    transform = identity

    for (i, scen) in enumerate(agg.scenarios)
        # Recursively read in scenario draws
        scen_draws, transform = read_scenario_output(m, scen, class, product, var_name)

        # Sample if desired
        agg_draws[i] = if agg.sample
            pct = agg.proportions[i]
            actual_ndraws = size(scen_draws, 1)
            desired_ndraws = convert(Int, round(pct * agg.total_draws))

            sampled_inds = if agg.replace
                sample(1:actual_ndraws, desired_ndraws, replace = true)
            else
                if desired_ndraws == 0
                    Int[]
                else
                    quotient  = convert(Int, floor(actual_ndraws / desired_ndraws))
                    remainder = actual_ndraws % desired_ndraws
                    vcat(repmat(1:actual_ndraws, quotient),
                         sample(1:actual_ndraws, remainder, replace = false))
                end
            end
            sort!(sampled_inds)
            scen_draws[sampled_inds, :]
        else
            # Record number of draws in this scenario
            n_scen_draws[i] = size(scen_draws, 1)
            scen_draws
        end
    end

    # Stack draws from all component scenarios
    fcast_series = cat(1, agg_draws...)

    # If not sampling, update `agg.proportions` and `agg.total_draws`
    if !agg.sample
        agg.total_draws = sum(n_scen_draws)
        agg.proportions = n_scen_draws ./ agg.total_draws
    end

    return fcast_series, transform
end

"""
```
read_scenario_mb(m, scen::AbstractScenario, output_var; directory = "")
```

Read in an alternative scenario `MeansBands` object.
"""
function read_scenario_mb(m::AbstractModel, scen::AbstractScenario, output_var::Symbol;
                          directory::String = "")
    filepath = get_scenario_mb_output_file(m, scen, output_var, directory = directory)
    read_mb(filepath)
end


text:
module FactorGraph

using TupleTools
using OMEinsum
using TensorOperations: optimaltree, Power
using TensorOperations

include("utils.jl")
include("contract.jl")
include("incidence_matrix.jl")
include("fgraph.jl")
include("random_graphs.jl")
include("Zygote_patch.jl")

end # module


text:
"""implementation of a standard prioritized planning

Ref:
- Silver, D. (2005). Cooperative Pathfinding. AIIDE.
- Erdmann, M., & Lozano-Perez, T. (1987). On multiple moving objects. Algorithmica.
"""
module LibPP
export PP

import ...MRMP: AbsState, Node, now, elapsed_sec
import ..Solvers:
    SearchNode,
    find_timed_path,
    convert_paths_to_configurations,
    get_distance_tables,
    gen_g_func,
    PRMs,
    PRMs!
import Printf: @sprintf
import Base: @kwdef
import Random: randperm

"""
    PP(
        config_init::Vector{State},
        config_goal::Vector{State},
        connect::Function,
        collide::Function,
        check_goal::Function;
        g_func::Function = gen_g_func(stay_penalty = 0.1),
        num_vertices::Int64 = 100,
        rad::Union{Nothing,Real} = nothing,
        rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(config_init)),
        roadmaps_growing_rate::Union{Nothing,Float64} = nothing,
        max_makespan::Int64 = 20,
        order_randomize::Bool = true,
        TIME_LIMIT::Union{Nothing,Real} = nothing,
        VERBOSE::Int64 = 0,
    )::Tuple{
        Union{Nothing,Vector{Vector{Node{State}}}},  # solution
        Vector{Vector{Node{State}}},  # roadmaps
    } where {State<:AbsState}

implementation of PP on PRM
"""
function PP(
    config_init::Vector{State},
    config_goal::Vector{State},
    connect::Function,
    collide::Function,
    check_goal::Function;
    g_func::Function = gen_g_func(stay_penalty = 0.1),
    num_vertices::Int64 = 100,
    rad::Union{Nothing,Real} = nothing,
    rads::Union{Vector{Nothing},Vector{Float64}} = fill(rad, length(config_init)),
    roadmaps_growing_rate::Union{Nothing,Float64} = nothing,
    max_makespan::Int64 = 20,
    order_randomize::Bool = true,
    TIME_LIMIT::Union{Nothing,Real} = nothing,
    VERBOSE::Int64 = 0,
)::Tuple{
    Union{Nothing,Vector{Vector{Node{State}}}},  # solution
    Vector{Vector{Node{State}}},  # roadmaps
} where {State<:AbsState}

    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    # get initial roadmaps
    roadmaps = PRMs(
        config_init,
        config_goal,
        connect,
        num_vertices;
        rads = rads,
        TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
    )
    VERBOSE > 0 && @info @sprintf(
        "\t%f sec, construct initial roadmaps: |V|=%d",
        elapsed(),
        num_vertices
    )

    if timeover()
        VERBOSE > 0 && @info @sprintf("\tsolution is not found within time limit")
        return (nothing, roadmaps)
    end

    # try to find solution
    solution, roadmaps = PP(
        roadmaps,
        collide,
        check_goal,
        g_func;
        max_makespan = max_makespan,
        order_randomize = order_randomize,
        VERBOSE = VERBOSE,
        TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
    )

    # fail to find solution -> increase roadmaps
    while isnothing(solution) && !isnothing(roadmaps_growing_rate) && !timeover()
        num_vertices = Int64(floor(num_vertices * roadmaps_growing_rate))
        VERBOSE > 0 &&
            @info @sprintf("\t%f sec, update roadmaps: |V|=%d", elapsed(), num_vertices)
        # update roadmap
        roadmaps = PRMs!(
            roadmaps,
            connect,
            num_vertices;
            rads = rads,
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        )
        timeover() && break
        # retry
        solution, roadmaps = PP(
            roadmaps,
            collide,
            check_goal,
            g_func;
            max_makespan = max_makespan,
            order_randomize = order_randomize,
            VERBOSE = VERBOSE,
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        )
    end

    if VERBOSE > 0
        if isnothing(solution)
            @info "\tfailed to find solution"
        else
            @info "\tfound solution"
        end
    end
    return (solution, roadmaps)
end

function PP(
    roadmaps::Vector{Vector{Node{State}}},
    collide::Function,
    check_goal::Function,
    g_func::Function;
    max_makespan::Int64 = 20,
    order_randomize::Bool = true,
    TIME_LIMIT::Union{Nothing,Real} = nothing,
    VERBOSE::Int64 = 0,
)::Tuple{
    Union{Nothing,Vector{Vector{Node{State}}}},  # solution
    Vector{Vector{Node{State}}},  # roadmaps
} where {State<:AbsState}

    # timeout
    t_s = now()
    elapsed() = elapsed_sec(t_s)
    timeover() = TIME_LIMIT != nothing && elapsed() > TIME_LIMIT

    N = length(roadmaps)
    distance_tables = get_distance_tables(roadmaps; g_func = g_func)
    paths = map(i -> Vector{Node{State}}(), 1:N)
    costs_table = Vector{Float64}([0.0])
    for (k, i) in enumerate(order_randomize ? randperm(N) : (1:N))
        timeover() && return (nothing, roadmaps)
        VERBOSE > 1 && @info @sprintf("\t\tstart planning for agent-%d, %d / %d", i, k, N)

        # setup functions
        invalid =
            (S_from::SearchNode{State}, S_to::SearchNode{State}) -> begin
                t = S_to.t

                # avoid infinite search
                if t > max_makespan
                    return true
                end

                # check collisions
                for j = 1:N
                    if j == i || isempty(paths[j])
                        continue
                    end
                    l = length(paths[j])
                    q_j_from = paths[j][min(t - 1, l)].q
                    q_j_to = paths[j][min(t, l)].q
                    if collide(S_from.v.q, S_to.v.q, q_j_from, q_j_to, i, j)
                        return true
                    end
                end

                return false
            end

        check_goal_i =
            (S::SearchNode{State}) -> begin
                if !check_goal(S.v, i)
                    return false
                end

                # check additional constraints
                for j = 1:N
                    if j == i || isempty(paths[j])
                        continue
                    end

                    for t = S.t+1:length(paths[j])
                        if collide(S.v.q, S.v.q, paths[j][t-1].q, paths[j][t].q, i, j)
                            return false
                        end
                    end
                end

                return true
            end
        h_func_i = (v::Node{State}) -> distance_tables[i][v.id]

        # accurate cost calculation
        g_func_i =
            (S::SearchNode{State}, q::State) -> begin
                cost_to_come_S =
                    max(S.g, isassigned(costs_table, S.t) ? costs_table[S.t] : 0)
                return max(
                    cost_to_come_S + g_func(S.v.q, q, i),
                    isassigned(costs_table, S.t + 1) ? costs_table[S.t+1] : 0,
                )
            end

        path = find_timed_path(
            roadmaps[i],
            invalid,
            check_goal_i,
            h_func_i,
            g_func_i;
            TIME_LIMIT = (isnothing(TIME_LIMIT) ? nothing : TIME_LIMIT - elapsed()),
        )

        # case failure
        if isnothing(path)
            VERBOSE > 1 && @info @sprintf("\t\tfailed planning for agent-%d", i)
            return (nothing, roadmaps)
        end

        paths[i] = path

        # update costs table
        for t = 2:length(path)
            c_new = max(
                isassigned(costs_table, t) ? costs_table[t] : 0,
                costs_table[t-1] + g_func(path[t-1].q, path[t].q),
            )
            if isassigned(costs_table, t)
                costs_table[t] = c_new
            else
                push!(costs_table, c_new)
            end
        end
    end

    return (convert_paths_to_configurations(paths), roadmaps)
end

end


text:
using QuasinormalModes
using RootsAndPoles

# ------------------------------------------------------------------
# 1. Setting up Schwarzschild black hole data structure for
#    computing the QNMs numerically.
# ------------------------------------------------------------------

struct NSchwarzschildData{N,T} <: NumericAIMProblem{N,T}
    nIter::N
    x0::T
    l::N
    s::N
end

function NSchwarzschildData(nIter::N, x0::T, l::N, s::N) where {N,T}
    return NSchwarzschildData{N,T}(nIter, x0, l, s)
end

QuasinormalModes.λ0(::NSchwarzschildData{N,T}) where {N,T} = (x,ω) -> (-1 + (2*im)*ω + x*(4 - 3*x + (4*im)*(-2 + x)*ω))/((-1 + x)^2*x)
QuasinormalModes.S0(d::NSchwarzschildData{N,T}) where {N,T} = (x,ω) -> (d.l + d.l^2 + (-1 + d.s^2)*(-1 + x) + (4*im)*(-1 + x)*ω + 4*(-2 + x)*ω^2)/((-1 + x)^2*x)

QuasinormalModes.get_niter(d::NSchwarzschildData{N,T}) where {N,T} = d.nIter
QuasinormalModes.get_x0(d::NSchwarzschildData{N,T}) where {N,T} = d.x0

# ------------------------------------------------------------------
# 2. Benchmark
# ------------------------------------------------------------------

"""
    bench_iter(iter_start, iter_end, x0, l, s, xt, ft, nls_iter, reference_ωr, reference_ωi, repeat)
Measures the time required to compute a pair of quasinormal frequencies as a function of the
number of AIM iterations performed. The result is also compared with a know result in order
to provide a error measurement.
# Input
- `iter_start`: First number of iterations to perform.
- `iter_end`: Last number of iterations.
- `x0`: The point around which AIM functions will be Taylor expanded.
- `l`: The l value of the sought mode.
- `s`: The s value of the sought mode.
- `xt`: The `xtol` value passed to NLSolve.
- `ft`: The `ftol` value passed to NLSolve.
- `nls_iter`: The number of iterations that NLSolve will perform.
- `reference_ωr`: The reference value for the real QNM frequency.
- `reference_ωi`: The reference value for the imaginary QNM frequency.
- `repeat`: The amount of times to repeat the operation.
# Output
    nothing
"""
function bench_iter(iter_start, iter_end, x0, l, s, xt, ft, nls_iter, reference_ωr, reference_ωi, repeat)

    for i in 1:repeat
        println("Benchmark iteration ", i, ":")
        file = open("run_$(i).dat", "w")
        println(file, "# 1:time(ns) 2:iter 3:x0 4:l 5:s 6:xt 7:ft 8:w_r 9:w_i 10:error in w_r 11:error in w_i")

        for iter in iter_start:iter_end
            println("    AIM iteration ", iter)
            p_num = NSchwarzschildData(iter, x0, l, s);
            c_num = AIMCache(p_num)

            t0 = time_ns();
            ev = computeEigenvalues(Threaded(), p_num, c_num, typeof(x0)(reference_ωr, reference_ωi), nls_xtol = xt, nls_ftol = ft, nls_iterations = nls_iter)
            elapsed_time = time_ns() - t0

            if(ev.x_converged || ev.f_converged)
                println(file,
                        elapsed_time,          "    ",
                        iter,                  "    ",
                        convert(Float64, x0),  "    ",
                        l,                     "    ",
                        s,                     "    ",
                        convert(Float64, xt),  "    ",
                        convert(Float64, ft),  "    ",
                        ev.zero[1],            "    ",
                        ev.zero[2],            "    ",
                        ev.zero[1] - reference_ωr, "    ",
                        ev.zero[2] - reference_ωi
                        )

                flush(file)
            end
        end
        close(file)
    end
end

# We call the function once and discard the results in order so that the compilation time does not get included in the benchmark.
bench_iter(0x00001, 0x00064, Complex(big"0.39", big"0.0"), 0x00000, 0x00000, big"1.0e-55", big"1.0e-55", 5000000, big"0.2209098781608393", big"-0.2097914341737619", 1)
bench_iter(0x00001, 0x00064, Complex(big"0.39", big"0.0"), 0x00000, 0x00000, big"1.0e-55", big"1.0e-55", 5000000, big"0.2209098781608393", big"-0.2097914341737619", 20)


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    isdefined(Comonicon, Symbol("##LeafCommand#11")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##LeafCommand#11")),
        String,
        Array{Comonicon.Types.Arg,1},
        Array{Comonicon.Types.Option,1},
        Array{Comonicon.Types.Flag,1},
        Comonicon.Types.CommandDoc,
        LineNumberNode,
        Type{Comonicon.Types.LeafCommand},
        typeof(identity),
    })
    isdefined(Comonicon, Symbol("##LeafCommand#11")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##LeafCommand#11")),
        String,
        Array{Comonicon.Types.Arg,1},
        Array{Comonicon.Types.Option,1},
        Array{Comonicon.Types.Flag,1},
        String,
        LineNumberNode,
        Type{Comonicon.Types.LeafCommand},
        typeof(identity),
    })
    isdefined(Comonicon, Symbol("##partition#12")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##partition#12")),
        Int64,
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Arg,
        String,
        String,
    })
    isdefined(Comonicon, Symbol("##partition#12")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##partition#12")),
        Int64,
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Flag,
        String,
        String,
    })
    isdefined(Comonicon, Symbol("##partition#12")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##partition#12")),
        Int64,
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.LeafCommand,
        String,
    })
    isdefined(Comonicon, Symbol("##partition#12")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##partition#12")),
        Int64,
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.NodeCommand,
        String,
    })
    isdefined(Comonicon, Symbol("##partition#12")) && precompile(Tuple{
        getfield(Comonicon.Types, Symbol("##partition#12")),
        Int64,
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Option,
        String,
        String,
    })
    isdefined(Comonicon, Symbol("##xcall#6")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#6")),
        Base.Iterators.Pairs{Symbol,QuoteNode,Tuple{Symbol},NamedTuple{(:line,),Tuple{QuoteNode}}},
        typeof(Comonicon.Parse.xcall),
        Module,
        typeof(identity),
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#7")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#7")),
        Base.Iterators.Pairs{Symbol,QuoteNode,Tuple{Symbol},NamedTuple{(:line,),Tuple{QuoteNode}}},
        typeof(Comonicon.Parse.xcall),
        Module,
        Symbol,
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#7")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#7")),
        Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},
        typeof(Comonicon.Parse.xcall),
        Module,
        Symbol,
        String,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#8")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#8")),
        Base.Iterators.Pairs{Symbol,QuoteNode,Tuple{Symbol},NamedTuple{(:line,),Tuple{QuoteNode}}},
        typeof(Comonicon.Parse.xcall),
        typeof(identity),
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#9")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#9")),
        Base.Iterators.Pairs{Symbol,QuoteNode,Tuple{Symbol},NamedTuple{(:line,),Tuple{QuoteNode}}},
        typeof(Comonicon.Parse.xcall),
        GlobalRef,
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#9")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#9")),
        Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},
        typeof(Comonicon.Parse.xcall),
        GlobalRef,
        GlobalRef,
        Int,
    })
    isdefined(Comonicon, Symbol("##xcall#9")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("##xcall#9")),
        Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},
        typeof(Comonicon.Parse.xcall),
        GlobalRef,
        String,
        Int,
    })
    isdefined(Comonicon, Symbol("#@cast")) &&
        precompile(Tuple{getfield(Comonicon.Parse, Symbol("#@cast")),LineNumberNode,Module,Int})
    isdefined(Comonicon, Symbol("#@main")) &&
        precompile(Tuple{getfield(Comonicon.Parse, Symbol("#@main")),LineNumberNode,Module,Int})
    isdefined(Comonicon, Symbol("#cmd_script##kw")) && precompile(Tuple{
        getfield(Comonicon, Symbol("#cmd_script##kw")),
        NamedTuple{
            (:exename, :project, :sysimg, :compile, :optimize),
            Tuple{String,String,Nothing,Nothing,Int64},
        },
        typeof(Comonicon.cmd_script),
        Module,
        String,
    })
    isdefined(Comonicon, Symbol("#command##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#command##kw")),
        NamedTuple{(:line,),Tuple{LineNumberNode}},
        typeof(Comonicon.Parse.command),
        typeof(identity),
        Array{Tuple{String,DataType,Bool},1},
        Array{Any,1},
    })
    isdefined(Comonicon, Symbol("#command##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#command##kw")),
        NamedTuple{(:line,),Tuple{LineNumberNode}},
        typeof(Comonicon.Parse.command),
        typeof(identity),
        Array{Tuple{String,DataType,Bool},1},
        Array{Tuple{String,DataType,Bool},1},
    })
    isdefined(Comonicon, Symbol("#install##kw")) && precompile(Tuple{
        getfield(Comonicon, Symbol("#install##kw")),
        NamedTuple{(:bin, :completion),Tuple{String,Bool}},
        typeof(Comonicon.install),
        Module,
    })
    isdefined(Comonicon, Symbol("#xcall##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#xcall##kw")),
        NamedTuple{(:line,),Tuple{QuoteNode}},
        typeof(Comonicon.Parse.xcall),
        GlobalRef,
        Symbol,
        Expr,
        Expr,
    })
    isdefined(Comonicon, Symbol("#xcall##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#xcall##kw")),
        NamedTuple{(:line,),Tuple{QuoteNode}},
        typeof(Comonicon.Parse.xcall),
        Module,
        Symbol,
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("#xcall##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#xcall##kw")),
        NamedTuple{(:line,),Tuple{QuoteNode}},
        typeof(Comonicon.Parse.xcall),
        Module,
        typeof(identity),
        Symbol,
        Int,
    })
    isdefined(Comonicon, Symbol("#xcall##kw")) && precompile(Tuple{
        getfield(Comonicon.Parse, Symbol("#xcall##kw")),
        NamedTuple{(:line,),Tuple{QuoteNode}},
        typeof(Comonicon.Parse.xcall),
        typeof(identity),
        Symbol,
        Expr,
        Expr,
    })
    precompile(Tuple{
        typeof(Comonicon.CodeGen.codegen),
        Comonicon.CodeGen.ASTCtx,
        Comonicon.Types.LeafCommand,
    })
    precompile(Tuple{
        typeof(Comonicon.CodeGen.codegen),
        Comonicon.CodeGen.ASTCtx,
        Comonicon.Types.NodeCommand,
    })
    precompile(Tuple{
        typeof(Comonicon.CodeGen.codegen),
        Comonicon.CodeGen.ZSHCompletionCtx,
        Comonicon.Types.EntryCommand,
    })
    precompile(Tuple{typeof(Comonicon.CodeGen.codegen),Comonicon.Types.EntryCommand})
    precompile(Tuple{typeof(Comonicon.CodeGen.prettify),Expr})
    precompile(Tuple{typeof(Comonicon.CodeGen.prettify),Int64})
    precompile(Tuple{typeof(Comonicon.CodeGen.prettify),Symbol})
    precompile(Tuple{typeof(Comonicon.CodeGen.rm_lineinfo),Int64})
    precompile(Tuple{typeof(Comonicon.CodeGen.rm_lineinfo),Symbol})
    precompile(Tuple{typeof(Comonicon.CodeGen.xparse_args),Type{Int},Expr})
    precompile(Tuple{typeof(Comonicon.CodeGen.xparse_args),Type{Int},Int64})
    precompile(Tuple{typeof(Comonicon.PATH.project),Module,String,String})
    precompile(Tuple{typeof(Comonicon.PATH.project),String,String,String})
    precompile(Tuple{typeof(Comonicon.PATH.project),String,String})
    precompile(Tuple{typeof(Comonicon.Parse.cast_m),Module,QuoteNode,Expr})
    precompile(Tuple{typeof(Comonicon.Parse.create_entry),Module,QuoteNode,Expr,Expr})
    precompile(Tuple{
        typeof(Comonicon.Parse.create_option),
        String,
        Type{Int},
        Base.Dict{String,Tuple{String,String,Bool}},
        LineNumberNode,
    })
    precompile(Tuple{typeof(Comonicon.Parse.disable_cache)})
    precompile(Tuple{typeof(Comonicon.Parse.main_m),Module,QuoteNode,Expr,Expr})
    precompile(Tuple{typeof(Comonicon.Parse.main_m),Module,QuoteNode,Expr})
    precompile(Tuple{typeof(Comonicon.Parse.read_doc),Markdown.MD})
    precompile(Tuple{
        typeof(Comonicon.Parse.set_cmd!),
        Base.Dict{String,Any},
        Comonicon.Types.EntryCommand,
        String,
    })
    precompile(Tuple{
        typeof(Comonicon.Parse.set_cmd!),
        Base.Dict{String,Any},
        Comonicon.Types.LeafCommand,
    })
    precompile(Tuple{typeof(Comonicon.Parse.to_argument),Base.Dict{Symbol,Any},Expr})
    precompile(Tuple{typeof(Comonicon.Parse.to_argument),Base.Dict{Symbol,Any},Symbol})
    precompile(Tuple{typeof(Comonicon.Parse.to_option_or_flag),Expr})
    precompile(Tuple{typeof(Comonicon.Parse.xcall),GlobalRef,String,QuoteNode,Int})
    precompile(Tuple{typeof(Comonicon.Parse.xcall),Module,Symbol,String,Int})
    precompile(Tuple{typeof(Comonicon.Types.cmd_doc),Comonicon.Types.LeafCommand})
    precompile(Tuple{typeof(Comonicon.Types.cmd_doc),Comonicon.Types.NodeCommand})
    precompile(Tuple{typeof(Comonicon.Types.cmd_name),Comonicon.Types.LeafCommand})
    precompile(Tuple{typeof(Comonicon.Types.cmd_name),Comonicon.Types.NodeCommand})
    precompile(Tuple{typeof(Comonicon.Types.ignore_type),Type{Int}})
    precompile(Tuple{
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Arg,
        String,
        String,
    })
    precompile(Tuple{
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Flag,
        String,
        String,
    })
    precompile(Tuple{
        typeof(Comonicon.Types.partition),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.Option,
        String,
        String,
    })
    precompile(Tuple{
        typeof(Comonicon.Types.print_cmd),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.EntryCommand,
    })
    precompile(Tuple{
        typeof(Comonicon.Types.print_cmd),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.LeafCommand,
    })
    precompile(Tuple{
        typeof(Comonicon.Types.print_cmd),
        Base.IOContext{Base.GenericIOBuffer{Array{UInt8,1}}},
        Comonicon.Types.NodeCommand,
    })
    precompile(Tuple{typeof(Comonicon.precompile_script),Module})
end


text:
function rotmat(u,θ)
    ux,uy,uz = u
    cos(θ)*eye(3) + sin(θ)*[[0 -uz uy];[uz 0 -ux];[-uy ux 0]] + (1-cos(θ))*
        [[ux^2 ux*uy ux*uz];[ux*uy uy^2 uy*uz];[ux*uz uy*uz uz^2]]
end

function randrotmat()
    u = rand(3)
    u ./= norm(u)
    θ = 2π*rand()
    rotmat(u,θ)
end

function rotate!(Ω::AngularGrid, R::Matrix)
    for i in 1:length(Ω)
        Ω[i] = R * Ω[i]
    end
end

max_stars_repo_name: johncsnyder/Element


text:
"""
    DistributionUtils

Parameter and probability distribution utils.
"""
module DistributionUtils

using Distributions
using JLD2
using EnsembleKalmanProcesses.ParameterDistributions

export construct_priors, deserialize_prior
export logmean_and_logstd, mean_and_std_from_ln
export flatten_config_dict
export flat_dict_keys_where, identity, above_eps

"""
    flatten_config_dict(param_dict::Dict{String, Vector{T}})

Given a dictionary of parameter names to parameter vectors of arbitrary length,
return a new dictionary that maps a unique parameter name to each element of the full
flattened vector of parameters.

Inputs:

 - param_dict :: Dictionary of parameter names to constraints.

Outputs:
 - u_names :: Vector{String} :: vector of parameter names
 - values :: Vector{Vector{T}} :: vector of single-valued vectors encapsulating parameter values.
"""
function flatten_config_dict(param_dict::Dict{String, Vector{T}}) where {T}

    u_names = Vector{String}()
    values = Vector{Vector{T}}()
    for (param, value) in param_dict
        if length(value) > 1
            for j in 1:length(value)
                push!(u_names, "$(param)_{$j}")
                push!(values, [value[j]])
            end
        else
            push!(u_names, param)
            push!(values, value)
        end
    end
    return (u_names, values)
end

"Boolean specifying whether the input value is significant to machine precision"
above_eps(x::FT) where {FT <: Real} = x > eps(FT)

"Identity function"
identity(x) = x

"""
    flat_dict_keys_where(dict::Dict{String, Vector{T}}, condition::Function = identity) where {T}

Flattens the values of a dictionary with parameter vectors as keys, and returns the indices of
entries in the flattened dictionary satisfying a given condition.

Inputs:
    param_dict :: Dictionary of parameter names to vectors.
    condition :: A condition function operating on each dictionary value.
Outputs:
    Indices of flattened entries satisfying the `condition`.
"""
function flat_dict_keys_where(param_dict::Dict{String, Vector{T}}, condition::Function = identity) where {T}
    # flatten global dict
    if any(1 .< [length(val) for val in collect(values(param_dict))])
        _, dict_values = flatten_config_dict(param_dict)
    else
        dict_values = collect(values(param_dict))
    end
    dict_values = vcat(dict_values...)

    return findall(x -> condition(x), dict_values)
end

"""
    construct_priors(
        params::Dict{String, Vector{Constraint}};
        unconstrained_σ::Float64 = 1.0,
        prior_mean::Union{Vector{Float64}, Nothing} = nothing,
        outdir_path::String = pwd(),
        to_file::Bool = true,
    )

Define a prior Gaussian ParameterDistribution in unconstrained space
from a dictionary of constraints.

This constructor assumes independent priors and the same unconstrained
standard deviation for each parameter. Note that the standard deviation
in unconstrained space is normalized with respect to the constrained
interval width, so it automatically takes into account parameter scales.

The constructor also allows passing a prior mean for each parameter in
constrained space.

Inputs:
 - params :: Dictionary of parameter names to constraints.
 - unconstrained_σ :: Standard deviation of the transformed gaussians (unconstrained space).
 - prior_mean :: The mean value of the prior in constrained space. If not given,
    the prior is selected to be 0 in unconstrained space.
 - outdir_path :: Output path.
 - to_file :: Whether to write the serialized prior to a JLD2 file.

Output:
 - The prior ParameterDistribution.
"""
function construct_priors(
    params::Dict{String, Vector{Constraint}};
    unconstrained_σ::Float64 = 1.0,
    prior_mean::Union{Dict{String, Vector{Float64}}, Nothing} = nothing,
    outdir_path::String = pwd(),
    to_file::Bool = true,
)
    # if parameter vectors found => flatten
    if any(1 .< [length(val) for val in collect(values(params))])
        u_names, constraints = flatten_config_dict(params)
    else
        u_names = collect(keys(params))
        constraints = collect(values(params))
    end

    n_param = length(u_names)

    # All vars are approximated as Gaussian in unconstrained space.
    if isnothing(prior_mean)
        distributions = repeat([Parameterized(Normal(0.0, unconstrained_σ))], n_param)
    else
        if any(1 .< [length(val) for val in collect(values(prior_mean))])
            u_names_mean, prior_mean = flatten_config_dict(prior_mean)
            @assert u_names_mean == u_names
        else
            prior_mean = collect(values(prior_mean))
        end
        @assert length(prior_mean) == n_param
        uncons_prior_mean =
            [c[1].constrained_to_unconstrained(μ_cons[1]) for (μ_cons, c) in zip(prior_mean, constraints)]
        distributions = [Parameterized(Normal(uncons_μ[1], unconstrained_σ)) for uncons_μ in uncons_prior_mean]
    end
    to_file ? jldsave(joinpath(outdir_path, "prior.jld2"); distributions, constraints, u_names) : nothing
    marginal_priors = ParameterDistribution.(distributions, constraints, u_names)
    return combine_distributions(marginal_priors)
end

"""
    deserialize_prior(prior_dict::Dict{String, Any})

Generates a prior ParameterDistribution from arguments stored
in a dictionary.
"""
function deserialize_prior(prior_dict::Dict{String, T}) where {T}
    marginal_priors =
        ParameterDistribution.(prior_dict["distributions"], prior_dict["constraints"], prior_dict["u_names"])
    return combine_distributions(marginal_priors)
end

"""
    logmean_and_logstd(μ, σ)

Returns the lognormal parameters μ and σ from the mean μ and std σ of the
lognormal distribution.
"""
function logmean_and_logstd(μ, σ)
    σ_log = sqrt(log(1.0 + σ^2 / μ^2))
    μ_log = log(μ / (sqrt(1.0 + σ^2 / μ^2)))
    return μ_log, σ_log
end

"""
    mean_and_std_from_ln(μ, σ)

Returns the mean and variance of the lognormal distribution
from the lognormal parameters μ and σ.
"""
function mean_and_std_from_ln(μ_log, σ_log)
    μ = exp(μ_log + σ_log^2 / 2)
    σ = sqrt((exp(σ_log^2) - 1) * exp(2 * μ_log + σ_log^2))
    return μ, σ
end

end # module


text:
export AbstractMoistureModel, DryModel, EquilMoist, NonEquilMoist

#### Moisture component in atmosphere model
abstract type AbstractMoistureModel end

vars_state(::AbstractMoistureModel, ::AbstractStateType, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::AbstractMoistureModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::AbstractMoistureModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

function compute_gradient_argument!(
    ::AbstractMoistureModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

internal_energy(atmos::AtmosModel, state::Vars, aux::Vars) =
    internal_energy(atmos, atmos.orientation, state, aux)

@inline function internal_energy(
    atmos::AtmosModel,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    Thermodynamics.internal_energy(
        density(atmos, state, aux),
        state.energy.ρe,
        state.ρu,
        gravitational_potential(orientation, aux),
    )
end

"""
    DryModel

Assumes the moisture components is in the dry limit.
"""
struct DryModel <: AbstractMoistureModel end

vars_state(::DryModel, ::Auxiliary, FT) = @vars(θ_v::FT, air_T::FT)
@inline function atmos_nodal_update_auxiliary_state!(
    moist::DryModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.air_T = air_temperature(ts)
    nothing
end

"""
    EquilMoist

Assumes the moisture components are computed via thermodynamic equilibrium.
"""
Base.@kwdef struct EquilMoist{FT, IT} <: AbstractMoistureModel
    maxiter::IT = nothing
    tolerance::FT = nothing
end

vars_state(::EquilMoist, ::Prognostic, FT) = @vars(ρq_tot::FT)
vars_state(::EquilMoist, ::Primitive, FT) = @vars(q_tot::FT)
vars_state(::EquilMoist, ::Gradient, FT) = @vars(q_tot::FT)
vars_state(::EquilMoist, ::GradientFlux, FT) = @vars(∇q_tot::SVector{3, FT})
vars_state(::EquilMoist, ::Auxiliary, FT) =
    @vars(temperature::FT, θ_v::FT, q_liq::FT, q_ice::FT)

@inline function atmos_nodal_update_auxiliary_state!(
    moist::EquilMoist,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.temperature = air_temperature(ts)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.q_liq = PhasePartition(ts).liq
    aux.moisture.q_ice = PhasePartition(ts).ice
    nothing
end

function compute_gradient_argument!(
    moist::EquilMoist,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
end

function compute_gradient_flux!(
    moist::EquilMoist,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    # diffusive flux of q_tot
    diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
end

"""
    NonEquilMoist

Does not assume that the moisture components are in equilibrium.
"""
struct NonEquilMoist <: AbstractMoistureModel end

vars_state(::NonEquilMoist, ::Prognostic, FT) =
    @vars(ρq_tot::FT, ρq_liq::FT, ρq_ice::FT)
vars_state(::NonEquilMoist, ::Primitive, FT) =
    @vars(q_tot::FT, q_liq::FT, q_ice::FT)
vars_state(::NonEquilMoist, ::Gradient, FT) =
    @vars(q_tot::FT, q_liq::FT, q_ice::FT)
vars_state(::NonEquilMoist, ::GradientFlux, FT) = @vars(
    ∇q_tot::SVector{3, FT},
    ∇q_liq::SVector{3, FT},
    ∇q_ice::SVector{3, FT}
)
vars_state(::NonEquilMoist, ::Auxiliary, FT) = @vars(temperature::FT, θ_v::FT)

@inline function atmos_nodal_update_auxiliary_state!(
    moist::NonEquilMoist,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ts = new_thermo_state(atmos, state, aux)
    aux.moisture.temperature = air_temperature(ts)
    aux.moisture.θ_v = virtual_pottemp(ts)
    nothing
end

function compute_gradient_argument!(
    moist::NonEquilMoist,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
    transform.moisture.q_liq = state.moisture.ρq_liq * ρinv
    transform.moisture.q_ice = state.moisture.ρq_ice * ρinv
end

function compute_gradient_flux!(
    moist::NonEquilMoist,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    # diffusive fluxes of moisture variables
    diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
    diffusive.moisture.∇q_liq = ∇transform.moisture.q_liq
    diffusive.moisture.∇q_ice = ∇transform.moisture.q_ice
end


text:
immutable Epanechnikov{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ::T

    Epanechnikov(μ::T, σ::T) = (@check_args(Epanechnikov, σ > zero(σ)); new(μ, σ))
end

Epanechnikov{T<:Real}(μ::T, σ::T) = Epanechnikov{T}(μ, σ)
Epanechnikov(μ::Real, σ::Real) = Epanechnikov(promote(μ, σ)...)
Epanechnikov(μ::Integer, σ::Integer) = Epanechnikov(Float64(μ), Float64(σ))
Epanechnikov(μ::Real) = Epanechnikov(μ, 1.0)
Epanechnikov() = Epanechnikov(0.0, 1.0)

@distr_support Epanechnikov d.μ - d.σ d.μ + d.σ

#### Conversions
function convert{T<:Real}(::Type{Epanechnikov{T}}, μ::Real, σ::Real)
    Epanechnikov(T(μ), T(σ))
end
function convert{T <: Real, S <: Real}(::Type{Epanechnikov{T}}, d::Epanechnikov{S})
    Epanechnikov(T(d.μ), T(d.σ))
end

## Parameters

location(d::Epanechnikov) = d.μ
scale(d::Epanechnikov) = d.σ
params(d::Epanechnikov) = (d.μ, d.σ)
@inline partype{T<:Real}(d::Epanechnikov{T}) = T

## Properties
mean(d::Epanechnikov) = d.μ
median(d::Epanechnikov) = d.μ
mode(d::Epanechnikov) = d.μ

var(d::Epanechnikov) = d.σ^2 / 5
skewness{T<:Real}(d::Epanechnikov{T}) = zero(T)
kurtosis{T<:Real}(d::Epanechnikov{T}) = -2.914285714285714*one(T)  # 3/35-3

## Functions
function pdf{T<:Real}(d::Epanechnikov{T}, x::Real)
    u = abs(x - d.μ) / d.σ
    u >= 1 ? zero(T) : 3 * (1 - u^2) / (4 * d.σ)
end

function cdf{T<:Real}(d::Epanechnikov{T}, x::Real)
    u = (x - d.μ) / d.σ
    u <= -1 ? one(T) :
    u >= 1 ? zero(T) :
    1//2 + u * (3//4 - u^2/4)
end

function ccdf{T<:Real}(d::Epanechnikov{T}, x::Real)
    u = (d.μ - x) / d.σ
    u <= -1 ? one(T) :
    u >= 1 ? zero(T) :
    1//2 + u * (3//4 - u^2/4)
end

@quantile_newton Epanechnikov

function mgf{T<:Real}(d::Epanechnikov{T}, t::Real)
    a = d.σ * t
    a == 0 ? one(T) :
    3exp(d.μ * t) * (cosh(a) - sinh(a) / a) / a^2
end

function cf{T<:Real}(d::Epanechnikov{T}, t::Real)
    a = d.σ * t
    a == 0 ? one(T)+zero(T)*im :
    -3exp(im * d.μ * t) * (cos(a) - sin(a) / a) / a^2
end


text:
import Pkg

function get_my_version()
    version_number = try
    projecttoml_filename = joinpath(dirname(dirname(@__FILE__)),"Project.toml")
    projecttoml_parsed = Pkg.TOML.parse(read(projecttoml_filename, String))
    VersionNumber(projecttoml_parsed["version"])
    catch e
        @warn(string("Ignoring error: ", e))
        VersionNumber(0)
    end
    return version_number
end

macro async_errhandle(ex)
    ex_string = string(ex)
    quote
        t = @async try
            $(esc(ex))
        catch err
            bt = catch_backtrace()
            @warn "@async error from @async $($ex_string)"
            showerror(stderr, err, bt)
        end
    end
end


text:
function residual(model::Model, mass_matrix, dynamics_bias, kinematics, kinematics_jacobians, z, θ, μ)
    # dimensions
    nq = model.nq
    nu = model.nu
    nw = model.nw
    nc = model.nc
    nf = length(friction_coefficients(model))

    # unpack data
    q0 = θ[1:nq]
    q1 = θ[nq .+ (1:nq)]
    u1 = θ[2nq .+ (1:nu)]
    w1 = θ[2nq + nu .+ (1:nw)]
    friction = θ[2nq + nu + nw .+ (1:nf)]
    h = θ[2nq + nu + nw + nf .+ (1:1)]

    # unpack variables
    q2 = z[1:nq]
    γ1 = z[nq .+ (1:nc)]
    sγ1 = z[nq + nc .+ (1:nc)]
    ψ1 = z[nq + nc + nc .+ (1:nc)]
    b1 = z[nq + nc + nc + nc .+ (1:nc)]
    sψ1 = z[nq + nc + nc + nc + nc .+ (1:nc)]
    sb1 = z[nq + nc + nc + nc + nc + nc .+ (1:nc)]

    # rotation matrix at surface
    rotation = [vector_rotation_matrix([0.0; 1.0]) for i = 1:nc] #TODO non-flat surfaces

    # tangential velocity
    v1 = (q2 - q1) / h[1]
    vT = [(rotation[i] * kinematics_jacobians[i](q2) * v1)[1] for i = 1:nc]

    # contact forces
    J = contact_jacobian(model, q1)
    λ1 = transpose(J) * vcat([transpose(rotation[i]) * [b1[i]; γ1[i]] for i = 1:nc]...)

    # residual
    [
    dynamics(model, mass_matrix, dynamics_bias, h, q0, q1, u1, w1, λ1, q2);
    signed_distance(model, q2) - sγ1;
    vcat([vT[i] - sb1[i] for i = 1:nc]...);
    vcat([friction[i] * γ1[i] - ψ1[i] for i = 1:nc]...);
    γ1 .* sγ1 .- μ[1];
    vcat([cone_product([ψ1[i]; b1[i]], [sψ1[i]; sb1[i]]) - [μ[1]; 0.0] for i = 1:nc]...);
    ]
end


text:
using DelimitedFiles
using Statistics
using PyPlot
using PyCall
@pyimport matplotlib.patches as patches

# define my colors -
negative_color = (1/255)*[255,51,51]
positive_color = (1/255)*[51,153,255]

function calculate_x_axis_tick_positions(data_array;box_width=1.0)

    # what is the size?
    (number_of_rows, number_of_cols) = size(data_array)

    # initialize -
    x_axis_points = Float64[]
    epsilon = 0.2;

    for col_index = 1:number_of_cols

        # compute origin point -
        origin_point = (col_index - 1)+(col_index - 1)*epsilon + 1
        push!(x_axis_points,origin_point+0.5*box_width)
    end

    return x_axis_points
end

function calculate_y_axis_tick_positions(data_array; box_height=1.0)

    # what is the size?
    (number_of_rows, number_of_cols) = size(data_array)

    # initialize -
    y_axis_points = Float64[]

    for col_index = 1:number_of_cols

        # get the data array -
        data_scaled = data_array[:,col_index]

        # how many patches per col?
        number_of_patches = length(data_scaled)
        epsilon = 0.2;
        for row_index = 1:number_of_patches

            # compute origin point -
            origin_point = [(col_index - 1)+(col_index - 1)*epsilon + 1,(row_index - 1)+(row_index - 1)*epsilon+1];

            push!(y_axis_points,origin_point[2]+0.5*box_height)
        end
    end

    return y_axis_points
end

function visualize(data_array,x_axis_ticks,y_axis_ticks,y_label_array)

    # what is the size?
    (number_of_rows, number_of_cols) = size(data_array)

    epsilon = 0.2;
    ax = gca()

    # set the axis tick postions -
    #ax.set_xticks(x_axis_points)
    ax.set_yticks(y_axis_ticks)
    ax.set_yticklabels(y_label_array, fontsize=10,fontname="Arial")
    subplots_adjust(left=0.25)

    # create x-label array -
    x_label_array = ["S$(x)" for x=1:number_of_cols]
    #ax.set_xticklabels(x_label_array, rotation=45,fontsize=8)
    ax.spines["right"].set_visible(false)
    ax.spines["top"].set_visible(false)
    ax.spines["bottom"].set_visible(false)
    #ax.spines["left"].set_visible(false)

    ax.get_xaxis().set_visible(false)

    for col_index = 1:number_of_cols

        # get the data array -
        data_scaled = data_array[:,col_index]

        # how many patches per col?
        number_of_patches = length(data_scaled)
        for row_index = 1:number_of_patches

            # compute origin point -
            origin_point = [(col_index - 1)+(col_index - 1)*epsilon + 1,(row_index - 1)+(row_index - 1)*epsilon+1];

            # compute alpha -
            data_value = data_array[row_index, col_index]
            if (data_value == 0.0)

                # compute a new color -
                patch_face_color = "#ffffff"

                # draw the square -
                ax.add_patch(

                           patches.Rectangle(origin_point,   # (x,y)
                               1.0,          # width
                               1.0,          # height
                               facecolor=patch_face_color,
                               edgecolor="grey",
                               linewidth=0.5,
                           )
                       )

            else

                #beta = (data_value - a)/(b - a)

                # compute a new color -
                #patch_face_color = (1-beta)*negative_color .+ beta*positive_color
                patch_face_color = "#ffffff"
                if (data_value == 1.0)
                    patch_face_color="#DDE3ED"
                elseif (data_value == 2.0)
                    patch_face_color="#C8D1E0"
                elseif (data_value == 3.0)
                    patch_face_color="#AFBACC"
                elseif (data_value == 4.0)
                    patch_face_color="#8E99AB"
                elseif (data_value == 5.0)
                    patch_face_color="#707A8A"
                elseif (data_value == 6.0)
                    patch_face_color="#58606E"
                elseif (data_value == 7.0)
                    patch_face_color="#434A54"
                elseif (data_value == 8.0)
                    patch_face_color="#333840"
                end

                # draw the square -
                ax.add_patch(

                           patches.Rectangle(origin_point,   # (x,y)
                               1.0,          # width
                               1.0,          # height
                               facecolor=patch_face_color,
                               edgecolor="grey",
                               linewidth=0.5,
                               alpha = 1.0
                           )
                       )
            end
        end
    end

    axis("square")
end

function reduce(path_to_sensitivity_array::String; epsilon::Float64=1e-2)

    # parameter name array =
    parameter_name_array = [
		"dG_deGFP_RNAP"	            ;	# 1
		"dG_deGFP_σ70"	            ;	# 2
		"n_deGFP_σ70"	            ;	# 3
		"K_deGFP_σ70"	            ;	# 4
        "tc mod mRNA_deGFP"         ;   # 5
        "tc mod deGFP"              ;   # 6
        "stability mRNA_deGFP"      ;   # 7
        "stability deGFP"           ;   # 8
        "stability σ70"                 ;   # 9
        "half-life translation"         ;   # 10
        "saturation translation"        ;   # 11
    ];

    # load -
    sensitivity_array = readdlm(path_to_sensitivity_array)

    # parameters to kepp -
    index_keep_p = Int64[]

    # the parameters are on the rows - are any values abpove a threshold -
    (NR,NC) = size(sensitivity_array)
    for row_index = 1:NR

        # get the row data -
        row_data = sensitivity_array[row_index, :]

        # check -
        idx_check = findall(x->x>=epsilon, row_data)
        if (isempty(idx_check) == false)
            push!(index_keep_p, row_index)
        end
    end

    # pull out rows -
    reduced_array = sensitivity_array[index_keep_p, :]

    # last thing - lets put things in groups of 0,1,2,3 (or none,low,medium,high)
    (NR,NC) = size(reduced_array)
    category_array = zeros(NR,NC)
    for row_index = 1:NR
        for col_index = 1:NC

            old_value = reduced_array[row_index, col_index]
            order_of_magnitude = floor(Int,log10(old_value))
            oom_eps = floor(Int,log10(epsilon))

            # rules -
            if (order_of_magnitude<-4)                                  # below -4
                category_array[row_index, col_index] = 0
            elseif (order_of_magnitude == -4)                           # O(-3)
                category_array[row_index, col_index] = 1
            elseif (order_of_magnitude == -3)                           # O(-2)
                category_array[row_index, col_index] = 2
            elseif (order_of_magnitude == -2)                           # O(-1)
                category_array[row_index, col_index] = 3
            elseif (order_of_magnitude == -1 )                           # O(0)
                category_array[row_index, col_index] = 4
            elseif (order_of_magnitude == 0)                            # O(1)
                category_array[row_index, col_index] = 5
            elseif (order_of_magnitude == 1)                            # O(2)
                category_array[row_index, col_index] = 6
            elseif (order_of_magnitude == 2)                            # O(3)
                category_array[row_index, col_index] = 7
            else
                category_array[row_index, col_index] = 8
            end
        end
    end

    # grab p labels -
    p_label_array = parameter_name_array[index_keep_p]

    # return -
    return (index_keep_p, p_label_array, category_array)
end

# setup path to sensitvity array -
path_to_sensitivity_array = "./Sensitivity-abs-MVP-T20.dat"

# execute -
(idx_p, p_label_array, rsa) = reduce(path_to_sensitivity_array; epsilon=1e-3)

# calculate the tick locations -
x_axis_points = calculate_x_axis_tick_positions(rsa)
y_axis_points = calculate_y_axis_tick_positions(rsa)

# call plotting code -
visualize(rsa,x_axis_points,y_axis_points, p_label_array)

max_stars_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_issues_repo_name: varnerlab/Biophysical-TXTL-Model-Code

max_forks_repo_name: varnerlab/Biophysical-TXTL-Model-Code


text:
# Truly random stuff going on here

function comp_choosing_card()
  _Cards ::Array = []
  @async GenCard(4, _Cards)
  thinking("Computer is thinking", "Computer chose its card. It's showtime.")
  compCardStats ::Dict = _Cards[rand(1:end)]
  return compCardStats
end

# Should probably make this better and more AI like
function comp_pick_move()
  moveNumber = rand(1:3)
  moveChoice = ""

  if moveNumber == 1
    moveChoice = "attack"
    return moveChoice
  elseif moveNumber == 2
    moveChoice = "dodge"
    return moveChoice
  else
    moveChoice = "defend"
    return moveChoice
  end
end

max_stars_repo_name: NarutoXY/austroporean-cards

max_issues_repo_name: NarutoXY/austroporean-cards

max_forks_repo_name: NarutoXY/austroporean-cards