text:
cdir = pwd(); ii = findlast("MDP", cdir); MDPpath = cdir[1:ii[end]] * "/";
include(MDPpath * "Installation/setup.jl");

# Read xyz file to obtain configurations
filename = "/Users/ngoccuongnguyen/Dropbox (MIT)/MDP/data/LJCluster/curated_lj_cluster.xyz";
species = ["Ar"];
config = readEXTXYZ(filename, species);

# shift atom positions to put them in the box [0,10]x[0,10]x[0,10]
config.x = config.x .+ 5.0

# cumalative sum of numbers of atoms
natom = [0; cumsum(config.natom[:])];

ci = 1; # choose a configuration to start MD simulation
x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
t = config.t[:,(natom[ci]+1):natom[ci+1]] # atom types of configuration ci
q = [];                                   # atom charges of configuration ci
v = 0.0*x # zero velocity

# 3 vectors define the simulation box
a = [10.0; 0.0; 0.0]
b = [0.0; 10.0; 0.0]
c = [0.0; 0.0; 10.0]

# periodic boundary conditions
pbc = [0; 0; 0]

# LJ parameters
epsilon = 1.0;
sigma = 1.0;
A = 2*epsilon*sigma^12;
B = 2*epsilon*sigma^6;
rcut = 3.2934420858948856;

# define LJ potential
potentials[1] = addpotential("pair", "nonbonded", [0], "ljpotential", [A; B], rcut);
include(potentials[1].potentialfunction * ".jl");  # include the potentials file

app.unitstyle = "lj"; # unit system
app.ensemblemode = "nve";
app.time = 0.0;             # initial time
app.dt = 0.5e-3;            # time step size
app.ntimesteps = 1000;       # time steps
app.globalfreq = 10;        # time frequency to print out global outputs
app.neighskin = 0.5;        # add a skin distance for neighborlist
app.atommasses = [1.0 1.0];

# run MD simulation using velocity Verlet algorithm
VelocityVerlet(x, v, q, t, a, b, c, pbc, app, potentials)

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# ---
# title: 1307. Verbal Arithmetic Puzzle
# id: problem1307
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Backtracking
# link: <https://leetcode.com/problems/verbal-arithmetic-puzzle/description/>
# hidden: true
# ---
#
# Given an equation, represented by `words` on left side and the `result` on
# right side.
#
# You need to check if the equation is solvable under the following rules:
#
#   * Each character is decoded as one digit (0 - 9).
#   * Every pair of different characters they must map to different digits.
#   * Each `words[i]` and `result` are decoded as one number **without** leading zeros.
#   * Sum of numbers on left side (`words`) will equal to the number on right side (`result`).
#
# Return `True` if the equation is solvable otherwise return `False`.
#
#
#
# **Example 1:**
#
#
#
#     Input: words = ["SEND","MORE"], result = "MONEY"
#     Output: true
#     Explanation: Map 'S'-> 9, 'E'->5, 'N'->6, 'D'->7, 'M'->1, 'O'->0, 'R'->8, 'Y'->'2'
#     Such that: "SEND" + "MORE" = "MONEY" ,  9567 + 1085 = 10652
#
# **Example 2:**
#
#
#
#     Input: words = ["SIX","SEVEN","SEVEN"], result = "TWENTY"
#     Output: true
#     Explanation: Map 'S'-> 6, 'I'->5, 'X'->0, 'E'->8, 'V'->7, 'N'->2, 'T'->1, 'W'->'3', 'Y'->4
#     Such that: "SIX" + "SEVEN" + "SEVEN" = "TWENTY" ,  650 + 68782 + 68782 = 138214
#
# **Example 3:**
#
#
#
#     Input: words = ["THIS","IS","TOO"], result = "FUNNY"
#     Output: true
#
#
# **Example 4:**
#
#
#
#     Input: words = ["LEET","CODE"], result = "POINT"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `2 <= words.length <= 5`
#   * `1 <= words[i].length, result.length <= 7`
#   * `words[i], result` contains only upper case English letters.
#   * Number of different characters used on the expression is at most 10.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end

ext: jl


text:
#=
P37 (**) Calculate Euler's totient function phi(m) (improved).
See problem P34 for the definition of Euler's totient function. If the list of the prime factors of a number m is known in the form of problem P36 then the function phi(m) can be efficiently calculated as follows: Let ((p1 m1) (p2 m2) (p3 m3) ...) be the list of prime factors (and their multiplicities) of a given number m. Then phi(m) can be calculated with the following formula:
phi(m) = (p1 - 1) * p1 ** (m1 - 1) + (p2 - 1) * p2 ** (m2 - 1) + (p3 - 1) * p3 ** (m3 - 1) + ...

Note that a ** b stands for the b'th power of a.
=#

function phi(prime_factors)
end

pf = Vector{}[[3, 2], [5, 1], [7, 1]]


text:
"""
Dimless collision between bounding boxes. NOTE: if they don't collide
then object don't.
"""
function bbox_collides{D}(b::Box{D}, B::Box{D})
    !any((b.Y[i] + EPS < B.X[i] || B.Y[i] < b.X[i] - EPS) for i in 1:D)
end

# Generate bbox collisions for other objects
for typeb in (:Segment, :Box)
    for typeB in (:Segment, :Box)
        if !(typeb == :Box && typeB == :Box)

            @eval begin
                function bbox_collides{D}(b::$(typeb){D}, B::$(typeB){D})
                    bbox_collides(b.box, B.box)
                end
            end
        end
    end
end

# Collision of point and segment is dimless
"""Is point contained in the segment"""
function in{D}(p::Point{D}, line::Segment{D})
    t = dot(p - line.O, line.t)/line.length
    -EPS < t < 1 + EPS
end

"""Is point contained in the box"""
function in{D}(p::Point{D}, box::Box{D})
    all(box.X[i] - EPS < p[i] < box.Y[i] + EPS for i in 1:D)
end

max_stars_repo_name: MiroK/fractional-dim


text:
using Test, MLToolkit
using Distributions: Poisson

@testset "Russian Roulette" begin
    function X(i)
        if i < 0
            return 0
        end
        if i > 10
            return 0
        end
        return i <= 5 ? i : 10 - i
    end

    S_true = sum(X.(1:10))

    dist = Poisson(5.0)

    n_mc = 100_000
    S_rr = mean([roll(X, dist) for _ = 1:n_mc])

    @test S_rr ≈ S_true atol=0.1
    @test roll(X, dist, n_mc) ≈ S_true rtol=0.01

    dist = Poisson(5.0)
    T(i) = pdf(dist, i) * X(i)

    S_true = sum(X.(1:10) .* pdf.(Ref(dist), 1:10))

    n_mc = 100_000
    S_rr_1 = mean([roll(T, dist) for _ = 1:n_mc])
    S_rr_2 = mean([roll_expectation(X, dist, Poisson(5.1)) for _ = 1:n_mc])
    S_rr_3 = mean([roll_expectation(X, dist, LogitNPD()) for _ = 1:n_mc])

    @test S_rr_1 ≈ S_true atol=0.1
    @test S_rr_2 ≈ S_true atol=0.1
    @test S_rr_3 ≈ S_true atol=0.1

    dist = LogitNPD()
    T(i) = pdf(dist, i) * X(i)

    S_true = sum(X.(1:10) .* pdf.(Ref(dist), 1:10))

    n_mc = 100_000
    S_rr_1 = mean([roll(T, dist) for _ = 1:n_mc])
    S_rr_2 = mean([roll_expectation(X, dist, Poisson(5.1)) for _ = 1:n_mc])
    S_rr_3 = mean([roll_expectation(X, dist, LogitNPD()) for _ = 1:n_mc])
    S_rr_4 = mean([roll_expectation(X, dist, dist) for _ = 1:n_mc])

    @test S_rr_1 ≈ S_true atol=0.1
    @test S_rr_2 ≈ S_true atol=0.1
    @test S_rr_3 ≈ S_true atol=0.1
    @test S_rr_4 ≈ S_true atol=0.1
end


text:
using RoboGrid
using Darwin
using CGP
import Base.isless
import JSON.json

CGP.Config.init("cfg/cgp.yaml")

function robo_eval(ind::Individual, nin::Int64, nout::Int64)
    c = PCGPChromo(ind.genes, nin, nout)
    f(x::Array{Float64}) = process(c, x)
    [ind.func(f; seed=ind.seed)]
end

mutable struct CGPInd <: Individual
    genes::Array{Float64}
    fitness::Array{Float64}
    func::Function
    seed::Int64
end

isless(i1::CGPInd, i2::CGPInd) = i1.fitness[1] < i2.fitness[1]

function CGPInd(cfg::Dict)
    c = PCGPChromo(cfg["nin"], cfg["nout"])
    CGPInd(c.genes, [-Inf], RoboGrid.coverage_fitness, 1)
end

function CGPInd(genes::Array{Float64}, fitness::Array{Float64})
    CGPInd(genes, fitness, RoboGrid.coverage_fitness, 1)
end

function json(i::CGPInd)
    json(Dict("genes"=>i.genes, "fitness"=>i.fitness))
end

function generation(e::Evolution)
  fits = [RoboGrid.coverage_fitness, RoboGrid.water_memorize_fitness,
            RoboGrid.cross_memorize_fitness, RoboGrid.cross_strategy_fitness]
    for i in e.population
        i.seed = floor(Int64, (e.gen - 1) / 10)
        func = mod(floor(Int64, (e.gen - 1) / e.cfg["goal_gen"]), length(fits)) + 1
        i.func = fits[func]
    end
end


text:
module TimeSeriesA3P6

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD,
  Gadfly, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  ParallelAccelerator, Measures, Formatting

importall CTIO, CTReg, CTStat

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const FRED_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

const DEFAULT_BURN = 1_000
const NUM_SIMS = 10_000
const DEFAULT_LAG = 5

#generates the VAR data
function getVARData!(T::Int,
    y::Vector{Float64} = Vector{Float64}(T),
    x::Vector{Float64} = Vector{Float64}(T))::NTuple{2, Vector{Float64}}

  #pregenerate RVs
  ε::Vector{Float64} = rand(Normal(), T)
  η::Vector{Float64} = rand(Normal(), T)

  y[1] = 0.0
  x[1] = 0.0

  @fastmath for i::Int ∈ 2:T
    y[i] = x[i-1] + ε[i]
    x[i] = x[i-1] + η[i]
  end

  return (y, x)
end

#hold the regression data
mutable struct VARReg
  N::Int
  T::Int
  lag::Int
  βs::Matrix{Float64}
  SEs::Matrix{Float64}
  Ts::Matrix{Float64}
end

function VARReg(N::Int, T::Int; lag::Int=DEFAULT_LAG, verbose::Bool = false,
    βs::Matrix{Float64} = Matrix{Float64}(N,2),
    SEs::Matrix{Float64} = Matrix{Float64}(N,2),
    Ts::Matrix{Float64} = Matrix{Float64}(N,2))::VARReg

  #preallocate regression components
  Z::Matrix{Float64} = ones(T,2)
  Y::Vector{Float64} = Vector{Float64}(T)
  z2::Vector{Float64} = Vector{Float64}(T) #holds VAR results for Z

  ZNames::Vector{Symbol} = [:intercept, :Z]
  YName::Symbol = :Y

  #container for holding the X data, even though we don't need it
  X::Vector{Float64} = Vector{Float64}(T) #holds VAR results for X

  for i::Int ∈ 1:N
    getVARData!(T, z2, X)
    Z[:, 2] .= z2
    getVARData!(T, Y, X)
    reg::CTLM = CTLM(Z, Y, XNames=ZNames, YName=YName)
    βs[i,:] .= reg.β
    SEs[i,:] .= (diag(getNeweyWest!(reg, lag))).^0.5
    Ts[i,:] .= βs[i,:] ./ SEs[i,:]

    if verbose
      println("NW: $(SEs[i,:]), NWSlow: $((diag(getNeweyWestSlow(reg, lag))).^0.5)")
    end
  end

  return VARReg(N, T, lag, βs, SEs, Ts)
end

#convenience method for the above
VARReg(;N::Int=NUM_SIMS, T::Int=50, lag=DEFAULT_LAG,
    βs::Matrix{Float64} = Matrix{Float64}(N,2),
    SEs::Matrix{Float64} = Matrix{Float64}(N,2),
    Ts::Matrix{Float64} = Matrix{Float64}(N,2))::VARReg =
  VARReg(N, T, lag=lag, βs=βs, SEs=SEs, Ts=Ts)

