text:
"""
Stores simulation data in memory
"""
mutable struct InMemorySimulationStore <: SimulationStore
    params::SimulationStoreParams
    dm_data::OrderedDict{Symbol, DecisionModelStore}
    em_data::EmulationModelStore
    container_key_lookup::Dict{String, OptimizationContainerKey}
end

function InMemorySimulationStore()
    return InMemorySimulationStore(
        SimulationStoreParams(),
        OrderedDict{Symbol, DecisionModelStore}(),
        EmulationModelStore(),
        Dict{String, OptimizationContainerKey}(),
    )
end

function open_store(
    func::Function,
    ::Type{InMemorySimulationStore},
    directory::AbstractString,  # Unused. Need to match the interface.
    mode=nothing,
    filename=nothing,
)
    store = InMemorySimulationStore()
    return func(store)
end

function Base.empty!(store::InMemorySimulationStore)
    for val in values(get_dm_data(store))
        empty!(val)
    end
    empty!(get_em_data(store))
    @debug "Emptied the store" _group = LOG_GROUP_SIMULATION_STORE
    return
end

Base.isopen(::InMemorySimulationStore) = true
Base.close(::InMemorySimulationStore) = nothing
Base.flush(::InMemorySimulationStore) = nothing
get_params(store::InMemorySimulationStore) = store.params

function get_decision_model_params(store::InMemorySimulationStore, model_name::Symbol)
    return get_params(store).decision_models_params[model_name]
end

get_container_key_lookup(store::InMemorySimulationStore) = store.container_key_lookup

list_decision_models(x::InMemorySimulationStore) = collect(keys(x.dm_data))
log_cache_hit_percentages(::InMemorySimulationStore) = nothing

function list_decision_model_keys(
    store::InMemorySimulationStore,
    model_name::Symbol,
    container_type::Symbol,
)
    return list_fields(_get_model_results(store, model_name), container_type)
end

function list_emulation_model_keys(store::InMemorySimulationStore, container_type::Symbol)
    return list_fields(store.em_data, container_type)
end

function write_optimizer_stats!(
    store::InMemorySimulationStore,
    model_name,
    stats::OptimizerStats,
    index::DecisionModelIndexType,
)
    write_optimizer_stats!(get_dm_data(store)[model_name], stats, index)
    return
end

function write_optimizer_stats!(
    store::InMemorySimulationStore,
    stats::OptimizerStats,
    index::EmulationModelIndexType,
)
    write_optimizer_stats!(get_em_data(store), stats, index)
    return
end

function write_result!(
    store::InMemorySimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
    update_timestamp::Dates.DateTime,
    array,
)
    write_result!(
        get_dm_data(store)[model_name],
        model_name,
        key,
        index,
        update_timestamp,
        array,
    )
    return
end

function write_result!(
    store::InMemorySimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
    update_timestamp::Dates.DateTime,
    array,
)
    write_result!(get_em_data(store), model_name, key, index, update_timestamp, array)
    return
end

function read_optimizer_stats(store::InMemorySimulationStore, model_name)
    # TODO EmulationModel: this interface is TBD
    return read_optimizer_stats(get_dm_data(store)[model_name])
end

function initialize_problem_storage!(
    store::InMemorySimulationStore,
    params::SimulationStoreParams,
    dm_problem_reqs::Dict{Symbol, SimulationModelStoreRequirements},
    em_problem_reqs::SimulationModelStoreRequirements,
    ::CacheFlushRules,
)
    store.params = params
    for problem in keys(store.params.decision_models_params)
        get_dm_data(store)[problem] = DecisionModelStore()
        for type in STORE_CONTAINERS
            for (key, reqs) in getfield(dm_problem_reqs[problem], type)
                container = getfield(get_dm_data(store)[problem], type)
                container[key] =
                    OrderedDict{Dates.DateTime, DatasetContainer{DataFrameDataset}}()
                store.container_key_lookup[encode_key_as_string(key)] = key
                @debug "Added $type $key in $problem" _group = LOG_GROUP_SIMULATION_STORE
            end
        end
    end

    for type in STORE_CONTAINERS
        for (key, reqs) in getfield(em_problem_reqs, type)
            container = get_data_field(get_em_data(store), type)
            container[key] = DataFrameDataset(
                DataFrames.DataFrame(
                    OrderedDict(c => fill(NaN, reqs["dims"][1]) for c in reqs["columns"]),
                ),
            )
            store.container_key_lookup[encode_key_as_string(key)] = key
            @debug "Added $type $key in emulation store" _group = LOG_GROUP_SIMULATION_STORE
        end
    end

    return
end

function read_result(
    ::Type{DataFrames.DataFrame},
    store::InMemorySimulationStore,
    model_name::Symbol,
    key::OptimizationContainerKey,
    index::DecisionModelIndexType,
)
    return read_results(get_dm_data(store)[model_name], key, index=index)
end

function read_result(
    ::Type{DataFrames.DataFrame},
    store::InMemorySimulationStore,
    ::Symbol,
    key::OptimizationContainerKey,
    index::EmulationModelIndexType,
)
    return read_results(get_em_data(store), key, index=index)
end

function read_results(
    store::InMemorySimulationStore,
    key::OptimizationContainerKey;
    index::EmulationModelIndexType=nothing,
    len::Int=nothing,
)
    return read_results(get_em_data(store), key, index=index, len=len)
end

function get_emulation_model_dataset_size(
    store::InMemorySimulationStore,
    key::OptimizationContainerKey,
)
    return get_dataset_size(get_em_data(store), key)
end

# Note that this function is not type-stable.
function _get_model_results(store::InMemorySimulationStore, model_name::Symbol)
    if model_name in keys(get_dm_data(store))
        results = get_dm_data(store)
    else
        # TODO EmulationModel: this interface is TBD
        error("model name $model_name is not stored")
    end

    return results[model_name]
end

function write_optimizer_stats!(
    store::InMemorySimulationStore,
    model::DecisionModel,
    index::DecisionModelIndexType,
)
    stats = get_optimizer_stats(model)
    dm_data = get_dm_data(store)
    write_optimizer_stats!(dm_data[get_name(model)], stats, index)
    return
end

function write_optimizer_stats!(
    store::InMemorySimulationStore,
    model::EmulationModel,
    index::EmulationModelIndexType,
)
    stats = get_optimizer_stats(model)
    em_data = get_em_data(store)
    write_optimizer_stats!(em_data, stats, index)
    return
end


text:
import Pluto
import Pluto.ExploreExpression: SymbolsState, compute_symbolreferences

"Calls `ExploreExpression.compute_symbolreferences` on the given `expr` and test the found SymbolsState against a given one, with convient syntax.

# Example

```jldoctest
julia> @test testee(:(
    begin
        a = b + 1
        f(x) = x / z
    end),
    [:b, :+], # 1st: expected references
    [:a, :f], # 2nd: expected definitions
    [:+],     # 3rd: expected function calls
    [
        :f => ([:z, :/], [], [:/], [])
    ])        # 4th: expected function definitions, with inner symstate using the same syntax
true
```
"
function testee(expr, expected_references, expected_definitions, expected_funccalls, expected_funcdefs; verbose::Bool=true)
    expected = easy_symstate(expected_references, expected_definitions, expected_funccalls, expected_funcdefs)
    result = compute_symbolreferences(expr)

    # Anonymous function are given a random name, which looks like anon67387237861123
    # To make testing easier, we rename all such functions to anon
    new_name(sym) = startswith(string(sym), "anon") ? :anon : sym

    result.assignments = Set(new_name.(result.assignments))
    result.funcdefs = let
        newfuncdefs = Dict()
        for (k, v) in result.funcdefs
            newfuncdefs[new_name.(k)] = v
        end
        newfuncdefs
    end

    if verbose && expected != result
        println()
        println("FAILED TEST")
        println(expr)
        println()
        dump(expr, maxdepth=20)
        println()
        @show expected
        resulted = result
        @show resulted
        println()
    end
    return expected == result
end

function easy_symstate(expected_references, expected_definitions, expected_funccalls, expected_funcdefs)
    new_expected_funccalls = map(expected_funccalls) do k
        new_k = k isa Symbol ? [k] : k
        return new_k
    end |> Set

    new_expected_funcdefs = map(expected_funcdefs) do (k, v)
        new_k = k isa Symbol ? [k] : k
        new_v = v isa SymbolsState ? v : easy_symstate(v...)
        return new_k => new_v
    end |> Dict

    SymbolsState(Set(expected_references), Set(expected_definitions), new_expected_funccalls, new_expected_funcdefs)
end

function occursinerror(needle, haystack::Pluto.Cell)
    return haystack.errored && occursin(needle, haystack.output_repr)
end

"Test notebook equality, ignoring cell UUIDs and such."
function notebook_inputs_equal(nbA, nbB)
    x = normpath(nbA.path) == normpath(nbB.path)

    to_compare(cell) = (cell.cell_id, cell.code)
    y = to_compare.(nbA.cells) == to_compare.(nbB.cells)

    x && y
end

"Whether the given .jl file can be run without any errors. While notebooks cells can be in arbitrary order, their order in the save file must be topological.

If `only_undefvar` is `true`, all errors other than an `UndefVarError` will be ignored."
function jl_is_runnable(path; only_undefvar=false)
    🔖 = Symbol("lab", time_ns())
    🏡 = Core.eval(Main, :(module $(🔖) end))
    try
        Core.eval(🏡, :(include($path)))
        true
    catch ex
        if (!only_undefvar) || ex isa UndefVarError || (ex isa LoadError && ex.error isa UndefVarError)
            println(stderr, "\n$(path) failed to run. File contents:")

            println(stderr, "\n\n\n")
            println.(enumerate(readlines(path; keep=true)))
            println(stderr, "\n\n\n")

            showerror(stderr, ex, stacktrace(catch_backtrace()))
            println(stderr)
            false
        else
            true
        end
    end
end

"Whether the `notebook` runs without errors."
function nb_is_runnable(notebook::Pluto.Notebook)
    run_reactive!(notebook, notebook.cells)
    errored = filter(c -> c.errored, notebook.cells)
    if !isempty(errored)
        @show errored
    end
    isempty(errored)
end

"The converse of Julia's `Base.sprint`."
function sread(f::Function, input::String, args...)
    io = IOBuffer(input)
    output = f(io, args...)
    close(io)
    return output
end

function num_backups_in(dir::AbstractString)
    count(readdir(dir)) do fn
        occursin("backup", fn)
    end
end

ext: jl


text:
using Manifolds, Manopt, Test, ManifoldsBase, Dates

struct TestProblem <: Problem end
struct TestOptions <: Options end

