text:
import Pkg
Pkg.activate(mktempdir())
Pkg.add(["HTTP","EzXML"])

using Downloads
using HTTP, EzXML
pyrhelio_dir = "/gpfs/group/ebf11/default/pipeline/data/neid_solar/pyrheliometer"
url = "https://neid.ipac.caltech.edu/pyrheliometer.php"

try

response = HTTP.request("GET",url)
body = String(response.body)
idx = findall(r"https\:\/\/\S+\/pyroheliotel\/neid_\w+_\S{4}_\d{8}\.tel",body)
for i in eachindex(idx)
   this_url = body[idx[i]]
   m = match(r"neid_l\S+\.tel",this_url)
   if isnothing(m) continue end
   filename = joinpath(pyrhelio_dir,string(m.match))
   if filesize(filename) > 0  continue  end
   println("# Need to download ", filename)
   Downloads.download(this_url, filename)
   if !(filesize(filename) > 0)
      println("# Somehow ", filename, " is still missing/empty.")
   end
end

catch ex
   @warn "Error: " * string(ex)
end


text:
module _TestIndirectImportsDownstream

using IndirectImports

@indirect import _TestIndirectImportsUpstream
@indirect _TestIndirectImportsUpstream.fun(x) = x + 1

dispatch(::typeof(_TestIndirectImportsUpstream.fun)) = :fun
dispatch(::typeof(_TestIndirectImportsUpstream.fun2)) = :fun2

# Test other importing syntax:
@indirect import _TestIndirectImportsUpstream: f1
@indirect import _TestIndirectImportsUpstream: f2, f3
@indirect import _TestIndirectImportsUpstream: f4, f5, f6

struct Config1 end
struct Config2 end

@indirect _TestIndirectImportsUpstream.op(::Config1, x, y) = x + y
@indirect _TestIndirectImportsUpstream.op(::Config2, x, y) = x - y

end # module

ext: jl


text:
######################################################################
# searchdomain.jl
# Contains search domain type
######################################################################
type SearchDomain
	cell_size::Int64		# side of one grid cell, in meters
	num_cells::Int64		# number of cells per side
	domain_size::Int64
	num_vehicles::Int64
	noise_sigma::Float64  # std deviation of noise, degrees
	d::Normal
	F						# fisher information stuff
	b::Belief
	theta::NTuple{2, Int}

	function SearchDomain(num_cells::Int64, theta_x::Int, theta_y::Int)
		m = new()
		m.cell_size = 1
		m.num_cells = num_cells
		m.domain_size = num_cells
		m.num_vehicles = 1
		m.noise_sigma = 15.0
		m.d = Normal(0, m.noise_sigma)
		m.F = fisher(m)
		m.b = ones(num_cells, num_cells) / (num_cells * num_cells)
		m.theta = (theta_x, theta_y)
		return m
	end
	function SearchDomain(num_cells::Int64)
		return SearchDomain(num_cells, num_cells-1, num_cells-1)
	end

end

# noise_sigma is in degrees
function set_noise(m::SearchDomain, noise_sigma::Float64)
	m.noise_sigma = noise_sigma
	m.d = Normal(0, m.noise_sigma)
end

"""
`initial_belief(m::SearchDomain)`
"""
initial_belief(m::SearchDomain) = ones(m.num_cells, m.num_cells) / (m.num_cells * m.num_cells)

"""
`update_belief(m::SearchDomain, X::VehicleSet, Z::ObsSet)`

Modifies the `b` field of `m`.
"""
function update_belief!(m::SearchDomain, X::VehicleSet, Z::ObsSet)

	# Loop over all possible jammer positions
	bp_sum = 0.0
	for theta_x = 1:m.num_cells
		for theta_y = 1:m.num_cells
			temp_val = m.b[theta_x, theta_y]
			for (i,xi) in enumerate(X)
				temp_val *= O(m, xi, (theta_x-1,theta_y-1), Z[i])
			end
			m.b[theta_x, theta_y] = temp_val
			bp_sum += temp_val
		end
	end

	# normalize
	for theta_x = 1:m.num_cells
		for theta_y = 1:m.num_cells
			m.b[theta_x, theta_y] /= bp_sum
		end
	end
end

# allows you to initialize vehicle set
# Assumes there is an odd number of cells
function Vehicle(m::SearchDomain)
	mid_cell = (m.num_cells-1) / 2
	return Vehicle(mid_cell, mid_cell)
end


text:
function lag_matrix(x::T, lag::Int) where {T<:AbstractMatrix}
    nobs, nx = size(x)
    xlag = zeros(nobs - lag, nx * lag)
    for i = 1:lag
        xlag[1:(nobs-lag), (nx*(i-1)+1):(nx*i)] = x[(lag+1-i):(nobs-i), 1:nx]
    end
    return xlag
end

function lag_matrix(x::T, lag::Int) where {T<:AbstractVector}
    x2 = reshape(x, length(x), 1)
    return lag_matrix(x2, lag)
end