function visualizeA3P6(VRs::Vector{VARReg})

  #convenience containers
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  #convenience method for number formatting
  n2s(x::T where T<:Real) = num2Str(x, 3, Ints=true)

  push!(plotNames, "Betas")
  push!(plots, Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.βs[:,2], Geom.histogram(bincount=25, density=true, position=:dodge),
        Guide.xlabel("β"),Guide.ylabel("Density"),
        Guide.title("T=$(n2s(VR.T)), " *
          "μ(β)=$(n2s(mean(VR.βs[:,2]))), σ(β)=$(n2s(std(VR.βs[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "SEs")
  push!(plots,Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.SEs[:,2], Geom.histogram(bincount=25, density=true, position=:dodge),
        Guide.xlabel("SE"),Guide.ylabel("Density"),Coord.Cartesian(xmin=0.0, xmax=1.0),
        Guide.title("T=$(n2s(VR.T)), " *
          "μ(SE)=$(n2s(mean(VR.SEs[:,2]))), σ(SE)=$(n2s(std(VR.SEs[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "Ts")
  push!(plots,Vector{Plot}())
  for VR::VARReg ∈ VRs
    push!(plots[end],
      plot(x=VR.Ts[:,2], Geom.histogram(bincount=25, density=true, position=:dodge),
        Guide.xlabel("t-stat (Newey-West Lag=5)"),Guide.ylabel("Density"),
        Guide.title("T=$(n2s(VR.T)), " *
          "μ(t)=$(n2s(mean(VR.Ts[:,2]))), σ(t)=$(n2s(std(VR.Ts[:,2])))")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
  end

end

#script for the problem
function A3P6Script()::Void
  y::Vector{Float64}, x::Vector{Float64} = getVARData!(10_000)
  println("Autocorrelation: $(autoρ(y,1))")
  println("Autocorrelation 1D: $(autoρ(y[2:end]-y[1:(end-1)],1))")
  #VARReg(100, 120, verbose=true)
  @time VRs::Vector{VARReg} = [ VARReg(N=10_000, T=120), VARReg(N=10_000, T=600)]
  visualizeA3P6(VRs)

  return nothing
end

@time begin
  A3P6Script()
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
include("utils.jl")

using HybridArrays, Random
using Manifolds: default_metric_dispatch
using StaticArrays: Dynamic

Random.seed!(42)

@testset "Power manifold" begin

    Ms = Sphere(2)
    Ms1 = PowerManifold(Ms, 5)
    @test power_dimensions(Ms1) == (5,)
    @test manifold_dimension(Ms1) == 10
    @test injectivity_radius(Ms1) == π
    Ms2 = PowerManifold(Ms, 5, 7)
    @test power_dimensions(Ms2) == (5,7)
    @test manifold_dimension(Ms2) == 70
    Mr = Manifolds.Rotations(3)
    Mr1 = PowerManifold(Mr, 5)
    Mrn1 = PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5)
    @test manifold_dimension(Mr1) == 15
    @test manifold_dimension(Mrn1) == 15
    Mr2 = PowerManifold(Mr, 5, 7)
    Mrn2 = PowerManifold(Mr, Manifolds.NestedPowerRepresentation(), 5, 7)
    @test manifold_dimension(Mr2) == 105
    @test manifold_dimension(Mrn2) == 105

    @test repr(Ms1) == "PowerManifold(Sphere(2), 5)"
    @test repr(Mrn1) == "PowerManifold(Rotations(3), NestedPowerRepresentation(), 5)"

    @test Manifolds.allocation_promotion_function(Ms, exp, ([1],)) == Manifolds.allocation_promotion_function(Ms1, exp, (1,))

    @test Ms^5 === Oblique(3,5)
    @test Ms^(5,) === Ms1
    @test Mr^(5, 7) === Mr2

    @test is_default_metric(Ms1, PowerMetric())
    @test default_metric_dispatch(Ms1, PowerMetric()) === Val{true}()
    types_s1 = [Array{Float64,2},
                HybridArray{Tuple{3,Dynamic()}, Float64, 2}]
    types_s2 = [Array{Float64,3},
                HybridArray{Tuple{3,Dynamic(),Dynamic()}, Float64, 3}]

    types_r1 = [Array{Float64,3},
                HybridArray{Tuple{3,3,Dynamic()}, Float64, 3}]

    types_rn1 = [Vector{Matrix{Float64}}, ]
    TEST_STATIC_SIZED && push!(types_rn1, Vector{MMatrix{3,3,Float64}})

    types_r2 = [Array{Float64,4},
                HybridArray{Tuple{3,3,Dynamic(),Dynamic()}, Float64, 4}]
    types_rn2 = [Matrix{Matrix{Float64}}]

    retraction_methods = [Manifolds.PowerRetraction(ManifoldsBase.ExponentialRetraction())]
    inverse_retraction_methods = [Manifolds.InversePowerRetraction(ManifoldsBase.LogarithmicInverseRetraction())]

    sphere_dist = Manifolds.uniform_distribution(Ms, @SVector [1.0, 0.0, 0.0])
    power_s1_pt_dist = Manifolds.PowerPointDistribution(Ms1, sphere_dist, randn(Float64, 3, 5))
    power_s2_pt_dist = Manifolds.PowerPointDistribution(Ms2, sphere_dist, randn(Float64, 3, 5, 7))
    sphere_tv_dist = Manifolds.normal_tvector_distribution(Ms, (@MVector [1.0, 0.0, 0.0]), 1.0)
    power_s1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms1), rand(power_s1_pt_dist), sphere_tv_dist)
    power_s2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Ms2), rand(power_s2_pt_dist), sphere_tv_dist)

    id_rot = @SMatrix [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0]
    rotations_dist = Manifolds.normal_rotation_distribution(Mr, id_rot, 1.0)
    power_r1_pt_dist = Manifolds.PowerPointDistribution(Mr1, rotations_dist, randn(Float64, 3, 3, 5))
    power_rn1_pt_dist = Manifolds.PowerPointDistribution(Mrn1, rotations_dist, [randn(Float64, 3, 3) for i in 1:5])
    power_r2_pt_dist = Manifolds.PowerPointDistribution(Mr2, rotations_dist, randn(Float64, 3, 3, 5, 7))
    power_rn2_pt_dist = Manifolds.PowerPointDistribution(Mrn2, rotations_dist, [randn(Float64, 3, 3) for i in 1:5, j in 1:7])
    rotations_tv_dist = Manifolds.normal_tvector_distribution(Mr, MMatrix(id_rot), 1.0)
    power_r1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr1), rand(power_r1_pt_dist), rotations_tv_dist)
    power_rn1_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mrn1), rand(power_rn1_pt_dist), rotations_tv_dist)
    power_r2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mr2), rand(power_r2_pt_dist), rotations_tv_dist)
    power_rn2_tv_dist = Manifolds.PowerFVectorDistribution(TangentBundleFibers(Mrn2), rand(power_rn2_pt_dist), rotations_tv_dist)

    trim(s::String) = s[1:min(length(s), 20)]

    basis_types = (DefaultOrthonormalBasis(),
        ProjectedOrthonormalBasis(:svd)
    )
    for T in types_s1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_s1_pt_dist)) for _ in 1:3]
            @test injectivity_radius(Ms1, pts1[1]) == π
            basis_diag = get_basis(
                Ms1,
                pts1[1],
                DiagonalizingOrthonormalBasis(log(Ms1, pts1[1], pts1[2])),
            )
            basis_arb = get_basis(Ms1, pts1[1], DefaultOrthonormalBasis())
            test_manifold(
                Ms1,
                pts1;
                test_reverse_diff = true,
                test_musical_isomorphisms = true,
                test_injectivity_radius = false,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_s1_pt_dist],
                tvector_distributions = [power_s1_tv_dist],
                basis_types_to_from = (basis_diag, basis_arb, basis_types...),
                rand_tvector_atol_multiplier = 6.0,
                retraction_atol_multiplier = 12.0,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end
    for T in types_s2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_s2_pt_dist)) for _ in 1:3]
            test_manifold(
                Ms2,
                pts2;
                test_reverse_diff = true,
                test_musical_isomorphisms = true,
                test_injectivity_radius = false,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_s2_pt_dist],
                tvector_distributions = [power_s2_tv_dist],
                rand_tvector_atol_multiplier = 6.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end

    for T in types_r1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_r1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr1,
                pts1;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_r1_pt_dist],
                tvector_distributions = [power_r1_tv_dist],
                basis_types_to_from = basis_types,
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0
            )
        end
    end

    for T in types_rn1
        @testset "Type $(trim(string(T)))..." begin
            pts1 = [convert(T, rand(power_rn1_pt_dist)) for _ in 1:3]
            test_manifold(
                Mrn1,
                pts1;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_rn1_pt_dist],
                tvector_distributions = [power_rn1_tv_dist],
                basis_types_to_from = basis_types,
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0
            )
        end
    end
    for T in types_r2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_r2_pt_dist)) for _ in 1:3]
            test_manifold(
                Mr2,
                pts2;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_r2_pt_dist],
                tvector_distributions = [power_r2_tv_dist],
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end
    for T in types_rn2
        @testset "Type $(trim(string(T)))..." begin
            pts2 = [convert(T, rand(power_rn2_pt_dist)) for _ in 1:3]
            test_manifold(
                Mrn2,
                pts2;
                test_reverse_diff = false,
                test_injectivity_radius = false,
                test_musical_isomorphisms = true,
                test_vee_hat = true,
                retraction_methods = retraction_methods,
                inverse_retraction_methods = inverse_retraction_methods,
                point_distributions = [power_rn2_pt_dist],
                tvector_distributions = [power_rn2_tv_dist],
                rand_tvector_atol_multiplier = 5.0,
                retraction_atol_multiplier = 12,
                is_tangent_atol_multiplier = 12.0,
            )
        end
    end

    @testset "Power manifold of Circle" begin
        pts_t = [[0.0, 1.0, 2.0], [1.0, 1.0, 2.4], [0.0, 2.0, 1.0]]
        MT = PowerManifold(Circle(), 3)
        @test representation_size(MT) == (3,)
        test_manifold(
            MT,
            pts_t;
            test_reverse_diff = false,
            test_forward_diff = false,
            test_injectivity_radius = false,
            test_musical_isomorphisms = true,
            test_vee_hat = true,
            retraction_methods = retraction_methods,
            inverse_retraction_methods = inverse_retraction_methods,
            rand_tvector_atol_multiplier = 5.0,
            retraction_atol_multiplier = 12,
            is_tangent_atol_multiplier = 12.0,
        )
    end

    @testset "Basis printing" begin
        p = hcat([[1.0, 0.0, 0.0] for i in 1:5]...)
        Bc = get_basis(Ms1, p, DefaultOrthonormalBasis())
        @test sprint(show, "text/plain", Bc) == """
        DefaultOrthonormalBasis(ℝ) for a power manifold
        Basis for component (1,):
        DefaultOrthonormalBasis(ℝ) with 2 basis vectors:
         E1 =
          3-element Array{Int64,1}:
           0
           1
           0
         E2 =
          3-element Array{Int64,1}:
           0
           0
           1
        Basis for component (2,):
        DefaultOrthonormalBasis(ℝ) with 2 basis vectors:
         E1 =
          3-element Array{Int64,1}:
           0
           1
           0
         E2 =
          3-element Array{Int64,1}:
           0
           0
           1
        Basis for component (3,):
        DefaultOrthonormalBasis(ℝ) with 2 basis vectors:
         E1 =
          3-element Array{Int64,1}:
           0
           1
           0
         E2 =
          3-element Array{Int64,1}:
           0
           0
           1
        Basis for component (4,):
        DefaultOrthonormalBasis(ℝ) with 2 basis vectors:
         E1 =
          3-element Array{Int64,1}:
           0
           1
           0
         E2 =
          3-element Array{Int64,1}:
           0
           0
           1
        Basis for component (5,):
        DefaultOrthonormalBasis(ℝ) with 2 basis vectors:
         E1 =
          3-element Array{Int64,1}:
           0
           1
           0
         E2 =
          3-element Array{Int64,1}:
           0
           0
           1
        """
    end

    @testset "Power manifold of Circle" begin
        pts_t = [[0.0, 1.0, 2.0], [1.0, 1.0, 2.4], [0.0, 2.0, 1.0]]
        MT = PowerManifold(Circle(), 3)
        @test representation_size(MT) == (3,)
        test_manifold(
            MT,
            pts_t;
            test_reverse_diff = false,
            test_forward_diff = false,
            test_injectivity_radius = false,
            test_musical_isomorphisms = true,
            retraction_methods = retraction_methods,
            inverse_retraction_methods = inverse_retraction_methods,
            rand_tvector_atol_multiplier = 5.0,
            retraction_atol_multiplier = 12,
            is_tangent_atol_multiplier = 12.0,
        )
    end