@testset "StoppingCriteria" begin
    struct myStoppingCriteriaSet <: StoppingCriterionSet end
    @test_throws ErrorException get_stopping_criteria(myStoppingCriteriaSet())

    s = StopWhenAll(StopAfterIteration(10),StopWhenChangeLess(0.1))
    s2 = StopWhenAll([StopAfterIteration(10),StopWhenChangeLess(0.1)])
    @test get_stopping_criteria(s)[1].maxIter == get_stopping_criteria(s2)[1].maxIter

    s3 = StopWhenCostLess(0.1)
    p = GradientProblem(Euclidean(1), x->x^2, x->2x)
    o = GradientDescentOptions(1.0)
    @test !s3(p,o,1)
    @test length(s3.reason) == 0
    o.x = 0.3
    @test s3(p,o,2)
    @test length(s3.reason) > 0
    # repack
    sn = StopWhenAny(StopAfterIteration(10), s3)
    # then s3 is the only active one
    @test get_active_stopping_criteria(sn) == [s3]
    @test get_active_stopping_criteria(s3) == [s3]
    @test get_active_stopping_criteria(StopAfterIteration(1)) == []
    sm = StopWhenAll(StopAfterIteration(10),s3)
    @test !sm(p,o,9)
    @test sm(p,o,11)
    an = sm.reason
    m = match(r"^((.*)\n)+",an)
    @test length(m.captures)==2 # both have to be active
end

@testset "Test StopAfter" begin
    p = TestProblem()
    o = TestOptions()
    s = StopAfter(Second(1))
    s(p,o,0) # Start
    @test s(p,o,1) == false
    sleep(1.02)
    @test s(p,o,2) == true
end


text:
export sym

"""
    sym(x, s)

Evaluate a symbol within the context of an object, `x`. By default uses
`getindex`, which allows for use with `Dict`s, `NamedTuple`s, `Tuple`s, and
`Array`s. `Pair`s are handled as a special case, calling `getfield` instead.
`Pair`s also include handling for use of `:1` and `:2` as shorthands for
`:first` and `:second`.

### Examples

```jdoctext
julia> sym(Dict(:x => 1), :x)
1

julia> sym(3 => 4, :1 + :second)
7

julia> sym((x = 1, y = 2), :x + :y)
3
"""
sym(x, s) = getindex(x, s)

# Pair
sym(x::Pair, s::Symbol) = getfield(x, s)
sym(x::Pair, s::Int) = getfield(x, [:first, :second][s])

# Dict
sym(x::Dict{String,}, s::Symbol) = getindex(x, string(s))

"""
"""
function match_syms(x, ss)
    s_in_x = syms_in_context(x, ss)
    ss[s_in_x], ss[.! s_in_x]
end

syms_in_context(x, ss) =
    repeat([true], length(ss))

syms_in_context(x::Union{Dict,NamedTuple,Base.Iterators.Pairs}, ss) =
    in.(ss, [keys(x)])

syms_in_context(x::Pair, ss) =
    in.(ss, [[:first, :second, :1, :2]])

syms_in_context(x::Array, ss) =
    indexin.(ss, [x]) .!== nothing

syms_in_context(x::Tuple, ss) =
    1 .<= ss .<= length(x)


text:
# [[file:../fresa.org::testminlp][testminlp]]
using Distributed
using Printf
@everywhere using Fresa
# define new type for mixed integer problems
# in general, this would be vectors of real and integer values
@everywhere struct MI
    x :: Float64
    y :: Int32
end
import Base
Base.show(io::IO, m::MI) = print(io, m.x, " ", m.y)
f = s -> (3s.y - 5s.x,
          max(2s.y + 3s.x - 24,
              3s.x - 2s.y - 8,
              2s.y^2 - 2*√s.y + 11s.y + 8s.x - 39 - 2*√s.x*s.y^2))
# bounds
domain = Fresa.Domain(x -> MI(1.0, 1),
                      x -> MI(6.0, 6))
# function to find a neighbouring solution for MI type decision points
function Fresa.neighbour(s :: MI,
                         a :: MI,
                         b :: MI,
                         f :: Float64) :: MI
    x = s.x + (b.x-a.x)*(1-f)*2*(rand()-0.5)
    x = x < a.x ? a.x : (x > b.x ? b.x : x)
    # for the integer variable, we move in one direction or the other
    # a random number of places depending on fitness
    positive = rand(Bool)
    r = rand()
    # @printf(": neighbour: f=%g r=%g\n", f, r)
    inc = ceil(f*r*(b.y-a.y)/2)
    # @printf(": neighbour: positive=%s inc=%d\n", positive, inc)
    y = s.y + (positive ? inc : -inc)
    y = y < a.y ? a.y : (y > b.y ? b.y : y)
    return MI(x,y)
end
# create the initial population consisting of a single MI point
p0 = [Fresa.createpoint(MI(1.0, 1),f)]
# now invoke Fresa to solve the problem
best, pop = Fresa.solve(f, p0, domain; ngen=100)
println("Population: $pop")
println("Best: f($(best.x)) = $(best.z), $(best.g)")
# testminlp ends here

# [[file:../fresa.org::testminlpsupplement][testminlpsupplement]]
println("#+plot: ind:3 deps:(2) with:\"linespoints pt 7\" set:nokey set:\"yrange [0:1]\"")
ancestor = best.ancestor;
while ancestor != Some(nothing) && typeof(ancestor) != Nothing
    global ancestor
    println("| $(ancestor.point.z) | $(ancestor.fitness) | $(ancestor.generation) |")
    ancestor = ancestor.point.ancestor
end
# testminlpsupplement ends here


text:
using LinearAlgebra

"""
    normal_hessenberg_matrix(Float64, vals::AbstractVector)

Generate a normal hessenberg matrix with eigenvalues `vals`
"""
function normal_hessenberg_matrix(T::Type{<:Number}, vals::AbstractVector)
    n = length(vals)
    Q, R = qr(randn(T, n, n))
    A = Q * Diagonal(vals) * Q'
    return triu(hessenberg!(A).factors, -1)
end