"""
    companion(A, p) -> Matrix

Construct the companion matrix for VAR coefficient matrix `A`, where `A` has
    dimensions (Kp × K). K is the number of endogenous variables, p is the
    number of lags.
"""
function companion(A, p)
    K = size(A, 2)
    subI = [Matrix(1.0I, K * (p - 1), K * (p - 1)) zeros(K * (p - 1), K)]
    Acomp = vcat(A', subI)
    return Acomp
end

function rhs_matrix(y, ylag, constant, trend, x)
    rhs = lag_matrix(y, ylag)
    if !isa(x, Nothing)
        rhs = hcat(x[ylag+1:end, :], rhs)
    end
    if trend
        ltrend = collect(1.0:size(rhs, 1))
        rhs = hcat(ltrend, rhs)
    end
    if constant
        rhs = hcat(ones(Float64, size(rhs, 1)), rhs)
    end
    return rhs
end


text:
using AbstractTrees

# Replace module in Expr tree
function expr_replace_module(expr, pat::Pair{Module, Module})
    expr_replace_module!(deepcopy(expr), pat)
end
function expr_replace_module!(expr, pat::Pair{Module, Module})
    treemap!(PreOrderDFS(expr)) do node
        if node isa GlobalRef && node.mod == pat.first
            GlobalRef(pat.second, node.name)
        else
            node
        end
    end
end

# Replace Module in CodeInfo
function code_replace_module(code::Core.CodeInfo, pat::Pair{Module, Module})
    code2 = copy(code)  # I guess this is sufficient...
    code_replace_module!(code2, pat)
end
function code_replace_module!(code::Core.CodeInfo, pat::Pair{Module, Module})
    for e in code.code
        expr_replace_module!(e, pat)
    end
    code
end

# Deepcopy Function

compressed_ast(m::Method, ci::Core.CodeInfo) =
    ccall(:jl_compress_ast, Any, (Any, Any), m, ci)

function deepcopy_function(destmodule::Module, f::Function)
    name = nameof(f)
    f2 = try
            Core.eval(destmodule, name)
        catch
            @eval destmodule function $name end
        end
    src_ms = methods(f).ms
    for m in src_ms
        deepcopy_method(destmodule, m, f2)
    end
end

_func_sig_params(t) = t isa UnionAll ? _func_sig_params(t.body) : t.parameters[2:end]
_func_sig_typeparams(t) = t isa UnionAll ? (t.var, _func_sig_typeparams(t.body)...) : ()
function deepcopy_method(destmodule, m, func::Function)
    # TODO: is it okay to use m.sig.parameters
    argdata = Core.svec(Core.svec(typeof(func), _func_sig_params(m.sig)...),
                        Core.svec(_func_sig_typeparams(m.sig)...))

    ci = Base.uncompressed_ast(m)
    ci = code_replace_module!(ci, m.module => destmodule)

    ccall(:jl_method_def, Cvoid, (Core.SimpleVector, Any, Ptr{Module}), argdata, ci, pointer_from_objref(destmodule))
end


text:
__precompile__()
module Gallium

    using ASTInterpreter2
    using ASTInterpreter2: @enter

    export @enter

end


text:
# explore these commands in the REPL:
using HTTP

## Request:
req = HTTP.Request(
   "GET",                        # 1 - Could be GET, POST, UPDATE etc
   "http://localhost:8081/search",    # 2 - URL -
   ["Content-Type" => "text/plain"],  # 3 - Header fields -
   "Hi there!"                        # 4 - Payload/body
)

# 1
req.method           # "GET"
# 2
req.target           # "http://localhost:8081/search"
# 3
req.headers
# 1-element Vector{Pair{SubString{String}, SubString{String}}}:
#  "Content-Type" => "text/plain"
# or more specific:
req["Content-Type"]   # "text/plain"
# same as: HTTP.header(req, "Content-Type")
# 4
HTTP.payload(req)
# 9-element Vector{UInt8}:
#  0x48
#  0x69
#  0x20
#  0x74
#  0x68
#  0x65
#  0x72
#  0x65
#  0x21
String(HTTP.payload(req))
# "Hi there!"
splitpath(req.target)[2]#  returns "search"

## Response:
resp = HTTP.Response(
   200,                       # 1 - Status code, 200 means success.
   ["Content-Type" => "text/plain"], # 2 - Header fields -
   body="Hi there!"         # 3 - Payload/body
)

# 1:
resp.status   # 200
# 2:
resp.headers
# 1-element Vector{Pair{SubString{String}, SubString{String}}}:
#  "Content-Type" => "text/plain"
# or more specific:
resp["Content-Type"]   # "text/plain"
# same as: HTTP.header(resp, "Content-Type")
# 3:
String(resp.body)          # "Hi there!"
# or:
String(HTTP.payload(resp)) # "Hi there!"


text:
module TestReparametrizedBautin
include("preamble.jl")

using ForwardDiff
using Setfield: compose
using StaticArrays: SVector, SMatrix

using Bifurcations: resolved_points, reparametrize
using Bifurcations.Codim2: first_lyapunov_coefficient
using Bifurcations.Codim2LimitCycle: FoldLimitCycleProblem
using Bifurcations.Continuations: get_prob_cache, get_u0,
    residual!, residual_jacobian
using Bifurcations.Examples.Bautin
using Bifurcations.Examples.Reparametrization: orig_p

# Putting @info so that Travis won't stop testing.
@testset for opts in [
        [],
        [:extra => SVector(-1.0)],
        [:extra => SVector(-1.0), :shift => 0.1, :seed => 1],
        ]
    @info "TestReparametrizedBautin:"
    @info "opts = $opts"

    prob = reparametrize(
        Bautin.make_prob(),
        ; opts...)
    solver1 = init(prob)
    solve!(solver1)

    codim1_points = resolved_points(solver1)
    @test length(codim1_points) == 1
    @test codim1_points[1].point_type === Codim1.PointTypes.hopf
    @test codim1_points[1].u[end] ≈ 0  atol=1e-6

    @info "Continuation of Hopf bifurcation..."
    hopf_prob = BifurcationProblem(
        codim1_points[1],
        solver1,
        compose(orig_p, (@lens _.β₂)),
        (-2.0, 2.0),
    )
    hopf_solver = init(
        hopf_prob;
        verbose = true,
    )
    @time solve!(hopf_solver)

    hopf_β₁ = [u[end-1] for sweep in hopf_solver.super.sol.sweeps for u in sweep.u]
    hopf_β₂ = [u[end]   for sweep in hopf_solver.super.sol.sweeps for u in sweep.u]
    @test all(@. abs(hopf_β₁) < 1e-6)
    @test maximum(hopf_β₂) > 2
    @test minimum(hopf_β₂) < -2

    l1s = [
        first_lyapunov_coefficient(hopf_solver.cache.super.prob_cache, u)
        for sweep in hopf_solver.super.sol.sweeps for u in sweep.u
    ]
    if !all(@. (l1s > 0) == (hopf_β₂ > 0))
        println("[hopf_β₂ l1s]' =")
        display([hopf_β₂ l1s]')
    end
    @test all(@. (l1s > 0) == (hopf_β₂ > 0))

    @info "Resolving Bautin point..."
    @time codim2_points = resolved_points(hopf_solver)
    @test length(codim2_points) == 1
    @test codim2_points[1].point_type === Codim2.PointTypes.bautin
    β_bautin = codim2_points[1].u[end-1:end]
    @test all(@. abs(β_bautin) < 1e-6)

    @info "Creating FoldLimitCycleProblem..."
    @time flc_prob = FoldLimitCycleProblem(
        codim2_points[1],
        hopf_solver;
        num_mesh = 20,
        degree = 3,
    )
    @show flc_prob.super.num_mesh
    @show flc_prob.super.degree
    @test flc_prob.t_domain == ([-2.0, -2.0], [2.0, 2.0])

    @info "Checking Jacobian..."
    @testset "Manual Jacobian vs ForwardDiff" begin
        prob_cache = get_prob_cache(flc_prob)
        u0 = get_u0(flc_prob)
        rng = MersenneTwister(0)
        for _ in 1:5
            # Trying to avoid state values whose abs is larger than 1;
            # fitting polynomial to random numbers is not nice thing
            # to try here...
            u_lc = 0.01 .* tanh.(randn(rng, eltype(u0), length(u0)))
            d = length(flc_prob.xs0)
            u_lc[d + 1] = abs(u_lc[d + 1])  # period > 0
            H_actual, J_actual = residual_jacobian(u_lc, prob_cache)

            H_desired = similar(@view u_lc[1:end-1])
            J_desired = similar(H_desired, (length(H_desired), length(u_lc)))
            ForwardDiff.jacobian!(
                J_desired,
                (y, x) -> residual!(y, x, prob_cache),
                H_desired,
                u_lc,
            )

            @test J_actual ≈ J_desired
            @test H_actual ≈ H_desired
        end
    end

    @info "Continuation of fold bifurcation of a limit cycle..."
    flc_solver = init(
        flc_prob;
        start_from_nearest_root = true,
        max_branches = 0,  # TODO: stop manually doing this
        verbose = true,
    )
    @time solve!(flc_solver)

    flc_β₁ = [u[end-1] for sweep in flc_solver.super.sol.sweeps for u in sweep.u]
    flc_β₂ = [u[end]   for sweep in flc_solver.super.sol.sweeps for u in sweep.u]
    @show maximum(@. abs(4 * flc_β₁ + flc_β₂^2))
    @test all(@. abs(4 * flc_β₁ + flc_β₂^2) < 5e-2)
    @test maximum(flc_β₂) > 2
    @test minimum(flc_β₂) > -5e-2
end

end  # module


text:
function distance(pnt1, pnt2, norm=2)
    sum = 0
    for i=1:length(pnt1)
        sum += (pnt1[i] - pnt2[i]) ^ norm
    end
    return sum ^ (1 / norm)
end

function angleBetweenVectors(pnt1, pnt2, pnt3)
    c = distance(pnt3, pnt1)
    a = distance(pnt1, pnt2)
    b = distance(pnt2, pnt3)

    try
        return acos((a ^ 2 + b ^ 2 - c ^ 2) /
         (a * b * 2))
    catch
        return pi
    end
end

function convertKinectPntToSimPnt(kinectPnt::Vector{NIFloat}, origin::Vector{NIFloat})
    for i=1:length(origin)
        kinectPnt[i] -= origin[i]
        #now scale down, first convert to inches, then to game points
        kinectPnt[i] *= 0.0393701 / .27
    end
    kinectPnt[2] += 70 #add to y? TODO way to automate
    kinectPnt[3] *= -1 #flip z
    return kinectPnt
end

function convertSkeletonToSim(jointPos::Matrix{NIFloat}, origin::Vector{NIFloat})
    for i=1:length(jointPos[:, 1])
        jointPos[i, :] = convertKinectPntToSimPnt(jointPos[i, :], origin)
    end
end

function convertSkeletonToMJ(sim::MJSimEnv, skeleton::NISkeleton, origin)
    jntLHandID = Int(OpenNI.NI_SKEL_RIGHT_HAND)
    jntLElbowID = Int(OpenNI.NI_SKEL_RIGHT_ELBOW)
    jntLShoulderID = Int(OpenNI.NI_SKEL_RIGHT_SHOULDER)
    jntRHandID = Int(OpenNI.NI_SKEL_LEFT_HAND)
    jntRElbowID = Int(OpenNI.NI_SKEL_LEFT_ELBOW)
    jntRShoulderID = Int(OpenNI.NI_SKEL_LEFT_SHOULDER)
    jntNeckID = Int(OpenNI.NI_SKEL_NECK)
    jntTorsoID = Int(OpenNI.NI_SKEL_TORSO)
    jntLHipID = Int(OpenNI.NI_SKEL_RIGHT_HIP)
    jntRHipID = Int(OpenNI.NI_SKEL_LEFT_HIP)
    jntLKneeID = Int(OpenNI.NI_SKEL_RIGHT_KNEE)
    jntRKneeID = Int(OpenNI.NI_SKEL_LEFT_KNEE)

    #get positions and offset to simplify computations
    LHandPos = skeleton.jointPos[jntLHandID, :] - origin
    LElbowPos = skeleton.jointPos[jntLElbowID, :] - origin
    LShoulderPos = skeleton.jointPos[jntLShoulderID, :] - origin
    RHandPos = skeleton.jointPos[jntRHandID, :] - origin
    RElbowPos = skeleton.jointPos[jntRElbowID, :] - origin
    RShoulderPos = skeleton.jointPos[jntRShoulderID, :] - origin
    neckPos = skeleton.jointPos[jntNeckID, :] - origin
    torsoPos = skeleton.jointPos[jntTorsoID, :] - origin
    hipMidPos = (skeleton.jointPos[jntLHipID, :] + skeleton.jointPos[jntRHipID, :]) / 2 - origin
    kneeMidPos = (skeleton.jointPos[jntLKneeID, :] + skeleton.jointPos[jntRKneeID, :]) / 2 - origin

    #compute the reference angles via law of cosines
    refAngleShoulder1 = angleBetweenVectors([0 -0.17 .12], [0 -0.17 .06], [.18 -0.35 -.12])
    refAngleShoulder2 = angleBetweenVectors([-.01 0 .06], [0 -0.17 .06], [.18 -0.35 -.12])
    refAngleElbow = angleBetweenVectors([0 -0.17 .06], [.18 -0.35 -.12], [.36 -.17 .06])
    refTorsoZ = pi / 2
    refTorsoY = pi

    #compute the actual angles via law of cosines
    thetaLShoulder1 = Float64(angleBetweenVectors(deepcopy(LShoulderPos) + [0, 0, 100], LShoulderPos, LElbowPos) - refAngleShoulder2)
    thetaLShoulder2 = -Float64(angleBetweenVectors(neckPos, LShoulderPos, LElbowPos) - refAngleShoulder1)
    thetaLElbow = -Float64(angleBetweenVectors(LShoulderPos, LElbowPos, LHandPos) - refAngleElbow)
    thetaRShoulder1 = -Float64(angleBetweenVectors(deepcopy(RShoulderPos) + [0, 0, 100], RShoulderPos, RElbowPos) - refAngleShoulder2)
    thetaRShoulder2 = -Float64(angleBetweenVectors(neckPos, RShoulderPos, RElbowPos) - refAngleShoulder1)
    thetaRElbow = -Float64(angleBetweenVectors(RShoulderPos, RElbowPos, RHandPos) - refAngleElbow)
    thetaTorsoZ = Float64(angleBetweenVectors(deepcopy(hipMidPos) + [0, 0, 100], hipMidPos, skeleton.jointPos[jntRHipID, :] - origin) - refTorsoZ)
    thetaTorsoY = -Float64(angleBetweenVectors(kneeMidPos, hipMidPos, torsoPos) - refTorsoY)

    #set the angles
    icxx"""
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "left_shoulder1")]] = $thetaLShoulder1;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "left_shoulder2")]] = $thetaLShoulder2;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "left_elbow")]] = $thetaLElbow;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "right_shoulder1")]] = $thetaRShoulder1;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "right_shoulder2")]] = $thetaRShoulder2;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "right_elbow")]] = $thetaRElbow;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "abdomen_z")]] = $thetaTorsoZ;
    $(sim.data.d)->qpos[$(sim.model.m)->jnt_dofadr[mj_name2id($(sim.model.m), mjOBJ_JOINT, "abdomen_y")]] = $thetaTorsoY;
    """
end

function delimVector(pnt)
    str = ""
    for i=1:length(pnt)
        str = "$str$(pnt[i]),"
    end
    return str
end


text:
# This file was generated, do not modify it. # hide
MLJ.fit!(mach, rows=train, verbosity=3)


text:
# This implementation is WIP. In particular, need to implement switching function that is typically employed within LAMMPS.
# See https://github.com/cesmix-mit/InteratomicPotentials.jl/issues/11

"""
    ZBL{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{()},NamedTuple{(:rcutoff,)}}

The Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion for describing high-energy collisions between atoms.
"""
struct ZBL{T<:AbstractFloat} <: EmpiricalPotential{NamedTuple{()},NamedTuple{(:rcutoff,)}}
    Z₁::Int
    Z₂::Int
    e::T
    rcutoff::T
    species::Tuple
    function ZBL(Z₁, Z₂, e, rcutoff, species)
        Z₁, Z₂, e, rcutoff = promote(Z₁, Z₂, austrip(e), austrip(rcutoff))
        new{typeof(rcutoff)}(Z₁, Z₂, e, rcutoff, Tuple(species))
    end
    function ZBL{T}(; Z₁::Int, Z₂::Int, e::T, rcutoff::T, species::Tuple) where {T}
        new{T}(Z₁, Z₂, e, rcutoff, species)
    end
end

get_rcutoff(zbl::ZBL) = zbl.rcutoff
get_species(zbl::ZBL) = zbl.species

const ϕ_coeffs = [0.18175, 0.50986, 0.28022, 0.02817]
const ϕ_exps = [3.19980, 0.94229, 0.40290, 0.20162]
ϕ(d::AbstractFloat, e::AbstractFloat) = sum(coeff * e^(-exp * d) for (coeff, exp) ∈ zip(ϕ_coeffs, ϕ_exps))
dϕdr(d::AbstractFloat, e::AbstractFloat) = -sum(exp * coeff * e^(-exp * d) for (coeff, exp) ∈ zip(ϕ_coeffs, ϕ_exps))

_zbl_coeff(zbl::ZBL) = kₑ * zbl.Z₁ * zbl.Z₂ * zbl.e^2
_zbl_d(R::AbstractFloat, zbl::ZBL) = R / (0.8854 * 0.529 / (zbl.Z₁^(0.23) + zbl.Z₂^(0.23)))

potential_energy(R::AbstractFloat, zbl::ZBL) = _zbl_coeff(zbl) * ϕ(_zbl_d(R, zbl), zbl.e) / R
force(R::AbstractFloat, r::SVector{3}, zbl::ZBL) = (_zbl_coeff(zbl) * (ϕ(_zbl_d(R, zbl), zbl.e) + _zbl_d(R, zbl) * dϕdr(_zbl_d(R, zbl), zbl.e)) / R^3)r


text:
@info "Loading StatsModels"
# using StatsModels
import StatsModels: @formula, FormulaTerm, ModelFrame, ModelMatrix

@info "Loading DataFrames"
# using DataFrames
import DataFrames: DataFrame

@info "Loading Distributions"
# using Distributions
import Distributions: TDist, Chisq, std, cov, median

@info "Loading Clustering"
# using Clustering
import Clustering: Hclust, hclust, cutree

@info "Loading StatsBase"
# using StatsBase
import StatsBase: quantile, standardize, ZScoreTransform, mean, sample

@info "Loading LinearAlgebra"
# using LinearAlgebra
import LinearAlgebra: inv, pinv, det, diag, eigen, norm, Diagonal, qr, rank

@info "Loading Plots"
# using Plots
import Plots: scatter, title!, xlabel!, ylabel!, hline!, vline!, RGB, plot, font, text

@info "Loading Optim"
# using Optim
import Optim: optimize, NelderMead

@info "Loading Combinatorics"
import Combinatorics: combinations

@info "Loading source"
include("../basis.jl")
include("../data.jl")
include("../diagnostics.jl")
include("../ols.jl")
include("../hs93.jl")
include("../ks89.jl")
include("../smr98.jl")
include("../asm2000.jl")
include("../lms.jl")
include("../lts.jl")
include("../mve.jl")
include("../mveltsplot.jl")
include("../bch.jl")
include("../py95.jl")
include("../satman2013.jl")
include("../satman2015.jl")
include("../lad.jl")
include("../lta.jl")
include("../hadi1992.jl")
include("../hadi1994.jl")
include("../dataimage.jl")
include("../cga.jl")
include("../ga.jl")
include("../gwlts.jl")
include("../summary.jl")
include("../ransac.jl")
include("../atkinson94.jl")
include("../atkinsonstalactiteplot.jl")
include("../imon2005.jl")
include("../ccf.jl")
include("../bacon.jl")

println("ready")

using BenchmarkTools

reg = createRegressionSetting(@formula(y ~ x1 + x2 + x3), hbk)

max_stars_repo_name: akadal/LinRegOutliers

max_issues_repo_name: akadal/LinRegOutliers

max_forks_repo_name: akadal/LinRegOutliers


text:
## test data available for use in PopGen

"""
    nancycats()
Returns a `PopObj` of corresponding "nancycats" dataset as featured in
the R package `adegenet`. This is microsatellite data of 9 loci in 237
individuals across 17 populations.

Example:

`nancy = nancycats()`
"""
function nancycats()
    filename = normpath(joinpath(@__DIR__,"..","data", "data", "nancycats.csv"))
    return csv(filename, marker = "msat")
end

"""
    gulfsharks()
Returns a `PopObj` corresponding the Blacknose shark dataset as used in
Dimens et al. 2019. This is a mid-sized SNP dataset of 2213 SNPs in 212
individuals, across 7 populations.

Example:

`sharks = gulfsharks()`
"""
function gulfsharks()
    filename = normpath(joinpath(@__DIR__,"..","data", "data", "gulfsharks.gen"))
    xloc = [-80.59928
            -80.59954
            -80.59958
            -80.42248
            -80.35779
            -80.35786
            -80.34815
            -80.34822
            -80.43027
            -80.24937
            -80.24964
            -80.24959
            -80.24976
            -80.24954
            -80.24944
            -80.49627
            -80.49682
            -80.4967
            -80.49664
            -80.49605
            -80.10169
            -81.30037
            -81.3217
            -81.42224
            -81.29336
            -81.29336
            -81.38756
            -81.38756
            -81.39475
            -81.36869
            -81.36879
            -81.4027
            -81.20213
            -81.39859
            -81.39798
            -81.15039
            -81.37432
            -81.38336
            -81.34095
            -81.34075
            -81.34147
            -81.34094
            -81.34122
            -81.3417
            -81.33707
            -81.43061
            -81.20041
            -81.43111
            -81.43081
            -81.43102
            -81.25828
            -79.81894
            -79.79036
            -79.79086
            -80.35753
            -80.35754
            -80.35783
            -80.35765
            -80.35712
            -80.35712
            -80.35799
            -80.357
            -80.44488
            -80.44511
            -79.81454
            -79.82027
            -79.80382
            -79.8041
            -79.80383
            -79.80408
            -79.80471
            -79.80422
            -79.80381
            -79.80425
            -79.80451
            -79.80383
            -79.80434
            -80.70928
            -79.84989
            -81.56379
            -81.56447
            -81.42769
            -81.41631
            -81.16435
            -81.16054
            -81.09332
            -80.73249
            -80.93939
            -80.76868
            -80.76853
            -80.74511
            -80.74504
            -80.80831
            -81.78288
            -81.78309
            -81.80621
            -81.80677
            -81.81822
            -81.81819
            -81.76506
            -81.75001
            -81.74955
            -81.7501
            -81.73252
            -81.72459
            -81.86906
            -82.07268
            -82.07234
            -82.07196
            -82.07248
            -82.07196
            -82.07217
            -82.07264
            -82.08975
            -82.05538
            -82.08971
            -82.0546
            -82.05489
            -82.0748
            -82.07446
            -82.07476
            -82.09188
            -82.09226
            -82.09236
            -82.0916
            -82.09153
            -82.09168
            -82.09171
            -82.09156
            -80.9965
            -81.04498
            -81.00124
            -81.00353
            -80.99296
            -81.00691
            -80.99782
            -81.00029
            -80.99829
            -80.99768
            -81.00383
            -81.00061
            -80.99762
            -81.75011
            -81.03748
            -82.84887
            -83.60492
            -82.72491
            -82.72514
            -82.72485
            -82.72501
            -82.61679
            -82.61668
            -82.56186
            -83.08778
            -82.49909
            -82.74918
            -82.7396
            -82.80298
            -82.66804
            -82.66804
            -82.73849
            -82.74948
            -82.74963
            -82.66823
            -82.80976
            -82.80976
            -82.74893
            -82.80714
            -82.74321
            -82.79686
            -82.77452
            -82.8552
            -82.79715
            -84.72682
            -84.13857
            -84.13877
            -84.13855
            -83.42212
            -84.49099
            -84.45066
            -84.44933
            -84.48637
            -84.48524
            -85.15517
            -85.15494
            -84.48488
            -84.50144
            -84.50057
            -84.47698
            -84.48949
            -84.50102
            -85.23794
            -87.71888
            -88.00071
            -88.04933
            -88.04904
            -88.16747
            -88.05461
            -88.2974
            -88.35997
            -88.39216
            -88.16801
            -86.03987
            -85.64937
            -87.44028
            -86.53741
            -86.53756
            -86.09053
            -87.36605
            -87.36617
            -85.71432
            -85.71432
            ]
    yloc = [28.30624
            28.30787
            28.30234
            28.61234
            27.86661
            27.86659
            27.86821
            27.8711
            28.48153
            28.24211
            28.24276
            28.2378
            28.24141
            28.23596
            28.23644
            28.39057
            28.39016
            28.38607
            28.38686
            28.38649
            27.57661
            31.16833
            31.15513
            30.82858
            30.78698
            30.78698
            30.83904
            30.83904
            30.84482
            30.96765
            30.96915
            30.84907
            31.11922
            30.85875
            30.86431
            31.3305
            30.95857
            30.40459
            31.01552
            31.01451
            31.0157
            31.01682
            31.01409
            31.01971
            31.00965
            30.58793
            31.09504
            30.58181
            30.58476
            30.59014
            31.22919
            32.72165
            32.66413
            32.66126
            32.44526
            32.45119
            32.45071
            32.44581
            32.4508
            32.45339
            32.44588
            32.45174
            32.45064
            32.45696
            32.73458
            32.73266
            32.72372
            32.73138
            32.72676
            32.73036
            32.7323
            32.72279
            32.7283
            32.72665
            32.72764
            32.73162
            32.72347
            32.31162
            32.73778
            24.76982
            24.77288
            24.79486
            24.82328
            24.59827
            24.60102
            24.61792
            24.75817
            24.71942
            24.83066
            24.83652
            24.84129
            24.83924
            24.80488
            24.45817
            24.45688
            24.46559
            24.46232
            24.44635
            24.44414
            24.45884
            24.46137
            24.46082
            24.45485
            24.47263
            24.45921
            24.44529
            24.43813
            24.44609
            24.44349
            24.44051
            24.44121
            24.44217
            24.44578
            24.44059
            24.4335
            24.44744
            24.43314
            24.43995
            24.44429
            24.43938
            24.44675
            24.44683
            24.45088
            24.45091
            24.44666
            24.44509
            24.44284
            24.44985
            24.44255
            25.0083
            25.10191
            25.01428
            25.01907
            25.01096
            25.01279
            25.00765
            25.01504
            25.01565
            25.01249
            25.01759
            25.02244
            25.00738
            24.4617
            25.10086
            27.45678
            27.35171
            26.24867
            26.2486
            26.25634
            26.25514
            25.4607
            25.45863
            25.92942
            25.51734
            25.13769
            24.83092
            24.92929
            24.92038
            24.66909
            24.67271
            24.88658
            24.82947
            24.82769
            24.66979
            28.33518
            28.33691
            28.34884
            28.32341
            28.44358
            28.77552
            28.31135
            28.40218
            28.77383
            28.61922
            29.9992
            29.99503
            29.99472
            29.60836
            29.88518
            29.85259
            29.86432
            29.84045
            29.8507
            29.50347
            29.51108
            29.85473
            29.83699
            29.82836
            29.88213
            29.83377
            29.83382
            29.04839
            29.89008
            30.19427
            30.00212
            30.00689
            29.8362
            29.50566
            30.14284
            30.20737
            30.11506
            29.8362
            29.94663
            29.99695
            29.6966
            30.06404
            30.0696
            29.90653
            30.05317
            30.05217
            29.82344
            29.82344
            ]
    data = genepop(filename)
    data.samples.latitude = yloc ; data.samples.longitude = xloc
    renames = Dict(
        1 => "Cape Canaveral",
        2 => "Georgia",
        3 => "South Carolina",
        4 => "Florida Keys",
        5 => "Mideast Gulf",
        6 => "Northeast Gulf",
        7 => "Southeast Gulf"
    )
    # manually rename things so as to not trigger @info prompt in populations!()
    for eachkey in keys(renames)
        replace!(data.samples.population, eachkey => renames[eachkey])
    end

    return data
end


text:
using DFControl

#This example goes through how one would create a job from scratch, using a .cif file for the Structure.

#First go to your favourite Crystal structure database and download the .cif you want to use.
#e.g. Si (F d -3 m :1) : http://www.crystallography.net/cod/9011998.cif
cif_file = Downloads.download("http://www.crystallography.net/cod/9011998.cif", "Si.cif")

# Construct your structure from it, and set the pseudos to one of the sets we defined in the previous example
str = Structure(cif_file; name = "Si")
set_pseudos!(str, :pbesol)

# Since we are going to run QuantumEspresso, we define the executables to be used
pw_execs = [Exec("mpirun", "", :np => 4), Exec("pw.x", "/opt/qe/bin/", :nk => 4)]

# Then we generate the first calculation for our job, we name it scf, pass the executables to be used
# and set some specific flags.
# Afterwards the kpoints can also be set
scf_calculation = Calculation{QE}("scf", pw_execs, :calculation => "scf")
set_kpoints!(scf_calculation, (6, 6, 6, 1, 1, 1))
# Or

scf_calculation = Calculation{QE}("scf", pw_execs, :calculation => "scf";
                                  data = [InputData(:k_points, :automatic,
                                                    (6, 6, 6, 1, 1, 1))])

# Using these we can now define our job, if we would have more calculations they would be added
# to the list [scf_calculation].
# The flag => value pairs will set the specified flags to that value for all calculations in the job
# that allow that flag to be set, so it's ideal for things like cutoffs and smearing etc.
# The header can contain any lines that will be pasted in the job script in front of
# the calculation lines.
job = DFJob("Si", str, [scf_calculation], :ecutwfc => 20, #these flags will be set_ on all calculations that are passed to the job
            :verbosity => "high", :conv_thr => 1e-6; header = ["export OMP_NUM_THREADS=1"])

# Now the job can be submitted to be ran.
submit(job)
# this first saves the job and it's calculation files to the `job.local_dir` then pushes the `job.tt` file and the calculations to the `job.server_dir` on `job.server`, and tries to first submit it through slurm, otherwise just runs `bash job.tt`
# You can check the job.local_dir to see the calculation files and `job.tt` script.

# If the job was submitted via slurm, one can check if it's running through
slurm_isrunning(job)

# Hopefully everything went according to plan and we can retrieve the outputs
out = outputdata(job)
# Or, to retrieve the output of a specific calculation,
out = outputdata(job["scf"])

# Indeed specific calculations of the job can be accessed through their name
# If for some reason we discover that smearing was needed, the flags related to that
# can be changed either job-wide (the flag will be set for every calculation that allows for that flag):
job[:occupations] = "smearing"
job[:degauss] = 0.001
# Or on a specific calculation
job["scf"][:occupations] = "smearing"
job["scf"][:degauss] = "smearing"

# Next step would be to run a bands calculation:
bands_in = gencalc_bands(job["scf"], high_symmetry_kpath(job.structure))
# Here we generate a bands calculation using the scf as the template, and the k-points
# generated from the structure symmetries.
# k-points can also be passed as e.g. :
bands_in = gencalc_bands(job["scf"], [(0.5, 0.5, 0.5, 10.0), (0.0, 0.0, 0.0, 1.0)])
# mimicking the usual structure in the bands pw.x calculation.

# Then we add the bands calculation to the job,
push!(job, bands_in)
# and since we already ran the scf we can do
set_flow!(job, "scf" => false, "bands" => true)

# The String => Bool pairs supplied to setflow will sequentially fuzzy match any name
# in which the string occurs. i.e.
set_flow!(job, "" => false, "bands" => true)
# would first set all calculations to not run, followed by setting the bands calculation to run
set_flow!(job, "" => false, "scf" => true)
# would run calculations with both scf and nscf names.

# Anyway now we submit the bands calculation
set_flow!(job, "" => false, "bands" => true)
submit(job)

#now the bandstructure can be plotted
bands = outputdata(job["bands"])[:bands]
#or
bands = readbands(job["bands"])
#or
bands = readbands(job)
# The last will search for a bands calculation in the job, otherwise it will use the output of an nscf calculation
# to generate the bands

# similarly
fermi = outputdata(job["scf"])[:fermi]
fermi = readfermi(job["scf"])
fermi = readfermi(job)

using Plots

# Now the bands are plotted
plot(bands; fermi = fermi)

# To access the dos, we can first generate an nscf to generate a uniform k-grid calculation
# followed by a projwfc

push!(job, gencalc_nscf(job["scf"], (6, 6, 6)))
push!(job, gencalc_projwfc(job["nscf"], fermi - 10, fermi + 10, 0.1))

setflow!(job, "" => false, "prowjfc" => true, "nscf" => true)
submit(job)

# Now we can do the following:
plot(job, -5, 5)
# which will plot the band structure colored by the dos between fermi-5 and fermi+5

# This concludes this example.


text:
@everywhere begin
    import AARMBEM.AARMBEMheat.heattransfer
end

heattransfer(ARGS)


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# void hts_srand48(long seed);
@defun hts_srand48(seed::Clong)::Cvoid

# double hts_erand48(unsigned short xseed[3]);
@defun hts_erand48(xseed::Ptr{Cushort})::Cdouble

# double hts_drand48(void);
@defun hts_drand48()::Cdouble

# long hts_lrand48(void);
@defun hts_lrand48()::Clong


text:
using ..Fed: STATS

"""
    QDiffPayload{T <: Unsigned}(data, minval, maxval)

Payload generated by the quantization compression scheme. Data are structured
in a different way than `QPayload` in order to diff-deduplicate the client-side
payloads by chunks.
"""
mutable struct QDiffPayload{T <: Unsigned}
    data::Patch{T}
    minval::Float32
    maxval::Float32
end

"""
    QDiffPayloadSerde{T <: Unsigned}(chunksize::Int, is_client::Bool) <: PayloadSerde

Apply both quantization and diff-deduplication to the weights before
serialization in order to compress the payload. Quantization is a lossy process,
thus a loss of information is to be expected after deserialization.
"""
mutable struct QDiffPayloadSerde{T <: Unsigned} <: PayloadSerde
    chunksize::Int
    is_client::Bool
    original_data::Vector{T}

    QDiffPayloadSerde{T}(chunksize, is_client) where T <: Unsigned = new(
        chunksize,
        is_client,
        Vector{Vector{T}}(undef, 0)
    )
end

"""
    serialize_payload(::QDiffPayloadSerde, weights)

Serialize `weights` with the `QuantizedDedupPayloadSerde` where quantization
and diff-deduplication of values are applied before serialization.
"""
function serialize_payload(
    p::QDiffPayloadSerde{T},
    weights::Vector{Float32}
)::Vector{UInt8} where T <: Unsigned
    # quantize weights
    q = Quantizer{T}(weights)
    qweights = [quantize(q, w) for w in weights]
    STATS.base.req_data = qweights

    if p.is_client
        patched_file = diff(p.original_data, qweights, sizeof(T))
    else
        p.original_data = qweights
        patched_file = Patch{T}(UnitRange[], qweights)
    end

    payload = QDiffPayload{T}(patched_file, q.minval, q.maxval)

    return pack(payload)
end

"""
    deserialize_payload(::QuantizedDedupPayloadSerde, data, from)

Deserialize `data` with the `QuantizedDedupPayloadSerde` where regeneration of
deduplicated values and dequantizatoin are applied after deserialization.
"""
function deserialize_payload(
    p::QDiffPayloadSerde{T},
    data::Vector{UInt8},
    from::String
)::Vector{Float32} where T <: Unsigned
    payload = unpack(data)

    if p.is_client
        qweights = payload.data.data
        p.original_data = qweights
    else
        delta_earned = sum(length.(payload.data.range)) * sizeof(T) - (8 * length(payload.data.range))
        STATS.network.delta_earned += delta_earned

        qweights = patch(p.original_data, payload.data)
        push!(STATS.base.res_data, qweights)
    end

    # dequantize weights
    q = Quantizer{T}(payload.minval, payload.maxval)
    weights = [dequantize(q, w) for w in qweights]

    return weights
end


text:
module DeprecationBotTest

function foobar(x, y::T) where T
    println(x, y)
end

if !(isdefined(Core, :String) && isdefined(Core, :AbstractString))
    function jlconvert(T::Union{Type{String},Type{String}}, ::JldFile, ptr::Ptr)
        strptr = unsafe_load(convert(Ptr{Ptr{UInt8}}, ptr))
        n = Int(ccall(:strlen, Csize_t, (Ptr{UInt8},), strptr))
        T(unsafe_wrap(Array, strptr, n, true))
    end
end

_foo(x) = foo(x)

bar(x) = _foo(x)

function foobar(a,
            b::T,
            c::T) where T
    (a, b, c)
end

end # module

max_stars_repo_name: nkottary/DeprecationBotTest

max_issues_repo_name: nkottary/DeprecationBotTest

max_forks_repo_name: nkottary/DeprecationBotTest


text:
using SIMD

# Rank one update
#   C := C + ap * bp^T
# with rank-1 updates computed using axpy with vector ap
#   (ap * bp^T)_j = bpj * ap
function mygemm_pji(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for p = 1:k
    for j = 1:n
      for i = 1:m
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

function mygemm_jpi(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for j = 1:n
    for p = 1:k
      for i = 1:m
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# block Matrix times block column vector
#   Cj := A * Bj + Cj
# with inner A * Cj computed using block dot products:
#   (A * Bj)_i = Ai^T Bj
function mygemm_jip_pji(C, A, B, NB=4, MB=4)
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for J = 1:NB:n
    J2 = min(n, J+NB-1)
    for I = 1:MB:m
      I2 = min(m, I+MB-1)
      @views mygemm_pji(C[I:I2, J:J2], A[I:I2, :], B[:, J:J2])
    end
  end
end

# block Matrix times block column vector
#   Cj := A * Bj + Cj
# with inner A * Cj computed using block dot products:
#   (A * Bj)_i = Ai^T Bj
function mygemm_ji_4x4(C, A, B)
  (m, n) = size(C)
  (~, k) = size(A)
  @assert mod(m, 4) == 0
  @assert mod(n, 4) == 0
  @inbounds for J = 1:4:n
    for I = 1:4:m
      @views mygemm_pji_4x4(k,
                            pointer(C, I + (J-1) * m),
                            pointer(A, I),
                            B,
                            m,
                            m,
                            J)
      # @views mygemm_pji_4x4(k,
      #                       pointer(C) + I-1 + (J-1)*m,
      #                       pointer(A) + I-1,
      #                       pointer(V) + B[:, J:J2], m, m, k)
    end
  end
end

function mygemm_pji_4x4(k, C, A, B, ldc, lda, j)
  T = eltype(A)

  c1 = vload(Vec{4, T}, C + 0ldc*sizeof(T))
  c2 = vload(Vec{4, T}, C + 1ldc*sizeof(T))
  c3 = vload(Vec{4, T}, C + 2ldc*sizeof(T))
  c4 = vload(Vec{4, T}, C + 3ldc*sizeof(T))

  @inbounds for p = 1:k
    αv = vload(Vec{4, T}, A + (p-1)*lda*sizeof(T))

    # β = Vec{4, T}(B[p, 1])
    β = B[p, j + 0]
    c1 += β * αv

    # β = Vec{4, T}(B[p, 2])
    β = B[p, j + 1]
    c2 += β * αv

    # β = Vec{4, T}(B[p, 3])
    β = B[p, j + 2]
    c3 += β * αv

    # β = Vec{4, T}(B[p, 4])
    β = B[p, j + 3]
    c4 += β * αv
  end
  vstore(c1, C + 0ldc*sizeof(T))
  vstore(c2, C + 1ldc*sizeof(T))
  vstore(c3, C + 2ldc*sizeof(T))
  vstore(c4, C + 3ldc*sizeof(T))
end

function mygemm_ji(C, A, B, ::Val{MB}, ::Val{NB}) where {NB, MB}
  (m, n) = size(C)
  (~, k) = size(A)
  @assert mod(m, MB) == 0
  @assert mod(n, NB) == 0
  @inbounds for J = 1:NB:n
    for I = 1:MB:m
      @views mygemm_p(k,
                      pointer(C, I + (J-1) * m),
                      pointer(A, I),
                      B,
                      m,
                      m,
                      J, Val(MB), Val(NB))
    end
  end
end

function mygemm_p(k, C, A, B, ldc, lda, j, ::Val{MB}, ::Val{NB}) where {MB, NB}
  T = eltype(A)

  c1 = vload(Vec{MB, T}, C + 0ldc*sizeof(T))
  c2 = vload(Vec{MB, T}, C + 1ldc*sizeof(T))
  c3 = vload(Vec{MB, T}, C + 2ldc*sizeof(T))
  c4 = vload(Vec{MB, T}, C + 3ldc*sizeof(T))

  @inbounds for p = 1:k
    αv = vload(Vec{MB, T}, A + (p-1)*lda*sizeof(T))

    β = B[p, j + 0]
    c1 = muladd(β, αv, c1)

    if NB > 1
      β = B[p, j + 1]
      c2 = muladd(β, αv, c2)
    end

    if NB > 2
      β = B[p, j + 2]
      c3 = muladd(β, αv, c3)
    end

    if NB > 3
      β = B[p, j + 3]
      c4 = muladd(β, αv, c4)
    end

    if NB > 4
      error("NB too big")
    end
  end
  vstore(c1, C + 0ldc*sizeof(T))
  vstore(c2, C + 1ldc*sizeof(T))
  vstore(c3, C + 2ldc*sizeof(T))
  vstore(c4, C + 3ldc*sizeof(T))
end

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP


text:
mutable struct ChernMat{T<:Integer}
    B::Matrix{T}
    m::Int64     # Number of constraints
    n::Int64     # Dimension
    rows::Int    # Current number of rows
end

function ChernMat(A::Matrix{T}) where T<:Integer
    m, n = size(A)
    ChernMat([I A'], m, n, n)
end

function remove_row(A::Matrix{T}, k::Int) where T<:Integer
    m, n = size(A)
    B = Array{T}(undef, m - 1, n)
    i = 1
    for j in 1:m
        if j != k
            B[i,:] = A[j,:]
            i+=1
        end
    end
    return B
end

function add_row(A::Matrix{T}, a::Matrix{T}) where T<:Integer
    # Here a is a row vector
    m, n = size(A)
    B = Array{T}(undef, m+1, n)
    B[1:m, 1:n] = A
    B[m+1, :] = a
    return B
end

function pprint(mat::ChernMat)
    @printf("# Rays: %7d\n", mat.rows)
    for i in 1:mat.rows
        for j in 1:mat.n
            @printf("%6d", mat.B[i,j])
        end
        print(" |")
        for j in mat.n+1:mat.n+mat.m
            @printf("%6d", mat.B[i,j])
        end
        print("\n")
    end
end

function chernikova(mat::ChernMat{T}, verbosity::Int = 0) where T<:Integer
    m, n = mat.m, mat.n
    counter = 1
    while true
        if verbosity > 0
            println("Iteration ", counter)
            println("---------------")
            pprint(mat)
        end
        (col, p) = leading_col(mat)
        if col == 0
            if verbosity > 0; println("Generators have been found!"); end;
            return norm_cols(LHS(mat)')
        elseif p == mat.rows
            if verbosity > 0; println("Zero is the unique solution"); end;
            return Array{T}(undef, mat.n, 0)
        end

        # Inspect matrix to find the maximum size of the next
        # if verbosity > 0; println("Using column $(col) which has $(p) negatives"); end;
        pos, nil, neg = get_pos_neg_indices(mat, col)

        if verbosity > 0
            println("Adding constraint $(col-mat.n) (column $(col)) which is currently violated by $(length(neg)) of $(mat.rows) extremal rays")
            println("New basis will consist of a maximum of $(mat.rows) - $(length(neg)) + $(length(pos)) * $(length(neg)) = $(mat.rows - length(neg) + length(pos)*length(neg)) extremal rays")
        end;

        max_rows = length(pos) + length(nil) + length(pos) * length(neg)
        new_B = Array{T}(undef, max_rows, mat.m + mat.n)

        # Handle case of two rows separately
        if mat.rows == 2 && max_rows == 2
            if verbosity > 0; println("Handling special two row case..."); end;
            new_B[1,:] = mat.B[pos[1],:]
            new_row = combine(mat.B, col, pos[1], neg[1])
            if verbosity > 1; println("Candidate row: ", new_row); end;
            if !collinear(vec(mat.B[pos[1],:]), vec(new_row))
                new_B[2,:] = new_row
                mat.rows = 2
                mat.B = new_B
            else
                mat.rows = 1
                mat.B = new_B[1,:]
            end
            counter += 1
            continue
        end

        # Add rows with non-negative intersections
        # with leading column to next matrix
        for (k,i) in enumerate([pos; nil]); new_B[k,:] = mat.B[i,:]; end;
        new_rows = length(pos) + length(nil)

        # Go through pairs of row with positive and negative
        # coefficients in leading column and combine if necessary
        for i1 in pos, i2 in neg
            if combinable(mat, i1, i2, verbosity - 1)
                new_rows += 1
                new_B[new_rows, :] = combine(mat.B, col, i1, i2)
            end
        end
        mat.B = new_B[1:new_rows,:]
        mat.rows = new_rows
        counter += 1

        if verbosity > 0
            @printf("Added %d / %d positive combinations of extremal rays\n",
                    new_rows - length(pos) - length(nil), length(pos) * length(neg))
            println("------------------------------------------------")
        end
    end
end

@doc """
# Description
Calculates the finite generators of the cone formed
by the intersection of a polyhedral cone and the positive quadrant;
that is the cone formed by the solutions
of the following problem:

   Ax >= 0
    x >= 0

This is solved by Chernikova's algorithm:
"Algorithm for finding a general formula for the non-negative solutions of linear inequalities"
by N.V. Chernikova (1964)

# Arguments
* `A::Matrix{Integer}`: Matrix of cone constraints (each column defines a constraint)
* `verbosity`: level of information printed to standard output

# Returns
* `rays::Matrix{Integer}`: a matrix whose columns are the required cone generators
""" ->
function chernikova(A::Matrix{T}, verbosity::Int = 0) where T<:Integer
    mat = ChernMat(A)
    chernikova(mat, verbosity)
end

function chernikova(A::Matrix{T}, verbosity::Int = 0) where T<:Rational
    C = intmat(A)
    chernikova(C, verbosity)
end

function leading_col(mat::ChernMat{T}) where T<:Integer
    p = mat.rows
    I = 0
    for i in mat.n+1:mat.n+mat.m
        r = sum(mat.B[:,i] .< 0.0)
        if r == mat.rows
            return (i, r)
        elseif (r > 0) && (r < p)
           I, p = i, r
        end
    end
    return (I, p)
end

LHS(mat::ChernMat) = mat.B[:, 1:mat.n]

function get_pos_neg_indices(mat::ChernMat{T}, col::Int64) where T<:Integer
    pos, nil, neg = Int[],Int[],Int[]
    for k in 1:mat.rows
        if mat.B[k,col] > 0.0; push!(pos, k)
        elseif mat.B[k,col] < 0.0; push!(neg, k)
        else push!(nil, k) end
    end
    pos, nil, neg
end

function combinable(mat::ChernMat{T}, i1::Int64, i2::Int64, verbosity::Int = 0) where T<:Integer
    # Find all columns for which the rows have a common zero
    # LHS
    zero_cols = Int[]
    for j in 1:mat.n
        if (mat.B[i1,j] == 0.0) && (mat.B[i2,j] == 0.0)
            push!(zero_cols, j)
        end
    end

    # RHS
    for j in mat.n+1:mat.n + mat.m
        if all(mat.B[:,j] .>= 0.0)
            if (mat.B[i1,j] == 0.0) && (mat.B[i2,j] == 0)
                push!(zero_cols, j)
            end
        end
    end

    if length(zero_cols) == 0
        if verbosity > 0; println("Rays $(i1) and $(i2) saturate no common constraint - they do not combine to form a new extremal ray") end
        return false
    end

    if verbosity > 0;
        println("Rays $(i1) and $(i2) have ", length(zero_cols), " are both saturated by constraints ", zero_cols);
    end;

    for i in 1:mat.rows
        if i == i1 || i == i2; continue; end

        if verbosity > 1; println("\tRay $i - values on saturated constraints: ",  mat.B[i,zero_cols]); end;
        if vec(mat.B[i,zero_cols]) == zeros(length(zero_cols))
            if verbosity > 0; println("\tRay $i saturates all saturated constraints of these rays so do not combine") end;
            return false
        end
    end
    if verbosity > 0;
        println("\tNo other ray saturates the common saturated constraints of these rays so combine to form a new extreme ray")
    end;
    return true
end

# function combine(mat::Matrix{T}, col::Int, i1::Int, i2::Int) where T<:Integer
#     a, b = 1.0, -mat[i1,col]/mat[i2,col]
#     return a*mat[i1,:] + b*mat[i2,:]
# end

function combine(mat::Matrix{T}, col::Int, i1::Int, i2::Int) where T<:Integer
    t = gcd(mat[i1, col], mat[i2, col])
    a1, a2 = div(mat[i1, col],t), div(mat[i2, col],t)
    return a1 * mat[i2,:] - a2 * mat[i1,:]
end

function combine(col::Int, vec1::Vector{T}, vec2::Vector{T}) where T<:Integer
    t = gcd(vec1[col], vec2[col])
    a1, a2 = div(vec1[col],t), div(vec2[col],t)
    return a1 * vec2 - a2 * vec1
end


text:
@testset "Misc" begin
  include("Misc/psibound.jl")
  include("Misc/Primes.jl")
  include("Misc/Places.jl")
  include("Misc/meataxe.jl")
  include("Misc/stable_subgroups.jl")
  include("Misc/NumberField.jl")
  include("Misc/Poly.jl")
  include("Misc/Localization.jl")
end


text:
# This script is called automatically upon package installation.
using PyCall
if PyCall.conda
    import Conda
    Conda.add("spglib"; channel="conda-forge")
    Conda.add("pymatgen"; channel="conda-forge")
end


text:
using Test
using Elemental
using DistributedArrays
using LinearAlgebra

A  = drandn(50,50)
Al = Matrix(A)
B  = drandn(50,10)
Bl = Matrix(B)

@test inv(Al) ≈ inv(A)
@test Al\Bl ≈ A\B
@test logdet(Al'Al) ≈ logdet(A'A)
@test svdvals(Al) ≈ Elemental.svdvals(A)


text:
using Test
using Complementarity

@testset "Complementarity Tests" begin

    include("../examples/mcp/transmcp.jl")
    include("transmcp_numeric.jl")

    @testset "Array{Variable,4} Test" begin
        m = MCPModel()
        @variable(m, x[1:2, 1:3, 1:4, 1:5] >=0 )
        @mapping(m, F[i=1:2, j=1:3, k=1:4, l=1:5], x[i,j,k,l])
        # for i=1:2, j=1:3, k=1:4, l=1:5
        #     @complementarity(m, F[i,j,k,l], x[i,j,k,l])
        # end
        @complementarity(m, F, x)
        status = solveMCP(m)
        @test status == :Solved
    end

    include("hhoeschle_mcp1.jl")

    include("runtests_mcp.jl")

    include("runtests_mcp_nlsolve.jl")

    include("runtests_mpec.jl")

    include("test_bounds.jl")
end


text:
using ..Cos
using ..Common

latin_charset_encoding() = load_data_file("latin-charset-encoding.txt")

function get_latin_charset_dict(col)
    dict = Dict{UInt8,CosName}()
    m = latin_charset_encoding()
    map(m[:,1], m[:,col]) do x, y
        if y != "-"
            v = to_uint8(y)
            dict[v] = CosName(strip(x))
        end
    end
    return dict
end

const STDEncoding_to_GlyphName = get_latin_charset_dict(2)
MACEncoding_to_GlyphName = get_latin_charset_dict(3)
WINEncoding_to_GlyphName = get_latin_charset_dict(4)
const PDFEncoding_to_Glyphname = get_latin_charset_dict(5)

#Special handling for few cases as per spec
WINEncoding_to_GlyphName[0xAD] = cn"sfthyphen"
WINEncoding_to_GlyphName[0xA0] = cn"colon"

MACEncoding_to_GlyphName[0xCA] = cn"colon"

function reverse_dict(dict)
    rdict = Dict()
    for (x, y) in dict
        rdict[y] = x
    end
    return rdict
end

const GlyphName_to_STDEncoding = reverse_dict(STDEncoding_to_GlyphName)
const GlyphName_to_MACEncoding = reverse_dict(MACEncoding_to_GlyphName)
const GlyphName_to_WINEncoding = reverse_dict(WINEncoding_to_GlyphName)
const Glyphname_to_PDFEncoding = reverse_dict(PDFEncoding_to_Glyphname)

macexpt_charset_encoding() = load_data_file("mac-expert.txt")
symbols_charset_encoding() = load_data_file("symbols-encoding.txt")
zapf_charset_encoding()    = load_data_file("zapfdingbats-encoding.txt")

function get_charset_dict(f::Function)
    dict = Dict{UInt8,CosName}()
    m = f()
    map(m[:,2], m[:,3]) do x, y
        v = to_uint8(y)
        dict[v] = CosName(strip(x))
    end
    return dict
end

const MEXEncoding_to_GlyphName = get_charset_dict(macexpt_charset_encoding)
const SYMEncoding_to_GlyphName = get_charset_dict(symbols_charset_encoding)
const ZAPEncoding_to_GlyphName = get_charset_dict(zapf_charset_encoding)

const Glyphname_to_MEXEncoding = reverse_dict(MEXEncoding_to_GlyphName)
const GlyphName_to_SYMEncoding = reverse_dict(SYMEncoding_to_GlyphName)
const GlyphName_to_ZAPEncoding = reverse_dict(ZAPEncoding_to_GlyphName)

using AdobeGlyphList

function agl_mapping_to_dict(m)
    dict = Dict{CosName, Char}()
    map(m[:,1], m[:,2]) do x, y
        dict[CosName(strip(x))] = y
    end
    return dict
end

const AGL_Glyph_to_Unicode = agl_mapping_to_dict(agl())
const AGL_ZAP_to_Unicode   = agl_mapping_to_dict(zapfdingbats())
const AGL_Unicode_to_Glyph = reverse_dict(AGL_Glyph_to_Unicode)
const AGL_Unicode_to_ZAP   = reverse_dict(AGL_ZAP_to_Unicode)

const STDEncoding_to_Unicode = dict_remap(STDEncoding_to_GlyphName, AGL_Glyph_to_Unicode)
const MACEncoding_to_Unicode = dict_remap(MACEncoding_to_GlyphName, AGL_Glyph_to_Unicode)
const WINEncoding_to_Unicode = dict_remap(WINEncoding_to_GlyphName, AGL_Glyph_to_Unicode)

const MEXEncoding_to_Unicode = dict_remap(MEXEncoding_to_GlyphName, AGL_Glyph_to_Unicode)
const SYMEncoding_to_Unicode = dict_remap(SYMEncoding_to_GlyphName, AGL_Glyph_to_Unicode)
const ZAPEncoding_to_Unicode = dict_remap(ZAPEncoding_to_GlyphName, AGL_ZAP_to_Unicode)


text:
function maketablea2(dist::DataFrame;
    panela::Bool = true,
    panelb::Bool = true,
    )::Nothing

  #adjust the dataframe for the Harzmark Solomon stats
  hs::DataFrame = makehs(dist)

  #########################Panel a
  TA2_PANEL_A_ROWKEYS::Vector{Symbol} = [
    :ann,
    :interim,
    :ex,
    :annex,
    :ex40]

  TA2_PANEL_A_ROWLABEL::Dict = Dict(
    :ann=>"Announcement day",
    :interim=>"Interim",
    :ex=>"Ex-day",
    :annex=>"Ann-to-ex",
    :ex40=>"40 days after"
    )
  TA2_PANEL_A_COLNAMES::Vector{Vector{String}} = [["Return"]]

  #####################Panel b#######
  TA2_PANEL_B_ROWKEYS::Vector{Symbol} = [
    :constant,
    :yield,
    :r2,
    :N,]

  TA2_PANEL_B_ROWLABEL::Dict = Dict(
    :constant=>"Constant",
    :yield=>"Div yield",
    :r2=>"\$R^2\$",
    :N=>"N",
    )

  TA2_PANEL_B_COLKEYS::Vector{Symbol} = [
    :ann,
    :interim,
    :ex,
    :annex,
    :ex40
  ]

  TA2_PANEL_B_COLNAMES::Vector{Vector{String}} = [[
    "Announcement day", "Interim", "Ex-day", "Ann-to-ex", "40 days after"]]

  TA2_PANEL_B_COLFILTER = Dict(
    :ann=>:isann,
    :interim=>:isinterim,
    :ex=>:isex,
    :annex=>:isannex,
    :ex40=>:isex40
  )

  TA2_PANEL_B_NAME::String = "ta2panelb-$(REPLICATION_TYPE[])-$(DIST_TYPE[])-$(OUT_SUFFIX[])"

  TA2_PANEL_B_Σ = clusteredΣ!

  panelb && tablea2panelb(hs,
    rowkeys=TA2_PANEL_B_ROWKEYS,
    rowlabel=TA2_PANEL_B_ROWLABEL,
    tablename=TA2_PANEL_B_NAME,
    colkeys = TA2_PANEL_B_COLKEYS,
    colfilter = TA2_PANEL_B_COLFILTER,
    colnames = TA2_PANEL_B_COLNAMES,
    Σfunc = TA2_PANEL_B_Σ
  )
end

function tablea2panelb(hs;
  rowkeys::Vector{Symbol} = error("rowkeys is required"),
  rowlabel::Dict = error("rowlabel is required"),
  colkeys::Vector{Symbol} = error("colkeys is required"),
  colnames::Vector{Vector{String}} = error("colnames is required"),
  colfilter::Dict = error("colfilter is required"),
  tablename::String = error("table name is required"),
  Fret::Symbol = DIST_RET,
  Σfunc::Function = error("Σfunc is required"),
  tablepath::String = TABLE_PATH,
  decimals::Int = 3,)

  #the colmns vary by the views
  local specdfs::Vector{SubDataFrame} = Vector{SubDataFrame}()

  #only changing the dependent variable with the regressions
  specs::FMSpecs = FMSpecs(length(colkeys))
  xnames::Vector{Symbol} = [:intercept, :yield]
  clusterspec::Vector{Symbol} = [:permno, :exdate]
  xspec::Symbol=:yield
  push!(specs, yspec=:ann,  xspec=xspec,  xnames=xnames,  clusterspec=clusterspec)
  push!(specs, yspec=:interim,  xspec=xspec,  xnames=xnames,  clusterspec=clusterspec)
  push!(specs, yspec=:ex,  xspec=xspec,  xnames=xnames,  clusterspec=clusterspec)
  push!(specs, yspec=:annex,  xspec=xspec,  xnames=xnames,  clusterspec=clusterspec)
  push!(specs, yspec=:ex40,  xspec=xspec,  xnames=xnames,  clusterspec=clusterspec)

  #println(describe(hs))

  specs.clusterspecs[3] = :permno
  #set the filters
  #=push!(specdfs, view(hs, hs[!, colfilter[colkeys[1]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[2]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[3]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[4]]], [Fret; :yield; :day; :permno; :eid]))
  push!(specdfs, view(hs, hs[!, colfilter[colkeys[5]]], [Fret; :yield; :day; :permno; :eid]))=#

  println(describe(hs))

  #generate the regression content (will display in Finometrics template)
  contentrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys[1:2])
  computeFMLMresults!(hs, specs#=, containsmissings=false=#)

  Ndescrows::Int = 2
  Ncols = length(colkeys)
  desccontent::Vector{Vector{String}} = (i->Vector{String}(undef, Ncols)).(1:Ndescrows)
  descrownames::Vector{String} = (s::Symbol->rowlabel[s]).(rowkeys[3:4])

  for c ∈ 1:Ncols
    desccontent[1][c] = "$(num2str(R²(specs.results[c]), decimals))"
    desccontent[2][c] = "$(specs.results[c].N)"
  end

  #=for (i,lin) ∈ enumerate(specs.results)
    println("Iteration $i: ")
    display(Σfunc(lin))
  end=#

  tt::String = textable(specs.results, Σfunc, [:intercept, :yield],
    colnames=colnames,
    scaling=[1.,1.],
    descrownames=descrownames,
    desccontent=desccontent,
    decimaldigits=3)

  writetextable(tt, path=tablepath, outname="$tablename.tex")

  return nothing
end

#creates the neccessary fields to make the hs dataframe
function makehs(dist::DataFrame; Fret= DIST_RET)
  local hs::DataFrame

  sdist::SubDataFrame = view(dist, dist.oneex .* dist.match,:)

  #use this to aggregate over the days where appropriate
  @inline function aggfuncret(v::AbstractVector{<:MFloat64})
    (sum((ismissing).(v)) == length(v)) && return missing

    return prod(skipmissing(1.0 .+ v)) - 1.0
  end

  #form the root of the aggregated dataframe
  aggnames::Vector{Symbol} = [:permno, :eid, :yield, :exdate, :distcd, :dclrdate]
  sdisttoagg::SubDataFrame = view(sdist, (!ismissing).(sdist.yield),aggnames)
  hs = aggregate(groupby(sdisttoagg, :eid), v->v[1])

  #fix the names
  @inline nofunction(s::Symbol) = Symbol(replace(string(s), "_function"=>""))
  rename!(hs, (s::Symbol-> s=>nofunction(s)).(names(hs)))

  #these will be the outcome variables
  N::Int =size(hs,1)
  hs.ann = Vector{MFloat64}(undef, N)
  hs.interim = Vector{MFloat64}(undef, N)
  hs.ex = Vector{MFloat64}(undef, N)
  hs.annex = Vector{MFloat64}(undef, N)
  hs.ex40 = Vector{MFloat64}(undef, N)

  #prepare for the aggregation and create tags for deleting or keeping rows
  hs.tokeep = trues(N)
  hsindex::Dict = Dict(r.eid=>r for r::DataFrameRow ∈ eachrow(hs))

  hs.tokeep = trues(N)
  @assert issorted(sdist, [:eid, :day])

  ssdists::GroupedDataFrame = groupby(sdist, :eid)
  @mpar for i ∈ 1:length(ssdists)
    ssdist::SubDataFrame = ssdists[i]
    r::DataFrameRow = hsindex[ssdist.eid[1]]

    #find the index of the exdate if possible
    dvec::Vector{Int} = digits(ssdist.distcd[1])[4:-1:1]
    exidx::MInt = searchsortedfirst(ssdist.day, 0)
    local preexprice::MFloat64
    if ssdist.date[exidx] == r.exdate
      preexprice = ssdist.day[exidx-1] == -1 ? ssdist.price[exidx-1] : -1.0
    else
      exidx = missing
      preexprice = missing
    end

    dclridx::MInt = searchsortedfirst(ssdist.date, r.dclrdate)
    (ssdist.date[dclridx] ≠ r.dclrdate) && (dclridx = missing)

    #get the distribution codes and filter using the HS criteria
    #also some basic integrity checks
    if ((dvec[1] ≠ 1) ||
      (dvec[2] ≠ 2) ||
      (dvec[3] > 3) || #allow quarterly dividends only
      ismissing(exidx) ||
      ismissing(preexprice) ||
      (preexprice < 5.0) ||
      ismissing(dclridx))

      r.tokeep = false
      continue
    end

    r.ann = ssdist[dclridx, Fret]
    r.ex = ssdist[exidx, Fret]

    if exidx - dclridx > 1
      r.interim = aggfuncret(ssdist[(dclridx+1):(exidx-1), Fret])
    end

    r.annex = aggfuncret(ssdist[(dclridx):(exidx), Fret])
    r.ex40 = aggfuncret(ssdist[(d::Int->(d>0) && (d≤40)).(ssdist.day), Fret])

  end

  hs = hs[hs.tokeep,:]

  #WARNING delete later
  filter!(r::DataFrameRow->year(r.exdate) ≤ 2011, hs)
  #println(describe(hs))

  return hs
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
function collect_data_sankey(datafile::String)
    data = Vector{Tuple{AbstractString, AbstractString, Int}}()

    open(datafile) do io
        while ! eof(io)
            n1, v, n2 = match(r"(.*)\s+\[(\d+)\]\s+(.*)", readline(io)).captures
            v = parse(Int, v)
            push!(data, (n1, n2, v))
        end
    end

    return data
end

function get_names_sankey(data::Vector{Tuple{AbstractString, AbstractString, Int}})
    names = Vector{AbstractString}()

    for line in data
        n1, n2, v = line
        push!(names, n1); push!(names, n2)
    end

    return unique(names)
end

function parse_sankey(datafile::String)
    data = collect_data_sankey(datafile)
    names = get_names_sankey(data)

    source = Vector{Union{Missing, Int}}()
    target = Vector{Union{Missing, Int}}()
    values = Vector{Union{Missing, Real}}()

    for line in data
        n1, n2, v = line
        push!(values, v)
        push!(source, findfirst(==(n1), names) - 1)
        push!(target, findfirst(==(n2), names) - 1)
    end

    return names, source, target, values
end

"""
    sankey(names, source, target, value)

Creates an `EChart` diagram displaying the path and frequencies between nodes.

## Methods
```julia
sankey(names::AbstractVector, source::AbstractVector{<:Union{Missing, Int}},
    target::AbstractVector{<:Union{Missing, Int}}, value::AbstractVector{<:Union{Missing, Real}};
sankey(datafile::AbstractString)
```

## Arguments
* `kwargs` : varargs to set any field of resulting `EChart` struct

## Notes

## Examples
```julia
name = ["Agricultural waste", "Bio-conversion", "Liquid", "Losses", "Solid", "Gas", "Biofuel imports",
    "Biomass imports", "Coal imports", "Coal"]
source = [0, 1, 1, 1, 1, 6, 7, 8, 10]
target = [1, 2, 3, 4, 5, 2, 4, 9, 9]
value = [124.729, 0.597, 26.862, 280.322, 81.144, 35, 35, 11.606, 63.965]
s = sankey(name, source, target, value)
```

Using the second method mentioned above--as per http://sankeymatic.com/build/'s convention--you can have a datafile (call it `input.txt`), where each line has the information
```
SOURCE [AMOUNT] TARGET
```
In this case.  `SOURCE` and `TARGET` are both `String`s, where `AMOUNT` is a `Real`.  For example:
```
Wages [2000] Budget
Interest [25] Budget
Budget [500] Taxes
Budget [450] Housing
Budget [310] Food
Budget [205] Transportation
Budget [400] Health Care
Budget [160] Other Necessities
```
And simply run
```julia
s = sankey("input.txt")
```
"""
function sankey(names::AbstractVector,
                source::AbstractVector{<:Union{Missing, Int}},
                target::AbstractVector{<:Union{Missing, Int}},
                value::AbstractVector{<:Union{Missing, Real}};
                kwargs...)

    ec = newplot(kwargs, ec_charttype = "sankey")
    ec.series = [SankeySeries(name = "Series 1", layout = "none", data = arrayofdicts(name = names),
        links = arrayofdicts(source = source, target = target, value = value))]

    return ec

end

sankey(datafile::AbstractString; kwargs...) = sankey(parse_sankey(datafile)...; kwargs...)


text:
##
import LinearAlgebra as linalg
import Plots as plt

include("KalmanFilter.jl")
import .KalmanFilter as kf

##
function get_data()
    # Observations
    x = [4000.0 4260.0 4550.0 4860.0 5110.0] # position [m]
    v = [ 280.0  282.0  285.0  286.0  290.0] # velocity [m/s]
    # x = [4000.0 4260.0 4550.0 4860.0 5110.0 5220.0 5115.0 4995.0 4871.0] # position [m]
    # v = [ 280.0  282.0  285.0  286.0  290.0  295.0  292.0  287.0  284.0] # velocity [m/s]
    X = [x; v]

    # Initial values
    a = 2.0 # acceleration
    Δt = 1.0 # Time step
    ΔP = [20.0; 5.0] # Process errors in process covariance matrix
    Δ = [25.0; 6.0] # Observation errors
    return X, a, Δ, ΔP, Δt
end

## Define the state space model

function define_model(X, a, Δ, ΔP, Δt)
    m, n = size(X)
    Q = zeros(m, m) # Process noise covariance matrix
    F = [1.0 Δt; 0.0 1.0] # Transition matrix
    B = [Δt^2/2.0 Δt]' # Control-input matrix
    u = 2.0 # # Control vector, acceleration
    R = linalg.Diagonal([Δ[1]^2 0.0; 0.0 Δ[2]^2]) # Sensor noise covariance matrix
    H = kf.eye(m) # Measurement matrix

    # Initial predicted process covariance matrix
    P = linalg.Diagonal([ΔP[1]^2 ΔP[1]*ΔP[2]; ΔP[2]*ΔP[1] ΔP[2]^2])
    return P, Q, R, F, H, B, u, m, n
end

## Data input
function data_input(X, n)
    X = [ [X[1,i], X[2,i]] for i in 1:n ]
    Xp = [ X[1] ] # Initialize predicted state
    Xk = copy(Xp) # Initialized new state matrix
    return Xp, Xk, X
end

## Smoothers
function run_smoothers_bts(Xk, F, Σ)
    k = length(Xk)
    sm = Vector{kf.RTSSmoother}(undef, k)
    sm[k] = kf.RTSSmoother(Xk[k], Σ[k])
    ksmoother = kf.RauchTungStriebel()
    while k > 1
        s = kf.kalman_smooth(ksmoother, sm[k].xk, Xk[k-1], Xk[k], Σ[k-1], Σ[k], F)
        k -= 1
        sm[k] = s
    end
    return sm
end

## UKF
function run_ukf(X, f, g, P, R, Q, α; filter=kf.CholeskySqrt())
    m, n = size(X)
    _, _, Xx = data_input(X, n)
    k = 1
    K = kf._kalman_gain(P, kf.eye(m), R)
    sukf = [ kf.UnscentedKF(
                Xx[k], Xx[k],
                P, K,
                zeros(2*m+1), zeros(2*m+1), zeros(m, 2*m+1),
                zeros((2*m)*2+1), zeros((2*m)*2+1), zeros(m, (2*m)*2+1),
                P, P,
                filter,
            )
    ]
    while k < n
        s = kf.unscented_kalman(
            Xx[k+1], sukf[k].xk, f, g, sukf[k].P, R, Q; sqrt_filter=filter, α=α,
        )
        push!(sukf, s)
        k += 1
    end
    return sukf
end

## Plots
function plot_kf(Xp, Xk, Xx, Δ)
    # Xp, Xk, Xx = reduce(hcat, Xp), reduce(hcat, Xk), reduce(hcat, Xx)
    p1 = plt.plot(Xp[1,:], m=:o, label="Predicted", legend=:topleft, ylabel="Position [m]", xlabel="Time [s]")
    plt.plot!(p1, Xx[1,:], m=:s, label="Measured", ribbon=Δ[1])
    plt.plot!(p1, Xk[1,:], m=:p, label="KF")
    p2 = plt.plot(Xp[2,:], m=:o, label="Predicted", legend=:topleft, ylabel="Velocity [m/s]", xlabel="Time [s]")
    plt.plot!(p2, Xx[2,:], m=:s, label="Measured", ribbon=Δ[2])
    plt.plot!(p2, Xk[2,:], m=:p, label="KF")
    p = plt.plot(p1, p2, size=(700,350))
    return p
end

function plot_kf_smooth(Xp, Xk, Xx, Xks, Δ)
    # Xp, Xk = reduce(hcat, Xp), reduce(hcat, Xk)
    # Xx, Xks = reduce(hcat, Xx), reduce(hcat, Xks)
    p1 = plt.plot(Xp[1,:], m=:o, label="Predicted", legend=:topleft, ylabel="Position [m]", xlabel="Time [s]")
    plt.plot!(p1, Xx[1,:], m=:s, label="Measured", ribbon=Δ[1])
    plt.plot!(p1, Xk[1,:], m=:p, label="KF")
    plt.plot!(p1, Xks[1,:], m=:h, label="Smooth")
    p2 = plt.plot(Xp[2,:], m=:o, label="Predicted", legend=:topleft, ylabel="Velocity [m/s]", xlabel="Time [s]")
    plt.plot!(p2, Xx[2,:], m=:s, label="Measured", ribbon=Δ[2])
    plt.plot!(p2, Xk[2,:], m=:p, label="KF")
    plt.plot!(p2, Xks[2,:], m=:h, label="Smooth")
    p = plt.plot(p1, p2, size=(700,350))
    return p
end

## Run KFs
X, a, Δ, ΔP, Δt = get_data()
P, Q, R, F, H, B, u, m, n = define_model(X, a, Δ, ΔP, Δt)

## Unscented Kalman Filter
f(x) = F*x .+ B .* u
g(x) = H*x
skf = run_ukf(X, f, g, P, R, Q, 1e-3; filter=kf.SVDSqrt())
Xp, Xk = kf.unpack(skf, [:xp, :xk])
p1 = plot_kf(Xp, Xk, X, Δ)

## Rauch Tung Striebel Smoother
Σ, Xkm = kf.unpack_matrix(skf, [:P, :xk])
skm = run_smoothers_bts(Xkm, F, Σ)
Xkm = kf.unpack(skf, :xk)
p1m = plot_kf_smooth(Xp, Xk, X, Xkm, Δ)


text:
using SpinDoctor
using Literate
using Documenter

DocMeta.setdocmeta!(SpinDoctor, :DocTestSetup, :(using SpinDoctor); recursive = true)

# Generate examples
examples = [
    "Solve BTPDE" => "solve_btpde",
    "Custom gradients" => "custom_gradients",
    "Compare ADCs" => "compare_adcs",
    "Matrix Formalism" => "matrix_formalism",
    "High Angular Resolution" => "hardi",
]
output = "generated"
for e ∈ examples
    e = joinpath(@__DIR__, "..", "examples", "$(e.second).jl")
    o = joinpath(@__DIR__, "src", output)
    Literate.markdown(e, o)
    # Literate.notebook(e, o)
    # Literate.script(e, o)
end

makedocs(;
    modules = [SpinDoctor],
    authors = "Syver Døving Agdestein and contributors",
    repo = "https://github.com/SpinDoctorMRI/SpinDoctor.jl/blob/{commit}{path}#{line}",
    sitename = "SpinDoctor.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://spindoctormri.github.io/SpinDoctor.jl",
        assets = String[],
    ),
    pages = [
        "Home" => "index.md",
        "Getting Started" => "getting_started.md",
        "Examples" => [e.first => joinpath(output, e.second * ".md") for e ∈ examples],
        "Theory" => [
            "Theory" => "theory/theory.md",
            "Bloch-Torrey PDE" => "theory/btpde.md",
            "Apparent Diffusion Coefficient" => "theory/adc.md",
            "Matrix Formalism" => "theory/matrix_formalism.md",
            "Discretization" => "theory/discretization.md",
        ],
        "Neurons" => "neurons.md",
        "API Reference" => "api.md",
    ],
)

deploydocs(; repo = "github.com/SpinDoctorMRI/SpinDoctor.jl", devbranch = "main")


text:
import FactCheck.roughly
roughly(A::Tuple; kvtols...) = (B::Tuple) -> begin
    length(A) != length(B) && return false
    return all(a->roughly(a[1]; kvtols...)(a[2]), zip(A,B) )
end


text:
using RecipesBase,Colors

RecipesBase.debug(false)

@recipe function f(::Type{T},C::T) where T <: AbstractCurve
    aspect_ratio --> 1.0
    plotdata(C)
end

@recipe function f(P::AbstractPath;vertices=false)
    delete!(plotattributes,:vertices)
    aspect_ratio --> 1.0

    @series begin
        vcat( [plotdata(c) for c in P]... )
    end

    if vertices
        @series begin
            label := ""
            markercolor --> :black
            markershape --> :circle
            seriestype := :scatter
            ComplexRegions.vertices(P)
        end
    end
end

@recipe function f(p::AbstractCircularPolygon)
    if isfinite(p)
        p.path
    else
        C = ComplexRegions.enclosing_circle(p,8)
        q = truncate(p,C)
        z,R = C.center,C.radius/3
        xlims --> (real(z)-R,real(z)+R)
        ylims --> (imag(z)-R,imag(z)+R)
        q.path
    end
end

@recipe function f(::Type{T},R::T) where T<:SimplyConnectedRegion
   if R isa ExteriorSimplyConnectedRegion
        seriestype := :shapecomplement
    else
        seriestype := :shape
    end

    C = R.boundary  # could be curve or path
    if C isa Line
        # need to fake with a polygon
        θ = angle(C)
        Polygon([C(0.5),(θ,θ+π)])
#    elseif C isa AbstractCurve
 #       ClosedPath(C)
    else
        C
    end
end

@recipe function f(R::ExteriorSimplyConnectedRegion)
    P = innerboundary(R)
    C = enclosing_circle(ClosedPath(P),8)
    zc = C.center
    r = 0.2*C.radius
    xlims --> [real(zc)-r,real(zc)+r]
    ylims --> [imag(zc)-r,imag(zc)+r]
    between(C,P)
end

@recipe function f(R::Union{ConnectedRegion,ExteriorRegion})
    p0 = outerboundary(R)
    p1 = innerboundary(R)
    z1 = [plotdata(p) for p in p1]
    if isnothing(p0)
        zc,rho = enclosing_circle(vcat(z1...),8)
        p0 = Circle(zc,rho)
        r = 0.2*rho
        xlims --> [real(zc)-r,real(zc)+r]
        ylims --> [imag(zc)-r,imag(zc)+r]
    end
    z0 = plotdata(p0)
    # This is not fast, but I don't see a shortcut...
    # find pairwise distances between components
    comp = [z1...,z0]
    n = length(comp)
    index = Array{Tuple}(undef,n,n)
    dist = fill(Inf,n,n)
    for i in 1:n
        for j in i+1:n
            ka,kb = argclosest(comp[i],comp[j])
            index[i,j] = (ka,kb)
            index[j,i] = (kb,ka)
            dist[i,j] = dist[j,i] = abs(comp[i][ka]-comp[j][kb])
        end
    end

    # find the hops between components
    unused = trues(length(z1))
    curr = n
    path = []
    while sum(unused) > 1
        u = findall(unused)
        k = argmin(dist[curr,u])
        next = u[k]
        push!(path,(curr,next))
        unused[next] = false
        curr = next
    end
    push!(path,(curr,findfirst(unused)))

    # accumulate into the last component
    p = path[1]
    idx = index[p...]
    data_in = z0[idx[1]]
    data_out = [ z0[idx[1]:-1:1]; z0[end:-1:idx[1]] ]
    for k = 1:length(path)-1
        a = idx[2]
        zc = comp[p[2]]
        p = path[k+1]
        idx = index[p...]
        b = idx[1]
        if a > b
            data_in = [data_in;zc[a:-1:b]]
            data_out = [data_out;zc[a:end];zc[1:b]]
        else
            data_in = [data_in;zc[a:-1:1];zc[end:-1:b]]
            data_out = [data_out;zc[a:b]]
        end
    end
    a = idx[2]
    zc = comp[p[2]]
    data = [ data_in; zc[a:-1:1]; zc[end:-1:a]; data_out[end:-1:1] ]

    aspectratio --> 1
    @series begin
        seriestype := :shape
        linealpha := 0
       data
    end

    seriestype := :path
    linecolor --> :black
    label := ""
    @series begin
        z0
    end
    for z in z1
        @series begin
            z
        end
    end
end


text:
R, x = RationalFunctionField(QQ, "x")
L = Localization(R, degree)

@testset "DegreeLocalization" begin

    @testset "Constructor" begin

      @test parent_type(KInftyElem{fmpq}) == KInftyRing{fmpq}
      @test elem_type(L) == KInftyElem{fmpq}
      @test function_field(L) == R

    end

    @testset "Canonicalisation" begin

      for i in 1:300
        a = rand(L, 0:10, -10:10)
        c = canonical_unit(a)
        @test (iszero(c) && !isunit(c)) || isunit(c)
      end
    end

    @testset "Valuation" begin

      @test valuation(L(x//(x^2 + 1))) == 1
      @test valuation(L(5//5)) == 0
      @test valuation(L((x^2 + 1)//(3//5*x^3))) == 1
      @test valuation(L((2//5)//(x^3 + 3x + 1))) == 3

      for i in 1:300
        a = rand(L, 0:10, -10:10)
        while iszero(a)
          a = rand(L, 0:10, -10:10)
        end
        @test valuation(a) >= 0
      end

    end

    @testset "Parent object call overloading" begin

      @test parent(L(3)) == L
      @test parent(L(1//x)) == L
      @test parent(L(fmpz(3))) == L
      @test parent(L(fmpq(3, 2))) == L
      @test parent(L()) == L
      @test L(6//3) == L(2)
      @test L(R(5)) == L(5)
      @test_throws ErrorException L(x)
    end

    @testset "GCDX" begin
      for i in 1:550
        a = rand(L, 0:10, -10:10)
        b = rand(L, 0:10, -10:10)
        (g, u, v) = gcdx(a, b)
        @test g == u*a + v*b
      end
    end

    @testset "GCD" begin
      for i in 1:550
        a = rand(L, 0:10, -10:10)
        b = rand(L, 0:10, -10:10)
        g = gcd(a, b)
        @test divides(a, g)[1] && divides(b, g)[1]
      end
    end

    @testset "Exact division" begin

      @test divides(L(1//x), L(1//x)) == (true, L(1))
      @test_throws DivideError divexact(L(1//(x + 1)), L())

      for i in 1:300
        a = rand(L, 0:10, -10:10)
        b = rand(L, 0:10, -10:10)
        d = divides(a, b)
        if d[1]
          @test a == d[2] * b
        end
      end
    end

    @testset "Inversion" begin

      @test inv(L((x + 1)//(x + 2))) == L((x + 2)//(x + 1))
      @test inv(L(23)) == L(1//23)
      @test_throws DivideError inv(L())
    end

    @testset "Binary operators" begin

      @test L((-1//3)//(x^2 - 3*x)) + L(-1//(x - 1//2)) == L((-x^2 + 8//3*x + 1//6)//(x^3 - 7//2*x^2 + 3//2*x))
      @test L((-1//3)//(x^2 - 3*x)) - L(-1//(x - 1//2)) == L((x^2 - 10//3*x + 1//6)//(x^3 - 7//2*x^2 + 3//2*x))
      @test L((-1//3)//(x - 3))*L(-1//(x - 1//2)) == (1//3)//(x^2 - 7//2*x + 3//2)

      @test L(18//2) + L(2//1) == L(11)
      @test L(18//3) - L(1) == L(5)
      @test L(32) * L(4) == L(128)
    end

    @testset "Basic manipulation" begin

      @test iszero(L()) == true
      @test iszero(L(0)) == true
      @test isone(L(1//x)) == false
      @test isunit(L((x + 1)//(x + 2))) == true

      @test iszero(L(4//3)) == false
      @test isone(L(13//13)) == true
      @test isunit(L()) == false
    end
end


text:
#####################

#TODO This code seems to assume logspace flows as floating point numbers. if so, enforca that on type F
function marginal_pass_up(circuit::UpFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    resize_flows(circuit, num_examples(data))
    cache = zeros(Float64, num_examples(data)) #TODO: fix type later
    marginal_pass_up_node(n::UpFlowΔNode, ::PlainXData) = ()

    function marginal_pass_up_node(n::UpFlowLiteral{O,F}, cache::Array{Float64}, data::PlainXData{E}) where {E <: eltype(F)} where {O,F}
        pass_up_node(n, data)
        # now override missing values by 1
        npr = pr(n)
        npr[feature_matrix(data)[:,variable(n)] .< zero(eltype(F))] .= 1
        npr .= log.( npr .+ 1e-300 )
        return nothing
    end

    function marginal_pass_up_node(n::UpFlow⋀Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 0
        for i=1:length(n.children)
            # pr(n) .+= pr(n.children[i])
            broadcast!(+, pr(n), pr(n), pr(n.children[i]))
        end
        return nothing
    end

    function marginal_pass_up_node(n::UpFlow⋁Cached, cache::Array{Float64}, ::PlainXData)
        pr(n) .= 1e-300
        for i=1:length(n.children)
            cache .= 0
            # broadcast reduced memory allocation, though accessing prob_origin(n).log_thetas[i] still allocates lots of extra memory,
            # it is proabably due to derefrencing the pointer
            broadcast!(+, cache, pr(n.children[i]), prob_origin(n).log_thetas[i])
            broadcast!(exp, cache, cache)
            broadcast!(+, pr(n), pr(n), cache)
        end
        broadcast!(log, pr(n), pr(n));
        return nothing
    end

    ## Pass Up on every node in order
    for n in circuit
        marginal_pass_up_node(n, cache, data)
    end
    return nothing
end

##### marginal_pass_down

function marginal_pass_down(circuit::DownFlowΔ{O,F}) where {O,F}
    resize_flows(circuit, flow_length(origin(circuit)))
    for n in circuit
        reset_downflow_in_progress(n)
    end
    for downflow in downflow_sinks(circuit[end])
        # initialize root flows to 1
        downflow.downflow .= one(eltype(F))
    end
    for n in Iterators.reverse(circuit)
        marginal_pass_down_node(n)
    end
end

marginal_pass_down_node(n::DownFlowΔNode) = () # do nothing
marginal_pass_down_node(n::DownFlowLeaf) = ()

function marginal_pass_down_node(n::DownFlow⋀Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
     # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for c in n.children
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n)
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n)
            end
        end
    end
end

function marginal_pass_down_node(n::DownFlow⋁Cached)
    # todo(pashak) might need some changes, not tested, also to convert to logexpsum later
    # downflow(n) = EF_n(e), the EF for edges or leaves are note stored
    for (ind, c) in enumerate(n.children)
        for sink in downflow_sinks(c)
            if !sink.in_progress
                sink.downflow .= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)) )
                sink.in_progress = true
            else
                sink.downflow .+= downflow(n) .* exp.(prob_origin(n).log_thetas[ind] .+ pr(origin(c)) .- pr(origin(n)))
            end
        end
    end
end

#### marginal_pass_up_down

function marginal_pass_up_down(circuit::DownFlowΔ{O,F}, data::XData{E}) where {E <: eltype(F)} where {O,F}
    @assert !(E isa Bool)
    marginal_pass_up(origin(circuit), data)
    marginal_pass_down(circuit)
end


text:
# bubble plots

#VegaLite
using VegaLite, VegaDatasets

dataset("cars") |>
@vlplot(:point, x=:Horsepower, y=:Miles_per_Gallon, size=:Acceleration)

# pyplot

using PyPlot

x = 100*rand(70)
y = 100*rand(70)
areas = 2000*rand(70)

fig = figure("pyplot_scatterplot",figsize=(10,10))
ax = PyPlot.axes()
PyPlot.scatter(x,y,s=areas,alpha=0.7, c=:pink)

xlabel("X")
ylabel("Y")
PyPlot.grid("on")

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
#!/usr/bin/env julia

# CUDAnative.jl version

using CUDAdrv, CUDAnative
using InteractiveUtils

kernel_dummy(ptr) = Base.pointerset(ptr, 0f0, Int(blockIdx().x), 8)

const len = 1000
const ITERATIONS = 5000

function benchmark(gpu_buf)
    @cuda threads=len kernel_dummy(Base.unsafe_convert(Ptr{Float32}, gpu_buf))
end

function main()
    cpu_time = Vector{Float64}(ITERATIONS)
    gpu_time = Vector{Float64}(ITERATIONS)

    gpu_buf = Mem.alloc(len*sizeof(Float32))
    @code_warntype benchmark(gpu_buf)
    for i in 1:ITERATIONS
        i == ITERATIONS-4 && CUDAdrv.Profile.start()

        gpu_tic, gpu_toc = CuEvent(), CuEvent()

        cpu_tic = time_ns()
        record(gpu_tic)
        benchmark(gpu_buf)
        record(gpu_toc)
        synchronize(gpu_toc)
        cpu_toc = time_ns()

        cpu_time[i] = (cpu_toc-cpu_tic)/1000
        gpu_time[i] = CUDAdrv.elapsed(gpu_tic, gpu_toc)*1000000
    end
    CUDAdrv.Profile.stop()
    Mem.free(gpu_buf)

    @printf("CPU time: %.2fus\n", median(cpu_time))
    @printf("GPU time: %.2fus\n", median(gpu_time))
end

main()


text:
include("odemethods.jl")

"""
given some ode model with the initial condition "ic" and the time span "tspan"

```
function ode(u, t)
    du = f(u,t)
    ...
return du
```
you can solve the ode using
```
sol = odesolve(ode, ic, tspan,  stepsize=1e-5; method="RK4")
```

sol.u will give you the numerical solution to the ode

"""
function odesolve(odefunction::Function, ic, tspan,  stepsize=1e-5; method="RK4")
    h = stepsize
    ts = collect(tspan[1]:h:tspan[2])
    t0= tspan[1]
    u0 = ic

    sol = [u0]
    tsol = [t0]
    if method=="RK4"
        for i in 1:round((tspan[2]-tspan[1])/h)
            u1 = RK4(odefunction, u0, h, t0)
            push!(sol, u1)
            push!(tsol, t0+h)
            u0 = u1
            t0 = t0+h
        end
    elseif method=="RK2"
        for i in 1:round((tspan[2]-tspan[1])/h)
            u1 = RK2(odefunction, u0, h, t0)
            push!(sol, u1)
            push!(tsol, t0+h)
            u0 = u1
            t0 = t0+h
        end
    else
        println("undefined method")
    end
    return (u=sol, t=tsol)

end


text:
using MCMCBenchmarks, Test, Random

@testset "LBA Tests " begin
    Random.seed!(2215)
    path = pathof(MCMCBenchmarks)
    include(joinpath(path, "../../Models/LBA/LBA_Models.jl"))
    include(joinpath(path, "../../Models/LBA/LinearBallisticAccumulator.jl"))
    Nd = 250
    v = [1.0,1.5,2.0]
    A = .8
    k = .2
    tau = .4
    Nreps = 1
    ProjDir = @__DIR__
    cd(ProjDir)
    samplers = (
      CmdStanNUTS(CmdStanConfig,ProjDir),
      AHMCNUTS(AHMClba,AHMCconfig),
      #DHMCNUTS(sampleDHMC)
      )
    options = (Nsamples=2000,Nadapt=1000,delta=.8,Nd=Nd)
    results = benchmark(samplers,simulateLBA,Nreps;options...)
    @test results[!,Symbol("v[1]_mean")][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ v[1] atol = .7
    @test results[!,Symbol("v[1]_mean")][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ v[1] atol = .7
    #@test results[!,Symbol("v[1]_mean")][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ v[1] atol = .7
    @test results[!,Symbol("v[2]_mean")][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ v[2] atol = .7
    @test results[!,Symbol("v[2]_mean")][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ v[2] atol = .7
    #@test results[!,Symbol("v[2]_mean")][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ v[2] atol = .7
    @test results[!,Symbol("v[3]_mean")][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ v[3] atol = .7
    @test results[!,Symbol("v[3]_mean")][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ v[3] atol = .7
    #@test results[!,Symbol("v[3]_mean")][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ v[3] atol = .7
    @test results[!,:A_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ A atol = .6
    @test results[!,:A_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ A atol = .6
    #@test results[!,:A_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ A atol = .6
    @test results[!,:k_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ k atol = .2
    @test results[!,:k_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ k atol = .2
    #@test results[!,:k_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ k atol = .2
    @test results[!,:tau_mean][results[!,:sampler] .== :AHMCNUTS,:][1] ≈ tau atol = .2
    @test results[!,:tau_mean][results[!,:sampler] .== :CmdStanNUTS,:][1] ≈ tau atol = .2
    #@test results[!,:tau_mean][results[!,:sampler] .== :DHMCNUTS,:][1] ≈ tau atol = .2
    # @test results[!,Symbol("v[1]_sampler_rhat")][1] ≈ 1 atol = .05
    # @test results[!,Symbol("v[2]_sampler_rhat")][1] ≈ 1 atol = .05
    # @test results[!,Symbol("v[3]_sampler_rhat")][1] ≈ 1 atol = .05
    # @test results[!,:A_sampler_rhat][1] ≈ 1 atol = .05
    # @test results[!,:k_sampler_rhat][1] ≈ 1 atol = .05
    # @test results[!,:tau_sampler_rhat][1] ≈ 1 atol = .05
    isdir("tmp") && rm("tmp", recursive=true)
end


text:
export
    Event,
    EventList,
    AllEventList,
    getevent,
    getlastevent,
    pushevent!,
    samplelocalpath,
    quadpathpoly,
    pathmean,
    tmax

"""
    AllowedVarType

Syntactic sugar to make it explicit what variable types are handled in the
local BPS (i.e, the type of the variables on each node of the factor graph).
See also, Event, EventList.
"""
const AllowedVarType = Union{Real, Vector{<:Real}}

"""
    Event{T <: AllowedVarType}

Encapsulation of a local event: a tuple made out of (position, velocity, time)
for a node in the local BPS algorithm.
See also: `EventList`.
"""
struct Event{T <: AllowedVarType}
    x::T
    v::T
    t::Real
end

"""
    EventList{T <: AllowedVarType}

Storage for all events that happened for a variable in the local BPS.
Each variable is stored explicitly (as opposed to storing a list of Event)
so that the times can be traversed efficiently. The xs, vs must all have the
same type within the list.
See also: `Event`, `AllEventList`.
"""
mutable struct EventList{T <: AllowedVarType}
    xs::Vector{T}     # coordinates
    vs::Vector{T}     # associated velocities
    ts::Vector{Real}  # associated times
end

function EventList(T::Type)
    @assert T <: AllowedVarType
    EventList(Vector{T}(undef, 0), Vector{T}(undef, 0), Vector{Real}(undef, 0))
end

"""
    AllEventList

Storage for all the `EventList` objects corresponding to a factor graph in a
local BPS setting. The element `evl` can be accessed by the variable index so
`a.evl[k]` will return the eventlist associated with variable `k`. Note that
the different `EventList` may have different types (in particular, different
dimensions). See also: `Event`, `EventList`.
"""
mutable struct AllEventList
    evl::Vector{EventList}
    types::Vector{DataType}
    function AllEventList(types::Vector{DataType})
        new([EventList(types[i]) for i ∈ 1:length(types)], types)
    end
end
AllEventList(T::DataType, nvars::Int) = AllEventList([T for i ∈ 1:nvars])

tmax(aev::AllEventList) = maximum(aev.evl[i].ts[end] for i ∈ 1:length(aev.evl))

"""
    getevent(evl, evidx)

Retrieve a specfic event (at index `evidx`) in an `EventList` object `evl` and
return it as an `Event` object.
See also: `Event`, `EventList`.
"""
function getevent(evl::EventList, evidx::Int)
    Event(evl.xs[evidx], evl.vs[evidx], evl.ts[evidx])
end

"""
    getlastevent(evl)

Retrieve the last `Event` in an `EventList` object `evl`.
See also: `getevent`, `EventList`.
"""
function getlastevent(evl::EventList)
    Event(evl.xs[end], evl.vs[end], evl.ts[end])
end

"""
    pushevent!(evl, ev)

Append an `Event` object `ev` to an `EventList` object `evl`. The `Event`
object is unpacked and pushed to the vectors encapsulated by `evl`.
See also: `EventList`.
"""
function pushevent!(evl::EventList, ev::Event)
    push!(evl.xs, ev.x)
    push!(evl.vs, ev.v)
    push!(evl.ts, ev.t)
    evl
end

"""
    getlocalsegment(evl, j)

Retrieve the `j`th event in an `EventList` object `evl` and return as a tuple
that event and the next event. If there is no next event (end of event list),
an "infinity event" is returned (in order to allow extrapolation).
See also: `samplelocalpath`.
"""
function getlocalsegment(evl::EventList, j::Int)
    if j == length(evl.ts)
        return (getevent(evl,j), Event(NaN, NaN, Inf))
    end
    (getevent(evl,j), getevent(evl,j+1))
end

"""
    samplelocalpath(evl, t)

Sample x_k(t) given the corresponding `EventList` object `evl`. The times `t`
must be a list of sorted times (or a single positive time).
See also: `getlocalsegment`.
"""
function samplelocalpath(evl::EventList, t::AllowedTimeType)
    @assert issorted(t) "list of times should be sorted"
    ti  = t[1]
    @assert ti >= 0.0 "the first time should be greater than 0.0"
    # find the starting local segment
    j = searchsortedlast(evl.ts, ti)
    # event and next event (describes the "segment")
    eva, evb = getlocalsegment(evl,j)
    # allocation of samples
    samples = Vector{typeof(eva.x)}(undef, length(t))
    # recursively work with local segments
    i=1
    while i<= length(t)
        ti = t[i]
        if ti <= evb.t
            # stay in the current segment
            samples[i] = eva.x + (ti - eva.t) * eva.v
            i += 1
        else
            # get the next segment
            j += searchsortedlast(evl.ts[j+1:end], ti)
            eva, evb = getlocalsegment(evl, j)
        end
    end
    samples
end
samplelocalpath(evl::EventList, t::Real) = samplelocalpath(evl, [t])[1]

"""
    quadpathpoly(evl, pol, T)

Integrate a polynomial on the elements along a path for a single node.
So if we are sampling in Rd and have phi(X)=Poly(xk), then we can integrate the
polynomial along the path. Note that we can't do the same thing if phi(X) mixes
2 or more nodes. See also syntactic sugar to perform the same operation on all
nodes.
"""
function quadpathpoly(evl::EventList, pol::Poly, T::Real)
    dim  = length(evl.xs[1])
    res  = zeros(dim)
    nseg = searchsortedlast(evl.ts,T)
    for segidx = 1:nseg-1
        # for each segment in the path
        eva, evb = getlocalsegment(evl, segidx)
        # segment time, initial point and velocity
        tau = evb.t - eva.t
        xa  = eva.x
        v   = (evb.x - xa) / tau
        # integrate along that segment
        for d = 1:dim
            res[d] += polyint(polyval(pol,Poly([xa[d],v[d]])))(tau)
        end
    end
    # last segment
    eva = getevent(evl, nseg)
    # characteristics
    tau = T - eva.t
    for d = 1:dim
        res[d] += polyint(polyval(pol,Poly([eva.x[d],eva.v[d]])))(tau)
    end
    return res / T
end
function quadpathpoly(aev::AllEventList, pol::Poly, T::Real)
    res = Vector{Vector{<:Real}}(undef, length(aev.evl))
    for (d, evl) in enumerate(aev.evl)
        res[d] = quadpathpoly(evl, pol, T)
    end
    return res
 end
pathmean(evl::EventList, T::Real) = quadpathpoly(evl, Poly([0.0,1.0]), T)
pathmean(aev::AllEventList, T::Real) = quadpathpoly(aev, Poly([0.0,1.0]), T)
pathmean(aev::AllEventList) = pathmean(aev, tmax(aev))