end


text:
using LazySets

function lorenz_specification()
    # initial set
    X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])

    # time horizon: 10 time units
    time_horizon = 10.0

    # specification
    O = Dict{Symbol, Any}(:T => time_horizon)

    return X0, O
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
function convert_second_order(model::M) where M <: Lombardo_Sutherland_Soln

  ssh       = copy(model.ssh)
  hx        = copy(model.hx)
  hxx       = copy(model.hxx)
  ssg       = copy(model.ssg)
  gx        = copy(model.gx)
  gxx       = copy(model.gxx)
  eta       = copy(model.eta)
  sigma     = copy(model.sigma)
  grc       = copy(model.grc)
  soln_type = model.soln_type

  (n1,n2) = size(hx)
  (n3,n4) = size(gx)

  new_hxx = reshape(vech_to_vec(n2)*hxx[1,:],n2,n2)
  for i = 2:n1
    new_hxx = [new_hxx; reshape(vech_to_vec(n2)*hxx[i,:],n2,n2)]
  end

  (n5,n6) = size(new_hxx)

  j = 0
  for i = 1:n5
    j += 1
    new_hxx[i,j] = 2*new_hxx[i,j]
    if j == n6
      j = 0
    end
  end

  new_hxx = new_hxx/2

  new_gxx = reshape(vech_to_vec(n4)*gxx[1,:],n4,n4)
  for i = 2:n3
    new_gxx = [new_gxx; reshape(vech_to_vec(n4)*gxx[i,:],n4,n4)]
  end

  (n7,n8) = size(new_gxx)

  j = 0
  for i = 1:n7
    j += 1
    new_gxx[i,j] = 2*new_gxx[i,j]
    if j == n8
      j = 0
    end
  end

  new_gxx = new_gxx/2

  soln = Gomme_Klein_Soln(ssh,hx,new_hxx,ssg,gx,new_gxx,eta,sigma,grc,soln_type)

  return soln

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:2] <= 1)

                     @variable(m, -0.9705716506753843 <= q <= 1.6199082189427347)

                     add_NL_constraint(m, :(sigmoid(0.27812717705086953 + -0.4967460082986306*sigmoid(-0.6740257057257448 + -0.0007917455797805673*sigmoid(-0.08584163188401472 + 0.6683752749724832*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*sigmoid(-0.7769463282105753 + 0.25325339074048747*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + -0.883304587160544*sigmoid(0.3709718635886192 + -0.47514710764751555*sigmoid(-0.08584163188401472 + 0.6683752749724832*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*sigmoid(-0.7769463282105753 + 0.25325339074048747*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) + sigmoid(0.48450040123340576 + -0.32897580116872405*sigmoid(-0.6740257057257448 + -0.0007917455797805673*sigmoid(-0.08584163188401472 + 0.6683752749724832*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.5330540695702357*sigmoid(-0.7769463282105753 + 0.25325339074048747*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2])))) + 0.43090905819531766*sigmoid(0.3709718635886192 + -0.47514710764751555*sigmoid(-0.08584163188401472 + 0.6683752749724832*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + -0.39230803394698555*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))) + 0.6627028336583534*sigmoid(-0.7769463282105753 + 0.25325339074048747*sigmoid(0.3679628458703208 + -0.7706264460934666*$(x[1]) + -0.03424767458492983*$(x[2])) + 0.5840452662855768*sigmoid(0.40575897127542815 + 0.13151803175485233*$(x[1]) + 0.5944153777408365*$(x[2]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
This script acts as a bridge between "einsteinpy.geodesic" and "KerrSolver.jl".
It solves the system, using "KerrSolver.jl", and saves the output to .csv files in
`/tmp/epy_geod_jl_temp/` (on *nix) or `%TEMP%/epy_geod_jl_temp/` (on Windows).
These are accessed by the `_julia_wrapper` utility function in EinsteinPy.

# Stand-alone Usage Example
```
> julia run.jl 10., 1.57, 0., 0., 0., 2.427, 0.99, 0.9, 0.5, 0.05
("Success", "C:\\Users\\UserName\\AppData\\Local\\Temp\\epy_geod_jl_temp")
```
where the argument order is q1, q2, q3, p1, p2, p3, a, E, end_lambda, step_size.
See "einsteinpy.geodesic.Geodesic" or "KerrSolver.jl" for information on the arguments.
The output is a 2-tuple of strings, containing the "return code" from the solver
and "folder-path", where the output is stored.

"""

import Base
import DelimitedFiles: writedlm
import Main: include

Main.include("KerrSolver.jl")

function run(q, p, params)
    # Solving
    retcode, lambdas, vecs = Main.KerrSolver.solveSystem(q, p, params, end_lambda, step_size)

    # Saving to file in 'temp' directory
    temp_dir_name = "epy_geod_jl_temp/"
    cd(tempdir())
    mkpath(temp_dir_name) # Does nothing, if folder exists
    cd(temp_dir_name)

    # Writing to csvs (Overwrites existing files)
    writedlm("lambdas.csv", lambdas, ",")
    writedlm("vecs.csv", vecs, ",")

    # Return retcode and folder-path as a 2-tuple of strings
    return retcode, pwd()
end

# Parsing input
q1, q2, q3, p1, p2, p3, a, E, end_lambda, step_size = map(x -> parse(Float64, x), Base.ARGS)

# Packaging arguments
q = [q1, q2, q3]
p = [p1, p2, p3]
params = [a, E]

# Final output
println(run(q, p, params))

max_stars_repo_name: einsteinpy/einsteinpy-geodesics


text:
# using Distributed
# addprocs(1)
 include("par_module.jl")

 @everywhere begin
    include("par_module.jl")
    push!(LOAD_PATH, "/Users/brandonbehring/Desktop/Leap_Frog_2019")
    using .par_module
    N = 1000
    E = ones(N) * .5
    A = 1:N;
    B = A + ones(N)
end

pmap(par_module.f, A, B, E)


text:
# 1D data
@recipe function f(metric::Metric{T, 1}; axes="nil", coordinates=nothing) where {T <: UltrasoundMetric}
    axes == "nil" && error("Keyword missing: axes=\"x\" or axes=\"y\"..")
    check_coordinates(coordinates)

    title --> get_metric_display_name(metric)
    xguide --> "$(axes)-axis"
    aspect_ratio --> :equal
    legend --> false

    target_dimension = 1
    axis_coordinates = get_axis_coordinates(axes, target_dimension, coordinates)

    return axis_coordinates, metric.val
end

function check_coordinates(coordinates)
    if isnothing(coordinates)
        error("You must set coordinates from scan in keyword arg:
              ;coordinates=scan.coordinates")
    end
end

function get_axis_coordinates(axis_name, target_dimension, coordinates)
    axis_name = Val(Symbol(axis_name))
    target_dimension = Val(target_dimension)
    axis_coordinates = get_coordinates(axis_name, get_dimension(target_dimension, coordinates))
    return round.(typeof(1.0mm), uconvert.(mm, axis_coordinates); digits=2)
end

get_dimension(::Val{1}, coordinates) = coordinates[:, :, 1, 1]
get_dimension(::Val{2}, coordinates) = coordinates[:, 1, :, 1]
get_dimension(::Val{3}, coordinates) = coordinates[:, 1, 1, :]

get_coordinates(::Val{:x}, coordinates) = coordinates[1, :]
get_coordinates(::Val{:y}, coordinates) = coordinates[2, :]
get_coordinates(::Val{:z}, coordinates) = coordinates[3, :]

# 2D data
@recipe function f(metric::Metric{T, 2}; axes="nil", coordinates=nothing) where {T <: UltrasoundMetric}
    axes == "nil" && error("Keyword missing: axes=\"xy\" or axes=\"yz\"..")
    check_coordinates(coordinates)

    seriestype := :heatmap
    title --> get_metric_display_name(metric)
    xguide --> "$(axes[1])-axis"
    yguide --> "$(axes[2])-axis"
    aspect_ratio --> :equal

    metric_values_2d = metric.val

    axis_1_coordinates = get_axis_coordinates(axes[1], 1, coordinates)
    axis_2_coordinates = get_axis_coordinates(axes[2], 2, coordinates)
    axis_1_step_size = diff(axis_1_coordinates)[1]
    axis_2_step_size = diff(axis_2_coordinates)[1]
    if axis_1_step_size < 0mm
        axis_1_coordinates = reverse(axis_1_coordinates)
        axis_1_step_size = diff(axis_1_coordinates)[1]
        metric_values_2d = reverse(metric_values_2d, dims=1)
    end
    if axis_2_step_size < 0mm
        axis_2_coordinates = reverse(axis_2_coordinates)
        axis_2_step_size = diff(axis_2_coordinates)[1]
        metric_values_2d = reverse(metric_values_2d, dims=2)
    end
    xlims --> (axis_1_coordinates[1]-axis_1_step_size/2, axis_1_coordinates[end]+axis_1_step_size/2)
    ylims --> (axis_2_coordinates[1]-axis_2_step_size/2, axis_2_coordinates[end]+axis_2_step_size/2)

    return axis_1_coordinates, axis_2_coordinates, transpose(metric_values_2d)
end

# 3D data
@recipe function f(metric::Metric{T, 3}; coordinates=nothing, axes="xyz") where {T <: UltrasoundMetric}
    check_coordinates(coordinates)

    # remove the keywords from the attribute dictionary to avoid potential errors
    xslice = pop!(plotattributes, :xslice)
    yslice = pop!(plotattributes, :yslice)
    zslice = pop!(plotattributes, :zslice)

    slice_position, slice_axis_name, slice_axis_dim = get_slice_info(xslice, yslice, zslice)
    slice_axis_coordinates = get_axis_coordinates(slice_axis_name, slice_axis_dim, coordinates)
    slice_idx, nearest_slice_position = get_slice_idx(slice_axis_coordinates, slice_position)
    metric_values_2d = get_slice_data(Val(slice_axis_name), slice_idx, metric.val)
    _, slice_image_dims = get_slice_image_axes(Val(slice_axis_name))

    seriestype := :heatmap
    title := get_metric_display_name(metric) * "\n($(slice_axis_name)=$nearest_slice_position)"
    xguide --> "$(axes[slice_image_dims[1]])-axis"
    yguide --> "$(axes[slice_image_dims[2]])-axis"
    aspect_ratio --> :equal

    axis_1_coordinates = get_axis_coordinates(axes[slice_image_dims[1]], slice_image_dims[1], coordinates)
    axis_2_coordinates = get_axis_coordinates(axes[slice_image_dims[2]], slice_image_dims[2], coordinates)
    axis_1_step_size = diff(axis_1_coordinates)[1]
    axis_2_step_size = diff(axis_2_coordinates)[1]
    if axis_1_step_size < 0mm
        axis_1_coordinates = reverse(axis_1_coordinates)
        axis_1_step_size = diff(axis_1_coordinates)[1]
        metric_values_2d = reverse(metric_values_2d, dims=1)
    end
    if axis_2_step_size < 0mm
        axis_2_coordinates = reverse(axis_2_coordinates)
        axis_2_step_size = diff(axis_2_coordinates)[1]
        metric_values_2d = reverse(metric_values_2d, dims=2)
    end
    xlims --> (axis_1_coordinates[1]-axis_1_step_size/2, axis_1_coordinates[end]+axis_1_step_size/2)
    ylims --> (axis_2_coordinates[1]-axis_2_step_size/2, axis_2_coordinates[end]+axis_2_step_size/2)

    return axis_1_coordinates, axis_2_coordinates, transpose(metric_values_2d)
end

function get_slice_info(xslice, yslice, zslice)
    if reduce(+, isnan.([xslice, yslice, zslice])) != 2
        error("Only one of the keyword arguments must be set:
              xslice, yslice, or zslice")
    end
    if !isnan(xslice)
        slice_position = xslice
        slice_axis_name = :x
        slice_axis_dim = 1
    elseif !isnan(yslice)
        slice_position = yslice
        slice_axis_name = :y
        slice_axis_dim = 2
    elseif !isnan(zslice)
        slice_position = zslice
        slice_axis_name = :z
        slice_axis_dim = 3
    end
    return slice_position, slice_axis_name, slice_axis_dim
end

function get_slice_idx(slice_axis_positions, slice_position)
    slice_idx = findnearest(slice_axis_positions, uconvert(u"mm", slice_position))
    nearest_slice_position = slice_axis_positions[slice_idx]
    return slice_idx, nearest_slice_position
end

findnearest(A::AbstractArray,t) = findmin(abs.(A.-t))[2]

get_slice_data(::Val{:x}, slice_idx, data) = data[slice_idx, :, :]
get_slice_data(::Val{:y}, slice_idx, data) = data[:, slice_idx, :]
get_slice_data(::Val{:z}, slice_idx, data) = data[:, :, slice_idx]

get_slice_image_axes(::Val{:x}) = "yz", [2,3]
get_slice_image_axes(::Val{:y}) = "xz", [1,3]
get_slice_image_axes(::Val{:z}) = "xy", [1,2]


text:
using Ranges
using Base.Test

include("linspace.jl")


text:
module EzXML

export
    # Types
    # -----
    # NOTE: These are not exported but considered as public APIs.
    # Node,
    # Document,
    # XMLError,
    # StreamReader,

    # Node constructors
    # -----------------
    XMLDocumentNode,
    HTMLDocumentNode,
    ElementNode,
    TextNode,
    CommentNode,
    CDataNode,
    AttributeNode,
    DTDNode,

    # Document constructors
    # ---------------------
    XMLDocument,
    HTMLDocument,

    # Functions
    # ---------
    hasparentnode,
    parentnode,
    hasparentelement,
    parentelement,
    hasnode,
    firstnode,
    lastnode,
    haselement,
    firstelement,
    lastelement,
    hasnextnode,
    nextnode,
    hasprevnode,
    prevnode,
    hasnextelement,
    nextelement,
    hasprevelement,
    prevelement,
    countnodes,
    countelements,
    countattributes,
    link!,
    linknext!,
    linkprev!,
    unlink!,
    addelement!,
    hasversion,
    version,
    hasencoding,
    encoding,
    hasroot,
    root,
    setroot!,
    hasdtd,
    dtd,
    setdtd!,
    nodetype,
    nodepath,
    iselement,
    isattribute,
    istext,
    iscdata,
    iscomment,
    isdtd,
    hasdocument,
    document,
    nodename,
    setnodename!,
    hasnodecontent,
    nodecontent,
    setnodecontent!,
    systemID,
    externalID,
    eachnode,
    nodes,
    eachelement,
    elements,
    eachattribute,
    attributes,
    namespace,
    namespaces,
    prettyprint,
    readxml,
    readhtml,
    parsexml,
    parsehtml,
    readdtd,
    validate,
    nodedepth,
    expandtree,
    hasnodevalue,
    nodevalue,
    hasnodeattributes,
    nodeattributes

using Libdl
using Printf: @printf

# Load libxml2.
const libxml2path = joinpath(dirname(@__FILE__), "..", "deps", "deps.jl")
if !isfile(libxml2path)
    error("EzXML.jl is not installed properly, run Pkg.build(\"EzXML\") and restart Julia.")
end
include(libxml2path)
check_deps()

include("error.jl")
include("node.jl")
include("document.jl")
include("buffer.jl")
include("xpath.jl")
include("streamreader.jl")

function __init__()
    init_error_handler()
end

# deprcated methods
import Base: findall, findfirst, findlast
@deprecate findall(doc::Document, xpath::AbstractString)   findall(xpath, doc)
@deprecate findfirst(doc::Document, xpath::AbstractString) findfirst(xpath, doc)
@deprecate findlast(doc::Document, xpath::AbstractString)  findlast(xpath, doc)
@deprecate findall(node::Node, xpath::AbstractString, ns=namespaces(node))   findall(xpath, node, ns)
@deprecate findfirst(node::Node, xpath::AbstractString, ns=namespaces(node)) findfirst(xpath, node, ns)
@deprecate findlast(node::Node, xpath::AbstractString, ns=namespaces(node))  findlast(xpath, node, ns)

end # module


text:
using Caesar
using IncrementalInference, RoME
using HDF5, JLD
using Base.Test

println("[TEST] CloudGraphs API calls...")
if false
  println("[TEST] with CloudGraphs with local DB data layer (multicore)...")
  include("fourdoortestcloudgraph.jl")
  println("[SUCCESS]")
else
  warn("[NOT TESTING] CloudGraphs interface, auth required -- you can enable it here at $(@__FILE__).")
end

# if false
#   # using VictoriaParkTypes
#
# else
#   warn("not running full Victoria park dataset now")
# end


text:
#Dot and whisker plots of coefficient values + confidence intervals

#For submodels
"""
    whiskerplot(m::UnmarkedSubmodel; level::Real = 0.95)

Generate a dot-and-whisker plot for all parameters in a given submodel `m`.
Error bars are drawn for given confidence `level`.
"""
function whiskerplot(m::UnmarkedSubmodel; level::Real = 0.95)

  #Build data frame of coefficients and upper/lower bounds
  df = DataFrame(names=coefnames(m), coef=coef(m),
                 se=stderror(m))
  zval = -quantile(Normal(), (1-level)/2)
  df[!,:lower] = df.coef - zval * df.se
  df[!,:upper] = df.coef + zval * df.se

  plot(df, x=:names, y=:coef, ymin=:lower, ymax=:upper,
       yintercept=[0.0],
       style(major_label_font_size=18pt,
             minor_label_font_size=16pt,
             point_size=5pt,line_width=2pt),
       Guide.xlabel("Parameter"), Guide.ylabel("Value"),
       Guide.title(string(m.name)),
       Geom.point, Geom.errorbar,
       Geom.hline(color="orange",style=:dash))

end

#For complete model
"""
    whiskerplot(m::UnmarkedModel; level::Real = 0.95)

Generate a dot-and-whisker plot for all parameters in a given model `m`.
Error bars are drawn for given confidence `level`.
"""
function whiskerplot(m::UnmarkedModel; level::Real = 0.95)
  sm = collect(m.submodels)
  vstack(whiskerplot.(sm, level=level))
end

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

#Partial effects plots for each parameter

#Get array filled with "baseline" value for variable (mean or reference level)
function get_var_baseline(x::Array{<:Number}, nrep::Int)
  return repeat([mean(x)], nrep)
end

function get_var_baseline(x::CategoricalArray, nrep::Int)
  out = categorical(repeat([levels(x)[1]], nrep))
  levels!(out, levels(x))
  return out
end

#Get sequence of values for focal parameter
function get_var_seq(x::Array{<:Number})
  nrow = 100
  out = collect(range(min(x...), stop=max(x...), length=nrow))
  return (nrow, out)
end

function get_var_seq(x::CategoricalArray)
  nrow = length(levels(x))
  out = categorical(levels(x))
  levels!(out, levels(x))
  return (nrow, out)
end

#For a single parameter
"""
    effectsplot(m::UnmarkedSubmodel, param::String; level::Real = 0.95)

Generate a partial effects plot for parameter `param` from submodel `m`.
The response variable is plotted on the original scale. A confidence envelope
is also plotted for the given `level`.
"""
function effectsplot(m::UnmarkedSubmodel, param::String; level::Real = 0.95)

  resp = string(m.formula.lhs.sym)
  psym = Symbol(param)
  dat = deepcopy(m.data)

  #Build newdata for predict
  nrow, var_seq = get_var_seq(dat[:, psym])
  nd = mapcols(col -> get_var_baseline(col, nrow), dat)
  nd[!, psym] = var_seq

  #Predict values
  pr = DataFrame(predict(m, nd, interval=true, level=level))
  pr[!, psym] = var_seq

  #Choose plot type depending on variable type
  estplot = Geom.line; errorplot = Geom.ribbon; lw = 1pt
  if typeof(var_seq) <: CategoricalArray
    estplot = Geom.point; errorplot = Geom.errorbar; lw = 2pt
  end

  plot(pr, x=psym, y=:prediction, ymin=:lower, ymax=:upper,
       estplot, errorplot, Guide.xlabel(param), Guide.ylabel(resp),
       style(major_label_font_size=18pt, minor_label_font_size=16pt,
             point_size=5pt, line_width=lw))

end

#For a submodel
"""
    effectsplot(m::UnmarkedSubmodel; level::Real = 0.95)

Partial effects plot for all parameters in a submodel `m`.
"""
function effectsplot(m::UnmarkedSubmodel; level::Real = 0.95)
  vars = string.(collect(m.formula.rhs))
  vstack(map(x -> effectsplot(m, x, level=level), vars))
end

#For a complete model
"""
    effectsplot(m::UnmarkedModel; level::Real = 0.95)

Partial effects plot for all parameters in a model `m`.
"""
function effectsplot(m::UnmarkedModel; level::Real = 0.95)
  cols = effectsplot.(collect(m.submodels), level=level)
  hstack(cols[1], cols[2])
end

max_stars_repo_name: kenkellner/julia-hme

max_issues_repo_name: kenkellner/julia-hme

max_forks_repo_name: kenkellner/julia-hme


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

using Markdown
using InteractiveUtils

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

# ╔═╡ f313719e-56c1-11eb-0151-7fdc92bc5635
using DrWatson

# ╔═╡ 5d891788-56c2-11eb-1723-8361ac5bd415
begin
	@quickactivate

	using Images
	using PlutoUI
	using Serialization
	using ImageFiltering
	using ImageTransformations

	pd = projectdir()

	md"""
	## Custom EDH Cards
	"""
end

# ╔═╡ 9891eb6c-88f8-11eb-261a-e5efa12c84ac
md"""
Mage.jl Deck directory path: $(@bind DECK_DIR TextField(default="/home/dusty/Shared/PlaymatProjects/PlaymatGames/Mage.jl/EDH/decks/Vannifar's Circus"))
"""

# ╔═╡ b28c3a68-88f8-11eb-0175-779fdb2da838
DECK_NAME = split(DECK_DIR, '/')[end]

# ╔═╡ 8a3f6140-56c4-11eb-1cea-ff327a21d57b
deck = deserialize("$DECK_DIR/$DECK_NAME.jls")

# ╔═╡ f78a72a0-5d48-11eb-198c-d742de310200
sz = size(deck[:CARD_BACK_IMG])

# ╔═╡ 0152ba30-5ba1-11eb-013e-8dc5191b3c17
md"""
## Custom PNGs, JPGs, etc.!
"""

# ╔═╡ 2f5a332c-56c4-11eb-259b-8b7d51af1b04
begin
	custom_card_faces = [ fn => load("$DECK_DIR/custom_images/$fn") for fn in readdir("$DECK_DIR/custom_images") if !occursin(split(fn,".")[begin], join(deck[:commander_names])) && (occursin("png", fn) || occursin("gif",fn))  ]

	card_names = [ k for (k,v) in custom_card_faces ]
	card_imgs = [ v for (k,v) in custom_card_faces ]

	for i in 1:length(card_imgs)
		csz = size(card_imgs[i])
		rat = csz[1] / csz[2]
		card_imgs[i] = imresize(card_imgs[i], ceil(Int, sz[2]*rat), sz[2])
	end
end;

# ╔═╡ 64e4a3a6-5d54-11eb-25f7-9dfff5507992
md"""
### Cards
"""

# ╔═╡ 70971440-5985-11eb-3d51-cdedac799904
md"""
card index: $(@bind card_index Slider(1:length(card_imgs), show_value=true))
"""

# ╔═╡ 898434cc-5ee6-11eb-0e99-6d69667336d9
card_imgs[card_index]

# ╔═╡ 08063488-5973-11eb-0fcd-97b6d199dd11
card_info = [ [v,k,i,size(v[begin]),size(v[begin])[1]/size(v[begin])[2] ] for (i,(k,v)) in enumerate(deck[:CARD_FACES]) if occursin(k, card_names[card_index]) ][begin]

# ╔═╡ 77ce8398-5989-11eb-2d3c-d3d0957ae270
md"""
"Replace card face with custom png / gif?" $(@bind swap_card_face CheckBox())
"""

# ╔═╡ e010eb3a-597a-11eb-19da-01375b1d8367
if swap_card_face
	deck[:CARD_FACES][ card_info[3] ] =
	card_names[card_index] => [ card_imgs[card_index] ]
end

# ╔═╡ cd27bfd0-5ba0-11eb-0bd4-559470d3907a
md"""
### Custom GIFs!
"""

# ╔═╡ 437c0424-56c5-11eb-29ce-7f0090186512
custom_gifs = [ fn=>LocalResource("$DECK_DIR/custom_images/$fn") for fn in readdir("$DECK_DIR/custom_images") if occursin("gif", fn) ];

# ╔═╡ 2d187afa-598b-11eb-1dec-7b3c22ea634d
md"""
custom gif index: $(@bind gif_index Slider(1:length(custom_gifs), show_value=true))
"""

# ╔═╡ dfa2d2e8-598a-11eb-2e3a-5f9c457e5cae
custom_gif_names = [ k for (k,v) in custom_gifs ];

# ╔═╡ 74abf0a2-5802-11eb-3bcf-79a66eabe3e5
custom_gifs[gif_index]

# ╔═╡ 3c25d5f2-5ba1-11eb-0229-0fbfc8dbbf44
md"""
save custom card face data? $(@bind save_data CheckBox())
"""

# ╔═╡ 72de84ba-598d-11eb-139d-975970c19cc0
if save_data
	serialize("$DECK_DIR/Vannifar's Circus.jls", deck)
end

# ╔═╡ Cell order:
# ╟─f313719e-56c1-11eb-0151-7fdc92bc5635
# ╟─5d891788-56c2-11eb-1723-8361ac5bd415
# ╟─9891eb6c-88f8-11eb-261a-e5efa12c84ac
# ╟─b28c3a68-88f8-11eb-0175-779fdb2da838
# ╟─8a3f6140-56c4-11eb-1cea-ff327a21d57b
# ╟─f78a72a0-5d48-11eb-198c-d742de310200
# ╟─0152ba30-5ba1-11eb-013e-8dc5191b3c17
# ╟─2f5a332c-56c4-11eb-259b-8b7d51af1b04
# ╟─64e4a3a6-5d54-11eb-25f7-9dfff5507992
# ╟─70971440-5985-11eb-3d51-cdedac799904
# ╟─898434cc-5ee6-11eb-0e99-6d69667336d9
# ╟─08063488-5973-11eb-0fcd-97b6d199dd11
# ╟─77ce8398-5989-11eb-2d3c-d3d0957ae270
# ╟─e010eb3a-597a-11eb-19da-01375b1d8367
# ╟─cd27bfd0-5ba0-11eb-0bd4-559470d3907a
# ╟─2d187afa-598b-11eb-1dec-7b3c22ea634d
# ╠═437c0424-56c5-11eb-29ce-7f0090186512
# ╠═dfa2d2e8-598a-11eb-2e3a-5f9c457e5cae
# ╠═74abf0a2-5802-11eb-3bcf-79a66eabe3e5
# ╟─3c25d5f2-5ba1-11eb-0229-0fbfc8dbbf44
# ╠═72de84ba-598d-11eb-139d-975970c19cc0


text:
# A `Table` presents itself as an `AbstractArray` of `NamedTuples`s

"""
    Table(name1 = array1, ...)

Create a column-storage-based `Table` with column names `name1`, etc, from arrays `array1`,
etc. The input arrays `array1`, etc, must share the same dimensionality and indices.

`Table` itself is an `AbstractArray` whose elements are `NamedTuple`s of the form
`(name1 = first(array1), ...)`, etc. Rows of the table are obtained via standard array
indexing `table[i]`, and columns via `table.name`.

`Table` differs from `FlexTable` in that the columns are immutable - you may add, remove,
rename and replace entire columns of a `FlexTable`, but not a `Table`. However, `Table` can
access and iterate rows in local scope with fast, fully type-inferred code while `FlexTable`
will be more efficient with a higher-order interface.
"""
struct Table{T <: NamedTuple, N, Data <: NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{<:Any,N}}}}} <: AbstractArray{T, N}
    data::Data

    # Inner constructor, to compare axes?
end

Table(;kwargs...) = Table(kwargs.data)
Table(nt::NamedTuple) = Table{_eltypes(nt), _ndims(nt), typeof(nt)}(nt)

function Table(x)
    cols = columns(x)
    if cols isa NamedTuple{<:Any, <:Tuple{Vararg{AbstractArray{N}} where N}}
        return Table(cols)
    else
        return Table(columntable(cols))
    end
end

Tables.AccessStyle(::Table) = Tables.ColumnAccess()
Tables.schema(::Table{T}) where {T} = T

"""
    columns(table::Table)

Convert a `Table` into a `NamedTuple` of it's columns.
"""
@inline columns(t::Table) = Core.getfield(t, :data)

@inline function columns(t::Table, names::NTuple{N, Symbol}) where N
	nt = NamedTuple{names}(map(name -> getproperty(t, name), names))
end

@inline function columns(t::Table, names::AbstractVector{<:Symbol})
    columns(t, tuple(names...))
end
@inline function columns(t::Table, name::Symbol)
    columns(t, tuple(name))
end
@inline function columns(t::Table, idxs::NTuple{N, Int}) where N
    all_names = columnnames(t)
    names = [all_names[idx] for idx in idxs]
    columns(t, names)
end
@inline function columns(t::Table, idxs::AbstractVector{<:Int})
    columns(t, tuple(idxs...))
end
@inline function columns(t::Table, idx::Int)
    columns(t, tuple(idx))
end

@inline rows(t::Table) = t

# Simple column access via `table.columnname`
@inline Base.getproperty(t::Table, name::Symbol) = getproperty(columns(t), name)

# Private fields are never exposed since they can conflict with column names
Base.propertynames(t::Table, private::Bool=false) = columnnames(t)

"""
    columnnames(table)

Return a tuple of the column names of a `Table`.
"""
columnnames(::AbstractArray{<:NamedTuple{names}}) where {names} = names

# Basic AbstractArray interface

@inline Base.size(t::Table) = size(first(columns(t)))
@inline Base.axes(t::Table) = axes(first(columns(t)))
@inline Base.IndexStyle(t::Table) = IndexStyle(first(columns(t)))

@inline ncolumns(t::Table) = length(columnnames(t))

function Base.checkbounds(::Type{Bool}, t::Table, i...)
    # Make sure we are in bounds for *every* column. Only safe to do
    # here because each column might mutate size independently of the others!
    all(col -> checkbounds(Bool, col, i...), columns(t))
end

@inline function Base.getindex(t::Table{T}, i::Int) where {T}
    @boundscheck checkbounds(t, i)
    # AbstractArray{T} expects a subtype of T, and NamedTuples are invariant.
    # For e.g. columns which are Union types we makes sure we emit a named tuple with Union
    # elements. The convert should be a no-op in the case of strongly-typed columns.

    # TODO find a way to make this faster than O(2^N) for N Union columns
    convert(T, map(col -> @inbounds(getindex(col, i)), columns(t)))::T
end

@inline function Base.getindex(t::Table{T}, i::Int...) where {T}
    @boundscheck checkbounds(t, i...)
    convert(T, map(col -> @inbounds(getindex(col, i...)), columns(t)))::T
end

@inline function Base.setindex!(t::Table{T}, v::T, i::Int) where {T}
    @boundscheck checkbounds(t, i)
    map((val, col) -> @inbounds(setindex!(col, val, i)), v, columns(t))
    return t
end
@inline Base.setindex!(t::Table{T}, v, i::Int) where {T} = setindex!(t, convert(T, v), i)

@inline function Base.setindex!(t::Table{T}, v::T, i::Int...) where {T}
    @boundscheck checkbounds(t, i...)
    map((val, col) -> @inbounds(setindex!(col, val, i...)), v, columns(t))
    return t
end
@inline Base.setindex!(t::Table{T}, v, i::Int...) where {T} = setindex!(t, convert(T, v), i...)

# similar

@inline Base.similar(t::Table, ::Type{NamedTuple{names, T}}, dims) where {names, T <: Tuple} = _similar(t, NamedTuple{names, T}, dims)

# Ambiguities...
@inline Base.similar(t::Table, ::Type{NamedTuple{names, T}}) where {names, T <: Tuple} = similar(t, NamedTuple{names, T}, axes(t))
@inline Base.similar(t::Table, ::Type{NamedTuple{names, T}}, dims::Union{Integer, AbstractUnitRange}...) where {names, T <: Tuple} = _similar(t, NamedTuple{names, T}, dims)
@inline Base.similar(t::Table, ::Type{NamedTuple{names, T}}, dims::Tuple{Vararg{Int64,N}}) where {names, T <: Tuple, N} = _similar(t, NamedTuple{names, T}, dims)
@inline Base.similar(t::Table, ::Type{NamedTuple{names, T}}, dims::Tuple{Union{Integer, OneTo},Vararg{Union{Integer, OneTo}}}) where {names, T <: Tuple} = _similar(t, NamedTuple{names, T}, dims)

@generated function _similar(t::Table{<:NamedTuple{oldnames}}, ::Type{NamedTuple{names, T}}, newaxes) where {oldnames, names, T <: Tuple}
    if isempty(names)
        return :(similar(first(columns(t)), NamedTuple{(),Tuple{}}, newaxes))
    end

    # Try and be clever - if the name is preserved, we preserve the similar type, otherwise
    # just be similar to the first one.
    exprs = [name ∈ oldnames ? :(similar(getproperty(t, $(QuoteNode(name))), $ElType, newaxes)) :
                               :(similar(first(columns(t)), $Eltype, newaxes)) for (name, ElType) in zip(names, T.parameters)]
    return :(Table(NamedTuple{names}(tuple($(exprs...)))))
end

# empty
Base.empty(t::Table) = empty(t, eltype(t))
@generated function Base.empty(t::Table{<:NamedTuple{oldnames}, 1}, ::Type{NamedTuple{names, T}}) where {oldnames, names, T <: Tuple}
    if isempty(names)
        return :(empty(first(columns(t)), NamedTuple{(),Tuple{}}))
    end

    # Try and be clever - if the name is preserved, we preserve the similar type, otherwise
    # just be similar to the first one.
    exprs = [name ∈ oldnames ? :(empty(getproperty(t, $(QuoteNode(name))), $ElType)) :
                               :(empty(first(columns(t)), $Eltype)) for (name, ElType) in zip(names, T.parameters)]
    return :(Table(NamedTuple{names}(tuple($(exprs...)))))
end

# Support Vector / deque interface (mutable-length vectors)

function Base.empty!(t::Table)
    map(empty!, columns(t))
    return t
end

function Base.pop!(t::Table)
    return map(pop!, columns(t))
end

function Base.push!(t::Table{<:NamedTuple{names}}, v::NamedTuple{names}) where {names}
    map(push!, columns(t), v)
    return t
end

function Base.append!(t::Table{<:NamedTuple{names}}, t2::Table{<:NamedTuple{names}}) where {names}
    map(append!, columns(t), columns(t2))
    return t
end

function Base.popfirst!(t::Table)
    return map(popfirst!, columns(t))
end

function Base.pushfirst!(t::Table{<:NamedTuple{names}}, v::NamedTuple{names}) where {names}
    map(pushfirst!, columns(t), v)
    return t
end

function Base.prepend!(t::Table{<:NamedTuple{names}}, t2::Table{<:NamedTuple{names}}) where {names}
    map(prepend!, columns(t), columns(t2))
    return t
end

function Base.deleteat!(t::Table, i)
    map(col -> deleteat!(col, i), columns(t))
    return t
end

function Base.insert!(t::Table{<:NamedTuple{names}}, i::Integer, v::NamedTuple{names}) where {names}
    map((col, val) -> insert!(col, i, val), columns(t), v)
    return t
end

function Base.splice!(t::Table, inds::Integer)
    return map(col -> splice!(col, inds), columns(t))
end

function Base.splice!(t::Table, inds::AbstractArray)
    return Table(map(col -> splice!(col, inds), columns(t)))
end

function Base.splice!(t::Table{<:NamedTuple{names}}, inds::Integer, ins::NamedTuple{names}) where {names}
    return map((col, vals) -> splice!(col, inds, vals), columns(t), ins)
end

function Base.splice!(t::Table{<:NamedTuple{names}}, inds::AbstractArray, ins::NamedTuple{names}) where {names}
    return Table(map((col, vals) -> splice!(col, inds, vals), columns(t), ins))
end

function Base.splice!(t::Table{<:NamedTuple{names}}, inds::Integer, ins::Table{<:NamedTuple{names}}) where {names}
    return map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins))
end

function Base.splice!(t::Table{<:NamedTuple{names}}, inds::AbstractArray, ins::Table{<:NamedTuple{names}}) where {names}
    return Table(map((col, vals) -> splice!(col, inds, vals), columns(t), columns(ins)))
end

# TODO splicing in an `AbstractArray{<:NamedTuple}` should be possible...

# Speedups for column-based storage

function Base.getindex(t::Table, inds::Union{AbstractArray, Colon}...)
    return Table(map(col -> getindex(col, inds...), columns(t)))
end

function Base.view(t::Table, inds::Union{AbstractArray, Colon}...)
    return Table(map(col -> view(col, inds...), columns(t)))
end

# Deprecated for .= syntax (via Base.Broadcast.materialize!)
# It seems `Ref` might be the new cool here. Could also consider `AbstractArray{<:NamedTuple, 0}`?
#function Base.setindex!(t::Table{<:NamedTuple{names}}, v::NamedTuple{names}, inds::Union{AbstractArray, Colon}...) where {names}
#    map((col, val) -> setindex!(col, val, inds...), columns(t), v)
#    return t
#end

function Base.setindex!(t::Table{<:NamedTuple{names}}, t2::Table{<:NamedTuple{names}}, inds::Union{AbstractArray, Colon}...) where {names}
    map((col, col2) -> setindex!(col, col2, inds...), columns(t), columns(t2))
    return t
end

function Base.vcat(t::Table{<:NamedTuple{names}}, t2::Table{<:NamedTuple{names}}) where {names}
    return Table(map(vcat, columns(t), columns(t2)))
end

function Base.hcat(t::Table{<:NamedTuple{names}}, t2::Table{<:NamedTuple{names}}) where {names}
    return Table(map(hcat, columns(t), columns(t2)))
end

function Base.hvcat(rows::Tuple{Vararg{Int}}, tables::Table{<:NamedTuple{names}}...) where {names}
    return Table(map((cols...,) -> hvcat(rows, cols...), map(columns, tables)...))
end

function Base.resize!(t::Table, len::Int)
	for c in columns(t)
		resize!(c, len)
	end
    t
end


text:
@testset "parameters.jl" begin
    @test χParam.allowed_min == -Inf
    @test χParam.allowed_max == Inf
    @test NParam.allowed_min == 0
    @test NParam.allowed_max == Inf
    @test fParam.allowed_min == 0.0
    @test fParam.allowed_max == 1.0
end


text:
function invertStreamfunction(
    lats :: AbstractArray{Float64, 1}, # in rad
    levs :: AbstractArray{Float64, 1},
    v    :: AbstractArray{Float64, 2};
    g0   :: AbstractFloat = 9.8,
    Re   :: AbstractFloat = 6371000.0,
)

    if (length(lats), length(levs)-1) != size(v)
        throw(ErrorException("Wrong dimension"))
    end

    dm = (levs[2:end] - levs[1:end-1]) / g0
    ψ  = zeros(Float64, length(lats), length(levs))

    for i = 1:length(lats)
        for k = size(v)[2]-1:-1:1
            ψ[i, k] = ψ[i, k+1] - v[i, k] * cos(lats[i]) * dm[k]
        end
    end

    ψ .*= 2.0*π * Re
    return ψ
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
@testset "Length" begin
    for len in [1, 16, 32, 64, 128]
        @test length(DNAKmer(random_dna_kmer(len))) == len
        @test length(RNAKmer(random_rna_kmer(len))) == len
        @test length(AAKmer(random_aa(len))) == len
    end
end


text:
"""
    posterior = bpca(X, Σ; [maxiter=1000], [S=300])

`X` is a D×N matrix of flux observations.

`σ` is a D×N matrix of error measurements.

Arguments in brackets are optional.
Default number of iterations is 1200.
Default number of samples for automatic variational inference is 300.

Returns posterior distribution of PCA parameters returned as a ```Distributions.MvNormal``` type.
"""
function bpca(X, σ; maxiter=1000, S=300)

    Σ = σ.^2

    #---------------------------------------------
    # Set constants and report message
    #---------------------------------------------

    Q = 1

    JITTER = 1e-8

    D, N = size(X)

    @printf("BPCA: %d number of data items of dimension %d, projected to %d components\n", N, D, Q)

    #---------------------------------------------
    function logPrior(W, μ)
    #---------------------------------------------

        -1e-10*(sum(W.^2) + sum(μ.^2))

    end

    #---------------------------------------------
    function marginalLogLikelihood(W, μ)
    #---------------------------------------------

       local ℓ = zero(eltype(W))

       local WWT = W*W'

       @inbounds for n in 1:N

           # check eqs. (12.36), (12.43) in Bishop
           local C = WWT + Diagonal(Σ[:,n]) + JITTER*I

           ℓ +=  - 0.5*logdet(C) - 0.5*sum((X[:,n].-μ)'*(C\(X[:,n].-μ)))

       end

       return ℓ

    end

    #---------------------------------------------
    function unpack(param)
    #---------------------------------------------

        local MARK = 0

        local W = reshape(param[MARK+1:MARK+D*Q], D, Q)

        MARK += D*Q

        local μ = reshape(param[MARK+1:MARK+D], D)

        MARK += D

        @assert(MARK == length(param))

        return W, μ

    end

    #---------------------------------------------
    function objective(param)
    #---------------------------------------------

        @assert(length(param) == D*Q + D)

        local W, μ = unpack(param)

        return marginalLogLikelihood(W, μ) + logPrior(W, μ)

    end

    #---------------------------------------------
    # Initialise with ppca estimates
    #---------------------------------------------

    function initsol()
        local indices = map(x -> ceil(Int,x), rand(N)*N)
        local W, μ, = @suppress ppca(X[:,indices])
        [vec(W);μ]
    end

    #---------------------------------------------
    # Run optimiser and return posterior
    #---------------------------------------------

    posterior, logevidence = VI(objective, [initsol() for i=1:5], S = S, Stest = 0*S, optimiser=LBFGS(), inititerations=10, iterations = maxiter, show_every = 10)

    return posterior

end


text:
__precompile__()
module Metaheuristics

export WOA, GOA, GSA, CGSA, SA, CMAES_AEP
export CFO, DE, ABC, PSO

# v2 items
export optimize, ECA, Options, State, Information, Engine, Problem
export MOEAD_DE

# v2.1 items
export convergence, minimizer, minimum, positions, fvals, nfes, get_position, fval, NSGA2
export hval, gval, hvals, gvals

export NSGA3, gen_ref_dirs

export SMS_EMOA, SPEA2, MCCGA

export PerformanceIndicators, pareto_front, nadir, ideal

include("externals.jl")

include("core/abstracts.jl")
include("core/information.jl")
include("core/options.jl")
include("core/state.jl")
include("core/structures.jl")

include("solutions/individual.jl")
include("solutions/constrained.jl")
include("solutions/display.jl")

include("common/multi-objective-functions.jl")
include("common/repair.jl")

include("common/stop.jl")
include("common/compare.jl")
include("common/crossover.jl")
include("common/mutation.jl")
include("common/selection.jl")
include("common/generic-reproduction.jl")
include("common/non-dominated-sorting.jl")

include("TestProblems/TestProblems.jl")
include("PerformanceIndicators/PerformanceIndicators.jl")

include("algorithms/algorithm.jl")

include("optimize.jl")

#######################################################
#                                                     #
#           A L G O R I T H M S                       #
#                                                     #
#######################################################

# template file
include("algorithms/template.jl")

# ECA algorithm
include("algorithms/ECA/ECA.jl")
include("algorithms/ECA/CECA.jl")

# Differential algorithm
include("algorithms/DE/DE.jl")

# PSO algorithm
include("algorithms/PSO/PSO.jl")

# The whale optimization algorithm
# S Mirjalili, A Lewis - Advances in Engineering Software, 2016
include("algorithms/WOA/WOA.jl")

# Mirjalili, Seyedali, and Amir H. Gandomi.
# "Chaotic gravitational constants for the gravitational search algorithm."
# Applied Soft Computing 53 (2017): 407-419.
include("algorithms/CGSA/CGSA.jl")

# SA: Simulated Annealing
# Kirkpatrick, S., Gelatt, C.D., & Vecchi, M.P. (1983). Optimization by
# Simulated Annealing. _Science, 220_, 671-680.
include("algorithms/SA/SA.jl")

# Aritifical Bee colony
include("algorithms/ABC/ABC.jl")

# MOEA decomposition based using Differential Evolution
include("algorithms/MOEAD_DE/MOEAD_DE.jl")

# Non-dominate sorting Genetic Algorithm
include("algorithms/NSGA2/NSGA2.jl")
include("algorithms/NSGA3/NSGA3.jl")
include("algorithms/SMS_EMOA/SMS_EMOA.jl")
include("algorithms/SPEA2/SPEA2.jl")

# genetic algorithm
include("algorithms/MCCGA/MCCGA.jl")

include("common/reproduction.jl")

#######################################################
#                                                     #
#           D E P R E C A T E D                       #
#                                                     #
#######################################################
# nothing to deprecate yet

end # module


text:
using .Tracker

(ev::Eval)(A::Tracker.TrackedArray, args...) = Tracker.track(ev, A, args...)
(ev::Eval)(A, B::Tracker.TrackedArray, args...) = Tracker.track(ev, A, B, args...)
(ev::Eval)(A::Tracker.TrackedArray, B::Tracker.TrackedArray, args...) = Tracker.track(ev, A, B, args...)

Tracker.@grad function (ev::Eval)(args...)
    Z = ev.fwd(Tracker.data.(args)...)
    Z, Δ -> begin
        isnothing(ev.rev) && error("no gradient definition here!")
        tuple(ev.rev(Tracker.data(Δ), Z, Tracker.data.(args)...)...)
    end
end


text:
working_directory = pwd()
println("The current working directory is", working_directory)

print(2+2)
"Hello", "World"

# Use the julia package tool
using Pkg

# Show what packages are currently installed
print(Pkg.installed())

print("Hello, World!")

print(3^2)
typeof(2)
typeof(5.23)
typeof("Howdy, Ya'll!")
typeof("23")

# How many secords are there in 42 minutes and 42 seconds
42*60 + 42

# How many miles are there in 10 kilometers
10 / 1.61

# To get pi, to \pi then TAB to autocomplete
π

x = 5
3x

parse(Int64, "32")
parse(Float64, "3.32")
string(45)
sin(45)
log10(123)

degrees = 45
x = sin(degrees / 360 * 2 * π)
x = exp(log(x+1))

sqrt(5)

function voidfunc()
    println("printing is a liar")
end

print(typeof(voidfunc()))

max_stars_repo_name: ryanorsinger/JuliaPlayground

max_issues_repo_name: ryanorsinger/JuliaPlayground

max_forks_repo_name: ryanorsinger/JuliaPlayground


text:
# Cairo backend for compose
import Cairo

import Cairo.CairoContext, Cairo.CairoSurface, Cairo.CairoARGBSurface,
       Cairo.CairoEPSSurface, Cairo.CairoPDFSurface, Cairo.CairoSVGSurface

export PNG, PDF, PS

abstract ImageBackend
abstract PNGBackend <: ImageBackend

abstract VectorImageBackend <: ImageBackend
abstract SVGBackend <: VectorImageBackend
abstract PDFBackend <: VectorImageBackend
abstract PSBackend  <: VectorImageBackend

# Native unit
type ImageMeasure{T <: ImageBackend} <: NativeMeasure
    value::Float64
end

copy{T}(u::ImageMeasure{T}) = ImageMeasure{T}(u.value)

convert(::Type{Float64}, u::ImageMeasure) = u.value
function convert{T}(::Type{ImageMeasure{T}}, u::Number)
    ImageMeasure{T}(convert(Float64, u))
end

function *{T}(u::Number, v::ImageMeasure{T})
    ImageMeasure{T}(convert(Float64, u) * v.value)
end

function /{T}(u::ImageMeasure{T}, v::Number)
    ImageMeasure{T}(u.value / convert(Float64, v))
end

function +{T}(u::ImageMeasure{T}, v::ImageMeasure{T})
    ImageMeasure{T}(u.value + v.value)
end

function -{T}(u::ImageMeasure{T}, v::ImageMeasure{T})
    ImageMeasure{T}(u.value - v.value)
end

type ImagePropertyState
    stroke::ColorOrNothing
    fill::ColorOrNothing
end

type Image{B <: ImageBackend} <: Backend
    out::IO
    surface::CairoSurface
    ctx::CairoContext

    # Current state
    stroke::ColorOrNothing
    fill::ColorOrNothing
    opacity::Float64 # in [0,1]
    visible::Bool

    # Keep track of property
    state_stack::Vector{ImagePropertyState}

    # Close the surface when finished
    owns_surface::Bool

    # Close the stream when finished
    close_stream::Bool

    # Emit on finish
    emit_on_finish::Bool

    function Image(surface::CairoSurface, ctx::CairoContext, out::IO)
        img = new()
        img.out = out
        img.surface = surface
        img.ctx = ctx
        img.stroke = RGB(0., 0., 0.)
        img.fill   = RGB(0., 0., 0.)
        img.opacity = 1.0
        img.visible = true
        img.state_stack = Array(ImagePropertyState, 0)
        img.owns_surface = false
        img.emit_on_finish = false
        img
    end

    function Image(surface::CairoSurface, ctx::CairoContext)
        Image{B}(surface, ctx, IOString())
    end

    Image(surface::CairoSurface) = Image{B}(surface, CairoContext(surface))

    function Image(out::IO,
                   width::MeasureOrNumber,
                   height::MeasureOrNumber,
                   emit_on_finish::Bool=true)

        w = native_measure(width, B).value
        h = native_measure(height, B).value

        local surface::CairoSurface
        if B == SVGBackend
            surface = CairoSVGSurface(out, w, h)
        elseif B == PNGBackend
            surface = CairoARGBSurface(round(w), round(h))
        elseif B == PDFBackend
            surface = CairoPDFSurface(out, w, h)
        elseif B == PSBackend
            surface = CairoEPSSurface(out, w, h)
        else
            error("Unkown Cairo backend.")
        end

        if Cairo.status(surface) != Cairo.CAIRO_STATUS_SUCCESS
            error("Unable to create cairo surface.")
        end

        img = Image{B}(surface, CairoContext(surface), out)
        img.owns_surface = true
        img.close_stream = false
        img.emit_on_finish = emit_on_finish
        img
    end

    function Image(filename::String,
                   width::MeasureOrNumber,
                   height::MeasureOrNumber)
        img = Image{B}(open(filename, "w"), width, height)
        img.close_stream = true
        img
    end
end

width{B}(img::Image{B}) = ImageMeasure{B}(Cairo.width(img.surface))
height{B}(img::Image{B}) = ImageMeasure{B}(Cairo.height(img.surface))

finish{B <: ImageBackend}(::Type{B}, img::Image) = nothing

function finish(::Type{PNGBackend}, img::Image)
    Cairo.write_to_png(img.surface, img.out)
end

function finish{B <: ImageBackend}(img::Image{B})
    if(img.owns_surface)
        Cairo.destroy(img.ctx)
    end
    finish(B, img)
    if(img.owns_surface)
        Cairo.destroy(img.surface)
    end

    if img.emit_on_finish && typeof(img.out) == IOString
        seek(img.out, 0)
        if B == PNGBackend
            mime = "image/png"
        elseif B == SVGBackend
            mime = "image/svg+xml"
        elseif B == "PSBackend"
            mime = "application/postscript"
        elseif B == "PDFBackend"
            mime = "application/pdf"
        end
        emit(Emitable(mime, readall(img.f)))
        close(img.out)
    end
end

typealias PNG Image{PNGBackend}
typealias PDF Image{PDFBackend}
typealias PS  Image{PSBackend}

# sizes

function root_box{B}(img::Image{B})
    NativeBoundingBox(
        ImageMeasure{B}(0.),
        ImageMeasure{B}(0.),
        width(img),
        height(img))
end

native_zero{T}(backend::Image{T}) = ImageMeasure{T}(0.0)

# PNG conversion to native units (i.e., pixels)

function native_measure{K <: ImageBackend}(u::Number, ::Type{K})
    ImageMeasure{K}(convert(Float64, u))
end

function native_measure{K <: ImageBackend}(u::SimpleMeasure{PixelUnit},
                                           ::Type{K})
    ImageMeasure{K}(u.value)
end

function native_measure{K <: ImageBackend}(u::SimpleMeasure{MillimeterUnit},
                                           ::Type{K})
    native_measure(convert(SimpleMeasure{PixelUnit}, u), K)
end

function native_measure{K <: ImageBackend}(u::Number, backend::Image{K})
    native_measure(u, K)
end

function native_measure{K <: ImageBackend}(u::SimpleMeasure{PixelUnit},
                        backend::Image{K})
    native_measure(u, K)
end

function native_measure{K <: ImageBackend}(u::SimpleMeasure{MillimeterUnit},
                        backend::Image{K})
    native_measure(u, K)
end

# SVG/PDF/PS conversion to native units (i.e., pts)

function native_measure{K <: VectorImageBackend}(u::Number, ::Type{K})
    ImageMeasure{K}(convert(Float64, u))
end

function native_measure{K <: VectorImageBackend}(u::SimpleMeasure{PixelUnit},
                                                 ::Type{K})
    native_measure(convert(SimpleMeasure{MillimeterUnit}, u), K)
end

function native_measure{K <: VectorImageBackend}(u::SimpleMeasure{MillimeterUnit},
                                                 ::Type{K})
    ImageMeasure{K}(u / pt)
end

function native_measure{K <: VectorImageBackend}(
        u::Number,
        backend::Image{K})
    native_measure(u, K)
end

function native_measure{K <: VectorImageBackend}(
        u::SimpleMeasure{PixelUnit},
        backend::Image{K})
    native_measure(u, K)
end

function native_measure{K <: VectorImageBackend}(
        u::SimpleMeasure{MillimeterUnit},
        backend::Image{K})
    native_measure(u, K)
end

# Drawing

move_to(img::Image, point::Point) = Cairo.move_to(img.ctx, point.x.value, point.y.value)
line_to(img::Image, point::Point) = Cairo.line_to(img.ctx, point.x.value, point.y.value)

# TODO: What is the equivalent Cairo call?
curve_to(img::Image, ctrl0::Point, ctrl1::Point, anchor1::Point) =
    curve_to(img.ctx, ctrl0.x.value, ctrl0.y.value, ctrl1.x.value, ctrl1.y.value,
            anchor1.x.value, anchor1.y.value)

close_path(img::Image) = Cairo.close_path(img.ctx)

function arc(img::Image, x::Float64, y::Float64, radius::Float64,
             angle1::Float64, angle2::Float64)
    Cairo.arc(img.ctx,x,y,radius,angle1,angle2)
end

translate(img::Image, tx::Float64, ty::Float64) = Cairo.translate(img.ctx, tx, ty)
scale(img::Image, sx::Float64, sy::Float64) = Cairo.scale(img.ctx, sx, sy)
rotate(img::Image, theta::Float64) = Cairo.rotate(img.ctx, theta)

function fillstroke(img::Image)
    if img.fill != nothing && img.opacity > 0.0 && img.visible
        rgb = convert(RGB, img.fill)
        Cairo.set_source_rgba(img.ctx, rgb.r, rgb.g, rgb.b, img.opacity)

        if img.stroke != nothing
            Cairo.fill_preserve(img.ctx)
        else
            Cairo.fill(img.ctx)
        end
    end

    if img.stroke != nothing
        rgb = convert(RGB, img.stroke)
        Cairo.set_source_rgb(img.ctx, rgb.r, rgb.g, rgb.b)

        Cairo.stroke(img.ctx)
    end
end

function save_state(img::Image)
    push!(img.state_stack, ImagePropertyState(img.stroke, img.fill))
    Cairo.save(img.ctx)
end

function restore_state(img::Image)
    state = pop!(img.state_stack)
    img.stroke = state.stroke
    img.fill = state.fill
    Cairo.restore(img.ctx)
end

function draw(img::Image, form::Lines)
    if isempty(form.points); return; end

    paths = make_paths(form.points)
    for path in paths
        move_to(img, path[1])
        for point in path[2:]
            line_to(img, point)
        end
        fillstroke(img)
    end
end

function draw(img::Image, form::Curve)
    move_to(img, form.anchor0)
    curve_to(img, form.ctrl0, form.ctrl1, form.anchor1)
end

function draw(img::Image, form::Polygon)
    if isempty(form.points); return; end

    move_to(img, form.points[1])
    for point in form.points[2:]
        line_to(img, point)
    end
    close_path(img)
    fillstroke(img)
end

function draw(img::Image, form::Ellipse)
    cx = form.center.x.value
    cy = form.center.y.value
    rx = sqrt((form.x_point.x.value - cx)^2 +
              (form.x_point.y.value - cy)^2)
    ry = sqrt((form.y_point.x.value - cx)^2 +
              (form.y_point.y.value - cy)^2)
    theta = atan2(form.x_point.y.value - cy,
                  form.x_point.x.value - cx)

    save_state(img)
    translate(img, cx, cy)
    rotate(img, theta)
    translate(img, -rx, -ry)
    scale(img, 2rx, 2ry)
    arc(img, 0.5, 0.5, 0.5, 0.0, 2pi)
    restore_state(img)
    fillstroke(img)
end

function draw(img::Image, form::Text)
    if !img.visible || img.opacity == 0.0 || img.fill === nothing
        return
    end

    pos = copy(form.pos)

    if form.halign != hleft || form.valign != vtop
        extents = Array(Float64, 6)
        Cairo.text_extents(img.ctx, bytestring(form.value), extents)

        width, height = extents[3], extents[4]

        if form.halign == hcenter
            pos.x.value -= width / 2
        elseif form.halign == hright
            pos.x.value -= width
        end

        if form.valign == vcenter
            pos.y.value += height / 2
        elseif form.valign == vtop
            pos.y.value += height
        end
    end

    rgb = convert(RGB, img.fill)
    Cairo.set_source_rgba(img.ctx, rgb.r, rgb.g, rgb.b, img.opacity)

    save_state(img)
    Cairo.translate(img.ctx, form.t.M[1,3], form.t.M[2,3])
    rotate(img, atan2(form.t.M[2,1], form.t.M[1,1]))
    move_to(img, pos)
    Cairo.show_text(img.ctx, form.value)
    restore_state(img)
end

# Applying properties

function push_property(img::Image, p::Property)
    save_state(img)
    while !is(p, empty_property)
        apply_property(img, p.primitive)
        p = p.next
    end
end

function pop_property(img::Image)
    restore_state(img)
end

# Nop catchall
function apply_property(img::Image, p::PropertyPrimitive)
end

function apply_property(img::Image, p::Stroke)
    img.stroke = p.value
end

function apply_property(img::Image, p::Fill)
    img.fill = p.value
end

function apply_property(img::Image, p::Opacity)
    img.opacity = p.value
end

function apply_property(img::Image, p::Visible)
    img.visible = p.value
end

function apply_property(img::Image, property::LineWidth)
    Cairo.set_line_width(img.ctx, property.value.value)
end

function apply_property(img::Image, property::Font)
    Cairo.select_font_face(img.ctx, property.family,
                           Cairo.CAIRO_FONT_SLANT_NORMAL,
                           Cairo.CAIRO_FONT_WEIGHT_NORMAL)
end

function apply_property(img::Image, property::FontSize)
    Cairo.set_font_size(img.ctx, property.value.value)
end

function apply_property(img::Image, property::Clip)
    if isempty(property.points); return; end
    move_to(img, property.points[1])
    for point in property.points[2:]
        line_to(img, point)
    end
    close_path(img)
    Cairo.clip(img.ctx)
end


text:
@testset "trades" begin

  @testset "trades" begin
    trades = CoinbasePro.trades("BTC-USD")
    @test length(trades) == 2
    @test isa(trades[1], DataFrame)
    @test length(trades[2]) == 2
  end

end


text:
using Thebes, Luxor

function main()

    Drawing(400, 400, "test2-pyramids.png")
    origin()
    background("white")
    setopacity(1)

    eyepoint(Point3D(-200, 150, 20))
    centerpoint(Point3D(0, 0, 10))
    uppoint(Point3D(0, 0, 20)) # relative to centerpoint
    perspective(500)

# The northernmost and oldest pyramid of the group was built
# for Khufu (Greek: Cheops), the second king of the 4th
# dynasty. Called the Great Pyramid, it is the largest of
# the three, the length of each side at the base averaging
# 755.75 feet (230 metres) and its original height being
# 481.4 feet (147 metres).

# The middle pyramid was built for Khafre (Greek:
# Chephren), the fourth of the eight kings of the 4th
# dynasty; the structure measures 707.75 feet (216
# metres) on each side and was originally 471 feet (143
# metres) high.

# The southernmost and last pyramid to be built was that
# of Menkaure (Greek: Mykerinus), the fifth king of the
# 4th dynasty; each side measures 356.5 feet (109
# metres), and the structure’s completed height was 218
# feet (66 metres).

    carpet(80)
    sethue("gold")

    khufu = make(Pyramid)
    khafre = make(Pyramid)
    menkaure = make(Pyramid)

    pyramids = [menkaure, khafre, khufu]

    scaleby!(khufu,   23.0/2, 23.0/2, 14.7)
    scaleby!(khafre,  21.6/2, 21.6/2, 14.3)
    scaleby!(menkaure, 10.9/2, 10.9/2, 6.6)

    moveby!(khufu,     30, 30, 0)
    moveby!(khafre,   -0, -0, 0)
    moveby!(menkaure, -30, -30, 0)

    # needs depth sorting of faces
    sortfaces!.(pyramids, eyepoint=eyepoint())
    pin(menkaure)
    pin(khafre)
    pin(khufu)
    finish()
end

main()


text:
struct Basis
    base_var_constrs::Vector{MOI.BasisStatusCode}
    base_lt_constrs::Vector{MOI.BasisStatusCode}
    base_gt_constrs::Vector{MOI.BasisStatusCode}
    base_et_constrs::Vector{MOI.BasisStatusCode}
    branch_lt_constrs::Vector{MOI.BasisStatusCode}
    branch_gt_constrs::Vector{MOI.BasisStatusCode}
end
Basis() = Basis([], [], [], [], [], [])

function Base.copy(src::Basis)
    return Basis(
        copy(src.base_var_constrs),
        copy(src.base_lt_constrs),
        copy(src.base_gt_constrs),
        copy(src.base_et_constrs),
        copy(src.branch_lt_constrs),
        copy(src.branch_gt_constrs),
    )
end

mutable struct PollingState
    next_polling_target_time_sec::Float64
    period_node_count::Int
    period_simplex_iters::Int
end
PollingState() = PollingState(0.0, 0, 0)

mutable struct BaseConstraintState
    var_constrs::Vector{CI{SV,IN}}
    lt_constrs::Vector{CI{SAF,LT}}
    gt_constrs::Vector{CI{SAF,GT}}
    et_constrs::Vector{CI{SAF,ET}}
end
BaseConstraintState() = BaseConstraintState([], [], [], [])

function Base.empty!(cs::BaseConstraintState)
    empty!(cs.var_constrs)
    empty!(cs.lt_constrs)
    empty!(cs.gt_constrs)
    empty!(cs.et_constrs)
    return nothing
end

mutable struct BranchConstraintState
    num_lt_branches::Int
    num_gt_branches::Int
    lt_general_constrs::Vector{CI{SAF,LT}}
    gt_general_constrs::Vector{CI{SAF,GT}}
end
BranchConstraintState() = BranchConstraintState(0, 0, [], [])

function Base.empty!(cs::BranchConstraintState)
    cs.num_lt_branches = 0
    cs.num_gt_branches = 0
    empty!(cs.lt_general_constrs)
    empty!(cs.gt_general_constrs)
    return nothing
end

mutable struct ConstraintState
    base_state::BaseConstraintState
    branch_state::BranchConstraintState
end
function ConstraintState()
    return ConstraintState(BaseConstraintState(), BranchConstraintState())
end

function Base.empty!(cs::ConstraintState)
    empty!(cs.base_state)
    empty!(cs.branch_state)
    return nothing
end

mutable struct CurrentState
    gurobi_env::Gurobi.Env
    gurobi_model::Gurobi.Optimizer
    rebuild_model::Bool
    tree::Tree
    on_a_dive::Bool
    # TODO: Hashing nodes might be more expensive than we'd like. Instead, just
    # attach an ID to each node and use Ints as keys.
    warm_starts::Dict{Node,Basis}
    primal_bound::Float64
    dual_bound::Float64
    current_solution::Vector{Float64}
    best_solution::Vector{Float64}
    starting_time::Float64
    total_elapsed_time_sec::Float64
    total_node_count::Int
    total_simplex_iters::Int
    total_model_builds::Int
    total_warm_starts::Int
    _variable_indices::Vector{VI}
    constraint_state::ConstraintState
    disjunction_state::Dict{DisjunctiveFormulater,Any}
    polling_state::PollingState

    function CurrentState(; primal_bound::Real = Inf)
        state = new()
        state.gurobi_env = Gurobi.Env()
        state.on_a_dive = true
        # gurobi_model is left undefined; build it before accessing.
        state.rebuild_model = true
        state.tree = Tree()
        push_node!(state.tree, Node())
        state.warm_starts = Dict{Node,Basis}()
        state.primal_bound = primal_bound
        state.dual_bound = -Inf
        state.starting_time = time()
        state.current_solution = Float64[]
        state.best_solution = Float64[]
        state.total_elapsed_time_sec = 0.0
        state.total_node_count = 0
        state.total_simplex_iters = 0
        state.total_model_builds = 0
        state.total_warm_starts = 0
        state._variable_indices = VI[]
        state.constraint_state = ConstraintState()
        state.disjunction_state = Dict()
        state.polling_state = PollingState()
        return state
    end
end

function reset_formulation_state!(state::CurrentState)
    empty!(state._variable_indices)
    empty!(state.constraint_state)
    empty!(state.disjunction_state)
    return nothing
end

function update_dual_bound!(state::CurrentState)
    if isempty(state.tree)
        # Tree is exhausted, so have proven optimality of best found solution.
        state.dual_bound = state.primal_bound
    else
        state.dual_bound =
            minimum(node.dual_bound for node in state.tree.open_nodes)
    end
    return nothing
end

function instantiate(cvi::CVI, state::CurrentState)
    return state._variable_indices[index(cvi)]
end

function instantiate(csaf::CSAF, state::CurrentState)
    return SAF(
        [
            SAT(coeff, instantiate(cvi, state)) for
            (coeff, cvi) in zip(csaf.coeffs, csaf.indices)
        ],
        csaf.constant,
    )
end

function attach_index!(state::CurrentState, vi::VI)
    push!(state._variable_indices, vi)
    return CVI(length(state._variable_indices))
end


text:
export SisPosteriorSampler, SisIexInsertDelete, SisIexInsertDeleteGibbs

abstract type SisPosteriorSampler end

struct SisIexInsertDeleteGibbs <: SisPosteriorSampler
    ν::Int
    β::Float64
    path_dist::PathDistribution
    K::Int
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions in Gibbs scan
    ind_add::Vector{Int} # Storage for indexing of additions in Gibbs scan
    vals::Vector{Int} # Storage for valuse to insert in Gibbs scan
    function SisIexInsertDeleteGibbs(
        path_dist::PathDistribution;
        K=100,
        ν=4, β=0.6,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν)
        ind_add = zeros(Int, ν)
        vals = zeros(Int, ν)
        par_info = Dict()
        par_info[:ν] = "(maximum number of edit operations in iMCMC-within-Gibbs conditional updates)"
        par_info[:path_dist] = "(path distribution for insertions)"
        par_info[:β] = "(probability of Gibbs scan)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν, β, path_dist, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals
            )
    end
end

struct SisIexInsertDelete<: SisPosteriorSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum increase/decrease in dimension
    β::Real  # Probability of trans-dimensional move
    α::Float64  # How much vertex proposal is informed by data (0.0 max and ∞ is uniform over vertice)
    len_dist::DiscreteUnivariateDistribution  # Distribution used to introduce new interactions
    ε::Float64 # Neighborhood for sampling γ
    aux_mcmc::SisMcmcSampler
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values where are deleted from interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    function SisIexInsertDelete(
        aux_mcmc::SisMcmcSampler;
        K=100,
        ν_ed=2, ν_td=1 , β=0.7, α=0.0, len_dist=TrGeometric(0.8,1,K),
        ε=0.05,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum change in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:α] = "(controls how much data informs entry insertion proposals)"
        par_info[:aux_mcmc] = "(mcmc sampler from auxiliary data)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"
        par_info[:ε] = "(neighborhood for γ proposals)"

        new(
            ν_ed, ν_td, β, α,
            len_dist,
            ε,
            aux_mcmc, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td_add, ind_td_del
            )
    end
end


text:
"""
    ParameterizedForcing(func, parameters)

Construct a forcing function with parameters. The forcing function, which is
applied at grid point `i, j, k`, is called with the signature

    `func(i, j, k, grid, clock, state, parameters)`

Example
=======

function cool_forcing_function(i, j, k, grid, clock, state, parameters) =
    return @inbounds - parameters.μ * exp(grid.zC[k] / parameters.λ) * state.velocities.u[i, j, k]
end

parameters = (μ=42, λ=π)

cool_forcing = ParameterizedForcing(cool_forcing_function, parameters)
"""
struct ParameterizedForcing{F, P}
    func :: F
    parameters :: P
end

@inline (F::ParameterizedForcing)(i, j, k, grid, clock, state) =
    F.func(i, j, k, grid, clock, state, F.parameters)


text:
contactStart(matPair::Shapes.NoContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, contactNormal::Frames.Vector3D, elasticContactReductionFactor::Float64) =
matPair

contactStart(matPair::Shapes.ObserverContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, contactNormal::Frames.Vector3D, elasticContactReductionFactor::Float64) = matPair

contactStart(matPair::Shapes.ImpulseContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, contactNormal::Frames.Vector3D, elasticContactReductionFactor::Float64) =
error("contactStart is not yet implemented for ImpulseContactPairMaterial.")

contactStart(matPair::Shapes.WheelRailContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, contactNormal::Frames.Vector3D, elasticContactReductionFactor::Float64) =
error("contactStart is not yet implemented for WheelRailContactPairMaterial.")

responseCalculation(material::Shapes.NoContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, e_n::Frames.Vector3D, s::Float64, time::Float64, file, sim) =
(Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D)

function responseCalculation(material::Shapes.ObserverContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, e_n::Frames.Vector3D, s::Float64, time::Float64, file, sim)
  if material.printAlarm && ModiaLang.isEvent(sim)
    println("At time event ", time, " two observer objects ", Modia3D.fullName(obj1), " and ", Modia3D.fullName(obj2), " are colliding.")
  end
  return (Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D)
end

responseCalculation(material::Shapes.ImpulseContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, e_n::Frames.Vector3D, s::Float64, time::Float64, file, sim) =
(Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D)

responseCalculation(material::Shapes.WheelRailContactPairMaterial, obj1::Object3D, obj2::Object3D, rContact::Frames.Vector3D, e_n::Frames.Vector3D, s::Float64, time::Float64, file, sim) =
(Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D)

contactEnd(mat::Shapes.NoContactPairMaterial,obj1,obj2)::Nothing        = nothing
contactEnd(mat::Shapes.ObserverContactPairMaterial,obj1,obj2)::Nothing  = nothing
contactEnd(mat::Shapes.ImpulseContactPairMaterial,obj1,obj2)::Nothing   = nothing
contactEnd(mat::Shapes.WheelRailContactPairMaterial,obj1,obj2)::Nothing = nothing