function normal_hessenberg_matrix(T::Type{<:Real}, vals::AbstractVector{<:Complex})
    n = length(vals)
    Q, R = qr(randn(T, n, n))
    D = zeros(T, n, n)
    i = 1
    while i ≤ n
        if imag(vals[i]) != 0
            D[i+0,i+0] = real(vals[i])
            D[i+1,i+0] = imag(vals[i])
            D[i+0,i+1] = -imag(vals[i])
            D[i+1,i+1] = real(vals[i])
            i += 2
        else
            D[i] = real(vals[i])
            i += 1
        end
    end
    return triu(hessenberg!(Q * D * Q').factors, -1)
end

"""
    realimag(1 + 3im) → (1, 3)

Split imaginary number into a tuple of real and imaginary part
"""
realimag(x) = (real(x), imag(x))

"""
    is_hessenberg(H) → bool

Test whether the sub-subdiagonals of H are zero.
"""
is_hessenberg(H) = norm(tril(H, -2)) == 0


text:
## All steps of batch processing
# Data from previous steps is always read from file, so any step can
# be run separately from previous steps
using EstimateSDE
using Random
using Serialization
using Base.Threads
using Pipe: @pipe
using DataFrames

## Experiment name —————————————
name = "experiment02"

## Generate settings —————————————
S = (
  data=BatchDataSettings(;
    rng=MersenneTwister(9845789),
    name=name,
    # proclabels=[:ICMLA_W,:ICMLA_A,:ICMLA_B,:ICMLA_C],
    proclabels=[:ICMLA_A,:ICMLA_B,:ICMLA_C],
    Tav=1.0,
    Tmin=0.01,
    Ns=[500],
    Nrun=25
  ),
  mle_trace=BatchMlEstimation(;
    p=8,
    Tkld=0.2,
    engines=[ElimQ,SdeRootsEngine],
    init_actuals=[false]
  )
)

results_dir = settings_resultsdir(name)
f = joinpath(results_dir,"settings.jld")
serialize(f,S); @info "Settings saved to $f"

## Load settings —————————————
results_dir = settings_resultsdir(name)
S = deserialize(joinpath(results_dir,"settings.jld"))

## Thread check —————————————
nt = nthreads()
if nt==1
  @warn "nthreads==1. Increase number of threads for faster processing"
else
  @info "Threads used: $nt"
end

## Data —————————————
printstyled("Data set generation settings:\n",color=:blue)
println(S.data)
output_data = @time batch_data(S.data)
nothing

## Maximum Likelihood trace —————————————
# Settings
printstyled("\nSettings MLE trace\n",color=:blue)
print(S.mle_trace)

df_data = deserialize(joinpath(results_dir,"data.jld")).df

@info "Maximum Likelihood trace"
output_mle_trace = @time batch_mle_trace(S.mle_trace,df_data,results_dir)

## Posterior samples —————————————
# Skipped

## Derived estimates —————————————
df_mle_trace = deserialize(joinpath(results_dir,"mle-trace.jld")).df

@info "Maximum Likelihood derived estimates"
df_deriv = derived_estimates(df_mle_trace,results_dir)

## Model evaluation —————————————
df_deriv = deserialize(joinpath(results_dir,"derived-estimates.jld")).df
df_data = deserialize(joinpath(results_dir,"data.jld")).df

@info "Model evaluation"
output_eval = @time batch_model_evaluation(df_data,S.mle_trace,df_deriv,results_dir)

## Aggregate —————————————
@info "Aggregation"
include("aggregate.jl")

## Report —————————————
@info "Report"
include("figures.jl")


text:
function as_sub(x::AbstractMatrix)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = subview(y, 2:2:size(y,1), 2:2:size(y,2))
    for j=1:size(x,2)
        for i=1:size(x,1)
            y[i,j] = x[i,j]
        end
    end
    y
end

as_sub(eye(2))


text:
using CUDA
using DLPack
using PyCall
using Test

@info "Installing JAX"
#
packages = "jax[cpu]"
run(PyCall.python_cmd(`-m pip install $packages`))

# Load JAX
jax = pyimport("jax")
np = pyimport("jax.numpy")
dlpack = pyimport("jax.dlpack")

@testset "DLPack.jl" begin
    v = np.asarray([1.0, 2.0, 3.0])

    dlv = DLVector{Float32}(@pycall dlpack.to_dlpack(v)::PyObject)
    tensor = dlv.manager.dl_tensor

    @test tensor.ndim == 1
    @test tensor.dtype == DLPack.jltypes_to_dtypes()[Float32]

    if DLPack.device_type(dlv) == DLPack.kDLCPU
        jv = unsafe_wrap(Array, dlv)
        jv[1] = 0  # mutate a jax's tensor
    elseif DLPack.device_type(dlv) == DLPack.kDLGPU
        jv = unsafe_wrap(CuArray, dlv)
        jv[1:1] .= 0  # mutate a jax's tensor
    end

    @test py"$np.all($v[:] == $np.asarray([0.0, 2.0, 3.0])).copy().item()"
end


text:
module DiffEqDevTools

using DiffEqBase, RecipesBase, RecursiveArrayTools, DiffEqNoiseProcess

using LinearAlgebra, Distributed

using Statistics

import Base: length

import DiffEqBase: AbstractODEProblem, AbstractDDEProblem,
       AbstractODESolution, AbstractRODEProblem, AbstractSDEProblem,
       AbstractSDDEProblem, AbstractEnsembleProblem,
       AbstractDAEProblem, @def, ConvergenceSetup, DEAlgorithm,
       ODERKTableau, AbstractTimeseriesSolution, ExplicitRKTableau,
       ImplicitRKTableau

import LinearAlgebra: norm, I

const TIMESERIES_ERRORS = Set([:l2,:l∞,:L2,:L∞])
const DENSE_ERRORS = Set([:L2,:L∞])
const WEAK_TIMESERIES_ERRORS = Set([:weak_l2,:weak_l∞])
const WEAK_DENSE_ERRORS = Set([:weak_L2,:weak_L∞])
const WEAK_ERRORS = union(Set([:weak_final]),
                          WEAK_TIMESERIES_ERRORS,WEAK_DENSE_ERRORS)

include("benchmark.jl")
include("convergence.jl")
include("plotrecipes.jl")
include("test_solution.jl")
include("ode_tableaus.jl")
include("tableau_info.jl")

export ConvergenceSimulation, Shootout, ShootoutSet, TestSolution

#Benchmark Functions
export Shootout, ShootoutSet, WorkPrecision, WorkPrecisionSet

export test_convergence, analyticless_test_convergence, appxtrue!, appxtrue

export get_sample_errors

#Tab Functions
export stability_region, residual_order_condition, check_tableau

#Tableaus
export deduce_Butcher_tableau
export constructEuler, constructKutta3, constructRK4, constructRK438Rule,
       constructImplicitEuler, constructMidpointRule, constructTrapezoidalRule,
       constructLobattoIIIA4, constructLobattoIIIB2, constructLobattoIIIB4,
       constructLobattoIIIC2, constructLobattoIIIC4, constructLobattoIIICStar2,
       constructLobattoIIICStar4, constructLobattoIIID2, constructLobattoIIID4,
       constructRadauIA3, constructRadauIA5,
       constructRadauIIA3, constructRadauIIA5,
       constructRalston, constructRalston4, constructHeun, constructRKF5, constructBogakiShampine3,
       constructCashKarp, constructRKF8, constructDormandPrince8,
       constructFeagin10, constructFeagin12, constructFeagin14,
       constructDormandPrince8_64bit, constructRKF5, constructRungeFirst5,
       constructCassity5, constructLawson5,
       constructLutherKonen5, constructLutherKonen52,
       constructLutherKonen53, constructPapakostasPapaGeorgiou5,
       constructPapakostasPapaGeorgiou52,
       constructTsitouras5, constructBogakiShampine5, constructSharpSmart5,
       constructButcher6, constructButcher7,
       constructDverk, constructClassicVerner6,
       constructClassicVerner7, constructClassicVerner8,
       constructVernerRobust7, constructEnrightVerner7,
       constructTanakaYamashitaStable7,
       constructTanakaYamashitaEfficient7, constructSharpSmart7,
       constructSharpVerner7,
       constructVernerEfficient7, constructCooperVerner8,
       constructCooperVerner82,
       constructTsitourasPapakostas8, constructdverk78, constructEnrightVerner8,
       constructCurtis8, constructVernerRobust9, constructVernerEfficient9,
       constructSharp9, constructTsitouras9,
       constructTsitouras92,
       constructOno12, constructCurtis10,
       constructOno10,
       constructCurtis10, constructBaker10,
       constructHairer10, constructButcher63,
       constructButcher6, constructButcher62,
       constructVerner6, constructDormandPrince6,
       constructSharpVerner6, constructVerner9162,
       constructVerner916, constructVernerRobust6,
       constructVernerEfficient6, constructPapakostas6, constructLawson6,
       constructTsitourasPapakostas6, constructDormandLockyerMcCorriganPrince6,
       constructTanakaKasugaYamashitaYazaki6D,
       constructTanakaKasugaYamashitaYazaki6C,
       constructTanakaKasugaYamashitaYazaki6B,
       constructTanakaKasugaYamashitaYazaki6A,
       constructMikkawyEisa, constructChummund6, constructChummund62,
       constructHuta62, constructHuta6, constructRKF4,
       constructVerner7, constructVerner8,
       constructVerner6, constructSSPRK22, constructSSPRK33,
       constructSSPRK43, constructSSPRK104,
       constructRKO65

end # module


text:
using ForwardDiff
using ReverseDiff: track, value, gradient, TrackedVector, TrackedMatrix, TrackedArray
using Test

@testset "fill" begin
    v = fill(track(1), 3)
    @test v isa TrackedArray
    @test value(v) == fill(1, 3)
    @test gradient(x -> sum(fill(x[1], 3)), rand(1)) == [3.0]
end

@testset "any & all" begin
    @test all(iszero, track(zeros(3)))
    @test !all(iszero, track([zeros(2); 1.0]))
    @test !any(iszero, track(ones(3)))
    @test any(iszero, track([ones(2); 0.0]))
end

function testcat(f, args::Tuple, type, kwargs=NamedTuple())
    x = f(track.(args)...; kwargs...)
    @test x isa type
    @test value(x) == f(args...; kwargs...)

    if length(args) == 1
        x = f(track(args[1]); kwargs...)
        @test x isa type
        @test value(x) == f(args...; kwargs...)
    else
        @assert length(args) == 2
        x = f(track(args[1]), args[2]; kwargs...)
        @test x isa type
        @test value(x) == f(args...; kwargs...)

        x = f(args[1], track(args[2]); kwargs...)
        @test x isa type
        @test value(x) == f(args...; kwargs...)
    end

    args = (args..., args...)
    x = f(track.(args)...; kwargs...)
    @test x isa type
    @test value(x) == f(args...; kwargs...)

    sizes = size.(args)
    F = vecx -> sum(f(unpack(sizes, vecx)...; kwargs...))
    X = pack(args)
    @test ForwardDiff.gradient(F, X) == gradient(F, X)
end
function pack(xs)
    return mapreduce(vcat, xs) do x
        x isa Number ? x : vec(x)
    end
end
function unpack(sizes, vecx)
    start = 0
    out = map(sizes) do s
        if s === ()
            x = vecx[start+1]
            start += 1
        else
            x = reshape(vecx[start+1:start+prod(s)], s)
            start += prod(s)
        end
    end
    return out
end

@testset "cat" begin
    v = rand(3)
    m = rand(3,3)
    a = rand(3,3,3)
    n = rand()

    testcat(cat, (n,), TrackedVector, (dims=1,))
    testcat(cat, (n, n), TrackedVector, (dims=1,))
    testcat(cat, (n, n), TrackedMatrix, (dims=2,))
    testcat(cat, (v, n), TrackedVector, (dims=1,))
    testcat(cat, (n, v), TrackedVector, (dims=1,))

    testcat(cat, (v, v), TrackedVector, (dims=1,))
    testcat(cat, (v, v), TrackedMatrix, (dims=2,))
    testcat(cat, (v, m), TrackedMatrix, (dims=2,))
    testcat(cat, (m, v), TrackedMatrix, (dims=2,))

    testcat(cat, (a, a), TrackedArray, (dims=1,))
    testcat(cat, (a, a), TrackedArray, (dims=2,))
    testcat(cat, (a, a), TrackedArray, (dims=3,))
    testcat(cat, (a, m), TrackedArray, (dims=3,))

    testcat(vcat, (n,), TrackedVector)
    testcat(vcat, (n, n), TrackedVector)
    testcat(vcat, (v, n), TrackedVector)
    testcat(vcat, (n, v), TrackedVector)
    testcat(vcat, (v, v), TrackedVector)

    testcat(hcat, (n,), TrackedMatrix)
    testcat(hcat, (n, n), TrackedMatrix)
    testcat(hcat, (v, v), TrackedMatrix)
    testcat(hcat, (v, m), TrackedMatrix)
    testcat(hcat, (m, v), TrackedMatrix)
end


text:
"""
```
function solve_slm_system(Mdes::AbstractMatrix{S}, rhs::AbstractVector{S},
                          Mreg::AbstractMatrix{S}, rhsreg::AbstractVector{S},
                          λ::S, Meq::AbstractMatrix{S}, rhseq::AbstractVector{S},
                          Mineq::AbstractMatrix{S}, rhsineq::AbstractVector{S};
                          method::Symbol = :ipopt, atol::S = NaN, rtol::S = NaN,
                          max_iter::Int = 0, max_time::S = 30.,
                          init::AbstractVector{S} = Vector{S}(undef, 0),
                          use_lls::Bool = true, use_sparse::Bool = true,
                          verbose::Symbol = :low) where {S <: Real}
```

estimates the coefficient of a least squares spline by solving the
constrained linear least squares problem

```
min ½ || Mfit * x - rhsfit ||²
```

subject to the constraints

```
Meq * x = rhseq,    Mineq * x ≤ rhsineq,
```

where `Mfit = vcat(Mdes, λ .* Mreg)` and `rhsfit = vcat(rhs, λ .* rhsreg)`.

### Inputs
- `Mdes`: design matrix in system of equations defining the spline
- `rhs`: RHS of the system of equations defining the spline
- `Mreg`: matrix in system of equations regularizing the spline's smoothness
- `rhsreg`: RHS of the system of equations regularizing the spline smoothnes
- `λ`: degree of smoothness (smaller ⇒ smoother)
- `Meq`: matrix for equality constraints
- `rhseq`: the RHS of equality constraints
- `Mineq`: matrix for inequality constraints
- `rhsineq`: upper bound for inequality constraints

### Keywords
- `method`: minimization algorithm, currently only the interior point optimizer is used (Ipopt).
    New methods are welcome! One algorithm that may be added is adapting the matrices to the trust-region method
    via the calculation of pseudo-inverses
- `atol`: absolute tolerance for minimization algorithm
- `rtol`: relative tolerance for minimization algorithm
- `max_iter`: maximum number of evaluations
- `max_time`: maximum time (in minutes) for optimization
- `init`: initial guess of coefficients for the spline
- `use_lls`: use an `ADNLSModel` to define the problem. Currently, using an `LLSModel` (or equivalent) does not work.
    See the `NLPModels.jl` package.
- `use_sparse`: use sparse matrices (assuming the inputs are not already sparse matrices)
- `verbose`: how much information will be printed. Can be `:low` or `:high`
"""
function solve_slm_system(Mdes::AbstractMatrix{S}, rhs::AbstractVector{S},
                          Mreg::AbstractMatrix{S}, rhsreg::AbstractVector{S},
                          λ::S, Meq::AbstractMatrix{S}, rhseq::AbstractVector{S},
                          Mineq::AbstractMatrix{S}, rhsineq::AbstractVector{S};
                          method::Symbol = :ipopt, atol::S = 1e-8, rtol::S = 1e-6,
                          max_iter::Int = 3000, max_time::S = 3600.,
                          init::AbstractVector{S} = Vector{S}(undef, 0),
                          use_lls::Bool = true, use_sparse::Bool = true,
                          verbose::Symbol = :low, kwargs...) where {S <: Real}

    if use_sparse && all([issparse(x) for x in [Mdes, Mreg, Meq, Mineq]])
        # Don't spend time converting already sparse matrices
        use_sparse = false
    end

    Mfit = use_sparse ? vcat(sparse(Mdes), sparse(Mreg) .* λ) : vcat(Mdes, Mreg .* λ)
    rhsfit = vcat(rhs, rhsreg .* λ)

    if verbose == :high
        println("Condition number of the regression: $(cond(Mdes))")
    end

    if isempty(Mineq) && isempty(Meq)
        # Backslash suffices in this case
        coef = Mfit \ rhsfit
    else
        # Set up NLPModel and constraints
        if isempty(Mineq)
            lcon = rhseq
            ucon = rhseq
            C = use_sparse ? sparse(Meq) : Meq
        else
            lcon = vcat(fill(-Inf, length(rhsineq)), rhseq) # -Inf for the inequalities since that's how lsqlin is written in MATLAB
            ucon = vcat(rhsineq, rhseq)
            C = use_sparse ? vcat(sparse(Mineq), sparse(Meq)) : vcat(Mineq, Meq)
        end

        if use_lls
            nls = LLSModel(Mfit, rhsfit; C = C, lcon = lcon, ucon = ucon)
        else
            if isempty(init)
                init = zeros(S, size(Mfit, 2))
            end

            nls = ADNLSModel(x -> Mfit * x - rhsfit, init, size(Mfit, 1), c = x -> C * x, lcon = lcon, ucon = ucon)
        end

        # Decide solver for minimizing the constrained linear least-squares problem
        if method == :ipopt
            if use_lls
                # Feasibility Form adds additional coefficient which we don't want, hence
                # the [1:size(Mdes, 2)] slicing.
                coef = ipopt(FeasibilityFormNLS(nls); print_level = verbose == :high ? 3 : 0,
                             tol = atol, acceptable_tol = rtol, max_iter = max_iter,
                             max_cpu_time = max_time, kwargs...).solution[1:size(Mdes, 2)]
            else
                #
                coef = ipopt(nls; print_level = verbose == :high ? 3 : 0, tol = atol, acceptable_tol = rtol,
                             max_iter = max_iter, max_cpu_time = max_time, kwargs...).solution
            end
        else
            error("Method $method not is not implemented.")
        end
    end

    if verbose == :high
        println("Solver employed: $(out_info[:solver])")
        if !isempty(Mineq)
            println(nls)
        end
    end

    return coef
end


text:
module SfCoriolis

using Unitful, UnitfulAngles, Documenter
import LinearAlgebra: cross

import ..SfGeometry: radius, unit_x, unit_y, unit_z
import ..SfUnits: AngularVelocity, AngularMomentum, Angle, to_angle

export coriolis_acceleration, coriolis_force, tangential_velocity, centrifugal_acceleration, centrifugal_force,
	angular_momentum, angular_velocity, radius

DocMeta.setdocmeta!(SfCoriolis, :DocTestSetup, :(using Unitful, UnitfulAstro, ..SfCoriolis, ..SfGeometry); recursive=true)

# this can't be a jldoctest example right now, as there's some regex/string comparison error that breaks it.
"""
    coriolis_acceleration(Ω::Vector, v::Vector)

The apparent acceleration due to coriolis effects for an object moving with velolcity
vector `v` in a rotating reference frame with angular velocity vector `Ω`.

# Example

```julia-repl
julia> coriolis_acceleration(0.001u"rad/s" * unit_y, 400u"m/s" * unit_x)
3-element Array{Quantity{Float64,𝐋  𝐓 ^-2,Unitful.FreeUnits{(m, s^-2),𝐋  𝐓 ^-2,nothing}},1}:
  0.0 m s^-2
  0.0 m s^-2
 -0.8 m s^-2
```
"""
coriolis_acceleration(Ω::Vector, v::Vector) = upreferred.(2cross(Ω, v))

"""
    coriolis_acceleration(Ω::Vector, v::Vector)

The strength of the coriolis force apparently acting on an object with mass `m` moving with velolcity
vector `v` in a rotating reference frame with angular velocity vector `Ω`.
"""
coriolis_force(m::Unitful.Mass, Ω::Vector, v::Vector) = m * coriolis_acceleration(Ω, v)

"""
    tangential_velocity(ω::AngularVelocity, r::Unitful.Length)

Calculate the tangential velocity for an object rotating about a point at a distance `r` with angular velocity `ω`.
"""
tangential_velocity(ω::AngularVelocity, r::Unitful.Length) = upreferred(ω * r)

"""
    tangential_velocity(ω::AngularVelocity, r::Unitful.Length, λ::Angle)

Calculate the tangential velocity of an object at latitude `λ` on a rotating sphere of radius `r` with angular velocity `ω`.
"""
tangential_velocity(ω::AngularVelocity, r::Unitful.Length, λ::Angle) = upreferred(ω * r * cos(λ))

"""
    centrifugal_acceleration(ω::AngularVelocity, r::Unitful.Length)

Calculate the centrifugal acceleration for an object rotating about a point at a distance `r` with angular velocity `ω`.
"""
centrifugal_acceleration(ω::AngularVelocity, r::Unitful.Length) = upreferred(ω^2 * r)

"""
    angular_velocity(ca::Unitful.Acceleration, r::Unitful.Length)

Calculate the angular velocity required to develop a centrifugal acceleration of `ca` at a radius of `r`.
"""
angular_velocity(ca::Unitful.Acceleration, r::Unitful.Length) = sqrt(ca / r) |> u"rad/s"

"""
    radius(ω::AngularVelocity, ca::Unitful.Acceleration)

Calculate the radius of rotation of a body with angular velocity `ω` to develop a centrifugal acceleration of `ca`.
"""
radius(ω::AngularVelocity, ca::Unitful.Acceleration) = upreferred(ca / ω^2)

"""
    centrifugal_acceleration(ω::AngularVelocity, r::Unitful.Length, λ::Angle)

Calculate the radial part of the centrifugal acceleration of object at latitude `λ` on a rotating sphere of radius `r` with angular velocity `ω`.
"""
centrifugal_acceleration(ω::AngularVelocity, r::Unitful.Length, λ::Angle) = upreferred(ω^2 * r * cos(λ))

"""
    centrifugal_force(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length)

Calculate the centrifugal force acting on an object of mass `m` rotating about a point at a distance `r` with angular velocity `ω`.
"""
centrifugal_force(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length) = upreferred(m * centrifugal_acceleration(ω, r))

"""
    angular_velocity(m::Unitful.Mass, cf::Unitful.Force, r::Unitful.Length)

Calculate the angular velocity required to develop a centrifugal force of `cf` on a body of mass `m` at a radius of `r`.
"""
angular_velocity(m::Unitful.Mass, cf::Unitful.Force, r::Unitful.Length) = angular_velocity(cf / m, r)

"""
    centrifugal_force(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length, λ::Angle)

Calculate the radial part of the centrifugal force acting of object of mass `m` at latitude `λ` on a rotating sphere of radius `r` with angular velocity `ω`.
"""
centrifugal_force(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length, λ::Angle) = upreferred(m * centrifugal_acceleration(ω, r, λ))

"""
    angular_momentum(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length)

Calculate the angular momentum of an object of mass `m` rotating about a point at a distance `r` with angular velocity `ω`.
"""
angular_momentum(m::Unitful.Mass, ω::AngularVelocity, r::Unitful.Length) = upreferred(m * r^2 * ω)

"""
	angular_velocity(m::Unitful.Mass, L::AngularMomentum, r::Unitful.Length)

Calculate the angular velocity of an object with mass `m` and angular momentum `L` rotating about a point at a distance of `r`.
"""
angular_velocity(m::Unitful.Mass, L::AngularMomentum, r::Unitful.Length) = L / (m * r^2) |> u"rad/s"

end

max_stars_repo_name: icantbelieveitsnotphysics/SfPhysics

max_issues_repo_name: icantbelieveitsnotphysics/SfPhysics

max_forks_repo_name: icantbelieveitsnotphysics/SfPhysics


text:
## SuiteSparseQR

## ordering options
const SPQR_ORDERING_FIXED   = int32(0)
const SPQR_ORDERING_NATURAL = int32(1)
const SPQR_ORDERING_COLAMD  = int32(2)
const SPQR_ORDERING_GIVEN   = int32(3) # only used for C/C++ interface
const SPQR_ORDERING_CHOLMOD = int32(4) # CHOLMOD best-effort (COLAMD, METIS,...)
const SPQR_ORDERING_AMD     = int32(5) # AMD(A'*A)
const SPQR_ORDERING_METIS   = int32(6) # metis(A'*A)
const SPQR_ORDERING_DEFAULT = int32(7) # SuiteSparseQR default ordering
const SPQR_ORDERING_BEST    = int32(8) # try COLAMD, AMD, and METIS; pick best
const SPQR_ORDERING_BESTAMD = int32(9) # try COLAMD and AMD; pick best

# Let [m n] = size of the matrix after pruning singletons.  The default
# ordering strategy is to use COLAMD if m <= 2*n.  Otherwise, AMD(A'A) is
# tried.  If there is a high fill-in with AMD then try METIS(A'A) and take
# the best of AMD and METIS.  METIS is not tried if it isn't installed.

## Operations in qmult
const SPQR_QTX = int32(0)        # Y = Q'*X
const SPQR_QX  = int32(1)        # Y = Q*X
const SPQR_XQT = int32(2)        # Y = X*Q'
const SPQR_XQ  = int32(3)        # Y = X*Q

## Types of systems to solve
const SPQR_RX_EQUALS_B    = int32(0)    # solve R*X=B      or X = R\B
const SPQR_RETX_EQUALS_B  = int32(1)    # solve R*E'*X=B   or X = E*(R\B)
const SPQR_RTX_EQUALS_B   = int32(2)    # solve R'*X=B     or X = R'\B
const SPQR_RTX_EQUALS_ETB = int32(3)    # solve R'*X=E'*B  or X = R'\(E'*B)

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
using Documenter, Actors

makedocs(
    modules = [Actors],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        collapselevel=1,
        ),
    sitename = "Actors.jl",
    authors  = "Paul Bayer",
    pages = [
        "Introduction" => "index.md",
        "Tutorials" => [
            "tutorial/install.md",
            "tutorial/intro.md",
            "tutorial/tabletennis.md",
            "tutorial/stack.md",
            "tutorial/dictsrv.md",
            "tutorial/supervise.md",
        ],
        "How to" => [
            "`spawn` actors"              => "howto/spawn.md",
            "communicate with actors"     => "howto/communicate.md",
            "get information from actors" => "howto/information.md",
            "(not) share variables"       => "howto/share.md",
            "deal with failures"          => "howto/failure.md",
            "`register` actors"           => "howto/register.md",
        ],
        "On Actors" => [
            "manual/basics.md",
            "manual/actors.md",
            "manual/behaviors.md",
            "manual/protocol.md",
            "Error Handling" => [
                "manual/errors.md",
                "manual/connections.md",
                "manual/monitors.md",
                "manual/supervisors.md",
                "manual/node_failures.md",
                "manual/checkpoints.md",
                "manual/fault_tolerance.md"
                ],
            "manual/infrastructure.md",
            ],
        "Reference" => [
            "API" => [
                "api/api.md",
                "api/types.md",
                "api/starting.md",
                "api/primitives.md",
                "api/comm.md",
                "api/user_api.md",
                "api/registry.md",
                "api/connect.md",
                "api/monitor.md",
                "api/supervision.md",
                "api/checkpointing.md",
                "api/utils.md",
                "api/diagnosis.md",
                ],
            "Examples" => [
                "examples/dining_phil.md",
                "examples/prod_cons.md",
                "examples/examples.md"
                ],
            "Internals" => [
                "reference/diag.md",
                "reference/internals.md",
                "reference/messages.md",
                "reference/interface.md",
            ],
            "reference/glossary.md",
        ]
    ]
)

deploydocs(
    repo   = "github.com/JuliaActors/Actors.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing,
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#", "dev" => "dev"]
)


text:
"""
Chaudhuri, A., Kakde, D., Jahja, M., Xiao, W., Kong, S., Jiang, H. and Percdriy, S., 2018, January.
Sampling method for fast training of support vector data description.
In 2018 Annual Reliability and Maintainability Symposium (RAMS) (pp. 1-7). IEEE.

Python implementation available here: https://github.com/samplesvdd/sample_svdd/blob/master/sample_svdd.py

We dub the method random subsample SVDD (RSSVDD)

Parameter guideline from paper and code:
subsample_size = 10 (benchmark data ∈ [3, 20])
max_iter = 1000
convergence_stop_tol = 10^-6
convergence_num_iter = 30
"""
mutable struct RSSVDD <: Sampler
    subsample_size::Int
    max_iter::Int
    convergence_stop_tol::Real
    convergence_num_iter::Int
    init_strategy::SVDD.InitializationStrategy
    solver::OptimizerFactory
    function RSSVDD(subsample_size::Int, max_iter::Int, convergence_stop_tol::Real, convergence_num_iter::Int,
                    init_strategy::SVDD.InitializationStrategy, solver::OptimizerFactory)
        subsample_size > 0 || throw(ArgumentError("Invalid size of sub sample $(subsample_size)."))
        max_iter > 0 || throw(ArgumentError("Invalid number of iterations $(max_iter)."))
        convergence_stop_tol > 0 || throw(ArgumentError("Invalid connvergence tolerance $(convergence_stop_tol)."))
        convergence_num_iter > 0 || throw(ArgumentError("Invalid number of iterations for convergence $(convergence_num_iter)."))
        return new(subsample_size, max_iter, convergence_stop_tol, convergence_num_iter, init_strategy, solver)
    end
end

function RSSVDD(subsample_size::Int, max_iter::Int,
                init_strategy::SVDD.InitializationStrategy, solver::OptimizerFactory,
                convergence_stop_tol::Real=10^-6, convergence_num_iter::Int=30)
    return RSSVDD(subsample_size, max_iter, convergence_stop_tol, convergence_num_iter, init_strategy, solver)
end

function sample(sampler::RSSVDD, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    n = size(data, 2)
    sampler.subsample_size = min(n, sampler.subsample_size)

    model = SVDD.VanillaSVDD(data)
    SVDD.initialize!(model, sampler.init_strategy)
    SVDD.set_adjust_K!(model, true)

    master_svs = []
    R = nothing
    num_stable_iter = 0
    for i in 1:sampler.max_iter
        try
            subsample_indices = randperm(n)[1:sampler.subsample_size]
            subsample_svs = train_and_find_support_vectors(model, data[:, subsample_indices], sampler.solver)
            master_svs = unique(vcat(master_svs, subsample_indices[subsample_svs]))
            new_master_svs = train_and_find_support_vectors(model, data[:, master_svs], sampler.solver)
            master_svs = master_svs[new_master_svs]
        catch e
            # Model fitting on subsample failed. Skipping subsample.
            continue
        end
        if R !== nothing && abs(model.R - R) <= sampler.convergence_stop_tol * abs(R)
            num_stable_iter += 1
        else
            num_stable_iter = 0
        end
        if num_stable_iter >= sampler.convergence_num_iter
            break
        end
        R = model.R
    end
    sample_mask = falses(n)
    sample_mask[master_svs] .= true
    return sample_mask
end


text:
export printTimes,
       reset!

mutable struct Timers
  contract_t::Float64
  contract_c::Int
  gemm_t::Float64
  gemm_c::Int
  permute_t::Float64
  permute_c::Int
  svd_t::Float64
  svd_c::Int
  svd_store_t::Float64
  svd_store_c::Int
  svd_store_svd_t::Float64
  svd_store_svd_c::Int
end

global_timer = Timers()

function reset!(t::Timers)
  t.contract_t = 0.0; t.contract_c = 0;
  t.gemm_t = 0.0; t.gemm_c = 0;
  t.permute_t = 0.0; t.permute_c = 0;
  t.svd_t = 0.0; t.svd_c = 0;
  t.svd_store_t = 0.0; t.svd_store_c = 0;
  t.svd_store_svd_t = 0.0; t.svd_store_svd_c = 0;
end

function printTimes(t::Timers)
  @printf "contract_t = %.6f (%d), Total = %.4f \n" (t.contract_t/t.contract_c) t.contract_c t.contract_t
  @printf "  gemm_t = %.6f (%d), Total = %.4f \n" (t.gemm_t/t.gemm_c) t.gemm_c t.gemm_t
  @printf "  permute_t = %.6f (%d), Total = %.4f \n" (t.permute_t/t.permute_c) t.permute_c t.permute_t
  @printf "svd_t = %.6f (%d), Total = %.4f \n" (t.svd_t/t.svd_c) t.svd_c t.svd_t
  @printf "  svd_store_t = %.6f (%d), Total = %.4f \n" (t.svd_store_t/t.svd_store_c) t.svd_store_c t.svd_store_t
  @printf "  svd_store_svd_t = %.6f (%d), Total = %.4f \n" (t.svd_store_svd_t/t.svd_store_svd_c) t.svd_store_svd_c t.svd_store_svd_t
end


text:
using JLD2;
using Distributed;
using Printf;
using IterTools;
using Distributions
@everywhere include("runit.jl");
@everywhere include("../utilities/thresholds.jl");
include("../utilities/experiment_helpers.jl");

function f(x)
    return 9cos(x)/(x^2+10)
end

dist = Gaussian();
L = 0.9; # Lipschitz constant
arms = [1.25,1.5,
        1.75,2,2.25,2.5,2.75,
        3,3.25,3.5,3.75,4,
        4.25,4.5,4.75,5,5.25,
        5.5,5.75,6];
μ = [f(k) for k in arms];
pep = LipschitzBestArm(dist, arms, μ, L);
# methods to be compared
srs = [FWSampling(), Menard(CTracking, 1/oracle(pep, μ)[1]), TrackAndStop(DTracking)];

δs =  (0.1,0.01); # confidence
βs = GK16.(δs);
repeats = 100;
seed = 1234;

println("arms=$arms, μ=$μ, repeats=$repeats");
# compute
@time data = pmap(
    ((sr, i),) -> runit(seed + i, sr, μ, pep, βs),
    Iterators.product(srs, 1:repeats),
);
dump_stats(pep, μ, δs, βs, srs, data, repeats);
# save
@save isempty(ARGS) ? "BAI1.dat" : ARGS[1] dist μ pep srs data δs βs repeats seed
# visualise by loading viz_bai.jl


text:
const test_int = 10
const test_sym = AsapSim.DMEM
const test_loc_array = [Loc(AsapSim.DMEM, 10)]

@testset "Testing Random Misc. Stuff" begin
    quot = Meta.quot

    # Testing my convenience macro.
    kwargs = []
    a = 1; b = 2; c = 3;
    AsapSim.@pairpush! kwargs a b c
    @test kwargs == [:a => a, :b => b, :c => c]

    # --- Extract Loc --- #

    extract_loc = AsapSim.extract_loc

    # Passing an integer directly to extract_loc should result in a immediate
    # style Loc.
    @test eval(extract_loc(1)) == Loc(AsapSim.IMMEDIATE, 1)

    # If a symbol is passed to extract_loc that is not a reserved symbol, we
    # should get the expression :(Loc(test_int)). Since test_int is a Int,
    # when that expression is Evaled, we get an immediate Loc with the value
    # of test_int.
    @test eval(extract_loc(:test_int)) == Loc(AsapSim.IMMEDIATE, 10)

    # If instead the symbol given to extract_loc resolved to a Symbol, the
    # constructor Loc(x::Symbol) will be called when evaled.
    @test eval(extract_loc(:test_sym)) == Loc(AsapSim.DMEM, -1)

    # If we pass a symbol that is reserved, make sure we get the corresponding
    # Loc correctly.
    @test eval(extract_loc(:ibuf0)) == Loc(AsapSim.IBUF, 0)

    # If we pass "extract_loc" an expression that resolves to an integer, it
    # should keep that expression intact so that when it is evaluated, the
    # correct computation takes place.
    @test eval(extract_loc(:(1 + 1 + 1))) == Loc(AsapSim.IMMEDIATE, 3)

    # If a valid Source / Destination is given with a literal index, make sure
    # that this is discovered and a valid Src/Dest enum is created.
    @test eval(extract_loc(:(dmem[0]))) == Loc(AsapSim.DMEM, 0)

    # If an expression for an index is given where the item being index is
    # NOT a valid Src/Dest pair, it will be passed straight through to the
    # constructor of Loc. In this case, test_loc_array[1] IS a Loc, so we should
    # just get the same thing back.
    @test eval(extract_loc(:(test_loc_array[1]))) == test_loc_array[1]

    # Finally, test that if we pass a valid reference without a literal index
    # that this index is correctly extracted when the expression is Evaled.
    @test eval(extract_loc(:(dmem[test_int]))) == Loc(AsapSim.DMEM, test_int)

    # --------------------------------------------------- #
    # --- Test some of the simpler macros for correctness #
    # --------------------------------------------------- #

    expr = MacroTools.@q function test(branch)
        if branch
            return 10
        end
    end

    # Add a program vector initialization and return statement.
    split_expr = splitdef(expr)
    function_body = split_expr[:body]
    program_vector_sym = gensym("program")
    # Add an initializer for an empty vector of InstructionLabelTarget
    AsapSim.startvector!(function_body, program_vector_sym)

    # Replace returns. The actual invocation of this function needs to be a call
    # to "assemble" - but that's kind of hard to test. Rather, just use the
    # "identity" function to get the empty vector directly.
    function_body = AsapSim.replace_returns(
        function_body,
        program_vector_sym;
        return_fn = identity
    )
    # Remake the function and "eval" it.
    split_expr[:body] = function_body
    eval(MacroTools.combinedef(split_expr))

    # The function "test" should now be available.
    # If we call it (regardless of what we give for "branch", we should recieve
    # an empty vector of InstructionLabelTarget.
    @test test(false) == AsapSim.InstructionLabelTarget[]
    @test test(true) == AsapSim.InstructionLabelTarget[]
end

# TODO: getkwargs
# TODO: replace_instructions

macro evalexpr(ex)
    expr = eval(ex)
    return esc(:($expr))
end

@testset "Testing \"replace_instructions\"" begin
    # First, try a simple move instructions. Do this by evaling the
    # returned expression because comparing expressions is REALLY hard.
    let program = AsapSim.InstructionLabelTarget[]
        @evalexpr AsapSim.replace_instructions(:(MOVE(dmem[1], 1)), :program)

        @test length(program) == 1
        # Get the instruction type.

        inst = first(program).instruction
        expected_instruction = AsapInstructionKeyword(
            op = AsapSim.MOVE,
            dest = Loc(AsapSim.DMEM, 1),
            src1 = Loc(AsapSim.IMMEDIATE, 1),
            optype = AsapSim.ALU_TYPE,
        )

        @test inst == expected_instruction
        @test first(program).label == nothing
        @test first(program).branch_target == nothing
    end

    # Test a little more complicated program.
    let program = AsapSim.InstructionLabelTarget[]
        @evalexpr AsapSim.replace_instructions(quote
            @label test
                ADDCS(output[0], ibuf0, dmem[10], cxs0)
                BR(test, j, Z, neg)

            # Make sure stuff like this is left alone.
            x = 10

            @label another
                BRL(back, nop3)
        end, :program)

        # Make sure the assignment to "x" is seen.
        @test x == 10

        # Should be 3 instructions.
        @test length(program) == 3

        # Manually check the label and branch target vectors.
        branch_targets = [i.branch_target for i in program]
        labels = [i.label for i in program]

        # Only branch instructions should have branch targets.
        @test branch_targets == [nothing, :test, :back]
        @test labels == [:test, :test, :another]

        # Test that the instructions were decoded correctly.
        expected_add = AsapInstructionKeyword(
            op = AsapSim.ADDCS,
            dest = Loc(AsapSim.OUTPUT, 0),
            src1 = Loc(AsapSim.IBUF, 0),
            src2 = Loc(AsapSim.DMEM, 10),
            cxflag = AsapSim.CX_SET,
            cxindex = 0,
            signed = true,
            saturate = true,
            optype = AsapSim.ALU_TYPE,
        )
        @test program[1].instruction == expected_add

        expected_br = AsapInstructionKeyword(
            op = AsapSim.BR,
            # Bits 1 and 12 should be set
            src1 = Loc(AsapSim.IMMEDIATE, (1 << 1) | (1 << 12)),
            jump = true,
            optype = AsapSim.BRANCH_TYPE,
        )
        @test program[2].instruction == expected_br

        expected_brl = AsapInstructionKeyword(
            op = AsapSim.BRL,
            # Unconditional, so bit 11 should be set.
            src1 = Loc(AsapSim.IMMEDIATE, 1 << 11),
            nops = 3,
            isreturn = true,
            optype = AsapSim.BRANCH_TYPE,
        )
        @test program[3].instruction == expected_brl
    end
end

@testset "Testing Assembler" begin
    # Note: This program is taken from the aes137_keysub.cpp file.
    returneight() = 8
    instructions = InstructionLabelTarget[]

    # Push a bunch of instructions to the program.
    @evalexpr AsapSim.replace_instructions(quote
        # Declare a local variable. Should not be interpreted as an instruction
        localvar = 255

    	MOVI( dmem[130], localvar, SW)          # I1

        @label main
    	RPT(3)                                  # I2
            MOVE(output[0],  ibuf[0])           # I3
        END_RPT()

        BR(GOBACK)                              # I4

        @label GOBACK
    	RPT(4)                                  # I5
            MOVE(output[0],  ibuf_next[0])      # I6
            AND(pointer[0], ibuf[0], dmem[131]) # I7
        END_RPT()

        BRL(main)                               # I8

        BRL(back)                               # I9
    end, :instructions)

    @test length(instructions) == 11

    # Get out the constituent parts of the InstructionLabelTarget vector
    program = [i.instruction for i in instructions]
    labels = [i.label for i in instructions]
    branch_targets = [i.branch_target for i in instructions]

    # Make sure instructions got decorded correctly.
    @test program[1].op == AsapSim.MOVI
    @test program[2].op == AsapSim.RPT
    @test program[3].op == AsapSim.MOVE
    @test program[4].op == AsapSim.END_RPT
    @test program[5].op == AsapSim.BR
    @test program[6].op == AsapSim.RPT
    @test program[7].op == AsapSim.MOVE
    @test program[8].op == AsapSim.AND
    @test program[9].op == AsapSim.END_RPT
    @test program[10].op == AsapSim.BRL
    @test program[11].op == AsapSim.BRL

    # Make sure the label dictionary is created correctly.
    labeldict = AsapSim.makelabeldict(labels)
    @test labeldict[:main] == 2
    @test labeldict[:GOBACK] == 6

    # Test that handlerpt! does what it's supposed to.
    AsapSim.handle_rpt!(instructions)

    program = [i.instruction for i in instructions]
    @test program[1].op == AsapSim.MOVI
    @test program[2].op == AsapSim.RPT
    @test program[3].op == AsapSim.MOVE
    @test program[4].op == AsapSim.BR
    @test program[5].op == AsapSim.RPT
    @test program[6].op == AsapSim.MOVE
    @test program[7].op == AsapSim.AND
    @test program[8].op == AsapSim.BRL
    @test program[9].op == AsapSim.BRL

    # Make sure the first RPT starts and ends at instruction3
    @test AsapSim.repeat_start(program[2]) == 3
    @test AsapSim.repeat_end(program[2]) == 3

    # Make sure the second repeat starts at instruction 6 and ends at 7
    @test AsapSim.repeat_start(program[5]) == 6
    @test AsapSim.repeat_end(program[5]) == 7

    # Feed this program into the "assemble" function. CHeck that branch
    # targets are resolved correctly.
    #
    # NOTE: handle_rpt! should be save to call multiple times, so we don't
    # really have to worry about that messing things up.
    asap_program = AsapSim.assemble(instructions)

    # Check that branches resolved correctly.
    @test AsapSim.branch_target(asap_program.instructions[4]) == 5
    @test AsapSim.branch_target(asap_program.instructions[8]) == 2
    # This is the return - it's branch target should not have been modified.
    @test AsapSim.branch_target(asap_program.instructions[9]) == -1
end


text:
function get_trigger_limits(config_paths)
    triggers = Dict()
    triggers["BUY"] = Dict()
    triggers["SELL"] = Dict()
    for (symbol, paths) in config_paths
        !haskey(triggers["BUY"],  symbol) ? triggers["BUY"][symbol]  = Dict() : ""
        !haskey(triggers["SELL"], symbol) ? triggers["SELL"][symbol] = Dict() : ""
        for bot_path in paths
            open_orders = JSON.parsefile(joinpath(bot_path, "open_orders.json"))
            for (k,v) in open_orders["BUY"]
                if !haskey(triggers["BUY"][symbol], bot_path)
                    if haskey(v, "price")
                        triggers["BUY"][symbol][bot_path] = parse(Float64, v["price"])
                    end
                else
                    if parse(Float64, v["price"]) < triggers["BUY"][symbol][bot_path]
                        if haskey(v, "price")
                            triggers["BUY"][symbol][bot_path] = parse(Float64, v["price"])
                        end
                    end
                end
            end
            for (k,v) in open_orders["SELL"]
                if !haskey(triggers["SELL"][symbol], bot_path)
                    if haskey(v, "price")
                        triggers["SELL"][symbol][bot_path] = parse(Float64, v["price"])
                    end
                else
                    if parse(Float64, v["price"]) < triggers["SELL"][symbol][bot_path]
                        if haskey(v, "price")
                            triggers["BUY"][symbol][bot_path] = parse(Float64, v["price"])
                        end
                    end
                end
            end
        end
    end
    return triggers
end
function place_triggers(config_paths, symbol, price)
    if typeof(price) == String
        price = parse(Float64, price)
    end
    triggers = get_trigger_limits(config_paths)
    if haskey(triggers["SELL"], symbol)
        for (k,v) in triggers["SELL"][symbol]
            if v * 0.9985 <= price
                if !isfile(joinpath(k, "trigger.json"))
                    open(joinpath(k, "trigger.json"), "w") do f
                        JSON.print(f, Dict(), 4)
                    end
                end
            end
        end
    end
    if haskey(triggers["BUY"], symbol)
        for (k,v) in triggers["BUY"][symbol]
            if v * 1.0015 >= price
                if !isfile(joinpath(k, "trigger.json"))
                    open(joinpath(k, "trigger.json"), "w") do f
                        JSON.print(f, Dict(), 4)
                    end
                end
            end
        end
    end
end


text:
function SQ._collect{R<:Relation}(_src::Tuple{R}, q::SQ.Node{SQ.Filter})
    src = first(_src)
    # @time res = similar(src)
    res = similar(src)
    d = first(SQ.dos(q))
    # @time grow!(res, src, d)
    grow!(res, src, d)
    return res
end

# function grow!{F}(res, src, f::F)::Void
function grow!{F}(res, src, d::SQ.Filter{F})::Void
    f = first(SQ.parts(d))::F
    for i in src
        v = f(i)
        if ifelse(isnull(v), false, unsafe_get(v))
            push!(res, i)
        end
    end
end


text:
# TESTS FOR ALL 2D UNITCELLS
@testset "2D Unitcells" begin

    #----------------------
    #   SQUARE
    #----------------------
    @testset "Square" begin

        # test Default Wrapper
        @test_nowarn getUnitcellSquare();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellSquare()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,2}, LB,B<:AbstractBond{LB,2}}

        # test Default Wrapper getUnitcell(...)
        id = :square
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellSquare(String);
        @test_nowarn getUnitcellSquare(Symbol);
        @test_nowarn getUnitcellSquare(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquare(1);
        @test_nowarn getUnitcellSquare(String, 1);
        @test_nowarn getUnitcellSquare(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquare(3);
        @test_nowarn getUnitcellSquare(String, 3);
        @test_nowarn getUnitcellSquare(String, Symbol, 3);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquare(-3);
        @test_nowarn getUnitcellSquare(String, -3);
        @test_nowarn getUnitcellSquare(String, Symbol, -3);

    end

    #----------------------
    #   TRIANGULAR
    #----------------------
    @testset "Triangular" begin

        # test Default Wrapper
        @test_nowarn getUnitcellTriangular();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellTriangular()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,2}, LB,B<:AbstractBond{LB,2}}

        # test Default Wrapper getUnitcell(...)
        id = :triangular
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellTriangular(String);
        @test_nowarn getUnitcellTriangular(Symbol);
        @test_nowarn getUnitcellTriangular(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellTriangular(1);
        @test_nowarn getUnitcellTriangular(String, 1);
        @test_nowarn getUnitcellTriangular(String, Symbol, 1);

    end

    #----------------------
    #   HONEYCOMB
    #----------------------
    @testset "Honeycomb" begin

        # test Default Wrapper
        @test_nowarn getUnitcellHoneycomb();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellHoneycomb()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,2}, LB,B<:AbstractBond{LB,2}}

        # test Default Wrapper getUnitcell(...)
        id = :honeycomb
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellHoneycomb(String);
        @test_nowarn getUnitcellHoneycomb(Symbol);
        @test_nowarn getUnitcellHoneycomb(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellHoneycomb(1);
        @test_nowarn getUnitcellHoneycomb(String, 1);
        @test_nowarn getUnitcellHoneycomb(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellHoneycomb(4);
        @test_nowarn getUnitcellHoneycomb(String, 4);
        @test_nowarn getUnitcellHoneycomb(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcellHoneycomb(4))))) == 3

    end

    #----------------------
    #   KAGOME
    #----------------------
    @testset "Kagome" begin

        # test Default Wrapper
        @test_nowarn getUnitcellKagome();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellKagome()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,2}, LB,B<:AbstractBond{LB,2}}

        # test Default Wrapper getUnitcell(...)
        id = :kagome
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellKagome(String);
        @test_nowarn getUnitcellKagome(Symbol);
        @test_nowarn getUnitcellKagome(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellKagome(1);
        @test_nowarn getUnitcellKagome(String, 1);
        @test_nowarn getUnitcellKagome(String, Symbol, 1);

    end

    #----------------------
    #   SQUARE-OCTAGON
    #----------------------
    @testset "Square-Octagon" begin

        # test Default Wrapper
        @test_nowarn getUnitcellSquareOctagon();

        # Test the unitcell has the correct subtype
        @test typeof(getUnitcellSquareOctagon()) <: AbstractUnitcell{S,B} where {LS,S<:AbstractSite{LS,2}, LB,B<:AbstractBond{LB,2}}

        # test Default Wrapper getUnitcell(...)
        id = :squareoctagon
        @test_nowarn getUnitcell(id);
        @test_nowarn getUnitcell(id, 1);
        @test_nowarn getUnitcell(id, Symbol);
        @test_nowarn getUnitcell(id, Symbol, 1);
        @test_nowarn getUnitcell(id, Int64, Symbol);
        @test_nowarn getUnitcell(id, Int64, Symbol, 1);

        # Test Wrapper with different label types
        @test_nowarn getUnitcellSquareOctagon(String);
        @test_nowarn getUnitcellSquareOctagon(Symbol);
        @test_nowarn getUnitcellSquareOctagon(Int64, Symbol);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquareOctagon(1);
        @test_nowarn getUnitcellSquareOctagon(String, 1);
        @test_nowarn getUnitcellSquareOctagon(String, Symbol, 1);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquareOctagon(2);
        @test_nowarn getUnitcellSquareOctagon(String, 2);
        @test_nowarn getUnitcellSquareOctagon(String, Symbol, 2);

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquareOctagon(4);
        @test_nowarn getUnitcellSquareOctagon(String, 4);
        @test_nowarn getUnitcellSquareOctagon(String, Symbol, 4);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcellSquareOctagon(4))))) == 3

        # Test Wrapper with different label types and specific version
        @test_nowarn getUnitcellSquareOctagon(5);
        @test_nowarn getUnitcellSquareOctagon(String, 5);
        @test_nowarn getUnitcellSquareOctagon(String, Symbol, 5);

        # test that there are enough Kitaev bonds
        @test length(unique(label.(bonds(getUnitcellSquareOctagon(5))))) == 3

    end

