text:
# Author: Gabrio Rizzuti, rizzuti.gabrio@gatech.edu
# Date: August, 2020
# Copyright: Georgia Institute of Technology, 2020
#
# Run FWI on Sleipner data

## Module loading

using DrWatson
@quickactivate "UncertaintyQuantificationAVP"
using UncertaintyQuantificationAVP
using LinearAlgebra, PyPlot, Statistics, Optim
using Random; seed = 1; Random.seed!(seed)

## Loading model and synthetics

# Data
data_dict = wload(datadir("Sleipner", "Sleipner_data.bson"))
dat_n = data_dict[:dat_n]
σ2_n = data_dict[:σ2_n]
src_pars = data_dict[:src_pars]
abc_geom = data_dict[:abc_geom]
acc = data_dict[:acc_mod]

# Model
model_dict = wload(projectdir()*"/scripts/Sleipner/Sleipner_model.bson")
Mtrue = model_dict[:Mtrue]
Mbg = model_dict[:Mbg]

## Loss function

# Pre- and post-processing for log-posterior
nz = Mbg.geom.nz
dz = Mbg.geom.dz
mmin = 1f0/3000f0^2f0
mmax = Inf32
function proj_bounds(m::Array{Float32}, mmin, mmax)
    m_ = copy(m)
    m_[m_ .< mmin] .= mmin
    m_[m_ .> mmax] .= mmax
    return m_
end
prec(m::Array{Float32, 3}) = Model(Mbg.geom, proj_bounds(reshape(m, nz, :), mmin, mmax))
post(g::Model) = reshape(g.m, (nz, 1, :))

# Log-likelihood
F = seismod_fun(dat_n.pars, src_pars; abc_geom = abc_geom, acc = acc)
negloglike_fun(m::Model; compute_grad::Bool = false) = negLogLikelihood_gauss(dat_n, m, σ2_n, F; compute_grad = compute_grad)

# Log-prior
D = derivative_linop(Mbg.geom)
Mpr = Model(Mbg.geom, zeros(Float32, nz, 1))
Dx = D.eval(Mtrue-Mpr)
ϵ2 = 1f-3*norm(Dx.m, Inf)^2f0
# σ2_pr = 1f-2*sum(sqrt.(Dx.m.^2f0.+ϵ2))[1]
# neglogpr_fun(m::Model; compute_grad::Bool = false) = negLogPrior_L1smooth(m, σ2_pr, ϵ2; m0 = Mpr, D = D, compute_grad = compute_grad)
σ2_pr = 1f-3*norm(Dx)[1]^2f0
neglogpr_fun(m::Model; compute_grad::Bool = false) = negLogPrior_gauss(m, σ2_pr; m0 = Mpr, D = D, compute_grad = compute_grad)

# Log-posterior
neglogpost_fun(m::Model; compute_grad::Bool = false) = negLogPost(m, negloglike_fun, neglogpr_fun; compute_grad = compute_grad)

# Overloading log-posterior for tensor-type input+pre-/post-conditioning
function neglogpost_fun(m::Array{Float32, 3}; compute_grad::Bool = false)
    nb = size(m, 3)
    if compute_grad
        f, g_ = neglogpost_fun(prec(m); compute_grad = compute_grad)
        return sum(f)/nb, post(g_)
    else
        f = neglogpost_fun(prec(m); compute_grad = compute_grad)
        return sum(f)/nb
    end
end
function neglogpost_fun!(F, G::Union{Nothing, Array{Float32, 3}}, m::Array{Float32, 3})
    if G == nothing
        f = neglogpost_fun(m; compute_grad = false)
    else
        f, g = neglogpost_fun(m; compute_grad = true)
        G .= g
    end
    return f
end

## Optimization

M0 = deepcopy(Mbg)
method = LBFGS()
niter = 100
optimopt = Optim.Options(iterations = niter, store_trace = true, show_trace = true, show_every = 1)
result = optimize(Optim.only_fg!(neglogpost_fun!), reshape(M0.m, :, 1, 1), method, optimopt)
Mmap = prec(Optim.minimizer(result))
fval = Optim.f_trace(result)

## Saving results

invparams_dict = Dict(:σ2_n => σ2_n, :mmin => mmin, :mmax => mmax, :ϵ2 => ϵ2, :σ2_n => σ2_n, :σ2_pr => σ2_pr, :Mpr => Mpr, :M0 => M0)
results_dict = Dict(:Mmap => Mmap, :fval => fval)
wsave(datadir("Sleipner", "Sleipner_invpars.bson"), invparams_dict)
wsave(datadir("Sleipner", "Sleipner_MAP.bson"), results_dict)


text:
function get_action_indices(graph, fluent, level, set, except_set, all_actions)
    actions = [act for act in graph.acts[level] if act.name != :NoOp]
    # actions = all_actions
    action_indices = []
    if set == :add && except_set == :pre
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_eff && !(fluent in act.pos_prec)
                push!(action_indices, i)
            end
        end
    elseif set == :del && except_set == :pre
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.neg_eff && !(fluent in act.pos_prec)
                push!(action_indices, i)
            end
        end
    elseif set == :pre && except_set == :del
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_prec && !(fluent in act.neg_eff)
                push!(action_indices, i)
            end
        end
    else
        for i = 1:length(actions)
            act = actions[i]
            if fluent in act.pos_prec && fluent in act.neg_eff
                push!(action_indices, i)
            end
        end
    end
    return action_indices
end

function get_fluent_terminal_indices(graph, max_nf; loc=:init)
    indices = []
    not_ind = []
    T = graph.num_levels - 1
    if loc == :init
        for i=1:length(graph.initprops)
            ind = findall(x->x == graph.initprops[i], graph.props[0])[1]
            push!(indices, ind)
        end
    else
        for i=1:length(graph.goalprops)
            ind = findall(x->x == graph.goalprops[i], graph.props[T])[1]
            push!(indices, ind)
        end
    end
    not_ind = [i for i=1:max_nf if !(i in indices)]
    return indices, not_ind
end

function get_fluent_index(fluent, level, graph)
    for i=1:length(graph.props[level])
        if graph.props[level][i] == fluent
            return i
        end
    end
    return nothing
end

function optiplan_solve(domain_name, problem_name; max_levels=10)
    domain = load_domain(domain_name)
    problem = load_problem(problem_name)
    graph = create_graph(domain, problem; max_levels=max_levels)
    all_actions = mutex.get_all_actions(domain, problem)

    model = Model(Gurobi.Optimizer)
    set_silent(model)

    T = graph.num_levels
    max_nf = max([length(graph.props[i-1]) for i=1:T]...)
    max_na = max([length([act for act in graph.acts[i] if act.name != :NoOp]) for i=1:T-1]...)
    # max_na = length(all_actions)
    xs = []
    pa, pd, ad, dl, mn = 1,2,3,4,5

    #=
        VARIABLES
    =#
    [push!(xs, @variable(model, [i=1:max_nf, j=1:5], Bin, start=0)) for _=1:T]
    @variable(model, y[i=1:max_na, j=1:T-1], Bin, start=0)

    #=
        CONSTRAINTS
    =#
    #Eq. 2
    inits, notinits = get_fluent_terminal_indices(graph, max_nf; loc=:init)
    t=1
    #println("\ninits: ",inits)
    #println("notinits: ",notinits)
    for i in inits
        @constraint(model, xs[t][i, ad] == 1)
        @constraint(model, xs[t][i, pa] == 0)
        @constraint(model, xs[t][i, pd] == 0)
        @constraint(model, xs[t][i, mn] == 0)
        @constraint(model, xs[t][i, dl] == 0)
    end

    #Eq. 3
    for i in notinits
        @constraint(model, xs[t][i, ad] == 0)
        @constraint(model, xs[t][i, mn] == 0)
        @constraint(model, xs[t][i, pa] == 0)
        @constraint(model, xs[t][i, dl] == 0)
        @constraint(model, xs[t][i, pd] == 0)
    end

    #Eq. 4
    goals, _ = get_fluent_terminal_indices(graph, max_nf; loc=:goal)
    [@constraint(model, xs[T][i, ad] + xs[T][i, mn] + xs[T][i, pa] >= 1) for i in goals]

    # println("\nAll actions: ",all_actions)
    for t=1:T-1
        # println("\n\nLevel ",t)
        # println("\nAll actions: ",all_actions)
        for f = 1:length(graph.props[t-1])
            fluent = graph.props[t-1][f]
            # println("\nFluent: ",fluent)

            # Eq. 5
            act_inds = get_action_indices(graph, fluent, t, :add, :pre, all_actions)
            # println("Eq. 5 indices: ",act_inds)
            # println("Eq.5 actions: ",[all_actions[i] for i in act_inds])
            if !isempty(act_inds)
                @constraint(model, sum([y[i, t] for i in act_inds]) >= xs[t][f, ad] )
                # Eq. 6
                [@constraint(model, y[i, t] <= xs[t][f, ad]) for i in act_inds]
            end

            #Eq. 7
            act_inds = get_action_indices(graph, fluent, t, :del, :pre, all_actions)
            # println("Eq. 7 indices: ",act_inds)
            # println("Eq.7 actions: ",[all_actions[i] for i in act_inds])
            if !isempty(act_inds)
                @constraint(model, sum([y[i, t] for i in act_inds]) >= xs[t][f, dl] )
                # Eq. 8
                [@constraint(model, y[i, t] <= xs[t][f, dl]) for i in act_inds]
            end

            #Eq. 9
            act_inds = get_action_indices(graph, fluent, t, :pre, :del, all_actions)
            # println("Eq. 9 indices: ",act_inds)
            # println("Eq.9 actions: ",[all_actions[i] for i in act_inds])
            if !isempty(act_inds)
                @constraint(model, sum([y[i, t] for i in act_inds]) >= xs[t+1][f, pa] )
                # Eq. 10
                [@constraint(model, y[i, t] <= xs[t][f, pa]) for i in act_inds]
            end

            #Eq. 11
            act_inds = get_action_indices(graph, fluent, t, :pred, :anddel, all_actions)

            # println("Eq. 11 indices: ",act_inds)
            # println("Eq.11 actions: ",[all_actions[i] for i in act_inds])
            if !isempty(act_inds)
                @constraint(model, sum([y[i, t] for i in act_inds]) == xs[t][f, pd] )
            end
        end
    end
    for t=1:T-1
        # println("\n\nLevel ",t)
        # println("\nAll actions: ",all_actions)
        for f = 1:length(graph.props[t])
            fluent = graph.props[t][f]

            #Eq. 12
            @constraint(model, xs[t+1][f, ad] + xs[t+1][f, mn] + xs[t+1][f, dl] + xs[t+1][f, pd] <= 1)

            #Eq 13
            @constraint(model, xs[t+1][f, pa] + xs[t+1][f, mn] + xs[t+1][f, dl] + xs[t+1][f, pd] <= 1)

            #Eq. 14
            fprev = get_fluent_index(fluent, t-1, graph)
            if !isnothing(fprev)
            @constraint(model, xs[t+1][f, pa] + xs[t+1][f, mn] + xs[t+1][f, pd] <= xs[t][fprev, pa] + xs[t][fprev, ad] + xs[t][fprev, mn])
            end
        end
        # num_actions = length([act for act in graph.acts[t] if act.name != :NoOp])
        # @constraint(model, sum(y[1:num_actions, t ]) == 1)

        na = length([act for act in graph.acts[t] if act.name != :NoOp])
        for i=na+1:max_na
            @constraint(model, y[i,t] == 0)
        end
    end
    @objective(model, Min, sum(y))
    optimize!(model)
    sol = value.(y)
    # render_action(sol, graph, all_actions)
end

function render_action(y, graph, all_actions)
    T = graph.num_levels-1
    actions = []
    for t=1:T
        try
        ind = findall(x->x == 1.0, y[1:end,t])[1]
        act  = graph.acts[t][ind]
        # act = all_actions[ind]
        push!(actions, (act.name, act.args))
        catch
        end
    end
    # println(actions)
    return actions

end

ext: jl


text:
using Test

@testset "Spider.jl" begin
  @testset "$filename" for filename in (
    "basic_generate.jl",
  )
    println("Running $filename")
    include(filename)
  end
end


text:
isdefined(Base, :__precompile__) && __precompile__(false)

module ModelConstructors

    using DataFrames, Dates, Distributed, Distributions
    using ForwardDiff, Nullables, Printf, Random
    using LinearAlgebra, OrderedCollections, SpecialFunctions

    import Base.isempty, Base.<, Base.min, Base.max, Base.length, Base.iterate, Base.size
    import Base.sum, Base.getindex, Base.-
    import Distributions.log2π, Distributions.params, Distributions.mean, Distributions.std
    import Distributions.pdf, Distributions.logpdf, Distributions.Distribution
    import Distributions.rand, Distributions.Matrixvariate, Distributions.LinearAlgebra
    import LinearAlgebra.rank
    import SpecialFunctions.gamma
    import UnPack: unpack
    export
        # distributions_ext.jl
        BetaAlt, GammaAlt, RootInverseGamma, DegenerateMvNormal, DegenerateDiagMvTDist,
        MatrixNormal, <=, logpdf, init_deg_mvnormal,

        # settings.jl
        Setting, get_setting,

        # observables.jl
        Observable, PseudoObservable, check_mnemonics,

        # defaults.jl
        default_settings!, default_test_settings!,

        # abstractmodel.jl
        AbstractModel, description,
        n_states, n_states_augmented, n_shocks_exogenous, n_shocks_expectational,
        n_equilibrium_conditions, n_observables, n_parameters, n_parameters_regime_switching, n_parameters_steady_state,
        n_parameters_free, n_pseudo_observables, get_dict, get_key,
        spec, subspec, saveroot, dataroot,
        data_vintage, data_id, cond_vintage,
        logpath, workpath, rawpath, tablespath, figurespath, inpath, filestring_base, workpath,
        rawpath, tablespath, figurespath, logpath, savepath, inpath, filestring,

        # parameters.jl
        parameter, parameter_ad, Transform, NullablePrior, AbstractParameter,
        Parameter, ParameterVector, ScaledParameter,
        UnscaledParameter, SteadyStateParameter, SteadyStateParameterGrid,
        transform_to_real_line, transform_to_model_space,
        differentiate_transform_to_real_line, differentiate_transform_to_model_space,
        update, update!, transform_to_model_space, transform_to_real_line, Interval,
        ParamBoundsError, Untransformed, SquareRoot, moments,
        prior, get_untransformed_values, parameters2namedtuple,

        # distributions_ext.jl
        Uniform, Exponential, Normal, BetaAlt, GammaAlt, RootInverseGamma, pdf, logpdf, rand,
        DegenerateMvNormal, rank, length, DgenerateDiagMvTDist, mean, std,
        MatrixNormal, size, params, truncmean,

        AbstractVectorParameter, VectorParameter, VectorParameterVector, ScaledVectorParameter,
        UnscaledVectorParameter, Untransformed,

        # regimes.jl
        set_regime_val!, regime_val, set_regime_prior!, regime_prior, set_regime_fixed!, regime_fixed,
        set_regime_valuebounds!, regime_valuebounds, toggle_regime!,

        # statistics.jl
        prior, posterior, posterior!,

        # util.jl
        info_print, warn_print, println, print, @test_matrix_approx_eq, @test_matrix_eq2,

        # Generic Model
        GenericModel,

        # Linear Regression
        LinearRegression

    	const VERBOSITY = Dict(:none => 0, :low => 1, :high => 2)

    include("parameters.jl")
    include("distributions_ext.jl")
    include("abstractmodel.jl")
    include("settings.jl")
    include("observables.jl")
    include("defaults.jl")
    include("statistics.jl")
    include("util.jl")
    include("regimes.jl")
    include("generic_model.jl")
    include("linear_regression.jl")