end


text:
const __bodyfunction__ = Dict{Method,Any}()

# Find keyword "body functions" (the function that contains the body
# as written by the developer, called after all missing keyword-arguments
# have been assigned values), in a manner that doesn't depend on
# gensymmed names.
# `mnokw` is the method that gets called when you invoke it without
# supplying any keywords.
function __lookup_kwbody__(mnokw::Method)
    function getsym(arg)
        isa(arg, Symbol) && return arg
        @assert isa(arg, GlobalRef)
        return arg.name
    end

    f = get(__bodyfunction__, mnokw, nothing)
    if f === nothing
        fmod = mnokw.module
        # The lowered code for `mnokw` should look like
        #   %1 = mkw(kwvalues..., #self#, args...)
        #        return %1
        # where `mkw` is the name of the "active" keyword body-function.
        ast = Base.uncompressed_ast(mnokw)
        if isa(ast, Core.CodeInfo) && length(ast.code) >= 2
            callexpr = ast.code[end-1]
            if isa(callexpr, Expr) && callexpr.head == :call
                fsym = callexpr.args[1]
                if isa(fsym, Symbol)
                    f = getfield(fmod, fsym)
                elseif isa(fsym, GlobalRef)
                    if fsym.mod === Core && fsym.name === :_apply
                        f = getfield(mnokw.module, getsym(callexpr.args[2]))
                    elseif fsym.mod === Core && fsym.name === :_apply_iterate
                        f = getfield(mnokw.module, getsym(callexpr.args[3]))
                    else
                        f = getfield(fsym.mod, fsym.name)
                    end
                else
                    f = missing
                end
            else
                f = missing
            end
        else
            f = missing
        end
        __bodyfunction__[mnokw] = f
    end
    return f