end

max_stars_repo_name: FRBNY-DSGE/ModelConstructors

max_issues_repo_name: FRBNY-DSGE/ModelConstructors

max_forks_repo_name: FRBNY-DSGE/ModelConstructors


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

using Markdown
using InteractiveUtils

# ╔═╡ 54fc6142-ff38-11ea-2884-ed7abed903ca
using Pkg, DrWatson

# ╔═╡ 54fc9980-ff38-11ea-134f-cffec3a049e6
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
end

# ╔═╡ 550fc4ee-ff38-11ea-0bc0-8db9a84657b2
include(projectdir("models", "06", "m6.7s.jl"))

# ╔═╡ d6c8e622-ff39-11ea-38bb-513629d16e91
include(projectdir("models", "06", "m6.8s.jl"))

# ╔═╡ 8377e266-ff34-11ea-2963-b7dce70656f5
md"## Clip-06-20s.jl"

# ╔═╡ 54fd1328-ff38-11ea-22ee-03b3cbfb8b64
begin
	N = 1000
	df = DataFrame(
	  :h0 => rand(Normal(10,2 ), N),
	  :treatment => vcat(zeros(Int, Int(N/2)), ones(Int, Int(N/2))),
	  :M => rand(Bernoulli(), N)
	);

	d(i) = Binomial(1, 0.5 - 0.4 * df[i, :treatment] + 0.4 * df[i, :M])
	df.fungus = [rand(d(i), 1)[1] for i in 1:N]
	df.h1 = [df[i, :h0] + rand(Normal(5 + 3 * df[i, :M]), 1)[1] for i in 1:N]
end

# ╔═╡ 5509d90a-ff38-11ea-0563-a9cfe86b4e1e
md"##### Execute m6.7s & m6.8s."

# ╔═╡ 55103e56-ff38-11ea-0825-8b25859b9e32
begin
	(s, p) = plot_model_coef([m6_7s, m6_8s], [:a, :bt, :bf])
	p
end

# ╔═╡ 551932c2-ff38-11ea-28cd-5bfa4ad04d57
s

# ╔═╡ 5519a194-ff38-11ea-26c3-0140d22578e4
md"## End of clip-06-20s.jl"

# ╔═╡ Cell order:
# ╟─8377e266-ff34-11ea-2963-b7dce70656f5
# ╠═54fc6142-ff38-11ea-2884-ed7abed903ca
# ╠═54fc9980-ff38-11ea-134f-cffec3a049e6
# ╠═54fd1328-ff38-11ea-22ee-03b3cbfb8b64
# ╟─5509d90a-ff38-11ea-0563-a9cfe86b4e1e
# ╠═550fc4ee-ff38-11ea-0bc0-8db9a84657b2
# ╠═d6c8e622-ff39-11ea-38bb-513629d16e91
# ╠═55103e56-ff38-11ea-0825-8b25859b9e32
# ╠═551932c2-ff38-11ea-28cd-5bfa4ad04d57
# ╟─5519a194-ff38-11ea-26c3-0140d22578e4


text:
#########################
#   Geometric objects   #
#########################

### Declarations
if VERSION < v"0.6.0"
    include("geomtypes-0.5.jl")
else
    include("geomtypes.jl")
end

const TSP = TruncatedSquarePyramid
const Point = Vec

### Constructors
for S in (:Ellipsoid, :EllipticCylinder, :HollowCylinder, :Parallelepiped,
          :TriangularToroid, :Vec)
    @eval begin
        $S(a, b, c) = (t = promote(a, b, c); $S{eltype(t)}(t...))
    end
end

Cylinder(r, c) = (t = promote(r, c); Cylinder{eltype(t)}(t...))

SphericalCap(a, c) = (t = promote(a, c); SphericalCap{eltype(t)}(t...))

SquarePyramid(a, b) = (t = promote(a, b); m = 2b / a;
    SquarePyramid{eltype(t),typeof(m)}(t..., m))

SquarePyramid{T}(a, b, m::T) = (t = promote(a, b);
    SquarePyramid{eltype(t),T}(t..., m))

RectangularPyramid(a, b, c) = (t = promote(a, b, c); ma = 2c / a; mb = 2c / b;
    RectangularPyramid{eltype(t),typeof(ma)}(t..., ma, mb))

RectangularPyramid{T}(a, b, c, ma::T, mb::T) = (t = promote(a, b, c);
    RectangularPyramid{eltype(t),T}(t..., ma, mb))

TSP{R}(a, b, r::R) = (t = promote(a, b); m = 2b / a;
    TSP{eltype(t),R,typeof(m)}(t..., r, m))

TSP{R,S}(a, b, r::R, m::S) = (t = promote(a, b);
    TSP{eltype(t),R,S}(t..., r, m))

#############################
#   Coordinate transforms   #
#############################

@compat abstract type ShapePointTransformation end

immutable FunctionTransformation{F,T<:ShapePointTransformation}
    f::F
    t::T
end

for S in (:Ring, :Sphere, :Torus)

    T = Symbol(S, :PT)

    @eval begin
        immutable $T{S} <: ShapePointTransformation
            s::S
        end
    end
end

for S in (:Cube, :Cylinder, :Ellipsoid, :EllipticCylinder, :Parallelepiped,
          :RectangularPyramid, :SquarePyramid, :TriangularToroid)

    T = Symbol(S, :PT)

    @eval begin
        immutable $T{S,W} <: ShapePointTransformation
            s::S
            w::W
        end
    end
end

for S in (:HollowCylinder, :TSP, :SphericalCap)

    T = Symbol(S, :PT)

    @eval begin
        immutable $T{S,A,W} <: ShapePointTransformation
            s::S
            a::A
            w::W
        end
    end
end

### Constructors
for S in (:Ellipsoid, :EllipticCylinder, :Parallelepiped, :RectangularPyramid)

    T = Symbol(S, :PT)

    @eval begin
        $T(s::$S) = $T(s, s.a * s.b * s.c)
    end
end

TSPPT(s::TSP) = TSPPT(s, s.b * s.r, s.a^2 * s.b * s.r)

CubePT(s::Cube) = CubePT(s, s.a^3)

CylinderPT(s::Cylinder) = CylinderPT(s, s.c * s.r)

SquarePyramidPT(s::SquarePyramid) = SquarePyramidPT(s, s.a^2 * s.b)

TriangularToroidPT(s::TriangularToroid) = TriangularToroidPT(s, s.b * s.c)

SphericalCapPT(s::SphericalCap) =
    SphericalCapPT(s, 2 * (s.a / 2s.c)^2 + 2, s.c^3)

HollowCylinderPT(s::HollowCylinder) =
    HollowCylinderPT(s, s.R - s.r, s.c * (s.R - s.r))


text:
const time_names_table = ["ns", "µs", "ms", "s", "m", "h", "d", "y"]
const time_step_table =     [1000, 1000, 1000, 60,  60,  24, 365]

using Formatting
using Logging
using LightGraphs, SimpleWeightedGraphs

function format_seconds_readable(seconds, digits=2)::String
    fspec = FormatSpec(".$(digits)f")
    if 0.1 <= seconds <= 120.0 || seconds == 0.0
        return "$(fmt(fspec, seconds))s"
    end

    if seconds > 60.0
        # Up
        minutes = seconds / 60.0
        if minutes > 120.0
            hours = minutes / 60.0
            return "$(fmt(fspec, hours))h"
        end
        return "$(fmt(fspec, minutes))m"
    else
        # Down
        milliseconds = seconds * 1000.0
        if milliseconds < 0.1
            microseconds = milliseconds * 1000.0
            if microseconds < 0.1
                nanoseconds = milliseconds * 1000.0
                return "$(fmt(fspec, nanoseconds))ns"
            end
            return "$(fmt(fspec, microseconds))µs"
        end
        return "$(fmt(fspec, milliseconds))ms"
    end
    return string(seconds)
end

function format_ratio_readable(amount::Union{Int, Float64}, max::Union{Int, Float64} = 1, digits=2)::String
    if amount == 0
        if max == 0
            return "100%"
        else
            return "0%"
        end
    end
    fspec = FormatSpec(".$(digits)f")
    perc::Union{Float64, Rational{Int}} = 0.0
    if typeof(amount) == Float64
        perc = amount/max
    else
        max::Int
        perc = amount//max
    end
    perc *= 100
    if 1//5 <= perc || perc == 0
        return "$(fmt(fspec, perc))%"
    end

    perc *= 100

    if 1//5 <= perc
        return "$(fmt(fspec, perc))‰"
    end

    perc *= 100

    if 1//5 <= perc
        return "$(fmt(fspec, perc))‱"
    end

    perc *= 100

    if 1//5 <= perc
        return "$(fmt(fspec, perc))pcm"
    end

    perc *= 100

    return "$(fmt(fspec, perc))ppm"

end

function print_weighted_graph(graph::SimpleWeightedGraph, level=nothing, numberofvertices=nv(graph), numberofedges=ne(graph))
    if level === nothing
        println(numberofvertices)
    else
        @logmsg level numberofvertices
    end
    if level === nothing
        println(numberofedges)
    else
        @logmsg level numberofedges
    end
    line::Int = 0
    for (i, e) in enumerate(edges(graph))
        msg::String = "$line $(src(e) - 1) $(dst(e) - 1) $(Int(round(weight(e))))"
        if level === nothing
            println(msg)
        else
            @logmsg level msg
        end
        line += 1
    end
end

function print_graph(graph::AbstractGraph, level=nothing, weight=nothing, numberofvertices=nv(graph), numberofedges=ne(graph))
    if level === nothing
        println(numberofvertices)
    else
        @logmsg level numberofvertices
    end
    if level === nothing
        println(numberofedges)
    else
        @logmsg level numberofedges
    end
    if weight !== nothing
        if level === nothing
            println(weight)
        else
            @logmsg level weight
        end
    end
    line::Int = 0
    for (i, e) in enumerate(edges(graph))
        msg::String = "$line $(src(e) - 1) $(dst(e) - 1)"
        if level === nothing
            println(msg)
        else
            @logmsg level msg
        end
        line += 1
    end
end

max_stars_repo_name: Entze/jkMST

max_issues_repo_name: Entze/jkMST

max_forks_repo_name: Entze/jkMST


text:
# Autogenerated wrapper script for mold_jll for aarch64-linux-gnu-cxx03
export mold

using Zlib_jll
using OpenSSL_jll
JLLWrappers.@generate_wrapper_header("mold")
JLLWrappers.@declare_executable_product(mold)
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, OpenSSL_jll)
    JLLWrappers.@init_executable_product(
        mold,
        "bin/mold",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using .MonteCarloMeasurements

import .MonteCarloMeasurements: Particles

using OrderedCollections

function convert_a3d(a3d_array, cnames, ::Val{:particles};
    start=1, kwargs...)

  df = convert_a3d(a3d_array, Symbol.(cnames), Val(:dataframe))
  d = OrderedDict{Symbol, typeof(Particles(size(df, 1), Normal(0.0, 1.0)))}()

  for var in Symbol.(names(df))
    mu = mean(df[:, var])
    sigma = std(df[:, var])
    d[var] = Particles(size(df, 1), Normal(mu, sigma))
  end

  (; d...)

end

function Particles(df::DataFrame)

  d = OrderedDict{Symbol, typeof(Particles(size(df, 1), Normal(0.0, 1.0)))}()

  for var in Symbol.(names(df))
    d[var] = Particles(df[:, var])
  end

  (;d...)

end

export
  Particles

max_stars_repo_name: StanJulia/StanSample

max_issues_repo_name: StanJulia/StanSample

max_forks_repo_name: StanJulia/StanSample


text:
module MAB

"""
MAB.jl is a package for experimenting with multi arm bandit algorithms. This package comes
with different MAB algorithms, various arm models and a module to reproduce/evaluate your
own algorithms in different scenarios/experiments which are already available in literature.

This package is divied into three main modules
- [Algorithms](@ref)
- [Arm Models](@ref)
- [Experiments](@ref)
"""
MAB

"""
MAB.Algorithms module provides implementations of various multi-arm bandit algorithm and an
unified API to access those algorithms.
"""
module Algorithms
    using Distributions
    #--------------------------- Import Algorithms Here ---------------------------#
    include( "Algorithms/BanditAlgorithmBase.jl" )
    include( "Algorithms/UniformStrategy.jl" )
    include( "Algorithms/UCB.jl" )
    include( "Algorithms/epsGreedy.jl" )
    include( "Algorithms/EXP.jl" )
    include( "Algorithms/TS.jl" )
    include( "Algorithms/KLMANB.jl" )
    include( "Algorithms/SoftMax.jl" )
    include( "Algorithms/GradientBandits.jl" )

    export
        # Methods
        get_arm_index, update_reward!, reset!, info_str,
        show,
        # Agents
        BanditAlgorithmBase,
        epsGreedy, epsNGreedy,
        EXP3, EXP31, REXP3, EXP3IX,
        GradientBandit,
        KLMANB,
        SoftMax,
        TS, DynamicTS, OTS, TSNormal, dTS, dOTS, RestartTS, TSGaussPrior,
        UCB1, UCBNormal, DUCB, SWUCB, UCBV, BayesUCB, KLUCB, MUCB,
        UniformStrategy
end

"""
MAB.Arms module provides models for simulating various arms (reward sturctures) as well a
unified API to access them.
"""
module Arms
    using Distributions
    #--------------------------- Import Arm Models Here ---------------------------#
    include( "ArmModels/BanditArmBase.jl" )
    include( "ArmModels/Bernoulli.jl" )
    include( "ArmModels/Exponential.jl" )
    include( "ArmModels/Normal.jl" )
    include( "ArmModels/Beta.jl" )
    include( "ArmModels/Sinusoidal.jl" )
    include( "ArmModels/Pulse.jl" )
    include( "ArmModels/Square.jl" )
    include( "ArmModels/Variational.jl" )

    # export
        # Export only methods, not Arm types to avoid collision with Distributions
        # pull!, tick!, reset!
end

"""
MAB.Experiments provides a way to reproduce results from publish/preprint articles well known
in MAB community. You will be able to find implementations of experiments from popular papers
and a way to evaluate your algorithms in those environments.
"""
module Experiments
    using Distributions
    using Requires
    # import Base: run
    # using Plots
    @require Plots using Plots
    # info( "Some functionalities of MAB.Experiments will not be available until you load Plots.jl",
    #         prefix = "MAB: " )
    using ..Algorithms
    using ..Algorithms: make_agents_with_k  # Util for making agents inside Experiments
    using ..Arms
    # import
    #--------------------------- Import Experiments Here --------------------------#
    include( "Experiments/BanditExpBase.jl" )
    # include( "Experiments/utils.jl" )
    include( "Experiments/Compare.jl" )
    include( "Experiments/Sutton2017/Sutton2017.jl" )

    export
        # Methods
        run,
        # Experiments
        Sutton2017
end

    using .Algorithms
    using .Arms
    using .Experiments
    export
        # Export algorithms
        Algorithms,
        # Methods
        get_arm_index, update_reward!, reset!, info_str, show,
        # Agents
        BanditAlgorithmBase,
        epsGreedy, epsNGreedy,
        EXP3, EXP31, REXP3, EXP3IX,
        GradientBandit,
        KLMANB,
        SoftMax,
        TS, DynamicTS, OTS, TSNormal, dTS, dOTS, RestartTS, TSGaussPrior,
        UCB1, UCBNormal, DUCB, SWUCB, UCBV, BayesUCB, KLUCB, MUCB,
        UniformStrategy,
        # Export Arm models
        Arms,
        # Methods
        # pull!, tick!, reset!,  # Causes ambiguation with Algorithms.reset!
        # Export Experiments
        Experiments,
        # Methods
        run,
        # Experiments
        Sutton2017

end # module


text:
using ColorTypes
using ColorVectorSpace
using Documenter
using FixedPointNumbers
using OffsetArrays
import AxisArrays: AxisArray # No, we don't want to import `AxisArrays.axes`
using Statistics
using Test
using TiffImages

if VERSION >= v"1.4"
    DocMeta.setdocmeta!(TiffImages, :DocTestSetup, :(using TiffImages); recursive=true)
    doctest(TiffImages)
end

_wrap(name) = "https://github.com/tlnagy/exampletiffs/blob/master/$name?raw=true"

if VERSION >= v"1.6.0"
    using Downloads
    get_example(x) = Downloads.download(_wrap(x))
else
    get_example(x) = download(_wrap(x))
end

@testset "Tag loading" begin
    include("tags.jl")
end

@testset "Gray with ExtraSample" begin
    filepath = get_example("house.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (512, 512)
    @test eltype(img) == GrayA{N0f8}
    @test img[50,50] == GrayA{N0f8}(0.804, 1.0) # value from ImageMagick.jl
    img[50:300, 50:150] .= 0.0
    @test img[50, 50] == GrayA{N0f8}(0.0, 1.0)

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "MRI stack" begin
    filepath = get_example("mri.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (128, 128, 27)
    @test eltype(img) == RGB{N0f16}

    # TODO: inefficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt == img.data
    @test img_cvt !== img.data

    @testset "Multiple planes" begin
        # test slice based access
        @test red(img[92, 21, 2]) == 0.74118N0f16

        # make sure IFDs are included if whole slice is grabbed
        @test length(img[:, :, 2].ifds) == 1
        @test length(img[:, :, 2:5].ifds) == 4

        # make sure that subslicing in XY drops IFD info
        @test typeof(img[:, 50:60, 27]) == Array{RGB{Normed{UInt16,16}},2}
    end
end

@testset "Floating point RGB image" begin
    filepath = get_example("spring.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (619, 858)
    @test eltype(img) == RGB{Float16}

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Packbits image" begin
    filepath = get_example("coffee.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (378, 504)
    @test eltype(img) == Gray{N0f8}

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Adobe Deflate image" begin
    filepath = get_example("underwater_bmx.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (773, 1076)
    @test eltype(img) == RGB{N0f8}

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data

    # Dirty patch the TIFF_COMPRESSION tag from
    # COMPRESSION_ADOBE_DEFLATE (8) to COMPRESSION_DEFLATE (32946).
    # Those are in fact exactly the same so nothing else needs to be
    # adjusted.
    data = read(filepath)
    @assert data[1063947:1063948] == [0x08, 0x00]
    data[1063947:1063948] .= reinterpret(UInt8, [UInt16(32946)])
    _, io = mktemp()
    write(io, data)
    seekstart(io)
    img = TiffImages.load(io)
    @test size(img) == (773, 1076)
    @test eltype(img) == RGB{N0f8}
end

@testset "Bilevel image" begin
    filepath = get_example("capitol.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (378, 504)
    @test eltype(img) == Gray{Bool}

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Striped bilevel image" begin
    filepath = get_example("capitol2.tif")
    img = TiffImages.load(filepath)
    @test size(img) == (378, 504)
    @test eltype(img) == Gray{Bool}

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Signed integer type" begin
    filepath = get_example("4D-series.ome.tif")
    img = TiffImages.load(filepath)

    @test size(img) == (167, 439, 35)
    expected_rng = reinterpret.(Q0f7, Int8.((-1, 96)))
    @test extrema(img[:, :, 1]) == expected_rng

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Discontiguous striped image, Issue #38" begin
    filepath = get_example("julia.tif")
    img = TiffImages.load(filepath)

    @test size(img) == (300, 500)
    # if discontiguous striping is broken then the garbage padding data will
    # leak into the actual image
    @test all(img[3, 1:50] .== RGB{N0f8}(1, 1, 1))

    # Efficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt === img.data
end

@testset "Issue #12" begin
    @testset "Big endian striped file" begin
        filepath = get_example("flagler.tif")
        img = TiffImages.load(filepath)

        # verify that strip offsets are correctly bswapped for endianness
        img_stripoffsets = Int.(img.ifds[1][TiffImages.STRIPOFFSETS].data)
        @test img_stripoffsets == [8, 129848, 259688, 389528]

        # Efficient convert method
        img_cvt = convert(Array{eltype(img), ndims(img)}, img)
        @test img_cvt === img.data
    end

    @testset "Little endian striped file" begin
        filepath = get_example("house.tif")
        img = TiffImages.load(filepath)

        # verify that strip offsets are correctly bswapped for endianness
        img_stripoffsets = Int.(img.ifds[1][TiffImages.STRIPOFFSETS].data)
        @test issorted(img_stripoffsets)

        # Efficient convert method
        img_cvt = convert(Array{eltype(img), ndims(img)}, img)
        @test img_cvt === img.data
    end
end

@testset "Mmap" begin
    filepath = get_example("julia.tif")
    img = TiffImages.load(filepath, mmap=true)
    @test size(img) == (300, 500, 1)
    @test all(img[3, 1:50] .== RGB{N0f8}(1, 1, 1))
    # force close the stream behind the file to see if it's properly reopened
    close(img.data.file.io)
    @test all(img[3, 1:50] .== RGB{N0f8}(1, 1, 1))

    # TODO: inefficient convert method
    img_cvt = convert(Array{eltype(img), ndims(img)}, img)
    @test img_cvt == img.data
    @test img_cvt !== img.data
end

@testset "Writing" begin
    include("writer.jl")
end

@testset "Interpreting IFD layouts" begin
    include("layouts.jl")
end


text:
# ---------------------------------------------------------------
function run_cmd(cmd::Cmd; ios = [stdout], detach = false)

    # run
    _out = Pipe()
    cmd = pipeline(Cmd(cmd; detach), stdout = _out, stderr = _out)
    proc = run(cmd, wait = false)
    pid = _try_getpid(proc)
    wait(proc)

    # out
    close(_out.in)
    out = read(_out, String)

    # print
    for io in ios
        _append(io, out)
    end

    return (;pid, out)
end

# ---------------------------------------------------------------
function _short_run(cmd::Cmd;
        stdout_log::AbstractString = tempname(), stderr_log::AbstractString = tempname(),
        stdout_tee_ios = [stdout], stderr_tee_ios = [stderr],
        append = false,
        kwargs... # to ignore foreing kwargs
    )

    # empty
    !append && _empty!.([stderr_log, stdout_log])

    # run
    cmd = pipeline(cmd, stdout=stdout_log, stderr=stderr_log; append)
    proc = run(cmd; wait = false)
    pid = _try_getpid(proc)
    wait(proc)

    # tee
    for (logfile, ios) in [
            (stdout_log, stdout_tee_ios),
            (stderr_log, stderr_tee_ios)
        ]
        for io in ios
            !append && _empty!(io)
            _append(io, read(logfile, String))
        end
    end

    return pid
end

# ---------------------------------------------------------------
function _long_run(cmd::Cmd;
        stderr_log = tempname(), stdout_log = tempname(),
        stdout_tee_ios = [stdout], stderr_tee_ios = [stderr],
        printlk = ReentrantLock(), timeout = time() + 1e9,
        savetime = 60.0, # To wait for flushing
        append = false,
        upfrec = 1.0,
        finish_time = Ref{Float64}(time() + timeout), # define tee tout
        kwargs... # to ignore foreing kwargs
    )

    # empty
    !append && _empty!.([stderr_log, stdout_log])

    @sync begin

        proc = run(pipeline(cmd, stdout=stdout_log, stderr=stderr_log; append), wait=false)
        pid = _try_getpid(proc)

        # tee
        stdout_tsk = @async tee_file(stdout_tee_ios, stdout_log; finish_time, printlk, upfrec)
        stderr_tsk = @async tee_file(stderr_tee_ios, stderr_log; finish_time, printlk, upfrec)

        # ending
        wait(proc)
        finish_time[] = time() + savetime
        wait(stdout_tsk)
        wait(stderr_tsk)

        return pid
    end

end

# ---------------------------------------------------------------
function tee_run(cmd::Cmd; long::Bool = false, runkwargs...)
    long ? _long_run(cmd; runkwargs...) : _short_run(cmd; runkwargs...)
end

# ---------------------------------------------------------------
function run_bash(
        src_lines::Vector{String};
        startup::Vector{String} = String[],
        buff_file = tempname(),
        run_fun = run_cmd,
        rm_buff = true,
        runkwargs...
    )

    touch(buff_file)
    chmod(buff_file, 0o755)

    src = join(filter(!isempty, [
        startup;
        rm_buff ? """rm -f '$(buff_file)';""" : "";
        src_lines
    ]), "\n")

    write(buff_file, src)

    run_fun(`bash -c $(buff_file)`; runkwargs...)
end

run_bash(src::String; run_fun = run_cmd, runkwargs...) = run_fun(Cmd(`bash -c $(src)`); runkwargs...)


text:
using DynamicalSystems
using Plots
plotly()

# Plot a Van der Pol oscillator: \frac{d^2 x}{dt^2} - \mu(1-x^2) \frac{dx}{dt} + x = 0
# BUT with time varying parameters

X0 = [1.0, 1.0]
T = 2000.0

# --------------------------------- Generate a parameter profile --------------------------------- #
μ = unitStep(1000.0).*20.0 # A step from 0 to 20 at t = 1000
# ------------------------------------------------------------------------------------------------ #

# ----------------------------- Define an EOM, and possibly a Jacobian --------------------------- #
function vanderpol(X, μ, t)
    dX1 = X[2]
    dX2 = μ(t).*(1.0-X[1].^2).*X[2] - X[1]
    return SVector{2}(dX1, dX2)
end
# ------------------------------------------------------------------------------------------------ #

# ----------------------------------- Create a dynamical system ---------------------------------- #
ds = ContinuousDynamicalSystem(vanderpol, X0, μ)
# ------------------------------------------------------------------------------------------------ #

# --------------------------------------------- Solve -------------------------------------------- #
@time data = trajectory(ds, T; dt=0.01, tstops=sort(collect(μ.d)))
# ------------------------------------------------------------------------------------------------ #

plot(data[:, 1], seriestype=:path)
plot!(title="Van der Pol, parameter steps from 0 to 20")


text:
### Methods ###

function CerebellumModel(config::Dict)
  numZones = get(config, "numZones", 4)
  szGr = get(config, "szGr", 32)
  szPk = get(config, "szPk", 8)
  szDc = get(config, "szDc", 4)
  szIo = get(config, "szIo", 4)

  agent = Agent()

  # TODO: Split out Granule neurons for sharing purposes
  # TODO: Delays
  for i = 1:numZones
    agent_Zone = Agent()

    # Granule neurons
    addnode!(agent_Zone, GenericNeurons(szGr); name="Gr")

    # Purkinje neurons
    addnode!(agent_Zone, GenericNeurons(szPk); name="Pk")

    # Deep Cerebellar neurons
    addnode!(agent_Zone, GenericNeurons(szDc); name="Dc")

    # Inferior Olive neurons
    addnode!(agent_Zone, GenericNeurons(szIo); name="Io")

    # Gc -> Pk
    # TODO: Masking
    # TODO: Increasing individual delays
    addnode!(agent_Zone, GenericSynapses(szGr, szPk); name="Gr_Pk")
    addedge!(agent_Zone, agent_Zone["Gr_Pk"], (
      (agent_Zone["Gr"], :input),
      (agent_Zone["Pk"], :output)
    ))

    # Pk -(Inh)> Dc
    addnode!(agent_Zone, GenericSynapses(szPk, szDc; outputMask=-1); name="Pk_Dc")
    addedge!(agent_Zone, agent["Pk_Dc"], (
      (agent_Zone["Pk"], :input),
      (agent_Zone["Dc"], :output)
    ))

    # Dc -(Inh)> Io
    addnode!(agent_Zone, GenericSynapses(szDc, szIo; outputMask=-1); name="Dc_Io")
    addedge!(agent_Zone, agent["Dc_Io"], (
      (agent_Zone["Dc"], :input),
      (agent_Zone["Io"], :output)
    ))

    # Io -> Dc
    addnode!(agent_Zone, GenericSynapses(szIo, szDc); name="Io_Dc")
    addedge!(agent_Zone, agent["Io_Dc"], (
      (agent_Zone["Io"], :input),
      (agent_Zone["Dc"], :output)
    ))

    # Io -> Io
    addnode!(agent_Zone, GenericSynapses(szIo, szIo); name="Io_Io")
    addedge!(agent_Zone, agent["Io_Io"], (
      (agent_Zone["Io"], :input),
      (agent_Zone["Io"], :output)
    ))

    agent["$i"] = agent_Zone
  end

  return agent
end


text:
using StatsFuns: normcdf

export Probit, ProbitMeta

struct Probit end

struct ProbitMeta
    p :: Int
end

getp(meta::ProbitMeta) = meta.p

ProbitMeta(; p = 32) = ProbitMeta(p)

@node Probit Stochastic [ out, in ]

default_meta(::Type{ Probit }) = ProbitMeta(32)

default_functional_dependencies_pipeline(::Type{ <: Probit }) = RequireInboundFunctionalDependencies((2, ), (vague(NormalMeanPrecision), ))

default_interface_local_constraint(::Type{ <: Probit }, edge::Val{ :in })  = MomentMatching()
default_interface_local_constraint(::Type{ <: Probit }, edge::Val{ :out }) = Marginalisation()

@average_energy Probit (q_out::Union{PointMass, Bernoulli}, q_in::UnivariateNormalDistributionsFamily, meta::ProbitMeta) = begin

    # extract parameters
    p = mean(q_out)
    m, v = mean_var(q_in)

    # specify function
    h = (x) -> -p*log(normcdf(x)) - (1-p)*log(normcdf(-x))

    # calculate average average energy (default of 32 points)
    gh_cubature = GaussHermiteCubature(getp(meta))
    U = 0.0
    tmp = sqrt(2*v)
    for k = 1:getp(meta)
        U += gh_cubature.witer[k] * h( gh_cubature.piter[k] * tmp + m)
    end
    U /= sqrt(pi)

    # return average energy
    return U

end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    descendants, ancestors, roots, leafs,
    topologicalsort, reversetopologicalsort,
    transitive_reduction

descendants(graph::DirectedGraph, idx::Int) = reachablenodes(graph, idx)
ancestors(graph::DirectedGraph, idx::Int) = reversereachablenodes(graph, idx)

roots(graph::OrderedDiGraph) = [i for i in 1:nodecount(graph) if indegree(graph, i) == 0]
roots(graph::DirectedGraph) = [i for i in nodeset(graph) if indegree(graph, i) == 0]

leafs(graph::OrderedDiGraph) = [i for i in 1:nodecount(graph) if outdegree(graph, i) == 0]
leafs(graph::DirectedGraph) = [i for i in nodeset(graph) if outdegree(graph, i) == 0]

"""
    topologicalsort(graph::DirectedGraph) -> Vector{Int}
    topologicalsort(graph::DirectedGraph, root) -> Vector{Int}

Run topological sort to obtain an array of sorted nodes.

If `root` is given, the method is applied to only the downstream nodes of `root`.
"""
function topologicalsort(graph::DirectedGraph)
    result = Int[]
    stack = [i for i in nodeset(graph) if indegree(graph, i) == 0]
    used = Set{Int}()
    while !isempty(stack)
        n = pop!(stack)
        push!(result, n)
        for (ninc, nadj) in outneighbors(graph, n)
            ninc in used && continue
            push!(used, ninc)
            if isempty(setdiff(in_incidences(graph, nadj), used))
                push!(stack, nadj)
            end
        end
    end
    edgecount(graph) == length(used) || throw(ErrorException("cycle found"))
    return result
end

function topologicalsort(graph::DirectedGraph, root)
    sub = nodesubgraph(graph, union(descendants(graph, root), root))
    return topologicalsort(sub)
end

function reversetopologicalsort(graph::DirectedGraph)
    result = Int[]
    stack = [i for i in nodeset(graph) if outdegree(graph, i) == 0]
    used = Set{Int}()
    while !isempty(stack)
        n = pop!(stack)
        push!(result, n)
        for (ninc, nadj) in inneighbors(graph, n)
            ninc in used && continue
            push!(used, ninc)
            if isempty(setdiff(outincidences(graph, nadj), used))
                push!(stack, nadj)
            end
        end
    end
    edgecount(graph) == length(used) || throw(ErrorException("cycle found"))
    return result
end

function reversetopologicalsort(graph::DirectedGraph, leaf)
    sub = nodesubgraph(graph, union(ancestors(graph, leaf), leaf))
    return topologicalsort(sub)
end

"""
    longestpathedges(graph::DirectedGraph, root) -> ValueIterator

Return edges belong to the longest path from the given root node.
"""
function longestpathedges(graph::DirectedGraph, root)
    sub = nodesubgraph(graph, union(descendants(graph, root), root))
    dist = Dict{Int,Int}()
    inedges = Dict{Int,Int}()
    for n in nodeset(sub)
        dist[n] = n == root ? 0 : typemin(Int)
    end
    for n in topologicalsort(sub)
        for (inc, adj) in outneighbors(sub, n)
            if dist[adj] < dist[n] + 1
                dist[adj] = dist[n] + 1
                inedges[adj] = inc
            end
        end
    end
    return values(inedges)
end

"""
    transitive_reduction(graph::DirectedGraph) -> Set{Int}

Return transitive reduction as a subset of the DAG edges.
"""
function transitive_reduction(graph::DirectedGraph)
    edges = Set{Int}()
    for v in nodeset(graph)
        union!(edges, longestpathedges(graph, v))
    end
    return edges
end


text:
using SpecialFunctions
@testset "Two Atom scattering" begin
	# ----------------------------------------
	# Exact result for scattering a Gaussian mode
  σ = 1 + rand() # smaller σ need more careful numerical treatment
  function twoAtomOut(x)
    norm = √(√(π)*σ)

    return (1 + σ^2) * exp(-x^2/(2σ^2)) / norm -
      (norm / √(2)) * (x + 2 + σ^2 / 2) * exp((4x + σ^2)/8) * erfc((2x + σ^2) / (2 * √(2) * σ))
  end

  problem = WaveguideProblem(
    TwoLevelChain(2),
    WavePacket(GaussMode(σ = σ), Fock(1)),
    Mode(x -> twoAtomOut(-x); compression = Exponential(σ = σ), reltol = 1e-10),
    (-5σ, 20σ)
  )
  sol = ptrace(solve(problem, reltol = 1e-10)[2][end], [1, 2, 3])

  # Transfer chance
  n̂ = number(basis(sol))
  @test real( 1 - expect(n̂, sol) ) < 1e-5
end

@testset "Fock State scattering - Kiilerich" begin
	# ----------------------------------------
	# Recalculate Figure 2 of PRL 123, 123604 (2019)
	function outputModeFunctionFig2(t, γ = 1.0, τ = 4.0)
		abs(t-τ) > 10 && return 0.0
	  return exp(-(t - τ)^2/2) -
	    γ * √(π/2) * exp(-γ * (t - τ) / 2 + γ^2 / 8) * erfc((-2(t - τ) + γ) / (2 * √(2)))
	end

	WP = WavePacket(GaussMode(τ = 4.0), Fock(1))
	outputMode = Mode(t -> outputModeFunctionFig2(t, 1.0, 4.0), compression = :algebraic)

	problemA = WaveguideProblem(TwoLevelChain(1), WP, outputMode, 13.0)
	solA     = solve(problemA, reltol = 1e-10)[2]
	b        = basis(solA[end])

	n̂ᵢ, n̂ₒ = embed(b, 1, number(b.bases[1])), embed(b, 3, number(b.bases[3]))
	σ⁺σ⁻   = embed(b, 2, transition(NLevelBasis(2), 2, 2))
	excitations = n̂ᵢ + n̂ₒ + σ⁺σ⁻

	# particle Number conservation and correct final state
	particleNumbers = expect(excitations, solA) .|> real
	@test isapprox.(particleNumbers, 1.0, atol = 1e-5) |> all
	@test isapprox(ptrace(solA[end], [1, 2]).data[end, end], 1.0, atol = 1e-3)

	# ----------------------------------------
	# Two photon scattering
	WP2 = WavePacket(GaussMode(τ = 4.0), Fock(2))

	problemB = WaveguideProblem(TwoLevelChain(1), WP2, WP2.mode, 13.0)
	solB     = solve(problemB, reltol = 1e-10)[2]

	@test isapprox(
		ptrace(solB[end], [1, 2]).data[end, end] |> real,
		0.639387^2, atol = 1e-3
	) # See `Two Photon Overlaps.nb` for magic numbers

	problemB = WaveguideProblem(TwoLevelChain(1), WP2, outputMode, 13.0)
	solB     = solve(problemB, reltol = 1e-10)[2]

	@test isapprox(
		ptrace(solB[end], [1, 2]).data[end, end] |> real,
		0.17789^2, atol = 1e-3
	) # See `Two Photon Overlaps.nb` for magic numbers
end

@testset "Perfect Absorption" begin
	# ----------------------------------------
	# With a suitable state perfect absorption is almost possible
  γd =  0.5
  t₀ = -10.0/(1 - γd)

  problemA = WaveguideProblem(
		DissipativeLambdaChain(1, γd = γd),
    WavePacket(ExpMode(γ = 1 - γd), Fock(1)),
    (t₀, 0.0)
  )
  solA = ptrace.(solve(problemA, reltol = 1e-10)[2], 1)[end]

  # Project on the Rydberg manifold
  σ_Ryd = transition(NLevelBasis(3), 2, 2) + transition(NLevelBasis(3), 3, 3)
  @test real( 1 - expect(σ_Ryd, solA) ) < 1e-4

	problemB = WaveguideProblem(
		DissipativeLambdaChain(1, γd = γd),
		WavePacket(ExpMode(γ = 1 - γd), Fock(1)),
		(t₀, 10.0)
	)
  solB = ptrace.(solve(problemB, reltol = 1e-10)[2], 1)[end]

  # Project on the groundstate
  σ_GG = transition(NLevelBasis(3), 1, 1)
  @test abs( (1 - γd) / (1 + γd) - expect(σ_GG, solB) ) < 1e-4

  # Project on the dark state
  σ_DD = transition(NLevelBasis(3), 3, 3)
  @test abs( 2γd / (1 + γd) - expect(σ_DD, solB) ) < 1e-4
end

@testset "Exact Coherent Evolution" begin
	# ----------------------------------------
	# Compare against the exact time-evolution of coherently driven two level atom

	function blochVectorToDensity(vec)
	  ρ = nlevelstate(NLevelBasis(2), 1) |> dm

	  ρ.data[1, 1] = vec[1]
	  ρ.data[1, 2] = vec[2] / 2
	  ρ.data[2, 1] = conj(vec[2]) / 2
	  ρ.data[2, 2] = vec[3]

	  return ρ
	end

	function exactDynamics(t, κ, α; transform = false)
	  Ω = 4 * √(κ) * α
	  ω = √(4*Ω^2 - κ^2)

	  ρSS = [Ω^2 + 4κ^2, 4im * κ * Ω, Ω^2] ./ (2Ω^2 + 4κ^2)
	  ρ1  = [-1, (1.0im * κ + ω) / Ω, 1]
	  ρ2  = [-1, (1.0im * κ - ω) / Ω, 1]

	  c1 = Ω^2 / (-ω^2 + 3im * ω * κ) * exp((-3κ - 1.0im * ω) * t / 4)
	  c2 = Ω^2 / (-ω^2 - 3im * ω * κ) * exp((-3κ + 1.0im * ω) * t / 4)

	  ρ = ρSS .+ c1 .* ρ1 .+ c2 .* ρ2 |> blochVectorToDensity

	  if transform
	    ρ.data[1,2] *=  1im
	    ρ.data[2,1] *= -1im
	  end
	  return ρ
	end

	# ----------------------------------------
	# Test time evolution and steady state of superatom
  id, σWG, σGW = one(NLevelBasis(2)), transition(NLevelBasis(2), 2, 1),
		transition(NLevelBasis(2), 1, 2)
	αs = range(0.2, 2.0, length = 19)
	for αi ∈ αs
		problem = WaveguideProblem(TwoLevelChain(1), ContinuousWave(αi), 20.0)
		ts, sol = solve(problem, reltol = 1e-10)

		ρAtomsExact = exactDynamics.(ts, 1.0, αi, transform = true)
		@test sum(abs, tracedistance.(sol, ρAtomsExact)) / length(ts) < 1e-7

		@test abs(expect((αi*id + σWG) * (αi*id + σGW), sol[end]) - αi^2) < 1e-6
	end
end


text:
@doc raw"""
`TableauEPRK`: Tableau of an Explicit Partitioned Runge-Kutta method
```math
\begin{align*}
V_{n,i} &= \hphantom{-} \dfrac{\partial H}{\partial p} (Q_{n,i}, P_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= - \dfrac{\partial H}{\partial q} (Q_{n,i}, P_{n,i}) , &
P_{n,i} &= p_{n} + h  \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} ,
\end{align*}
```
usually satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + b_{j} a_{ji} &= b_{i} b_{j} , &
\bar{b}_i &= b_i .
\end{align*}
```
"""
struct TableauEPRK{T} <: AbstractTableauPRK{T}
    @HeaderTableau

    q::CoefficientsRK{T}
    p::CoefficientsRK{T}

    function TableauEPRK{T}(name, o, q, p) where {T}
        @assert q.s==p.s
        # TODO check that both tableaus are lower triangular and that only one element
        #      a_q[i,i] or a_p[i,i] is non-zero for all i.
        new(name, o, q.s, q, p)
    end
end

function TableauEPRK(name::Symbol, order::Int, q::CoefficientsRK{T}, p::CoefficientsRK{T}) where {T}
    TableauEPRK{T}(name, order, q, p)
end

function TableauEPRK(name::Symbol, order::Int, q::CoefficientsRK{T}) where {T}
    TableauEPRK{T}(name, order, q, q)
end

# TODO function readAbstractTableauPRKFromFile(dir::AbstractString, name::AbstractString)

"Explicit partitioned Runge-Kutta integrator."
struct IntegratorEPRK{DT,TT,VT,FT} <: DeterministicIntegrator{DT,TT}
    equation::PODE{DT,TT,VT,FT}
    tableau::TableauEPRK{TT}
    Δt::TT

    function IntegratorEPRK{DT,TT,VT,FT}(equation, tableau, Δt) where {DT,TT,VT,FT}
        new(equation, tableau, Δt)
    end
end

function IntegratorEPRK(equation::PODE{DT,TT,VT,FT}, tableau::TableauEPRK{TT}, Δt::TT) where {DT,TT,VT,FT}
    IntegratorEPRK{DT,TT,VT,FT}(equation, tableau, Δt)
end

equation(int::IntegratorEPRK) = int.equation
timestep(int::IntegratorEPRK) = int.Δt

"Explicit Runge-Kutta integrator cache."
mutable struct IntegratorCacheEPRK{DT,TT,D,S} <: PODEIntegratorCache{DT,D}
    n::Int
    t::TT
    t̅::TT

    q::Vector{TwicePrecision{DT}}
    q̅::Vector{TwicePrecision{DT}}
    p::Vector{TwicePrecision{DT}}
    p̅::Vector{TwicePrecision{DT}}

    Q::OffsetArray{Array{DT,1},1,Array{Array{DT,1},1}}
    P::OffsetArray{Array{DT,1},1,Array{Array{DT,1},1}}

    Y::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    F::Vector{Vector{DT}}

    function IntegratorCacheEPRK{DT,TT,D,S}() where {DT,TT,D,S}
        q = zeros(TwicePrecision{DT}, D)
        q̅ = zeros(TwicePrecision{DT}, D)
        p = zeros(TwicePrecision{DT}, D)
        p̅ = zeros(TwicePrecision{DT}, D)

        Q = create_internal_stage_vector_with_zero(DT, D, S)
        P = create_internal_stage_vector_with_zero(DT, D, S)

        Y = create_internal_stage_vector(DT, D, S)
        Z = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)

        new(0, zero(TT), zero(TT), q, q̅, p, p̅, Q, P, Y, Z, V, F)
    end
end

function create_integrator_cache(int::IntegratorEPRK{DT,TT}) where {DT,TT}
    IntegratorCacheEPRK{DT, TT, ndims(int), int.tableau.s}()
end

function CommonFunctions.reset!(cache::IntegratorCacheEPRK{DT,TT}, Δt::TT) where {DT,TT}
    cache.t̅  = cache.t
    cache.q̅ .= cache.q
    cache.p̅ .= cache.p
    cache.t += Δt
    cache.n += 1
end

function CommonFunctions.set_solution!(cache::IntegratorCacheEPRK, sol, n=0)
    t, q, p = sol
    cache.n  = n
    cache.t  = t
    cache.q .= q
    cache.p .= p
end

"Compute Q stages of explicit partitioned Runge-Kutta methods."
function computeStageQ!(int::IntegratorEPRK{DT,TT}, cache::IntegratorCacheEPRK{DT,TT}, i::Int, jmax::Int, t) where {DT,TT}
    fill!(cache.Y[i], zero(DT))
    for j in 1:jmax
        for k in 1:int.equation.d
            cache.Y[i][k] += int.tableau.q.a[i,j] * cache.V[j][k]
        end
    end
    for k in 1:int.equation.d
        cache.Q[i][k] = cache.q[k] + int.Δt * cache.Y[i][k]
    end
    int.equation.f(t, cache.Q[i], cache.P[jmax], cache.F[i])
end

"Compute P stages of explicit partitioned Runge-Kutta methods."
function computeStageP!(int::IntegratorEPRK{DT,TT}, cache::IntegratorCacheEPRK{DT,TT}, i::Int, jmax::Int, t) where {DT,TT}
    fill!(cache.Z[i], zero(DT))
    for j in 1:jmax
        for k in 1:int.equation.d
            cache.Z[i][k] += int.tableau.p.a[i,j] * cache.F[j][k]
        end
    end
    for k in 1:int.equation.d
        cache.P[i][k] = cache.p[k] + int.Δt * cache.Z[i][k]
    end
    int.equation.v(t, cache.Q[jmax], cache.P[i], cache.V[i])
end

"Integrate partitioned ODE with explicit partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorEPRK{DT,TT}, cache::IntegratorCacheEPRK{DT,TT}) where {DT,TT}
    local tqᵢ::TT
    local tpᵢ::TT

    # store previous solution
    cache.Q[0] .= cache.q
    cache.P[0] .= cache.p

    # reset cache
    reset!(cache, timestep(int))

    # compute internal stages
    for i in 1:int.tableau.s
        tqᵢ = cache.t̅ + int.Δt * int.tableau.q.c[i]
        tpᵢ = cache.t̅ + int.Δt * int.tableau.p.c[i]

        if int.tableau.q.a[i,i] ≠ zero(TT) && int.tableau.p.a[i,i] ≠ zero(TT)
            error("This is an implicit method!")
        elseif int.tableau.q.a[i,i] ≠ zero(TT)
            computeStageP!(int, cache, i, i-1, tpᵢ)
            computeStageQ!(int, cache, i, i, tqᵢ)
        elseif int.tableau.p.a[i,i] ≠ zero(TT)
            computeStageQ!(int, cache, i, i-1, tqᵢ)
            computeStageP!(int, cache, i, i, tpᵢ)
        else
            computeStageQ!(int, cache, i, i-1, tqᵢ)
            computeStageP!(int, cache, i, i-1, tpᵢ)
        end
    end

    # compute final update
    update_solution!(cache.q, cache.V, int.tableau.q.b, int.Δt)
    update_solution!(cache.p, cache.F, int.tableau.p.b, int.Δt)

    # take care of periodic solutions
    cut_periodic_solution!(cache.q, int.equation.periodicity)
end


text:
"""
Determines how to react to an exceeded rate limit.

- `ORL_THROW`: Throw a [`RateLimitedError`](@ref).
- `ORL_WAIT`: Block and wait for the rate limit to expire.
"""
@enum OnRateLimit ORL_THROW ORL_WAIT

"""
A generic rate limiter using the `[X-]RateLimit-Remaining` and `[X-]RateLimit-Reset` response headers.
The reset header is assumed to be a Unix timestamp in seconds.
"""
mutable struct RateLimiter
    remaining::Int
    reset::Int

    RateLimiter() = new(1, 0)
end

rate_limit_check(rl::RateLimiter) =
    rl.remaining == 0 && rate_limit_period(rl) > Millisecond(0)

rate_limit_wait(rl::RateLimiter) = sleep(rate_limit_period(rl))

rate_limit_period(rl::RateLimiter) = max(Millisecond(0), unix2datetime(rl.reset) - now(UTC))

function rate_limit_update!(rl::RateLimiter, r::HTTP.Response)
    remaining = tryheader(i -> parse(Int, i), r, "RateLimit-Remaining")
    reset = tryheader(i -> parse(Int, i), r, "RateLimit-Reset")
    if remaining === nothing  || reset === nothing
        @warn "Parsing rate limit headers failed"
    else
        rl.remaining = remaining
        rl.reset = reset
    end
end

function tryheader(f::Function, r::HTTP.Response, header::AStr)
    for h in [header, "X-$header"]
        HTTP.hasheader(r, h) && return f(HTTP.header(r, h))
    end
    return nothing
end


text:
using SchumakerSpline
from = 0.0
to = 10
x1 = collect(range(from, stop=to, length=40))
y1 = (x1).^2
s1 = Schumaker(x1,y1; extrapolation = (Constant, Curve))
#plot(s1)
abs(s1(x1[1]) - s1(x1[1]-0.5)) < eps()
abs(s1(x1[40]) - s1(x1[40]+0.5)) > 0.01

s2 = Schumaker(x1,y1; extrapolation = (Constant, Constant))
abs(s2(x1[40]) - s2(x1[40]+0.5)) < eps()
plot(s2)

s3 = Schumaker(x1,y1; extrapolation = (Constant, Linear))
abs((s3(x1[40]+0.5) - s3(x1[40])) - (s3(x1[40]+1.5) - s3(x1[40]+1.0))) < 2e-14
TopX = x1[40]
abs(s3(TopX - 10*eps()) + s3(TopX - 10*eps(),1)*0.5 - s3(TopX +0.5) ) < 1000*eps()


text:
struct Component{T <: Number}
    x::T
    y::T
    z::T
    rx::T
    ry::T
    rz::T
end

"""
    PunchType

Abstract supertype of PunchModalData and PunchFrequencyData
"""
abstract type PunchType end

"""
    PunchModalData

Nastran punch file modal data format.
"""
mutable struct PunchModalData{T <: Real} <: PunchType
    title::String
    subtitle::String
    label::String
    datatype::String
    outputtype::String
    subcaseid::Int
    eigenvalue::Float64
    mode::Int
    eigenvector::Vector{Component{T}} # only support for real modes so far
    PunchModalData{T}() where T = (x = new(); x.eigenvector = Vector{Component{T}}(undef, 0); x)
end

"""
    PunchFrequencyData

Nastran punch file modal data format.
"""
mutable struct PunchFrequencyData{T <: Real} <: PunchType
    title::String
    subtitle::String
    label::String
    datatype::String
    outputtype::String
    subcaseid::Int
    point::Int
    response::Vector{Component{Complex{T}}}
    PunchFrequencyData{T}() where T = (x = new(); x.response = Vector{Component{T}}(undef, 0); x)
end

"""
    PMData

Nastran punch file modal data load format.
"""
struct PMData{T <: Real}
    modal::Vector{PunchModalData{T}}
    node::Vector{Int}
end

"""
    PData

Nastran punch file response data load format.
"""
struct PFData{T <: Real}
    response::Vector{PunchFrequencyData{T}}
    frequency::Vector{T}
end

"""
    PunchData

Punch data union holding all other punch data types.
"""
PunchData{T} = Union{PunchModalData{T}, PunchFrequencyData{T}}

"""
    readpch(file::String, type::PunchData)

Read punch file.
"""
readpch(file::String, type::PunchData{T}) where T <: Real

"""
    readpch(file::String, ::Type{PunchModalData})

Read SOL103 real modal data from Nastran punch file according to SORT2.

This should be avoided, in favour of HDF5, but is included for existing result files.

# Reference
MSC Nastran 2017.1 Reference Manual
"""
function readpch(file::String, type::Type{PunchModalData{T}}) where T <: Real
    data_modal = Vector{type}(undef, 0)
    node = Vector{Int}(undef, 0)

    #try
        lines = readlines(file)
    #catch e
    #    println("Error $e. Can't read file")
    #end

    cnew = false
    data = PunchModalData{T}()

    next = iterate(lines)
    while next !== nothing
        (line, state) = next

        line = line[1:PCHLL]
        isdollar = line[1] == DOLLAR

        if isdollar && cnew
            push!(data_modal, data)
            data = PunchModalData{T}()
            cnew = false
        end

        if isdollar
            readmeta!(data, line)
        else
            state = readcomponents!(data, node, line, lines, state)
            cnew = true
        end

        next = iterate(lines, state)
    end
    push!(data_modal, data)

    return PMData(data_modal, node)
end

"""
    readpch(file::String, ::Type{PunchFrequencyData})

Read SOL111/108 frequency response data from Nastran punch file according to SORT2 and IMAG.

This should be avoided, in favour of HDF5, but is included for existing result files.

# Reference
MSC Nastran 2017.1 Reference Manual
"""
function readpch(file::String, type::Type{PunchFrequencyData{T}}) where T <: Real
    data_frequency = Vector{type}(undef, 0)
    frequency = Vector{T}(undef, 0)

    #try
        lines = readlines(file)
    #catch e
    #    println("Error $e. Can't read file")
    #end

    cnew = false
    data = PunchFrequencyData{T}()

    next = iterate(lines)
    while next !== nothing
        (line, state) = next

        line = length(line) < PCHLL ? line[:] : line[1:PCHLL]
        isdollar = line[1] == DOLLAR

        # can this be more general?
        if length(line) > 8
            if contains(line[1:8], "SET") || line[1:8] == "        "
                next = iterate(lines, state)
                continue
            end
        end

        if isdollar && cnew
            push!(data_frequency, data)
            data = PunchFrequencyData{T}()
            cnew = false
        end

        if isdollar
            readmeta!(data, line)
        else
            state = readcomponents!(data, frequency, line, lines, state)
            cnew = true
        end

        next = iterate(lines, state)
    end
    push!(data_frequency, data)

    return PFData(data_frequency, frequency)
end

"""
    readmeta!(data::PunchModalData, line)

Read modal metadata.
"""
function readmeta!(data::PunchModalData{T}, line) where T <: Real
    meta = split(line[2:end], EQUAL)
    meta1 = rstrip(meta[1])
    meta2 = length(meta) == 1 ? rstrip(lstrip(meta[1])) : rstrip(lstrip(meta[2]))
    #meta2 = meta2 === nothing ? "" : meta

    if meta1 == TITLE
        data.title = meta2
    elseif meta1 == SUBTITLE
        data.subtitle = meta2
    elseif meta1 == LABEL
        data.label = meta2
    elseif meta1 == SUBCASE_ID
        data.subcaseid = parse(Int64, meta2)
    elseif meta1 == EIGENVALUE
        data.eigenvalue = parse(Float64, split(meta2)[1])
        data.mode = parse(Int64, strip(meta[3]))
    else
        if !isdefined(data, :datatype)
            data.datatype = meta2
        else
            data.outputtype = meta2
        end
    end

    return nothing
end

"""
    reameta!(data::PunchFrequencyData, frequency, line, lines, state)

Read frequency metadata.
"""
function readmeta!(data::PunchFrequencyData{T}, line) where T <: Real
    meta = split(line[2:end], EQUAL)
    meta1 = rstrip(meta[1])
    meta2 = length(meta) == 1 ? rstrip(lstrip(meta[1])) : rstrip(lstrip(meta[2]))
    #meta2 = meta2 === nothing ? "" : meta

    if meta1 == TITLE
        data.title = meta2
    elseif meta1 == SUBTITLE
        data.subtitle = meta2
    elseif meta1 == LABEL
        data.label = meta2
    elseif meta1 == SUBCASE_ID
        data.subcaseid = parse(Int64, meta2)
    elseif meta1 == POINT
        data.point = parse(Int64, strip(split(meta[2])[1]))
    else
        if !isdefined(data, :datatype)
            data.datatype = meta2
        else
            data.outputtype = meta2
        end
    end

    return nothing
end

"""
    readcomponents!(data::PunchModalData, node, line, lines, state)

Read modal components.
"""
function readcomponents!(data::PunchModalData{T}, node, line, lines, state) where T <: Real
    c1 = strip(line[1:PCHL - 1]) # skip last character
    cmp(c1, CONT) == 0 ? nothing : push!(node, parse(Int, c1))

    x = parse(T, strip(line[PCHL + 1:2PCHL]))
    y = parse(T, strip(line[2PCHL + 1:3PCHL]))
    z = parse(T, strip(line[3PCHL + 1:4PCHL]))

    (line, state) = iterate(lines, state)
    rx = parse(T, strip(line[PCHL + 1:2PCHL]))
    ry = parse(T, strip(line[2PCHL + 1:3PCHL]))
    rz = parse(T, strip(line[3PCHL + 1:4PCHL]))

    push!(data.eigenvector, Component{T}(x, y, z, rx, ry, rz))

    return state
end

"""
    readcomponents!(data::PunchFrequencyData, frequency, line, lines, state)

Read real/imaginary translation and rotation components.
"""
function readcomponents!(data::PunchFrequencyData{T}, frequency, line, lines, state) where T <: Real
    # for now I think this is clearest
    c1 = strip(line[1:PCHL - 1]) # skip last character
    cmp(c1, CONT) == 0 ? nothing : push!(frequency, parse(Float64, c1))

    x = parse(T, strip(line[PCHL + 1:2PCHL]))
    y = parse(T, strip(line[2PCHL + 1:3PCHL]))
    z = parse(T, strip(line[3PCHL + 1:4PCHL]))

    (line, state) = iterate(lines, state)
    rx = parse(T, strip(line[PCHL + 1:2PCHL]))
    ry = parse(T, strip(line[2PCHL + 1:3PCHL]))
    rz = parse(T, strip(line[3PCHL + 1:4PCHL]))

    (line, state) = iterate(lines, state)
    x = x + im*parse(T, strip(line[PCHL + 1:2PCHL]))
    y = y + im*parse(T, strip(line[2PCHL + 1:3PCHL]))
    z = z + im*parse(T, strip(line[3PCHL + 1:4PCHL]))

    (line, state) = iterate(lines, state)
    rx = rx + im*parse(T, strip(line[PCHL + 1:2PCHL]))
    ry = ry + im*parse(T, strip(line[2PCHL + 1:3PCHL]))
    rz = rz + im*parse(T, strip(line[3PCHL + 1:4PCHL]))

    push!(data.response, Component{Complex{T}}(x, y, z, rx, ry, rz))

    return state
end

"""
    collect(data::PMData)

Collect all modal data into a matrix .
"""
function collect(data::PMData{T}) where T <: Real
    n, m = length(data.modal[1].eigenvector), length(data.modal)
    d = Matrix{Float64}(undef, n*DOF, m)
    for i in 1:m
        d[:, i] = collect(data.modal[i].eigenvector)
    end

    return d
end

"""
    collect(data::PFData)

Collect all response data into a matrix.
"""
function collect(data::PFData{T}) where T <: Real
    n, m = length(data.response[1].response), length(data.response)
    d = Matrix{Complex{T}}(undef, n*DOF, m)
    for i in 1:m
        d[:, i] = collect(data.response[i].response)
    end

    return d
end

"""
    collect(data::Component)

Collect components in x, y, z, rx, ry, and rz order
"""
function collect(data::Component{T}) where T <: Number
    return [data.x, data.y, data.z, data.rx, data.ry, data.rz]
end

"""
    collect(data::Vector{Component})

Collect components in x, y, z, rx, ry, and rz order
"""
function collect(data::Vector{Component{T}}) where T <: Number
    return collect(Iterators.flatten([collect(x) for x in data]))
end

max_stars_repo_name: mgcth/SDNastran

max_issues_repo_name: mgcth/SDNastran

max_forks_repo_name: mgcth/SDNastran


text:
struct IsingSampler <: Sampleable{Multivariate,Discrete}
    p::Int64                 # number of vertices
    θ::Vector{Float64}       # parameter matrix
    nbs::Vector{Vector{Int64}}
    state::BitVector
    thin::Int64
end

function IsingSampler(θ::Vector{Float64}; burn::Int64=3000, thin::Int64=1000)
    m = length(θ)
    p = convert(Int, ceil((1. + sqrt(1. + 8. * m)) / 2.))

    # Find neighbors
    nbs = Vector{Vector{Int64}}(undef, p)
    for node=1:p
        nbs[node] = Vector{Int64}(undef, 0)
    end
    for node=1:p
        for nb=1:p
            if node==nb
                continue
            end
            if !iszero(θ[trimap(node, nb)])
                push!(nbs[node], nb)
            end
        end
    end

    # initial state
    state = rand(Bernoulli(), p) .== 1.
    spl = IsingSampler(p, θ, nbs, state, thin)
    for j=1:burn
        _sample_one!(spl)
    end
    spl
end

Base.length(s::IsingSampler) = s.p
Base.eltype(s::IsingSampler) = Bool

Distributions.rand(s::IsingSampler) = Distributions._rand!(s, BitVector(undef, length(s)))
Distributions.rand(s::IsingSampler, n::Int) = Distributions._rand!(s, BitMatrix(undef, length(s), n))

function Distributions._rand!(s::IsingSampler, x::BitMatrix)
    n = size(x, 2)
    for i=1:n
        for j=1:s.thin
            _sample_one!(s)
        end
        _sample_one!(s)
        x[:, i] .= s.state
    end
    x
end

function Distributions._rand!(s::IsingSampler, x::BitVector)
    for j=1:s.thin
        _sample_one!(s)
    end
    _sample_one!(s)
    x .= s.state
end

function _sample_one!(spl::IsingSampler)
    for node=1:length(spl)
        _sample_one_coordinate!(spl, node)
    end
    spl
end

function _sample_one_coordinate!(
    spl::IsingSampler,
    node::Int64)

    θ, nbs, y = spl.θ, spl.nbs, spl.state

    v = 0.
    for nb in spl.nbs[node]
        k = trimap(node, nb)
        v += y[nb] ? θ[k] : -θ[k]
    end
    Q = exp(2. * v)
    y[node] = rand() <= Q / (Q + 1.)
end

### graph construction

# lenC --- length of a chain
# the graph consists of p / lenC chains
# sgn +1, 0, -1 --- whether edges are positive, mixed, or negative
# lb and ub are assumed positive
function chain(p::Integer; lenC::Integer=10, lb::Float64=0.1, ub::Float64=0.3, sgn::Integer=1)
    mod(p, lenC) == 0 || throw(ArgumentError("p=$p should be divisible by lenC=$lenC"))
    m = div(p*(p-1), 2)
    θ = zeros(Float64, m)

    du = Uniform(lb, ub)
    for ic=1:div(p, lenC)
        for j=1:lenC-1
            col = (ic - 1) * lenC + j
            row = col + 1

            v = rand(du)
            if sgn == -1
                v *= -1.
            elseif sgn == 0
                v *= 2. * rand(Bernoulli()) - 1.
            end
            θ[trimap(row, col)] = v
        end
    end
    θ
end

# k --- # of children
# the graph is a k-ary tree on p nodes
# sgn +1, 0, -1 --- whether edges are positive, mixed, or negative
# lb and ub are assumed positive
function kary_tree(p::Integer; k::Integer=3, lb::Float64=0.1, ub::Float64=0.3, sgn::Integer=1)
    k > 1 || throw(ArgumentError("k=$k should be at least 2"))
    m = div(p * (p - 1), 2)
    θ = zeros(Float64, m)

    du = Uniform(lb, ub)

    depth = ceil(Int64, log(k, (k - 1) * p + 1) -1)

    youngest = cumsum(k.^(0:depth))
    youngest[end] = p

    children = (youngest[1] + 1):youngest[2]
    for j in children
        v = rand(du)
        if sgn == -1
            v *= -1.
        elseif sgn == 0
            v *= 2. * rand(Bernoulli()) - 1.
        end
        θ[trimap(1, j)] = v
    end
    for d = 2:depth
        children = (youngest[d] + 1):youngest[d + 1]
        for j in children
            i = fld(j - youngest[d] - 1, k) + youngest[d-1] + 1
            v = rand(du)
            if sgn == -1
                v *= -1.
            elseif sgn == 0
                v *= 2. * rand(Bernoulli()) - 1.
            end
            θ[trimap(i, j)] = v
        end
    end
    θ
end

function removeEdges!(θ, numChanges::Integer=4)
    ind_change = sample(findall(!iszero, θ), numChanges; replace=false)
    for i in ind_change
        θ[i] = 0.
    end
    θ
end

function addEdges!(θ, numChanges::Integer=4, lb::Float64=0.2, ub::Float64=0.4, sgn::Integer=1)
    ind_change = sample(findall(iszero, θ), numChanges; replace=false)
    du = Uniform(lb, ub)
    for i in ind_change
        v = rand(du)
        if sgn == -1
            v *= -1.
        elseif sgn == 0
            v *= 2. * rand(Bernoulli()) - 1.
        end
        θ[i] = v
    end
    θ
end


text:
using Statistics, StatsBase
using PyCall
xr = pyimport("xarray")
using PyPlot
using Elefridge
using JLD
using BFloat16s

path = "/Users/milan/cams/gridded/"
filelist = filter(x->endswith(x,".grib"),readdir(path))

n = length(filelist)
varnames = fill("",n)

for (i,file) in enumerate(filelist)
    varname = split(split(file,"0_")[end],".")[1]
    varnames[i] = varname
    println("---")
    println("Reading $varname")
    X = xr.open_dataarray(joinpath(path,file),engine="cfgrib").data

    N = prod(size(X))
    nzero = sum(X .== zero(eltype(X)))
    nneg = sum(X .< zero(eltype(X)))

    # println("$nzero, $(nzero/N*100)% 0 entries.")
    # println("$nneg, $(nneg/N*100)% negative entries.")

    H = bitentropy(X)
    println("$H bit Float32 entropy.")
end


text:
"""
    showdistribution([io], [mime], d)

Show a UnicodePlots.barplot representation of a distribution.

# Keyword Arguments

- `title::String=string(typeof(d))*" distribution"`: title for the barplot.
"""
function showdistribution(io::IO, mime::MIME"text/plain", d; title=string(typeof(d))*" distribution")
    limited = get(io, :limit, false)
    strings = String[]
    probs = Float64[]

    rows = first(get(io, :displaysize, displaysize(io)))
    rows -= 6 # Yuck! This magic number is also in Base.print_matrix

    if limited && rows > 1 && length(support(d)) >= rows
        for (x,p) in Iterators.take(weighted_iterator(d), rows-1)
            push!(strings, sprint(show, x)) # maybe this should have conext=:compact=>true
            push!(probs, p)
        end

        push!(strings, "<everything else>")
        push!(probs, 1.0-sum(probs))
    else
        for (x,p) in weighted_iterator(d)
            push!(strings, sprint(show, x))
            push!(probs, p)
        end
    end
    show(io, mime, barplot(strings, probs, title=title))
end

showdistribution(io::IO, d; kwargs...) = showdistribution(io, MIME("text/plain"), d; kwargs...)
showdistribution(d; kwargs...) = showdistribution(stdout, d; kwargs...)


text:
function has_pandigital_product(n::Int)

	sqrtn=floor(Int,sqrt(Float64(n)))
	for i=1:sqrtn
		if  n%i==0
			str=string(n)*string(i)*string(Int(n/i))
			if Set(str)==Set("123456789")
				return true
			end
		end
	end

	return false
end

function compute()

	ans=sum( i for i in 1:10000 if has_pandigital_product(i) )
	return ans

end

println(compute())


text:
# Submodule for imaging functions

#######################################
# Function to resize an array representation
# of an image to a certain number of pixeles
function resize_image(f, N)
    if size(f, 1) > size(f, 2)
		f = reshape(f, size(f, 2), size(f, 1), size(f, 3));
    end
	P = size(f, 1);
	# add 1 pixel boundary
	g = f;
	g = cat(2, g, reshape(g[:,1,:], size(g, 1), 1, size(g, 3)));
    g = cat(1, g, reshape(g[1,:,:], 1, size(g, 2), size(g, 3)));
	# interpolate
	t = linspace(1, P, N);
	ti = round.(Int64, floor.(t)) ; tj = round.(Int64, ceil.(t));
	fi = round.(Int64, t - floor.(t)); fj = 1 - fi;
	h = zeros(N, N, size(f, 3));
	for s in 1:size(f, 3)
	    h[:,:,s] = g[ti,ti,s] .* (fj * fj') + g[tj,tj,s] .* (fi * fi') + g[ti,tj,s] .* (fj * fi') + g[tj,ti,s] .* (fi * fj');
	end
	return h;
end # resize_image

#####################################
# Function that rescales date in a,b
"""
...
rescale(x::Array,a=0,b=1) rescales data in
[a,b]
...
"""
function rescale(x, a=0, b=1)
	m = minimum(x[:]);
	M = maximum(x[:]);

	if M - m < 1e-10
	    y = x;
	else
	    y = (b - a) * (x - m) / (M - m) + a;
	end
	return y;

end ## rescale

#####################################################
# Function that rescales a wavelet transform in an array A
function rescaleWav(A)
    v = maximum(abs.(A[:]));
    B = copy(A)
    if v > 0
        B = .5 + .5 * A / v;
    end
    return B;
end


text:
function plothistory(problem, p = plot())
    hist = problem.model.history
    h = hist[:p]
    nonnan = !isnan.(hist[:p].values)
    val = hist[:p].values[nonnan]
    iter = hist[:p].iterations[nonnan]
    plot!(p, iter, val, yscale=:log10)
end


text:
const LLVec=[
"Aggregate",
"AutoPipeline",
"BaselineClassifier",
"BaselineRegressor",
"Batching",
"Both",
"ConcatFeatures",
"GridSearchCV",
"GroupBy",
"HalvingGridSearchCV",
"Hyperopt",
"IdentityWrapper",
"Join",
"Map",
"NoOp",
"Observing",
"Project",
"Relational",
"SMAC",
"SampleBasedVoting",
"Scan",
"TopKVotingClassifier"]


text:
# Broken code
# https://github.com/JuliaIO/SerialPorts.jl/issues/23

using SerialPorts
using Test

@testset "hello_serialports" begin

    function _write_end_of_command(ser)
        for i in 0:2
            write(ser, 0xff)
        end
    end

    function _execute_command(ser, cmd)
        println(cmd)
        write(ser, cmd)
        _write_end_of_command(ser)
    end

    function main()
        portname = "/dev/ttyUSB0"
        baudrate = 9600
        ser = SerialPort(portname, baudrate)
        _execute_command(ser, "page 0")
        sleep(2)
        _execute_command(ser, "page 1")
        sleep(2)
        _execute_command(ser, "t0.txt=\"Hello\"")
        close(ser)
    end

    main()
    @test true

end


text:
using FastGroupBy,BenchmarkTools

@time id6c = copy(id6)
@time res = sumby(id6c, v1) #17,20
@time sumby_sorted(id6c, v1)
@time sumby_sorted(id6c, v1)
@time sumby_sorted2(id6c, v1)
@time sumby_sorted2(id6c, v1)

@time id6c = copy(id6)
@time res = sumby(id6c, v1) #17,11
@time res1 = sumby_sorted(id6c, v1)
@time res1 = sumby_sorted(id6c, v1)

pool1 = [@sprintf "id%010d" k for k in 1:(N/K)]
id3 = rand(pool1, N)

function randstrarray1(pool, N)
    K = length(pool)
    PooledArray(PooledArrays.RefArray(rand(1:K, N)), pool)
end

@time id3 = rand(pool1, N)

@time a = sumby(hash.(id3), v1);
@time a = sumby(hash.(id3), v1);

@time hid3 = hash.(id3)
@time hid3 = hash.(id3)

@time a = sumby(hid3, v1);

@time a = sumby_dict(id3, v1)


text:
# Parses policy xml file and returns alpha vectors and alpha actions
# Should handle any policy in the .policy format
# Should handle policies written by POMDPs.jl or APPL
#
# alpha_vectors is a matrix containing the alpha vectors
#   Each row corresponds to a different alpha vector
#
# alpha_actions is a vector containing the list of action indices
#   Each index corresponds to action associated with the alpha vector of this row
#   These are 0-indexed... 0 means it is the first action
#
# TODO: Check that the input file exists and handle the case that it doesn't
# TODO: Handle sparse vectors
#

function read_momdp(filename::String)

    # Parse the xml file
    # TODO: Check that the file exists and handle the case that it doesn't
    xdoc = parse_file(filename)

    # Get the root of the document (the Policy tag in this case)
    policy_tag = root(xdoc)
    #println(name(policy_tag))      # print the name of this tag

    # Determine expected number of vectors and their length
    alphavector_tag = get_elements_by_tagname(policy_tag, "AlphaVector")[1]     #length of 1 anyway
    num_vectors = int(attribute(alphavector_tag, "numVectors"))
    vector_length = int(attribute(alphavector_tag, "vectorLength"))
    num_full_obs_states = int(attribute(alphavector_tag, "numObsValue"))

    # For debugging purposes...
    #println("AlphaVector tag: # vectors, vector length: $(num_vectors), $(vector_length)")

    # Arrays with vector and sparse vector tags
    vector_tags = get_elements_by_tagname(alphavector_tag, "Vector")
    sparsevector_tags = get_elements_by_tagname(alphavector_tag, "SparseVector")

    num_vectors_check = length(vector_tags) + length(sparsevector_tags)         # should be same as num_vectors

    # Initialize the gamma matrix. This is basically a matrix with the alpha
    #   vectors as rows.
    #alpha_vectors = Array(Float64, num_vectors, vector_length)
    alpha_vectors = Array{Float64}(vector_length, num_vectors)
    alpha_actions = Array{String}(num_vectors)
    observable_states = Array{String}(num_vectors)
    gammarow = 1

    # Fill in gamma
    for vector in vector_tags
        alpha = float(split(content(vector)))
        #alpha_vectors[gammarow, :] = alpha
        alpha_vectors[:,gammarow] = alpha
        alpha_actions[gammarow] = attribute(vector, "action")
        observable_states[gammarow] = attribute(vector, "obsValue")
        gammarow += 1
    end

    # TODO: Handle sparse vectors
    for vector in sparsevector_tags
        # Turn these into vectors as well
    end

    # Return alpha vectors and indices of actions
    return alpha_vectors, int(alpha_actions), int(observable_states)
end

function read_pomdp(filename::String)

    # Parse the xml file
    # TODO: Check that the file exists and handle the case that it doesn't
    xdoc = parse_file(filename)

    # Get the root of the document (the Policy tag in this case)
    policy_tag = root(xdoc)
    #println(name(policy_tag))      # print the name of this tag

    # Determine expected number of vectors and their length
    alphavector_tag = get_elements_by_tagname(policy_tag, "AlphaVector")[1]     #length of 1 anyway
    num_vectors = parse(Int64, attribute(alphavector_tag, "numVectors"))
    vector_length = parse(Int64, attribute(alphavector_tag, "vectorLength"))
    num_full_obs_states = parse(Int64, attribute(alphavector_tag, "numObsValue"))

    # For debugging purposes...
    #println("AlphaVector tag: # vectors, vector length: $(num_vectors), $(vector_length)")

    # Arrays with vector and sparse vector tags
    vector_tags = get_elements_by_tagname(alphavector_tag, "Vector")
    sparsevector_tags = get_elements_by_tagname(alphavector_tag, "SparseVector")

    num_vectors_check = length(vector_tags) + length(sparsevector_tags)         # should be same as num_vectors

    # Initialize the gamma matrix. This is basically a matrix with the alpha
    #   vectors as columns.
    #alpha_vectors = Array(Float64, num_vectors, vector_length)
    alpha_vectors = Array{Float64}(undef, vector_length, num_vectors)
    alpha_actions = Array{String}(undef, num_vectors)
    observable_states = Array{String}(undef, num_vectors)
    gammarow = 1

    # Fill in gamma
    for vector in vector_tags
        alpha = parse.(Float64, split(content(vector)))
        #alpha_vectors[gammarow, :] = alpha
        alpha_vectors[:,gammarow] = alpha
        alpha_actions[gammarow] = attribute(vector, "action")
        observable_states[gammarow] = attribute(vector, "obsValue")
        gammarow += 1
    end

    # TODO: Handle sparse vectors
    for vector in sparsevector_tags
        # Turn these into vectors as well
    end

    # Return alpha vectors and indices of actions
    return alpha_vectors, [parse(Int64,s) for s in alpha_actions]
end


text:
"""
    set_pixels!(o, v)

Internal function to set the number of pixels used for a color-range.
"""
set_pixels!(o::Colorbar, v::Int) = (o.pixels = v;)

function set_position!(o::Colorbar, v::String)
    v_lc = lowercase(v)
    v_lc ∈ ("left", "right", "top", "bottom") || throw(
        ArgumentError("Unknown position value $v")
    )
    o.position = v_lc
    return
end

set_ticks!(o::Colorbar, v::Vector{F64}) =
    (o.ticks.places = v; o.ticks.labels = TicksLabels();)

function set_labels!(o::Colorbar, v::Vector{String})
    isempty(o.ticks.places) && throw(
        ArgumentError("""
            You must specify the ticks position before the labels.
            """)
        )
    length(o.ticks.places) == length(v) || throw(
        ArgumentError("""
            The dimensions of ticks places and labels don't match.
            """)
        )
    o.ticks.labels.names = v
    return
end


text:
using DifferentialEquations, Statistics, Random

function Diffusion(μ,κ,α,σ,ρ,u0,tspan)
    f = function (du,u,p,t)
        du[1] = μ # drift in log prices
        du[2] = κ.*(α.-u[2]) # mean reversion in shocks
    end
    g = function (du,u,p,t)
        du[1] = exp(u[2]/2.0)
        du[2] = σ
    end
    Γ = [1.0 ρ;ρ 1.0] # Covariance Matrix
    noise = CorrelatedWienerProcess!(Γ,tspan[1],zeros(2),zeros(2))
    sde_f = SDEFunction{true}(f,g)
    SDEProblem(sde_f,g,u0,tspan,noise=noise)
end

@views function JDmodel(θ, burnin, rndseed=1234)
    Random.seed!(rndseed)
    TradingDays = burnin+1000 # the sample is 1000 days, also need initial burnin and burnin between samples
    Days = TradingDays + Int(TradingDays/5*2) # add weekends
    MinPerDay = 1440 # minutes per day
    MinPerTic = 10 # minutes between tics, lower for better accuracy
    tics = Int(MinPerDay/MinPerTic) # number of tics in day
    dt = 1/tics # divisions per day
    closing = Int(round(6.5*60/MinPerTic)) # tic at closing
    # solve the diffusion
    μ, κ, α, σ, ρ, λ0, λ1, τ = θ
    u0 = [μ; α]
    prob = Diffusion(μ, κ, α, σ, ρ, u0, (0.0,Days))
    ## jump in log price
    rate(u,p,t) = λ0.*(λ0>0.0) # the prior allows for negative rate, to allow an accumulation at zero
    # jump is random sign times  λ1 times current st. dev.
    affect1!(integrator) = (integrator.u[1] = integrator.u[1].+rand([-1.0,1.0]).*λ1.*exp(integrator.u[2]./2.0))
    jump = ConstantRateJump(rate,affect1!)
    jump_prob = JumpProblem(prob,Direct(), jump)
    # do the simulation
    sol = solve(jump_prob,SRIW1(), dt=dt, adaptive=false, seed=rndseed)
    # get log price, with measurement error
    lnPs = [sol(t)[1] .+ τ.*randn()   for t in dt:dt:Days]
    # get log price at end of trading days. We will compute lag, so loose first
    lnPtrading = zeros(TradingDays+1)
    RV = zeros(TradingDays+1)
    BV = zeros(TradingDays+1)
    DayofWeek = 0 # counter for day of week
    TradingDay = 0 # counter for trading days
    Day = 0
    lnPlag = 0.0
    @inbounds while TradingDay < TradingDays
        Day +=1
        DayofWeek +=1
        # set day of week, and record if it's a trading day
        if DayofWeek<6
            TradingDay +=1 # advance trading day
            # compute realized measures
            t1 = 0.0
            t2 = 0.0
            for tic = 1:closing
                lnP = lnPs[(Day-1)*tics+tic]
                ret = lnP - lnPlag
                lnPlag = lnP # update lag, this is inter-day for the first
                t2 = t1 # one lag
                t1 = abs(ret) # current
                # RV measures
                if tic > 1
                    RV[TradingDay] += ret*ret
                    BV[TradingDay] += t1*t2
                end
            end
            lnPtrading[TradingDay] = lnPs[(Day-1)*tics+closing]
        end
        if DayofWeek==7 # restart the week if Sunday
            DayofWeek = 0
        end
    end
    rets = lnPtrading[burnin+1:burnin+1000] .- lnPtrading[burnin:burnin+1000-1] # inter-day returns
    RV = RV[burnin+1:burnin+1000]
    BV = (pi/2.0) .* BV[burnin+1:burnin+1000]
    [rets RV BV]
end

# auxstats, using simulated data
@views function auxstat(θ, reps)
    auxstat.(JDmodel(θ, 50, rand(1:Int64(1e12))) for i = 1:reps)
end

# auxstats, given data
@views function auxstat(data)
    rets = data[:,1]
    RV = log.(data[:,2])
    BV = log.(data[:,3])
    jump = RV .> (1.5 .* BV)
    nojump = jump .== false
    n = size(data,1)
    # ensure variation
    jump[1:2] .= true
    nojump[1:2] .= true
    # jump stats
    jumpsize = mean(RV[jump]) - mean(BV[jump])
    jumpsize2 = std(rets[jump]) - std(rets[nojump])
    njumps = mean(jump[3:end])
    # ρ
    X = [ones(n-1) BV[2:end] BV[1:end-1]]
    y = rets[2:end]
    βrets = X\y
    ϵrets = y-X*βrets
    σrets = std(ϵrets)
    κrets = std(ϵrets.^2.0)
    # normal volatility: κ, α and σ
    X = [ones(n-2) BV[1:end-2] BV[2:end-1]]
    y = BV[3:end]
    βvol = X\y
    ϵvol = y-X*βvol
    σvol = std(ϵvol)
    κvol = std(ϵvol.^2.0)
    # jump size
    X = [ones(n) jump BV jump.*BV]
    y = RV
    βjump = X\y
    ϵjump = y-X*βjump
    σjump = std(ϵjump)
    κjump = std(ϵjump.^2.0)
    # jump frequency
    qs = quantile(abs.(rets),[0.5, 0.9])
    qs2 = quantile(RV,[0.5, 0.9])
    qs3 = quantile(BV,[0.5, 0.9])
    return sqrt(n)*vcat(βrets, βvol, βjump,σrets, σvol, σjump,κrets, κvol, κjump, mean(RV) - mean(BV), jumpsize, jumpsize2, qs[2]/qs[1], qs2[2]/qs2[1], qs3[2]./qs3[1], qs2 ./ qs3, njumps)
    # brets 1:3
    # bvol 4:6
    # bjump 7:10
    # sigrets 11
    # sigvol 12
    # sigjump 13
    # krets 14
    # kvol 15
    # kjump 16
    # mean RV - mean BV 17
    # jumpsize 18
    # jumpsize2 19
    # qs1 20
    # qs2 21
    # qs3 22
    # qs2/qs3 23-24
    # njumps 25
end

function TrueParameters()
    μ = 0.02
    κ = 0.2
    α = 0.3
    σ = 0.7
    ρ = -0.7
    λ0 = 0.005 # jump rate per day
    λ1 = 4.0 # scaling factor st. dev. of jumps
    τ = 0.005
    return [μ, κ, α, σ, ρ, λ0, λ1, τ]
end

function PriorSupport()
    lb = [-0.1, 0.001, -6.0, 0.5, -0.99, -0.02,  2.0, -0.02]
    ub = [0.1,  0.2, -2.0, 1.5,  -0.5, 0.05, 5.0, 0.05]
    lb,ub
end

# check if parameter is in support.
function InSupport(θ)
    lb,ub = PriorSupport()
    all(θ .>= lb) & all(θ .<= ub)
end

# prior checks that we're in the bounds, and that the unconditional std. dev. of log vol is not too high
# returns 1 if this is true, zero otherwise. Value is not important, as it's constant
function Prior(θ)
    InSupport(θ) ? 1.0 : 0.0
end

function PriorDraw()
    lb, ub = PriorSupport()
    θ = (ub-lb).*rand(size(lb,1)) + lb
end

function PriorMean()
    lb,ub = PriorSupport()
    (ub + lb) ./ 2.0
end

max_stars_repo_name: mcreel/ILSNM

max_issues_repo_name: mcreel/ILSNM

max_forks_repo_name: mcreel/ILSNM


text:
function ei(::Type{T}, i::Integer, n::Integer)::Vector{T} where T
    v = zeros(T, n)
    v[i] = one(T)
    return v
end
ei(i::Integer, n::Integer)::Vector{Float64} = ei(Float64, i, n)

function pivot!(M::Array{T,2}, i::Int, j::Int, Binv::Array{T, 2} = Array{T, 2}(undef, 0, 0); freelines::BitVector = trues(size(M, 1) - 1))::Matrix{T} where T
    updBinv = !iszero(length(Binv))
    m,n = size(M)
    @assert M[i, j] != 0 "pivot  illegal"
    vMi = @view M[i, :]
    updBinv && (vBinv = @view Binv[i, :])

    M[i, :] = vMi / M[i, j]
    updBinv && (Binv[i, :] = vBinv/M[i, j])

    linenottodo = findall(iszero, @view M[:, j])
    guilylines = findall(!, freelines)
    linenottodo[:] = linenottodo
    for (iter, k) in enumerate(setdiff(1:m, i, linenottodo, guilylines))
        #@show M
        #@show Binv
        M[k, :] -= M[k, j] * vMi
        (k != m) && updBinv && (Binv[iter, :] -= M[k, j] * vBinv)
    end
    return M
end

abignumber(::Type{Float64}) = reinterpret(Float64, (reinterpret(UInt64, typemax(Float64)) - UInt64(1)))/100
abignumber(::Type{Float32}) = reinterpret(Float32, (reinterpret(UInt32, typemax(Float32)) - UInt32(1)))/100
#abignumber(::Type{Float16}) = reinterpret(Float16, (reinterpret(UInt16, typemax(Float16)) - UInt16(1)))/100
abignumber(::Type{Rational{Int128}}) = div(typemax(Int128), Int128(1000))//Int128(1)
abignumber(::Type{Rational{Int64}}) = div(typemax(Int64), Int64(1000))//Int64(1)
abignumber(::Type{Rational{Int32}}) = div(typemax(Int32), Int32(1000))//Int32(1)

tol(::Type{T}) where T <: AbstractFloat = 10*eps(T)
tol(::Type{T}) where T <: Rational = T(0)

function isacanonicalcol(v::AbstractVector)::Bool
    m = length(v)
    return (sum(iszero, v) == (m - 1)) && (sum(isone, v) == 1)
end
function iscanonical(A::AbstractMatrix)::Tuple{Bool, Vector{Int}}
    m, n = size(A)
    basis = zeros(Int, m)
    for i in 1:n
        Ai = @view A[:, i]
        if isacanonicalcol(Ai)
            idx = findfirst(Bool.(Ai))
            basis[idx] = i
        end
    end
    return all(!iszero, basis), basis
end

max_stars_repo_name: fbastin/Jasmin

max_issues_repo_name: fbastin/Jasmin

max_forks_repo_name: fbastin/Jasmin


text:
mutable struct State
    a::Int
    b::Int
end

function part1(input)
    program = parse_line.(eachline(input))
    return run(program, a = 0)
end

function part2(input)
    program = parse_line.(eachline(input))
    return run(program, a = 1)
end

function parse_line(line)
    parts = split(replace(line, "," => ""))
    instruction = Symbol(first(parts))
    instruction == :jmp && return (instruction, parse(Int, parts[2]))
    if instruction == :jie || instruction == :jio
        return (instruction, Symbol(parts[2]), parse(Int, parts[3]))
    end
    return (instruction, Symbol(parts[2]))
end

function run(program; a)
    state = State(a, 0)
    n = 1
    while n <= length(program)
        instruction, arguments... = program[n]
        if instruction == :hlf
            setfield!(state, first(arguments), getfield(state, first(arguments)) >> 1)
            n += 1
        elseif instruction == :tpl
            setfield!(state, first(arguments), getfield(state, first(arguments)) * 3)
            n += 1
        elseif instruction == :inc
            setfield!(state, first(arguments), getfield(state, first(arguments)) + 1)
            n += 1
        elseif instruction == :jmp
            n += first(arguments)
        elseif instruction == :jie
            n += getfield(state, first(arguments)) % 2 == 0 ? last(arguments) : 1
        elseif instruction == :jio
            n += getfield(state, first(arguments)) == 1 ? last(arguments) : 1
        end
    end
    return state.b
end


text:
push!(LOAD_PATH, "../../Exercism/julia/primes/src/")
using YA_PRIMES

function consecutive_prime_sum(n::T) where {T <: Integer}
  limit = T(log(n) ÷ log(10)) ÷ 2 + 2 # heuristic ?
  primes = Primes{T}(10^limit)
  #
  # compute all the sum starting from 1st prime number
  # as the sum and the number of terms need to be maximize start scanning by
  # last candidate prime sum (downto)
  #
  max_sum = zero(T)
  sum_primes = [zero(T)]

  for p ∈ primes
    push!(sum_primes, sum_primes[end] + p)
    sum_primes[end] ≥ n && break
  end

  m = T(length(sum_primes))
  start_ix = one(T)
  nterms = one(T)

  for ix ∈ one(T):m
    for jx ∈ m:T(-1):(ix + nterms)
      psum = sum_primes[jx] - sum_primes[ix]

      if (jx - ix) > nterms && isprime(psum)  #
        psum > n && continue
        max_sum, start_ix, nterms = psum, ix, jx - ix
      end
    end
  end

  (max_sum, v(primes)[start_ix:start_ix + nterms - 1], nterms)
end

# ---------------------- A first benchmark - too slow!
# julia> @time consective_prime_sum(10_000)
#   0.008910 seconds (49.92 k allocations: 15.661 MiB)
# (9521, [3, 5, 7, 11, 13, 17, 19, 23, 29, 31  …  269, 271, 277, 281, 283, 293, 307, 311, 313, 317], 65)

# julia> @time consecutive_prime_sum(100_000)
#   0.163134 seconds (490.22 k allocations: 339.731 MiB, 9.41% gc time)
# (92951, [3, 5, 7, 11, 13, 17, 19, 23, 29, 31  …  1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097], 183)

# julia> @time consecutive_prime_sum(1_000_000)
#   6.825734 seconds (5.42 M allocations: 12.378 GiB, 4.85% gc time)
# (997651, [7, 11, 13, 17, 19, 23, 29, 31, 37, 41  …  3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931], 543)

# julia> @time consecutive_prime_sum(2_000_000)
#  21.232955 seconds (11.46 M allocations: 36.335 GiB, 3.30% gc time)
# (1999219, [11, 13, 17, 19, 23, 29, 31, 37, 41, 43  …  5653, 5657, 5659, 5669, 5683, 5689, 5693, 5701, 5711, 5717], 749)

# julia> @time consecutive_prime_sum(5_000_000)
#  88.024431 seconds (30.19 M allocations: 140.041 GiB, 2.52% gc time)
# (4975457, [2, 3, 5, 7, 11, 13, 17, 19, 23, 29  …  9203, 9209, 9221, 9227, 9239, 9241, 9257, 9277, 9281, 9283], 1150)

# ---------------------- A second benchmark - still too slow!
# julia> @time consecutive_prime_sum(10_000)
#   0.000316 seconds (1.52 k allocations: 848.375 KiB)
# (9521, [3, 5, 7, 11, 13, 17, 19, 23, 29, 31  …  269, 271, 277, 281, 283, 293, 307, 311, 313, 317], 65)

# julia> @time consecutive_prime_sum(100_000)
#   0.005229 seconds (10.62 k allocations: 15.791 MiB)
# (92951, [3, 5, 7, 11, 13, 17, 19, 23, 29, 31  …  1039, 1049, 1051, 1061, 1063, 1069, 1087, 1091, 1093, 1097], 183)

# julia> @time consecutive_prime_sum(1_000_000)
#   0.062776 seconds (81.37 k allocations: 342.527 MiB, 22.20% gc time)
# (997651, [7, 11, 13, 17, 19, 23, 29, 31, 37, 41  …  3877, 3881, 3889, 3907, 3911, 3917, 3919, 3923, 3929, 3931], 543)

# julia> @time consecutive_prime_sum(10_000_000)
#   0.664440 seconds (675.23 k allocations: 8.030 GiB, 13.65% gc time)
# (9951191, [5, 7, 11, 13, 17, 19, 23, 29, 31, 37  …  13309, 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399], 1587)

# julia> @time consecutive_prime_sum(100_000_000)
#  20.069132 seconds (11.56 M allocations: 202.315 GiB, 7.09% gc time)
# (99819619, [7, 11, 13, 17, 19, 23, 29, 31, 37, 41  …  45061, 45077, 45083, 45119, 45121, 45127, 45131, 45137, 45139, 45161], 4685)

# julia> @time consecutive_prime_sum(10_000_000)
#   0.644315 seconds (675.16 k allocations: 8.029 GiB, 14.15% gc time)
# (9951191, [5, 7, 11, 13, 17, 19, 23, 29, 31, 37  …  13309, 13313, 13327, 13331, 13337, 13339, 13367, 13381, 13397, 13399], 1587)

# julia> @time consecutive_prime_sum(100_000_000)
#  18.991867 seconds (11.56 M allocations: 202.313 GiB, 6.90% gc time)
# (99819619, [7, 11, 13, 17, 19, 23, 29, 31, 37, 41  …  45061, 45077, 45083, 45119, 45121, 45127, 45131, 45137, 45139, 45161], 4685)

# ---------------------- A better benchmark (more pruning)

# cf. README for better benchmark

max_stars_repo_name: pascal-p/Euler

max_issues_repo_name: pascal-p/Euler

max_forks_repo_name: pascal-p/Euler


text:
#= = # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Description
#
#   This file contains helpers related to space indices.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # = =#
# struct SPI{T} end # TODO
const SPI = Val
# import Base.length
for sym in [:F10, :F10obs, :F10adj, :F10M, :F10Mobs, :F10Madj,
            :Kp, :Ap, :Kp_vect, :Ap_vect,
            :Dst, :Dst_vect,
            :S10, :S81a, :M10, :M81a, :Y10, :Y81a, :DstΔTc]

    qsym = Meta.quot(sym)

    @eval begin
        export $sym
        $sym() = SPI{$qsym}()
        # MethodError: no method matching length(::HTsi.SPI
        # length(x::SPI{$qsym}) = 1
        # MethodError: no method matching iterate(::HTsi.SPI{:F10obs})
        # iterate(x::SPI{$qsym}) = x
    end
end


text:
#=
shape.jl
Types for describing 2D and 3D shapes, like ellipses and ellipsoids.
=#

export AbstractShape, AbstractShape2, AbstractShape3

"""
    AbstractShape
Generic shape type for `ImagePhantoms`.
"""
abstract type AbstractShape end

"""
    AbstractShape2 <: AbstractShape
Generic 2D shape type (with subtypes Ellipse...)
"""
abstract type AbstractShape2 <: AbstractShape end

"""
    AbstractShape2 <: AbstractShape
Generic 3D shape type (with subtypes Ellipsoid...)
"""
abstract type AbstractShape3 <: AbstractShape end