end

function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    Base.precompile(Tuple{Core.kwftype(typeof(fd2html)),NamedTuple{(:dir,), Tuple{String}},typeof(fd2html),String})   # time: 0.61810833
    Base.precompile(Tuple{Core.kwftype(typeof(serve)),NamedTuple{(:single,), Tuple{Bool}},typeof(serve)})   # time: 0.37084788
    Base.precompile(Tuple{typeof(lx_show),LxCom,Vector{LxDef}})   # time: 0.32274863
    Base.precompile(Tuple{typeof(fd2html),String})   # time: 0.23594038
    Base.precompile(Tuple{Core.kwftype(typeof(convert_md)),NamedTuple{(:isinternal,), Tuple{Bool}},typeof(convert_md),String})   # time: 0.22428645
    Base.precompile(Tuple{Core.kwftype(typeof(serve)),NamedTuple{(:clear, :single, :cleanup, :nomess), NTuple{4, Bool}},typeof(serve)})   # time: 0.14878392
    let fbody = try __lookup_kwbody__(which(optimize, ())) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Bool,Bool,Bool,String,Bool,Function,Bool,Bool,Bool,typeof(optimize),))
        end
    end   # time: 0.07003818
    Base.precompile(Tuple{Core.kwftype(typeof(optimize)),NamedTuple{(:prerender,), Tuple{Bool}},typeof(optimize)})   # time: 0.05545673
    Base.precompile(Tuple{typeof(find_tokens),String,LittleDict{Char, Vector{Pair{Tuple{Int64, Bool, Function, Union{Nothing, Bool, Function}}, Symbol}}, Vector{Char}, Vector{Vector{Pair{Tuple{Int64, Bool, Function, Union{Nothing, Bool, Function}}, Symbol}}}},LittleDict{Char, Symbol, Vector{Char}, Vector{Symbol}}})   # time: 0.040931374
    Base.precompile(Tuple{typeof(hfun_taglist)})   # time: 0.030687159
    Base.precompile(Tuple{typeof(add_rss_item)})   # time: 0.029599395
    Base.precompile(Tuple{typeof(hfun_toc),Vector{String}})   # time: 0.02444929
    Base.precompile(Tuple{typeof(resolve_code_block),SubString{String}})   # time: 0.01858709
    Base.precompile(Tuple{typeof(hfun_redirect),Vector{String}})   # time: 0.015353491
    Base.precompile(Tuple{typeof(lx_textoutput),LxCom,Vector{LxDef}})   # time: 0.013971993
    Base.precompile(Tuple{typeof(reprocess),String,Vector{LxDef{SubString{String}}}})   # time: 0.012774653
    Base.precompile(Tuple{typeof(lx_textinput),LxCom,Vector{LxDef}})   # time: 0.01045118
    Base.precompile(Tuple{typeof(convert_block),OCBlock,Vector{LxDef}})   # time: 0.009655179
    Base.precompile(Tuple{typeof(generate_tag_pages)})   # time: 0.009324851
    Base.precompile(Tuple{typeof(process_config)})   # time: 0.00678984
    isdefined(Franklin, Symbol("#180#186")) && Base.precompile(Tuple{getfield(Franklin, Symbol("#180#186")),SubString{String}})   # time: 0.006691009
    Base.precompile(Tuple{typeof(convert_block),Token,Vector{LxDef}})   # time: 0.006101899
    Base.precompile(Tuple{typeof(write_tag_page),String})   # time: 0.004991966
    Base.precompile(Tuple{typeof(lx_literate),LxCom,Vector{LxDef}})   # time: 0.004831332
    Base.precompile(Tuple{typeof(html_code),SubString{String},SubString{String}})   # time: 0.004530656
    Base.precompile(Tuple{typeof(hfun_insert),Vector{String}})   # time: 0.003900784
    Base.precompile(Tuple{Type{LxDef},String,Int64,Pair{String, String}})   # time: 0.003559113
    Base.precompile(Tuple{Core.kwftype(typeof(convert_md)),NamedTuple{(:isrecursive,), Tuple{Bool}},typeof(convert_md),String})   # time: 0.003385631
    Base.precompile(Tuple{typeof(hfun_paginate),Vector{String}})   # time: 0.003321519
    Base.precompile(Tuple{Core.kwftype(typeof(serve)),NamedTuple{(:clear, :single, :nomess), Tuple{Bool, Bool, Bool}},typeof(serve)})   # time: 0.002933776
    Base.precompile(Tuple{typeof(locvar),Symbol})   # time: 0.002928755
    isdefined(Franklin, Symbol("#λ#23")) && Base.precompile(Tuple{getfield(Franklin, Symbol("#λ#23")),Int64,Char})   # time: 0.002772743
    Base.precompile(Tuple{Core.kwftype(typeof(fd2html)),NamedTuple{(:dir, :internal), Tuple{String, Bool}},typeof(fd2html),String})   # time: 0.002679923
    Base.precompile(Tuple{Type{LxDef},String,Int64,SubString{String}})   # time: 0.002655788
    Base.precompile(Tuple{Core.kwftype(typeof(fd2html)),NamedTuple{(:nop,), Tuple{Bool}},typeof(fd2html),String})   # time: 0.002494217
    Base.precompile(Tuple{Core.kwftype(typeof(optimize)),NamedTuple{(:minify, :prerender, :prepath), Tuple{Bool, Bool, String}},typeof(optimize)})   # time: 0.002491155
    isdefined(Franklin, Symbol("#λ#22")) && Base.precompile(Tuple{getfield(Franklin, Symbol("#λ#22")),SubString{String},Bool})   # time: 0.002363617
    Base.precompile(Tuple{Core.kwftype(typeof(serve)),NamedTuple{(:single, :clear, :cleanup), Tuple{Bool, Bool, Bool}},typeof(serve)})   # time: 0.002335907
    let fbody = try __lookup_kwbody__(which(convert_md, (SubString{String},Vector{LxDef},))) catch missing end
        if !ismissing(fbody)
            precompile(fbody, (Base.Iterators.Pairs{Symbol, Bool, Tuple{Symbol, Symbol}, NamedTuple{(:isrecursive, :has_mddefs), Tuple{Bool, Bool}}},typeof(convert_md),SubString{String},Vector{LxDef},))
        end
    end   # time: 0.002099546
    Base.precompile(Tuple{typeof(convert_block),LxEnv,Vector{LxDef}})   # time: 0.002070447
    Base.precompile(Tuple{typeof(convert_md),String})   # time: 0.002023048
    Base.precompile(Tuple{Type{HFun},SubString{String},SubString{String},Vector{SubString{String}}})   # time: 0.001731392
    Base.precompile(Tuple{typeof(hfun_href),Vector{String}})   # time: 0.001713784
    Base.precompile(Tuple{Core.kwftype(typeof(fd2html)),NamedTuple{(:internal,), Tuple{Bool}},typeof(fd2html),String})   # time: 0.001555199
    isdefined(Franklin, Symbol("#λ#22")) && Base.precompile(Tuple{getfield(Franklin, Symbol("#λ#22")),String,Bool})   # time: 0.001472338
    Base.precompile(Tuple{Core.kwftype(typeof(convert_md)),NamedTuple{(:isrecursive, :has_mddefs, :nostripp), Tuple{Bool, Bool, Bool}},typeof(convert_md),SubString{String},Vector{LxDef}})   # time: 0.001085887
end


text:
abstract type AbstractSelectionRule end

function reset!(::AbstractSelectionRule)
    return nothing
end

"""
    SelectUniform(n::Integer)

Incoming cuts are placed into aggregates uniformly, so that each aggregate has at most `n` cuts. Behaves as [`PartialAggregation`](@ref).

"""
mutable struct SelectUniform <: AbstractSelectionRule
    n::Integer

    function SelectUniform(n::Integer)
        return new(max(1,n))
    end
end

function select(rule::SelectUniform, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    if num_subproblems(aggregates[1]) == rule.n - 1
        return 1, true
    end
    return 1, false
end

function str(::SelectUniform)
    return "uniform selection"
end

"""
    SelectDecaying(T₀::Integer, T̲::Integer = 1, γ::T)

Behaves like [`SelectUniform`](@ref), but the uniform aggregate size decays by `γ` each iteration, starting from `T₀`. `T̲` is an optional lower bound on the aggregate size.

"""
mutable struct SelectDecaying <: AbstractSelectionRule
    n::Int
    T̲::Int
    γ::Float64
    idx::Int

    function SelectDecaying(T₀::Integer, γ::Float64)
        return new(max(1,T₀), 1, γ, 1)
    end

    function SelectDecaying(T₀::Integer, T̲::Integer, γ::Float64)
        return new(max(T₀), max(1,T̲), γ, 1)
    end
end

function select(rule::SelectDecaying, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    if num_subproblems(aggregates[1]) == rule.n - 1
        return 1, true
    end
    return 1, false
end

function reset!(rule::SelectDecaying)
    rule.n = max(rule.T̲, round(Int, rule.γ*rule.n))
    rule.idx = 1
    return nothing
end

function str(::SelectDecaying)
    return "decaying aggregation level selection"
end

"""
    SelectRandom(max = Inf)

Incoming cuts are placed into aggregates randomly. An optional maximum number of cuts `max` can be specified.

"""
struct SelectRandom <: AbstractSelectionRule
    max::Float64
end

SelectRandom(; max = Inf) = SelectRandom(max)

function select(rule::SelectRandom, aggregates::Vector{<:AggregatedOptimalityCut}, ::HyperPlane{OptimalityCut})
    idx = rand(eachindex(aggregates))
    while num_subproblems(aggregates[idx]) > rule.max
        idx = rand(eachindex(aggregates))
    end
    return idx, (num_subproblems(aggregates[idx]) == rule.max-1)
end

function str(rule::SelectRandom)
    if rule.max < Inf
        return "random selection of at most $(rule.max)"
    end
    return "random selection"
end

"""
    SelectClosest(τ::AbstractFloat; distance::Function = absolute_distance)

Incoming cuts are placed into the closest aggregate, according the supplied `distance` function. An empty aggregate is chosen if no aggregate is within the tolerance `τ`

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
struct SelectClosest <: AbstractSelectionRule
    τ::Float64
    distance::Function
end

function SelectClosest(τ::AbstractFloat; distance::Function = absolute_distance)
    return SelectClosest(τ, distance)
end

function select(rule::SelectClosest, aggregates::Vector{<:AggregatedOptimalityCut}, cut::HyperPlane{OptimalityCut})
    separations = map(c -> rule.distance(c, cut), aggregates)
    (dist, idx) = findmin(separations)
    if dist > rule.τ
        zero_idx = findfirst(iszero, aggregates)
         if zero_idx == nothing
            return idx, true
        end
        return zero_idx, false
    end
    return idx, false
end

function str(::SelectClosest)
    return "distance based selection"
end

"""
    SelectClosestToReference(τ::AbstractFloat; distance::Function = absolute_distance)

Incoming cuts are placed into an aggregate based on the distance to a reference cut, according the supplied `distance` function. Behaves as [`SelectClosest`](@ref) if not withing the tolerance `τ` to the reference cut.

The following distance measures are available
- [`absolute_distance`](@ref)
- [`angular_distance`](@ref)
- [`spatioangular_distance`](@ref

"""
mutable struct SelectClosestToReference{T <: AbstractFloat} <: AbstractSelectionRule
    τ::T
    distance::Function
    reference::AggregatedOptimalityCut{T}
    buffer::AggregatedOptimalityCut{T}
end

function SelectClosestToReference(τ::AbstractFloat; distance::Function = absolute_distance)
    T = typeof(τ)
    return SelectClosestToReference{T}(τ, distance, zero(AggregatedOptimalityCut{T}), zero(AggregatedOptimalityCut{T}))
end

function select(rule::SelectClosestToReference, aggregates::Vector{<:AggregatedOptimalityCut}, cut::HyperPlane{OptimalityCut})
    # Fill buffer with total aggregate each iteration
    rule.buffer += cut
    if iszero(rule.reference)
        # Multicut before reference has been calculated
        return 1, true
    elseif rule.distance(cut, rule.reference) <= rule.τ
        # Distance to reference within tolerance
        return 1, false
    else
        for i = 2:length(aggregates)
            agg = aggregates[i]
            dist = rule.distance(cut, agg)
            if iszero(agg) || dist <= rule.τ / (num_subproblems(rule.reference) -
                                                num_subproblems(agg) + 1)
                if i == length(aggregates)
                    # Last position operates as multicut
                    return i, true
                end
                # Aggregate at i if within number adjusted tolerance
                return i, false
            end
        end
    end
end

function reset!(rule::SelectClosestToReference{T}) where T <: AbstractFloat
    # Swap in buffer and reset when all cuts have been seen
    rule.reference = zero(AggregatedOptimalityCut{T})
    rule.reference += rule.buffer
    rule.buffer = zero(AggregatedOptimalityCut{T})
    return nothing
end

function str(::SelectClosestToReference)
    return "distance to reference based selection"
end


text:
# Test if tf listner works correctly
using RobotOS
@rosimport geometry_msgs.msg: Wrench, Vector3
rostypegen(@__MODULE__)
import .geometry_msgs.msg: Wrench, Vector3

using ROSTransforms.tf
using Test

# Initialize node
init_node("listner", anonymous=true) # delete

# Declare an arbitrary transform
tf1 = Transform([1., 0, 0], [0, 0, 0, 1.])
wrench1 = Wrench(Vector3(0, 0, -1), Vector3(0, 1, 0))
wrench1_correct = Wrench(Vector3(0, 0, -1), Vector3(0, 2, 0))

# Transform that flips x/y to -x/-y
tf2 = Transform([1., 1, 0], [0, 0, 1, 0.])
wrench2 = Wrench(Vector3(0, 1, -1), Vector3(0, 0, 0))
wrench2_correct = Wrench(Vector3(0, -1, -1), Vector3(1, -1, 1))

@testset "Applying Transform" begin
    wrench_comp = transformWrench(tf1, wrench1)
    @test all([wrench_comp.force.x == wrench1_correct.force.x,
               wrench_comp.force.y == wrench1_correct.force.y,
               wrench_comp.force.z == wrench1_correct.force.z,
               wrench_comp.torque.x == wrench1_correct.torque.x,
               wrench_comp.torque.y == wrench1_correct.torque.y,
               wrench_comp.torque.z == wrench1_correct.torque.z])

    wrench_comp = transformWrench(tf2, wrench2)
    @test all([wrench_comp.force.x == wrench2_correct.force.x,
               wrench_comp.force.y == wrench2_correct.force.y,
               wrench_comp.force.z == wrench2_correct.force.z,
               wrench_comp.torque.x == wrench2_correct.torque.x,
               wrench_comp.torque.y == wrench2_correct.torque.y,
               wrench_comp.torque.z == wrench2_correct.torque.z])
end


text:
using Rankers.Storages
using Rankers.ProTeams
using LOLTools: SummonerV4, LeagueV4, http_action
using Rankers: Chart

api_key = get(ENV, "RIOT_TOKEN", "")
region = "kr"

read_from_file = !("u" in ARGS) # true # false
queue = "RANKED_SOLO_5x5"

get_ranking(::typeof(LeagueV4.challengerleagues), idx) = idx
get_ranking(::typeof(LeagueV4.grandmasterleagues), idx) = 300+idx
get_ranking(::typeof(LeagueV4.masterleagues), idx) = 1000+idx

rankers = Dict{String,Tuple{Union{Int,String},Int}}()
for league_func in (LeagueV4.challengerleagues, LeagueV4.grandmasterleagues, LeagueV4.masterleagues)
    league = league_func(api_key, region, queue; action=update_with(region, http_action, read_from_file=read_from_file, debug=true))
    for (idx, entry) in enumerate(Chart.entries(league))
        rankers[entry.summonerName] = (get_ranking(league_func, idx), entry.leaguePoints)
    end
end

chart = []
for player in players([DRX, FA], [Top, Jungle, Middle, ADC, Support])
    for summoner in summonersof(player)
        summonerName = summoner.name
        if !haskey(rankers, summonerName)
            summonerdto = SummonerV4.summoner_by_name(api_key, summoner.region, summoner.name; action=store_with(region, http_action, debug=true))
            encryptedSummonerId = summonerdto.id
            entries = LeagueV4.entries_by_summoner_id(api_key, summoner.region, encryptedSummonerId; action=update_with(summoner.region, http_action, read_from_file=read_from_file, debug=true))
            for entry in entries
                rankers[entry.summonerName] = (string(entry.tier, ' ', entry.rank), entry.leaguePoints)
                summonerName = entry.summonerName
            end
        end
        ranker = rankers[summonerName]
        push!(chart, (ranker..., summonerName))
    end
end
for entry in sort(chart, by = x -> x[2], rev=true)
    @info entry
end


text:
to_int(s) = parse(Int, s)

diff(x, y) = abs(x-y)

parse_file(filename) = to_int.(split(read(filename, String), ","))

function solve(filename, cost_fn)
    xs = parse_file(filename)
    costs = [sum([cost_fn(i, x) for x in xs]) for i in min(xs...):max(xs...)]
    min(costs...)
end

# Part 1
solve("7eg.txt", diff) == 37
solve("7.txt", diff) == 336701

cost(x,y) = sum(1:diff(x,y))

# Part 2
solve("7eg.txt", cost)
solve("7.txt", cost)


text:
"""
	path = choose_output_dir(pathinit)

Avoid overwriting (potentially important) existing results by appending to the directory name\n
Checks if the suggested output directory already exists. While yes, it appends _1, _2, etc till an unused name is found
"""
function choose_output_dir(pathinit::String)
	path = pathinit
	counter = 1
	while isdir(path)
		path = string(pathinit, "_", counter)
		counter += 1
	end
	return path
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
function rand_poisson(lambda)
	x::Int64 = 0
	t::Float64 = 0.0

	while t < lambda
		x = x + 1
		t = t + randexp()
	end

	return x - 1
end


text:
if basename(pwd()) == "aoc"
    cd("2016/25")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code