text:
"""
Sample a value For each parameter uniformly at random from the candidate vectors. Log-uniform sampling available by providing a log-spaced candidate vector.
"""
struct RandomSampler <: Sampler end

function (s::RandomSampler)(ho, iter)
    [list[rand(HO_RNG[threadid()], 1:length(list))] for list in ho.candidates]
end

# Latin Hypercube Sampler ======================================================

"""
Sample from a latin hypercube
"""
Base.@kwdef mutable struct LHSampler <: Sampler
    samples = zeros(0,0)
    iters = -1
end
function init!(s::LHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    ndims = length(ho.candidates)
    all(length(c) == ho.iterations for c in ho.candidates) ||
        throw(ArgumentError("Latin hypercube sampling requires all candidate vectors to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates))))) with $(ho.iterations) iterations"))
        s.iters == -1 && (s.iters = (1000*100*2)÷ho.iterations÷ndims)
    X, fit = LHCoptim(ho.iterations,ndims,s.iters)
    s.samples = copy(X')
end

function (s::LHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

"""
    CLHSampler(dims=[Continuous(), Categorical(2), ...])
Sample from a categorical/continuous latin hypercube. All continuous variables must have the same length of the candidate vectors.
"""
Base.@kwdef mutable struct CLHSampler <: Sampler
    samples = zeros(0,0)
    dims = []
end
function init!(s::CLHSampler, ho)
    s.samples != zeros(0,0) && return # Already initialized
    all(zip(s.dims, ho.candidates)) do (d,c)
        d isa Categorical || length(c) == ho.iterations
    end || throw(ArgumentError("Latin hypercube sampling requires all candidate vectors for Continuous variables to have the same length as the number of iterations, got lengths $(repr.(collect(zip(ho.params, length.(ho.candidates)))))"))
    ndims = length(ho.candidates)
    initialSample = randomLHC(ho.iterations,s.dims)
    X,_ = LHCoptim!(initialSample, 500, dims=s.dims)
    s.samples = copy(X')
end

function (s::CLHSampler)(ho, iter)
    [list[s.samples[dim,iter]] for (dim,list) in enumerate(ho.candidates)]
end

# GaussianProcesses sampler ====================================================

"""
Sample using Bayesian optimization. `GPSampler(Min)/GPSampler(Max)` fits a Gaussian process to the data and tries to use this model to figure out where the best point to sample next is (using expected improvement). Underneath, the package [BayesianOptimization.jl](https://github.com/jbrea/BayesianOptimization.jl/) is used. We try to provide reasonable defaults for the underlying model and optimizer and we do not provide any customization options for this sampler. If you want advanced control, use BayesianOptimization.jl directly.
"""
mutable struct GPSampler <: Sampler
    sense
    model
    # opt
    modeloptimizer
    logdims
    candidates
end

GPSampler() = error("You must specify GPSampler(Min)/GPSampler(Max) for minimization or maximization.")
GPSampler(sense) = GPSampler(sense,nothing,nothing,nothing,nothing)

function islogspace(x)
    all(x->x > 0, x) || return false
    length(x) > 4    || return false
    std(diff(log.(x))) < sqrt(eps(eltype(x)))
end

function to_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? log10.(x) : x
    end
end

function from_logspace(x,logdims)
    map(x,logdims) do x,l
        l ? exp10.(x) : x
    end
end

function init!(s::GPSampler, ho)
    # s.model === nothing || return # Already initialized
    ndims                = length(ho.candidates)
    logdims              = islogspace.(ho.candidates)
    candidates           = to_logspace(ho.candidates, logdims)
    lower_bounds         = [minimum.(candidates)...]
    upper_bounds         = [maximum.(candidates)...]
    widths               = upper_bounds - lower_bounds
    kernel_widths        = widths/10
    log_function_noise   = 0.
    hypertuning_interval = max(9, max(ho.iterations, length(ho.history)) ÷ 9)
    model = ElasticGPE(ndims, mean = MeanConst(0.),
                       kernel = SEArd(log.(kernel_widths), log_function_noise), logNoise = 0., capacity=max(ho.iterations, length(ho.history))+1)
    # set_priors!(model.mean, [GaussianProcesses.Normal(0, 100)])

    modeloptimizer = MAPGPOptimizer(every = hypertuning_interval, noisebounds = [-4, 3], # log bounds on the function noise?
                                    maxeval = 40, kernel_bounds=fill([0,0,0], length(kernel_widths))) # max iters for optimization of the GP hyperparams

    s.model = model
    s.modeloptimizer = modeloptimizer
    s.logdims = logdims
    s.candidates = candidates
end

function train_model!(s, ho)
    xv = [reshape(to_logspace(float.(ho.history[j]), s.logdims), :, 1) for j in 1:length(ho.history)]
    input = reduce(hcat, xv)
    targets = Int(s.sense)*ho.results
    BayesianOptimization.update!(s.model, input, targets)
    BayesianOptimization.optimizemodel!(s.modeloptimizer, s.model) # This determines whether to run or not internally?
end

function Base.iterate(ho::Hyperoptimizer{<:GPSampler}, iter=1)
    iter > ho.iterations && return nothing
    s = ho.sampler
    init!(s, ho)
    th = max(4, length(ho.candidates)+1)
    if length(ho.history) < th
        samples = [rand(HO_RNG[threadid()], list) for (dim,list) in enumerate(ho.candidates)]
        nt = (; Pair.((:i, ho.params...), (iter, samples...))...)
        # res = ho.objective(iter, samples...)
        # push!(ho.results, res)
        push!(ho.history, samples)
        return nt, iter+1
    else
        try
            train_model!(s, ho)
        catch ex
            @warn("BayesianOptimization failed with optimizemodel! at iter $iter: error: ", ex, maxlog=1)
        end
    end

    # We now optimize the acq function using the random sampler. This could potentially be improved upon
    @assert length(s.model.y) >= th "Hyperoptimizer does not contain enough stored samples. If you iterate the optimizer manually, do not forget to push the results to the internal storage, see https://github.com/baggepinnen/Hyperopt.jl#full-example for an example."
    acqfunc = BayesianOptimization.acquisitionfunction(ExpectedImprovement(maximum(s.model.y)), s.model)
    # plot(s) |> display
    # plot!(x->acqfunc([x]), 1, 5, sp=2)
    # sleep(0.2)
    iters = min(80, prod(length, s.candidates))
    ho2 = Hyperoptimizer(iterations=iters, params=ho.params, candidates=s.candidates, objective=acqfunc)
    for params in ho2
        params2 = collect(params)[2:end]
        res = -Inf
        try
            res = acqfunc(params2)
            if isnan(res)
                init!(s, ho); train_model!(s, ho)
                acqfunc = BayesianOptimization.acquisitionfunction(ExpectedImprovement(maximum(s.model.y)), s.model)
            end
        catch ex
            @warn("BayesianOptimization acqfunc failed at iter $iter: error: ", ex, maxlog=1)
        end
        push!(ho2.results, res)
    end

    # @show ho2.maximizer
    # xpl = reduce(vcat,ho2.history)
    # @show length(xpl)
    # scatter!(xpl, zeros(3000), sp=2) |> display

    samples = from_logspace(ho2.maximizer, s.logdims)
    push!(ho.history, samples)
    # push!(ho.results, res)
    nt = (; Pair.((:i, ho.params...), (iter, samples...))...)
    return nt, iter+1
end

# Hyperband ====================================================================

"""

"""
Base.@kwdef mutable struct Hyperband <: Sampler
    R
    η::Int = 3
    minimum = (Inf,)
    inner::Sampler = RandomSampler()
end
Hyperband(R) = Hyperband(R=R)

function hyperband_costfun(ex, params, candidates, sampler, ho_, objective)
    quote
        iters = $(esc(candidates[1]))
        if $(esc(sampler)).inner isa LHSampler
            smax = floor(Int, log($(esc(sampler)).η,$(esc(sampler)).R))
            B = (smax + 1)*$(esc(sampler)).R
            iters = floor(Int,$(esc(sampler)).R*$(esc(sampler)).η^smax)
            ss = string($(esc(sampler)).inner)
            @info "Starting Hyperband with inner sampler $(ss). Setting the number of iterations to R*η^log(η,R)=$(iters), make sure all candidate vectors have this length as well!"
        end
        costfun = $(objective)
        (::$typeof(costfun))($(esc(params[1])), $(esc(:state))) = $(esc(ex.args[2]))
        costfun, iters
    end
end

function hyperband(ho::Hyperoptimizer{Hyperband})
    hb = ho.sampler
    R, η = hb.R, hb.η
    hb.minimum = (Inf,)
    smax = floor(Int, log(η,R))
    B = (smax + 1)*R
    # ho.iterations >= R*η^smax || error("The number of iterations must be larger than R*η^log(η,R) = $(R*η^smax)")
    Juno.progress() do id
        for s in smax:-1:0
            n = ceil(Int, (B/R)*((η^s)/(s+1)))
            r = R / (η^s)
            minᵢ = successive_halving(ho, n, r, s)
            if minᵢ[1] < hb.minimum[1]
                hb.minimum = minᵢ
            end
            Base.CoreLogging.@logmsg -1 "Hyperband" progress=(smax-s)+1/(smax+1)  _id=id
        end
    end
    return hb.minimum
end

function successive_halving(ho, n, r=1, s=round(Int, log(hb.η, n)))
    hb = ho.sampler
    costfun = ho.objective
    η = hb.η
    minimum = Inf
    T = [ hb.inner(ho, i) for i=1:n ]
    # append!(ho.history, T)
    Juno.progress() do id
        for i in 0:s
            nᵢ = floor(Int,n/(η^i))
            rᵢ = floor(Int,r*(η^i))
            if i == 0
                LTend = [ costfun(rᵢ, t...) for t in T ]
            else
                LTend = [ costfun(rᵢ, t) for t in T ]
            end
            L, T = first.(LTend), last.(LTend)
            # if i == 0
            #     append!(ho.results, L)
            # else
            #     for t in eachindex(T) # Update results for those that were continued
            #         ho.results[findlast(x->x==T[t], ho.history)] = L[t]
            #     end
            # end
            append!(ho.history, T)
            append!(ho.results, L)
            if hb.inner isa BOHB
                update_observations(ho, rᵢ, T, L)
            end
                perm = sortperm(L)
            besti = perm[1]
            if L[besti] < minimum[1]
                minimum = (L[besti], rᵢ, T[besti])
            end
            T = T[perm[1:floor(Int,nᵢ/η)]]
            # T, minimum
            # T, minimum = top_k(hb,T,L,nᵢ,minimum)
            Base.CoreLogging.@logmsg -1 "successive_halving" progress=i/s  _id=id
        end
    end
    return minimum
end

# BOHB ====================================================================
# Acknowledgement: Code structure refers to official implementation of BOHB in ['HpBandSter'](https://github.com/automl/HpBandSter)
#
# Copyright of HpBandSter:
#
#     Copyright (c) 2017-2018, ML4AAD
#     All rights reserved.

# struct to record BOHB Observation
mutable struct ObservationsRecord
    dim::Int
    observation::Union{Vector, Tuple}
    loss::Real
end
function ObservationsRecord(observation, loss)
    ObservationsRecord(length(observation), observation, loss)
end

"""
    BOHB samplers
All variable names refer symbols in paper [`https://arxiv.org/pdf/1807.01774v1.pdf`]
 - `ρ`: Fraction of random samples
 - `q`: Fraction of best observations to build l and g
 - `N_s`: Sample batch number
 - `N_min`: Minimum number of points to build a model
 - `bw_factor`: Bandwidth factor
 - `D`: Evaluated observations
 - `max_valid_budget`: Maximum budget i that |D_{i}| is big enough to fit a model
 - `N_b`: |D_{max_valid_budget}|
 - `KDE_good`: KDE consists of "good observations", see BOHB paper
 - `KDE_bad`: KDE consists of "bad observations", see BOHB paper
"""
Base.@kwdef mutable struct BOHB <: Sampler
    dims::Union{Vector{DimensionType}, Nothing}=nothing
    # hyperparameters for BOHB
    N_min::Union{Int, Nothing} = nothing
    ρ::AbstractFloat = 1/3
    q::AbstractFloat = 0.15
    N_s::Int = 64
    bw_factor::Real = 3
    ## minimum bandwidth: this parameter doesn't occur in paper but used in the official implementation
    min_bandwidth::Real = 1e-3
    # Random sampler used for random sampling in BOHB algorithm
    random_sampler::RandomSampler = RandomSampler()
    # Context data
    ## Current observations, stored in a Dict, in which key is budget, value is a observation array to fit KDEs
    ## key of D: A real number represents budget
    ## value of D: An vector of ObservationsRecord, all the records of corresponding budget
    D::Dict{Real, Vector{ObservationsRecord}} = Dict{Real, Vector{ObservationsRecord}}()
    ## Current maximum budget that |D_{b}| > N_{min}+2, means it is valid for fit KDEs
    max_valid_budget::Union{Number, Nothing} = nothing
    ## |D| of max_valid_budget
    N_b::Union{Int, Nothing} = nothing
    ## Good and bad kernel density estimator
    KDE_good::Union{MultiKDE.KDEMulti, Nothing} = nothing
    KDE_bad::Union{MultiKDE.KDEMulti, Nothing} = nothing
end

# object call of BOHB sampler
function (s::BOHB)(ho, iter)
    # with probability ρ, return random sampled observations.
    # If max_valid_budget is nothing, which means currently we don't have enough sample for TPE, random sample as well.
    if rand() < s.ρ || s.max_valid_budget === nothing
        return s.random_sampler(ho, iter)
    end
    potential_samples = [sample_potential_hyperparam(s.KDE_good, s.min_bandwidth, s.bw_factor) for _ in 1:s.N_s]
    scores = [score(sample, s.KDE_good, s.KDE_bad) for sample in potential_samples]
    _, best_idx = findmax(scores)
    potential_samples[best_idx]
end

# Sample score l(x)/g(x), refers to line 6 of Algorithm2 in paper
function score(sample::Vector, KDE_good::MultiKDE.KDEMulti, KDE_bad::MultiKDE.KDEMulti)
    pdf(KDE_good, sample) / pdf(KDE_bad, sample)
end

# Update budget observations in BOHB
function update_observations(ho::Hyperoptimizer{Hyperband}, rᵢ, observations, losses)
    # history passed from hyperband is reversed, can not used for update
    observations = reverse.(observations)
    bohb = ho.sampler.inner
    if !haskey(bohb.D, rᵢ)
        bohb.D[rᵢ] = []
    end
    for (c, l) in zip(observations, losses)
        push!(bohb.D[rᵢ], ObservationsRecord(c, l))
    end
    D_length = length(bohb.D[rᵢ])
    if bohb.N_min === nothing
        bohb.N_min = length(ho.candidates)+1
    end
    if D_length > bohb.N_min+2 && (bohb.max_valid_budget===nothing || rᵢ >= bohb.max_valid_budget)
        bohb.max_valid_budget, bohb.N_b = rᵢ, D_length
        update_KDEs(ho)
    end
end

function update_KDEs(ho::Hyperoptimizer{Hyperband})
    bohb = ho.sampler.inner
    records = bohb.D[bohb.max_valid_budget]
    # fit KDEs according to Eqs. (2) and (3) in paper
    N_bl = Int(max(bohb.N_min, floor(bohb.q*bohb.N_b)))
    N_bg = max(bohb.N_min, bohb.N_b-N_bl)
    sort_idx = sortperm(records, by=d->d.loss)
    idx_N_bl = sort_idx[begin:N_bl]
    idx_N_bg = reverse(sort_idx)[N_bg:end]
    bohb.KDE_good = KDEMulti(bohb.dims, records[idx_N_bl], bohb.min_bandwidth, ho.candidates)
    bohb.KDE_bad = KDEMulti(bohb.dims, records[idx_N_bg], bohb.min_bandwidth, ho.candidates)
end

# sample from KDEMulti
function sample_potential_hyperparam(kde::MultiKDE.KDEMulti, min_bandwidth, bw_factor)
    idx = rand(1: size(kde.mat_observations)[2])
    param = [kde.observations[i][idx] for i in 1:length(kde.observations)]
    sample = Vector()
    for (_i, _param, dim_type, _kde) in zip(1:length(kde.dims), param, kde.dims, kde.KDEs)
        bw = max(_kde.bandwidth, min_bandwidth)
        local ele
        if dim_type isa MultiKDE.ContinuousDim
            bw = bw*bw_factor
            ele = rand(TruncatedNormal(_param, bw, -_param/bw, (1-_param)/bw))
        elseif dim_type isa MultiKDE.CategoricalDim
            if rand() < (1-bw)
                ele = _param
            else
                ele = rand(1:dim_type.levels)
            end
        elseif dim_type isa MultiKDE.UnorderedCategoricalDim
            if rand() < (1-bw)
                ele = _param
            else
                ele = rand(1:dim_type.levels)
            end
        else
            error(string("Dim type ", string(dim_type), " not supported. "))
        end
        if kde.mapped[_i]
            ele = kde.index_to_unordered[_kde][ele]
        end
        push!(sample, ele)
    end
    sample
end

# Constructor extensions and adapters for MultiKDE.jl
const DIMENSION_TYPE = Dict(Categorical=>MultiKDE.CategoricalDim, Continuous=>MultiKDE.ContinuousDim, UnorderedCategorical=>MultiKDE.UnorderedCategoricalDim)

function MultiKDE.KDEMulti(dim_types::Vector{DimensionType}, records::Vector{ObservationsRecord}, min_bandwidth::Real, candidates::Tuple)
    # Get KDEMulti with min_bandwidth
    dim = records[1].dim
    observations = Vector{Vector}()
    for record in records
        @assert record.dim == dim "All observations need to be same dimension. "
        _observations = record.observation
        if _observations isa Tuple
            _observations = [_obs for _obs in _observations]
        end
        push!(observations, _observations)
    end
    # bws = [max(min_bandwidth, KernelDensity.default_bandwidth(observations[i, :])) for i in 1:size(observations)[1]]
    # candidates = get_dict_candidates(dim_types, candidates)
    multi_kde = KDEMulti(dim_types, observations, candidates)
    for i in 1:length(multi_kde.KDEs)
        if multi_kde.KDEs[i].bandwidth < min_bandwidth
            multi_kde.KDEs[i].bandwidth = min_bandwidth
        end
    end
    multi_kde
end

function MultiKDE.KDEMulti(dims::Vector{DimensionType}, observations::Vector, candidates::Tuple)
    dims = Vector{MultiKDE.DimensionType}([DIMENSION_TYPE[typeof(dim)] === MultiKDE.ContinuousDim ? DIMENSION_TYPE[typeof(dim)]() :
                                            DIMENSION_TYPE[typeof(dim)](dim.levels) for dim in dims])
    MultiKDE.KDEMulti(dims, observations, candidates)
end

function Base.getproperty(dim_type::Union{MultiKDE.CategoricalDim, MultiKDE.UnorderedCategoricalDim}, v::Symbol)
    if v == :levels
        getfield(dim_type, :level)
    else
        getfield(dim_type, v)
    end
end


text:
module RMachine

using RCall
using DataFrames
using Random

using AMLPipelineBase
using AMLPipelineBase.AbsTypes
using AMLPipelineBase.Utils

import AMLPipelineBase.AbsTypes: fit!, transform!

const RLearners=["rf","svm"]

function __init__()
   #packages = ["caret","e1071","gam","randomForest",
   #            "nnet","kernlab","grid","MASS","pls"]
   rpackages = ["caret","e1071","gam","randomForest"]
   for pk in rpackages
      rcall(:library,pk,"lib=.libPaths()")
   end
end

export CRTLearner,fit!,transform!
# CARET wrapper that provides access to all learners.
#
# Options for the specific CARET learner is to be passed
# in `options[:impl_options]` dictionary.
mutable struct CRTLearner <: Learner
   name
   model
   function CRTLearner(args::Dict{Symbol,<:Any} = Dict())
      fitControl="trainControl(method = 'cv',number = 5,repeats = 5)"
      fitControl="trainControl(method = 'none')"
      default_args = Dict{Symbol,Any}(
                      :name => "rlearner",
                      :learner => "rf",
                      :fitControl => fitControl,
                      :impl_args => Dict{Symbol,Any}()
                     )
      cargs = nested_dict_merge(default_args, args)
      cargs[:name] = cargs[:name]*"_"*randstring(3)
      rl = cargs[:learner]
      if !(rl in RLearners)
         println("$rl is not supported.")
         println()
         throw(ArgumentError("Argument keyword error"))
      end
      new(cargs[:name],cargs)
   end
end

function fit!(crt::CRTLearner,x::DataFrame,y::Vector)
   rmodel = rcall(:train,x,y,method=crt.model[:learner],trControl = reval(crt.model[:fitControl]))
   crt.model[:rmodel] = rmodel
end

function transform!(crt::CRTLearner,x::DataFrame)
   res = rcall(:predict,crt.model[:rmodel],x)
   return rcopy(res) |> Array
end

end # module

ext: jl

max_stars_repo_name: ppalmes/AMLPipelineR

max_issues_repo_name: ppalmes/AMLPipelineR

max_forks_repo_name: ppalmes/AMLPipelineR


text:
"""
    single_variable_objective(com::CS.CoM, var_idx::Int, val::Int)

Return the best objective if `var_idx` is set to `val`
"""
function single_variable_objective(com::CS.CoM, var_idx::Int, val::Int)
    if com.objective.index != var_idx
        if com.sense == MOI.MIN_SENSE
            return com.search_space[com.objective.index].min
        else # MAX
            return com.search_space[com.objective.index].max
        end
    else
        return val
    end
end


text:
using Documenter
using Documenter: Selectors, Expanders
using Markdown

# Include external code inside documentation
struct LiteralInclude <: Expanders.ExpanderPipeline end

Selectors.order(::Type{LiteralInclude}) = 0.5
Selectors.matcher(::Type{LiteralInclude}, node::Markdown.Code, page, doc) = occursin("@literalinclude", node.language)
Selectors.matcher(::Type{LiteralInclude}, node, page, doc) = false

function Selectors.runner(::Type{LiteralInclude}, node, page, doc)
    if node.code != ""
        error("content is not supported in @literalinclude")
    end

    matched = match(r"@literalinclude (.*) (\d+)-(\d*)", node.language)
    if matched === nothing
        error("@literalinclude should look like '@literalinclude path/to/file.jl 3-5'")
    end

    path = matched[1]
    start = parse(Int, matched[2])
    if matched[3] == ""
        stop = nothing
    else
        stop = parse(Int, matched[3])
        if stop < start
            error("@literalinclude last line should be larger than start line")
        end
    end

    content = open(path) do fd
        lines = readlines(fd)
        if stop !== nothing
            join(lines[start:stop], '\n')
        else
            join(lines[start:end], '\n')
        end
    end
    page.mapping[node] = Markdown.Code("julia", content)
end

using Pkg
Pkg.activate(joinpath(@__DIR__, ".."))
Pkg.instantiate()
using Chemfiles

makedocs(
    sitename="Chemfiles.jl",
    modules=[Chemfiles],
    strict=true,
    checkdocs=:all,
    format=Documenter.HTML(
       prettyurls=get(ENV, "CI", nothing) == "true",
    ),
    pages=[
        "index.md",
        "tutorials.md",
        "reference/trajectory.md",
        "reference/frame.md",
        "reference/atom.md",
        "reference/cell.md",
        "reference/topology.md",
        "reference/residue.md",
        "reference/selection.md",
        "reference/misc.md",
    ]
)


text:
module TestMacrosSelect

import Test: @testset, @test
import Volcanito: @select, materialize
import DataFrames: DataFrame

@testset "@select" begin
    df = DataFrame(
        a = [1, 2, 3, 4],
        b = [2, 3, 4, 5],
    )

    plan = @select(df, c = a + b)
    res = materialize(plan)
    @test res == DataFrame(c = [3, 5, 7, 9])
end

end


text:
add_console_command(r::Regex, f::Function) =
    push!(GtkREPL.console_commands, ConsoleCommand(r, f, :normal))

add_console_command(r::Regex, f::Function, c::Symbol) =
    push!(GtkREPL.console_commands, ConsoleCommand(r, f, c))

#first try to match line number
function init_console_commands()

    add_console_command(r"^edit (.*):(\d+)", (m, c) -> begin
        try
            line = parse(Int, m.captures[2])
            file = normpath(strip(m.captures[1]))
            isempty(file) && return nothing

            file = file[1] == '~' ? remotecall_fetch(expanduser, worker(c), file) : file
            file = isabspath(file) ? file : joinpath(pwd(c), file)

            open_tab(file, _editor(c), line=line)
        catch err
            println("Invalid line number: $(m.captures[2])")
        end
        nothing
    end, :file)

    add_console_command(r"^edit (.*)", (m, c) -> begin
        file = normpath(strip(m.captures[1]))
        isempty(file) && return nothing

        file = file[1] == '~' ? remotecall_fetch(expanduser, worker(c), file) : file
        file = isabspath(file) ? file : joinpath(pwd(c), file)

        open_tab(file, _editor(c))
        nothing
    end, :file)

end


text:
# one EIF neuron with constant input

push!(LOAD_PATH, abspath(@__DIR__,".."))
using LinearAlgebra,Statistics,StatsBase
using Plots,NamedColors ; theme(:dark)
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using ProgressMeter

function onesparsemat(w::Real)
  mat=Matrix{Float64}(undef,1,1) ; mat[1,1]=w
  return sparse(mat)
end

## Time parameters
dt = 1E-4
Ttot = 30.
# EIF neuron parameters
# τ::Float64 # time constant (membrane capacitance)
# g_l::Float64 # conductance leak
# v_expt::Float64 # exp term threshold
# steep_exp::Float64 # steepness of exponential term
# v_threshold::Float64 # spiking threshold
# v_reset::Float64 # reset after spike
# v_leak::Float64 # reversal potential for leak term
# τ_refractory::Float64 # refractory time

τ = 0.5
g_l = 1.0
v_expt=10.0
steep_exp = 1.0
v_th = 20.
v_r = -5.0
τrefr = 0.0
v_leak = 0.0

myinput = 11.0
mynt = S.NTEIF(τ,g_l,v_expt,steep_exp,v_th,v_r,v_leak,τrefr)
myps = S.PSEIF(mynt,1)

## one static input above threshold
in_state = S.PSSimpleInput(S.InputSimpleOffset(myinput))
mypop = S.Population(myps,(S.FakeConnection(),),(in_state,))

## that's it, let's make the network
myntw = S.RecurrentNetwork(dt,(mypop,))

## Recorder
myrec1 = S.RecStateNow(myps,10,dt,Ttot)
myrec2 = S.RecCountSpikes(myps,dt)

## Run!
times = (0:myntw.dt:Ttot)
nt = length(times)
# initial conditions
myps.state_now[1] = v_r

S.reset!.([myrec1,myrec2])

for (k,t) in enumerate(times)
  S.dynamics_step!(t,myntw)
  myrec1(t,k,myntw)
  myrec2(t,k,myntw)
end
r0 = S.get_mean_rates(myrec2,dt,Ttot)[1]
##
plot(myrec1.times,myrec1.state_now[1,:])
##

# now, use Fokker-Planck theory and Richardson method
# it is hardcoded for now... maybe later I will convert it to
# a function
dv = 1E-3
v_lb = v_r - 0.1
psifun(V) = steep_exp * exp((V-v_expt)/steep_exp)
ndv = round(Int64,(v_th - v_lb)/dv)
k_re = round(Int64,(v_r - v_lb)/dv)
ks = collect(0:ndv)
Vks = @. v_lb + ks * dv
js = fill(NaN,ndv+1)
ps = fill(NaN,ndv+1)

ksrev = collect(ndv+1:-1:2)

js[ndv+1]=1.0
ps[ndv+1]=0.0

for k in ksrev
  krond = k==k_re ? 1.0 : 0.0
  js[k-1] = js[k] - krond
  ps[k-1] = js[k-1] * mynt.τ / (myinput-Vks[k]+ psifun(Vks[k]))
end

r0 = inv(dv*sum(ps))

##
histogram(myrec1.state_now[:];nbins=100,normalize=true)
plot!(Vks,r0.*ps;color=:white,linewidth=2)

## can I use this for the f-I curve?

function f_lif_num(h_in;dt=1E-4,Ttot=10.0)
  in_state = S.PSSimpleInput(S.InputSimpleOffset(h_in))
  mypop = S.Population(myps,(S.FakeConnection(),),(in_state,))
  myntw = S.RecurrentNetwork(dt,(mypop,))
  rec = S.RecCountSpikes(myps,dt)
  ## Run!
  times = (0:myntw.dt:Ttot)
  # initial conditions
  myps.state_now[1] = v_r
  S.reset!(rec)
  for (k,t) in enumerate(times)
    S.dynamics_step!(t,myntw)
    rec(t,k,myntw)
  end
  r0 = S.get_mean_rates(rec,dt,Ttot)[1]
  return r0
end

function f_lif_semianalitic(h_in;dv=1E-3)
  v_lb = v_r - 0.1
  psifun(V) = steep_exp * exp((V-v_expt)/steep_exp)
  ndv = round(Int64,(v_th - v_lb)/dv)
  k_re = round(Int64,(v_r - v_lb)/dv)
  ks = collect(0:ndv)
  Vks = @. v_lb + ks * dv
  js = fill(NaN,ndv+1)
  ps = fill(NaN,ndv+1)

  ksrev = collect(ndv+1:-1:2)

  js[ndv+1]=1.0
  ps[ndv+1]=0.0

  for k in ksrev
    krond = k==k_re ? 1.0 : 0.0
    js[k-1] = js[k] - krond
    ps[k-1] = js[k-1] * mynt.τ / (h_in-Vks[k]+ psifun(Vks[k]))
  end

  return inv(dv*sum(ps))
end

##
_ = let Ttot = 50.0
  nin = 50
  inputs = range(10.0,80.0;length=nin)
  y1 = similar(inputs)
  y2 = similar(inputs)
  @showprogress for i in eachindex(inputs)
    y1[i]=f_lif_num(inputs[i];Ttot=Ttot)
    y2[i]=f_lif_semianalitic(inputs[i])
  end
  plot(inputs,[y1 y2];
     label=["simulation" "Fokker-Planck"],
     linewidth=2,marker=:star,leg=:bottomright,
     xlabel="input current",
     ylabel="rate", title="exp IF model")
end

savefig("/tmp/tmp.png")


text:
function assemble_coherent_interface_condition!(
    sysmatrix,
    basis,
    interfacequads,
    stiffness,
    mesh,
    penalty,
    eta,
)

    check_eta(eta)
    dim = dimension(mesh)
    ncells = number_of_cells(mesh)
    jac = jacobian(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)
        if cellsign == 0
            quad1 = interfacequads[-1, cellid]
            quad2 = interfacequads[+1, cellid]

            normals = interface_normals(interfacequads, cellid)
            scaleareas = interface_scale_areas(interfacequads, cellid)

            top = coherent_traction_operators(
                basis,
                quad1,
                quad2,
                normals,
                stiffness[-1],
                stiffness[+1],
                dim,
                scaleareas,
                jac,
                vectosymmconverter,
                eta,
            )
            mop = coherent_mass_operators(
                basis,
                quad1,
                quad2,
                dim,
                penalty * scaleareas,
            )

            nodeids1 = nodal_connectivity(mesh, -1, cellid)
            nodeids2 = nodal_connectivity(mesh, +1, cellid)

            assemble_interface_condition!(
                sysmatrix,
                nodeids1,
                nodeids2,
                dim,
                top,
                mop,
            )
        end
    end
end

function assemble_incoherent_interface_condition!(
    sysmatrix,
    basis,
    interfacequads,
    stiffness,
    mesh,
    penalty,
    eta,
)

    check_eta(eta)
    dim = dimension(mesh)
    ncells = number_of_cells(mesh)
    jac = jacobian(mesh)
    vectosymmconverter = vector_to_symmetric_matrix_converter()

    for cellid = 1:ncells
        cellsign = cell_sign(mesh, cellid)
        check_cellsign(cellsign)
        if cellsign == 0
            quad1 = interfacequads[-1, cellid]
            quad2 = interfacequads[+1, cellid]

            normals = interface_normals(interfacequads, cellid)
            components = normals
            scaleareas = interface_scale_areas(interfacequads, cellid)

            top = component_traction_operators(
                basis,
                quad1,
                quad2,
                components,
                normals,
                stiffness[-1],
                stiffness[+1],
                dim,
                scaleareas,
                jac,
                vectosymmconverter,
                eta,
            )
            mop = component_mass_operators(
                basis,
                quad1,
                quad2,
                components,
                dim,
                penalty * scaleareas,
            )

            nodeids1 = nodal_connectivity(mesh, -1, cellid)
            nodeids2 = nodal_connectivity(mesh, +1, cellid)

            assemble_interface_condition!(
                sysmatrix,
                nodeids1,
                nodeids2,
                dim,
                top,
                mop,
            )
        end
    end
end


text:
# SimLynx/src/activity.jl
# Licensed under the MIT License. See LICENSE.md file in the project root for
# full license information.

# This is a placeholder for the future definition of an activity.

"""
    Activity

Represents an activity in a simulation.

# Fields
- `name::String`: the name of the activity (used in tracing)
- `enter::Union{Event, Nothing}`: the enter event for the activity
- `duration::Function': the function to compute the activity duration
- `enter::Union{Event, Nothing}`: the enter event for the activity
"""
struct Activity <: ActiveElement
    name::String
    enter::Union{Event, Nothing}
    duration::Function
    exit::Union{Event, Nothing}
end

"Prints the activity name to the IO stream."
Base.show(io::IO, activity::Activity) =
    print(io, "Activity $(activity.name)")

"""
    @activity <sig> begin
        <body>
    end

Define a simulation event with the specified signature and implemented by the
given body. The exact format is TBD.
"""
macro activity(sig, body)
end


text:
jl_bfr = memio()
c_bfr = memio()

abstract StructArg
  abstract CXUnsavedFile <: StructArg
  abstract CXSourceLocation <: StructArg
  abstract CXSourceRange <: StructArg
  abstract CXCursor <: StructArg
  abstract CXType <: StructArg
  abstract CXToken <: StructArg
  abstract CXString <: StructArg

abstract EnumArg
  abstract CXTypeKind <: EnumArg
  abstract CXCursorKind <: EnumArg
  abstract CXLinkageKind <: EnumArg
  abstract CXAvailabilityKind <: EnumArg
  abstract CXLanguageKind <: EnumArg

abstract PointerArg
  abstract CXTranslationUnit <: PointerArg
  abstract CXFile <: PointerArg

c_types = {
  Int32 => "int",
  Uint32 => "unsigned int",
  Int64 => "long long",
  Uint64 => "unsigned long long",
  CXUnsavedFile => "char*",
  CXSourceLocation => "char*",
  CXSourceRange => "char*",
  CXCursor => "char*",
  CXType => "char*",
  CXToken => "char*",
  CXString => "char*",
  CXFile => "char*",
  Void => "void",
  CXTypeKind => "CXTypeKind",
  CXCursorKind => "CXCursorKind",
  CXLinkageKind => "CXLinkageKind",
  CXAvailabilityKind => "CXAvailabilityKind",
  CXLanguageKind => "CXLanguageKind",
  CXTranslationUnit => "CXTranslationUnit",
  }

symbols = {}

is_struct(x) = (super(x) == StructArg ? true : false)
is_enum(x) = (super(x) == EnumArg ? true : false)

type Arg
  cname::String
  ctype::String
end

function flatten(f...)
  ret = {}
  upret(x) = push(ret, x)
  flt(k) = (typeof(k) == Array{Any,1}) ? (for x in k upret(x) end) : upret(k)
  map(flt, f)
  return ret
end

function to_c(io::IOStream, rtype, fname, args)
  c_ret = is_struct(rtype) ? "void " : c_types[rtype]
  c_args = map(x->ref(c_types,x), is_struct(rtype) ? flatten(args,{rtype}) : args )

  named_args = { Arg("a$i", a) for (i,a) in enumerate(c_args) }
  declargs = join( [x.ctype*" "*x.cname*"," for x in named_args])[1:end-1]

  println(io, "$c_ret wci_$fname($declargs) {")

  cl_args = {}
  for (i,a) in enumerate(args)
    if is_struct(a)
      println(io, "  $a l$i = wci_get_$a(a$i);")
      push(cl_args, "l$i")
    else
      push(cl_args, "a$i")
    end
  end
  cl_fargs = join( ["$a," for a in cl_args])[1:end-1]

  if is_struct(rtype)
    println(io, "  $rtype rx = clang_$fname($cl_fargs);")
    println(io, "  wci_save_$rtype(rx,a", string(length(args)+1), ");")
  else
    println(io, "  return clang_$fname($cl_fargs);")
  end
  println(io, "}")
end

function to_jl(io::IOStream, rtype, fname, args)
  push(symbols, "wci_$fname")
  arg_ids = ["a"*string(i) for i=1:length(args)]
  f_args = join([string(arg_ids[i])*"::"string(args[i])*"," for i in 1:length(args)])

  println(io, "function $fname($f_args)")
  f_rtype = rtype
  ret_data = ""
  ret_ptr = ""
  if( is_struct(rtype) )
    println(io, "  c = $rtype()")
    ret_data = "c.data,"
    ret_ptr = "Ptr{Void},"
    f_rtype = "Void"
  end
  cc_atype(x) = (
    if (is_struct(x)) return "Ptr{Void}"
    else return x end)
  cc_rtype(x,y) = (
    if (is_struct(x)) return "$y.data"
    else return y end)
  ccall_args1 = join([string(i)*"," for i in map(cc_atype, args)])
  ccall_args2 = join([i*"," for i in map(cc_rtype, args, arg_ids)])
  println(io, "  ccall( (:wci_$fname, libwci), $f_rtype, ($ccall_args1$ret_ptr), $ccall_args2 $ret_data)")
  if(rtype == CXString)
    println(io, "  return get_string(c)")
  elseif ( is_struct(rtype) )
    println(io, "  return c")
  end
  println(io, "end")
end
to_jl(rtype, fname, args) = to_jl(stdout_stream, rtype, fname, args)

function write_output()
  flush(jl_bfr)
  flush(c_bfr)
  seek(jl_bfr,0)
  seek(c_bfr,0)

  # Write Julia wrapper functions
  f_base = open("../src/cindex_base.jl", "w")

  for line in EachLine(jl_bfr)
    print(f_base, line)
  end

  # Write C wrapper functions
  f_hdr = open("../lib/wrapcindex.h", "w")
  for line in EachLine(c_bfr)
    print(f_hdr, line)
  end

  close(f_hdr)
  close(f_base)
end

macro cx(rtype, fname, args)
  r,f,a = eval(rtype),string(fname),eval(args)
  to_c(c_bfr, r, f, a)
  to_jl(jl_bfr, r,f,a)
end

# Four spaces: probably won't wrap.
# Two  spaces: might wrap
    #287 Pointer{"Char_S"} clang_getCString(CXString)
    #288 Void clang_disposeString(CXString)
    #289 Typedef{"Pointer CXIndex"} clang_createIndex(int, int)
    #290 Void clang_disposeIndex(CXIndex)
    #293 Void clang_CXIndex_setGlobalOptions(CXIndex, unsigned int)
    #294 UInt clang_CXIndex_getGlobalOptions(CXIndex)

#296 Typedef{"Record CXString"} clang_getFileName(CXFile)
@cx CXString getFileName {CXFile}

    #297 Typedef{"Typedef time_t"} clang_getFileTime(CXFile)

#TODO: enable
  #298 UInt clang_isFileMultipleIncludeGuarded(CXTranslationUnit, CXFile)
  #@cx Uint32 isFileMultipleIncludeGuarded {CXTranslationUnit, CXFile}

    #299 Typedef{"Pointer CXFile"} clang_getFile(CXTranslationUnit, const char *)

#304 Typedef{"Record CXSourceLocation"} clang_getNullLocation()
@cx CXSourceLocation getNullLocation {}
#305 UInt clang_equalLocations(CXSourceLocation, CXSourceLocation)
@cx Uint32 equalLocations {CXSourceLocation, CXSourceLocation}

    #306 Typedef{"Record CXSourceLocation"} clang_getLocation(CXTranslationUnit, CXFile, unsigned int, unsigned int)
    #307 Typedef{"Record CXSourceLocation"} clang_getLocationForOffset(CXTranslationUnit, CXFile, unsigned int)

#308 Typedef{"Record CXSourceRange"} clang_getNullRange()
@cx CXSourceRange getNullRange {}
#309 Typedef{"Record CXSourceRange"} clang_getRange(CXSourceLocation, CXSourceLocation)
@cx CXSourceRange getRange {CXSourceLocation, CXSourceLocation}
#310 UInt clang_equalRanges(CXSourceRange, CXSourceRange)
@cx Uint32 equalRanges {CXSourceRange, CXSourceRange}
#311 Int clang_Range_isNull(CXSourceRange)
@cx Int32 Range_isNull {CXSourceRange}

    #312 Void clang_getExpansionLocation(CXSourceLocation, CXFile *, unsigned int *, unsigned int *, unsigned int *)
    #313 Void clang_getPresumedLocation(CXSourceLocation, CXString *, unsigned int *, unsigned int *)
    #314 Void clang_getInstantiationLocation(CXSourceLocation, CXFile *, unsigned int *, unsigned int *, unsigned int *)
    #315 Void clang_getSpellingLocation(CXSourceLocation, CXFile *, unsigned int *, unsigned int *, unsigned int *)

    #316 Typedef{"Record CXSourceLocation"} clang_getRangeStart(CXSourceRange)
    #317 Typedef{"Record CXSourceLocation"} clang_getRangeEnd(CXSourceRange)
    #321 UInt clang_getNumDiagnosticsInSet(CXDiagnosticSet)
    #322 Typedef{"Pointer CXDiagnostic"} clang_getDiagnosticInSet(CXDiagnosticSet, unsigned int)
    #324 Typedef{"Pointer CXDiagnosticSet"} clang_loadDiagnostics(const char *, enum CXLoadDiag_Error *, CXString *)
    #325 Void clang_disposeDiagnosticSet(CXDiagnosticSet)
    #326 Typedef{"Pointer CXDiagnosticSet"} clang_getChildDiagnostics(CXDiagnostic)
    #327 UInt clang_getNumDiagnostics(CXTranslationUnit)
    #328 Typedef{"Pointer CXDiagnostic"} clang_getDiagnostic(CXTranslationUnit, unsigned int)
    #329 Typedef{"Pointer CXDiagnosticSet"} clang_getDiagnosticSetFromTU(CXTranslationUnit)
    #330 Void clang_disposeDiagnostic(CXDiagnostic)
    #332 Typedef{"Record CXString"} clang_formatDiagnostic(CXDiagnostic, unsigned int)
    #333 UInt clang_defaultDiagnosticDisplayOptions()
    #334 Enum{"CXDiagnosticSeverity"} clang_getDiagnosticSeverity(CXDiagnostic)
    #335 Typedef{"Record CXSourceLocation"} clang_getDiagnosticLocation(CXDiagnostic)
    #336 Typedef{"Record CXString"} clang_getDiagnosticSpelling(CXDiagnostic)
    #337 Typedef{"Record CXString"} clang_getDiagnosticOption(CXDiagnostic, CXString *)
    #338 UInt clang_getDiagnosticCategory(CXDiagnostic)
    #339 Typedef{"Record CXString"} clang_getDiagnosticCategoryName(unsigned int)
    #340 Typedef{"Record CXString"} clang_getDiagnosticCategoryText(CXDiagnostic)
    #341 UInt clang_getDiagnosticNumRanges(CXDiagnostic)
    #342 Typedef{"Record CXSourceRange"} clang_getDiagnosticRange(CXDiagnostic, unsigned int)
    #343 UInt clang_getDiagnosticNumFixIts(CXDiagnostic)
    #344 Typedef{"Record CXString"} clang_getDiagnosticFixIt(CXDiagnostic, unsigned int, CXSourceRange *)
    #345 Typedef{"Record CXString"} clang_getTranslationUnitSpelling(CXTranslationUnit)
    #346 Typedef{"Pointer CXTranslationUnit"} clang_createTranslationUnitFromSourceFile(CXIndex, const char *, int, const char *const *, unsigned int, struct CXUnsavedFile *)
# TODO: enable here, move all to jl
  #347 Typedef{"Pointer CXTranslationUnit"} clang_createTranslationUnit(CXIndex, const char *)
  #349 UInt clang_defaultEditingTranslationUnitOptions()
  #350 Typedef{"Pointer CXTranslationUnit"} clang_parseTranslationUnit(CXIndex, const char *, const char *const *, int, struct CXUnsavedFile *, unsigned int, unsigned int)
  #352 UInt clang_defaultSaveOptions(CXTranslationUnit)
  #354 Int clang_saveTranslationUnit(CXTranslationUnit, const char *, unsigned int)
  #355 Void clang_disposeTranslationUnit(CXTranslationUnit)
  #357 UInt clang_defaultReparseOptions(CXTranslationUnit)
  #358 Int clang_reparseTranslationUnit(CXTranslationUnit, unsigned int, struct CXUnsavedFile *, unsigned int)
  #360 Pointer{"Char_S"} clang_getTUResourceUsageName(enum CXTUResourceUsageKind)
  #365 Typedef{"Record CXTUResourceUsage"} clang_getCXTUResourceUsage(CXTranslationUnit)
  #366 Void clang_disposeCXTUResourceUsage(CXTUResourceUsage)

#370 Typedef{"Record CXCursor"} clang_getNullCursor()
@cx CXCursor getNullCursor {}
#371 Typedef{"Record CXCursor"} clang_getTranslationUnitCursor(CXTranslationUnit)
@cx CXCursor getTranslationUnitCursor {CXTranslationUnit}

#372 UInt clang_equalCursors(CXCursor, CXCursor)
@cx Uint32 equalCursors {CXCursor, CXCursor}
#373 Int clang_Cursor_isNull(CXCursor)
@cx Int32 Cursor_isNull {CXCursor}
#374 UInt clang_hashCursor(CXCursor)
@cx Uint32 hashCursor {CXCursor}
#375 Enum{"CXCursorKind"} clang_getCursorKind(CXCursor)
@cx Uint32 getCursorKind {CXCursor}
#376 UInt clang_isDeclaration(enum CXCursorKind)
@cx Uint32 isDeclaration {CXCursorKind}
#377 UInt clang_isReference(enum CXCursorKind)
@cx Uint32 isReference {CXCursorKind}

    # these should be called directly...
    #378 UInt clang_isExpression(enum CXCursorKind)
    #379 UInt clang_isStatement(enum CXCursorKind)
    #380 UInt clang_isAttribute(enum CXCursorKind)
    #381 UInt clang_isInvalid(enum CXCursorKind)
    #382 UInt clang_isTranslationUnit(enum CXCursorKind)
    #383 UInt clang_isPreprocessing(enum CXCursorKind)
    #384 UInt clang_isUnexposed(enum CXCursorKind)

#386 Enum{"CXLinkageKind"} clang_getCursorLinkage(CXCursor)
@cx CXLinkageKind getCursorLinkage {CXCursor}
#387 Enum{"CXAvailabilityKind"} clang_getCursorAvailability(CXCursor)
@cx CXAvailabilityKind getCursorAvailability {CXCursor}
#389 Enum{"CXLanguageKind"} clang_getCursorLanguage(CXCursor)
@cx CXLanguageKind getCursorLanguage {CXCursor}
  # TODO: enable
  #390 Typedef{"Pointer CXTranslationUnit"} clang_Cursor_getTranslationUnit(CXCursor)
  #393 Typedef{"Pointer CXCursorSet"} clang_createCXCursorSet()
  #394 Void clang_disposeCXCursorSet(CXCursorSet)

    # don't think this is useful...
    #395 UInt clang_CXCursorSet_contains(CXCursorSet, CXCursor)
    #396 UInt clang_CXCursorSet_insert(CXCursorSet, CXCursor)

#397 Typedef{"Record CXCursor"} clang_getCursorSemanticParent(CXCursor)
@cx CXCursor getCursorSemanticParent {CXCursor}
#398 Typedef{"Record CXCursor"} clang_getCursorLexicalParent(CXCursor)
@cx CXCursor getCursorLexicalParent {CXCursor}

  #TODO: wrap this manually
  #399 Void clang_getOverriddenCursors(CXCursor, CXCursor **, unsigned int *)
  #400 Void clang_disposeOverriddenCursors(CXCursor *)
  # maybes..
  #401 Typedef{"Pointer CXFile"} clang_getIncludedFile(CXCursor)
  #402 Typedef{"Record CXCursor"} clang_getCursor(CXTranslationUnit, CXSourceLocation)

    # probably not
    #403 Typedef{"Record CXSourceLocation"} clang_getCursorLocation(CXCursor)
    #404 Typedef{"Record CXSourceRange"} clang_getCursorExtent(CXCursor)

#409 Typedef{"Record CXType"} clang_getCursorType(CXCursor)
@cx CXType getCursorType {CXCursor}
#410 Typedef{"Record CXType"} clang_getTypedefDeclUnderlyingType(CXCursor)
@cx CXType getTypedefDeclUnderlyingType {CXCursor}
#411 Typedef{"Record CXType"} clang_getEnumDeclIntegerType(CXCursor)
@cx CXType getEnumDeclIntegerType {CXCursor}
#412 LongLong clang_getEnumConstantDeclValue(CXCursor)
@cx Int64 getEnumConstantDeclValue {CXCursor}
#413 ULongLong clang_getEnumConstantDeclUnsignedValue(CXCursor)
@cx Uint64 getEnumConstantDeclUnsignedValue {CXCursor}
#414 Int clang_Cursor_getNumArguments(CXCursor)
@cx Int32 Cursor_getNumArguments {CXCursor}
#415 Typedef{"Record CXCursor"} clang_Cursor_getArgument(CXCursor, unsigned int)
@cx CXCursor Cursor_getArgument {CXCursor, Int32}
#416 UInt clang_equalTypes(CXType, CXType)
@cx Uint32 equalTypes {CXType, CXType}
#417 Typedef{"Record CXType"} clang_getCanonicalType(CXType)
@cx CXType getCanonicalType {CXType}
#418 Uint32 clang_isConstQualifiedType(CXType)
@cx Uint32 isConstQualifiedType {CXType}
#419 UInt clang_isVolatileQualifiedType(CXType)
@cx Uint32 isVolatileQualifiedType {CXType}
#420 UInt clang_isRestrictQualifiedType(CXType)
@cx Uint32 isRestrictQualifiedType {CXType}
#421 Typedef{"Record CXType"} clang_getPointeeType(CXType)
@cx CXType getPointeeType {CXType}
#422 Typedef{"Record CXCursor"} clang_getTypeDeclaration(CXType)
@cx CXCursor getTypeDeclaration {CXType}
#423 Typedef{"Record CXString"} clang_getDeclObjCTypeEncoding(CXCursor)
@cx CXString getDeclObjCTypeEncoding {CXCursor}
#424 Typedef{"Record CXString"} clang_getTypeKindSpelling(enum CXTypeKind)
@cx CXString getTypeKindSpelling {CXTypeKind}
#425 Enum{"CXCallingConv"} clang_getFunctionTypeCallingConv(CXType)
@cx Int32 getFunctionTypeCallingConv {CXType}
#426 Typedef{"Record CXType"} clang_getResultType(CXType)
@cx CXType getResultType {CXType}
#427 Int clang_getNumArgTypes(CXType)
@cx Int32 getNumArgTypes {CXType}
#428 Typedef{"Record CXType"} clang_getArgType(CXType, unsigned int)
@cx CXType getArgType {CXType, Uint}
#429 UInt clang_isFunctionTypeVariadic(CXType)
@cx Uint32 isFunctionTypeVariadic {CXType}
#430 Typedef{"Record CXType"} clang_getCursorResultType(CXCursor)
@cx CXType getCursorResultType {CXCursor}
#431 UInt clang_isPODType(CXType)
@cx Uint32 isPODType {CXType}
#432 Typedef{"Record CXType"} clang_getElementType(CXType)
@cx CXType getElementType {CXType}
#433 LongLong clang_getNumElements(CXType)
@cx Int64 getNumElements {CXType}
#434 Typedef{"Record CXType"} clang_getArrayElementType(CXType)
@cx CXType getArrayElementType {CXType}
#435 LongLong clang_getArraySize(CXType)
@cx Int64 getArraySize {CXType}
#436 UInt clang_isVirtualBase(CXCursor)
@cx Uint32 isVirtualBase {CXCursor}
#438 Enum{"CX_CXXAccessSpecifier"} clang_getCXXAccessSpecifier(CXCursor)
@cx Int getCXXAccessSpecifier {CXCursor}
#439 UInt clang_getNumOverloadedDecls(CXCursor)
@cx Uint32 getNumOverloadedDecls {CXCursor}
#440 Typedef{"Record CXCursor"} clang_getOverloadedDecl(CXCursor, unsigned int)
@cx CXCursor getOverloadedDecl {CXCursor, Uint32}
    # what is this?
    #441 Typedef{"Record CXType"} clang_getIBOutletCollectionType(CXCursor)
    # Do not wrap.
    #444 UInt clang_visitChildren(CXCursor, CXCursorVisitor, CXClientData)
#445 Typedef{"Record CXString"} clang_getCursorUSR(CXCursor)
@cx CXString getCursorUSR {CXCursor}

    #446 Typedef{"Record CXString"} clang_constructUSR_ObjCClass(const char *)
    #447 Typedef{"Record CXString"} clang_constructUSR_ObjCCategory(const char *, const char *)
    #448 Typedef{"Record CXString"} clang_constructUSR_ObjCProtocol(const char *)
    #449 Typedef{"Record CXString"} clang_constructUSR_ObjCIvar(const char *, CXString)
    #450 Typedef{"Record CXString"} clang_constructUSR_ObjCMethod(const char *, unsigned int, CXString)
    #451 Typedef{"Record CXString"} clang_constructUSR_ObjCProperty(const char *, CXString)

#452 Typedef{"Record CXString"} clang_getCursorSpelling(CXCursor)
@cx CXString getCursorSpelling {CXCursor}
  # maybe
  #453 Typedef{"Record CXSourceRange"} clang_Cursor_getSpellingNameRange(CXCursor, unsigned int, unsigned int)
#454 Typedef{"Record CXString"} clang_getCursorDisplayName(CXCursor)
@cx CXString getCursorDisplayName {CXCursor}
#455 Typedef{"Record CXCursor"} clang_getCursorReferenced(CXCursor)
@cx CXCursor getCursorReferenced {CXCursor}
#456 Typedef{"Record CXCursor"} clang_getCursorDefinition(CXCursor)
@cx CXCursor getCursorDefinition {CXCursor}
#457 UInt clang_isCursorDefinition(CXCursor)
@cx Uint32 isCursorDefinition {CXCursor}
#458 Typedef{"Record CXCursor"} clang_getCanonicalCursor(CXCursor)
@cx CXCursor getCanonicalCursor {CXCursor}
    #459 Int clang_Cursor_getObjCSelectorIndex(CXCursor)
#460 UInt clang_CXXMethod_isStatic(CXCursor)
@cx Uint32 CXXMethod_isStatic {CXCursor}
#461 UInt clang_CXXMethod_isVirtual(CXCursor)
@cx Uint32 CXXMethod_isVirtual {CXCursor}
#462 Enum{"CXCursorKind"} clang_getTemplateCursorKind(CXCursor)
@cx Int32 getTemplateCursorKind {CXCursor}
#463 Typedef{"Record CXCursor"} clang_getSpecializedCursorTemplate(CXCursor)
@cx CXCursor getSpecializedCursorTemplate {CXCursor}

    #464 Typedef{"Record CXSourceRange"} clang_getCursorReferenceNameRange(CXCursor, unsigned int, unsigned int)
    # there is no way to get a CXToken..
    # also, cindex.py does not expose.
    #470 Typedef{"Enum CXTokenKind"} clang_getTokenKind(CXToken)
    @cx Int32 getTokenKind {CXToken}
    #471 Typedef{"Record CXString"} clang_getTokenSpelling(CXTranslationUnit, CXToken)
    @cx CXString getTokenSpelling {CXTranslationUnit, CXToken}

    #472 Typedef{"Record CXSourceLocation"} clang_getTokenLocation(CXTranslationUnit, CXToken)
    #473 Typedef{"Record CXSourceRange"} clang_getTokenExtent(CXTranslationUnit, CXToken)
    #474 Void clang_tokenize(CXTranslationUnit, CXSourceRange, CXToken **, unsigned int *)
    #475 Void clang_annotateTokens(CXTranslationUnit, CXToken *, unsigned int, CXCursor *)
    #476 Void clang_disposeTokens(CXTranslationUnit, CXToken *, unsigned int)

#477 Typedef{"Record CXString"} clang_getCursorKindSpelling(enum CXCursorKind)
@cx CXString getCursorKindSpelling {CXCursorKind}

    #478 Void clang_getDefinitionSpellingAndExtent(CXCursor, const char **, const char **, unsigned int *, unsigned int *, unsigned int *, unsigned int *)

  # these sound useful..
  #479 Void clang_enableStackTraces()
  #480 Void clang_executeOnThread(void (*)(void *), void *, unsigned int)

    # writing a C++ editor in Julia? let me be the first to salute you.
    #485 Enum{"CXCompletionChunkKind"} clang_getCompletionChunkKind(CXCompletionString, unsigned int)
    #486 Typedef{"Record CXString"} clang_getCompletionChunkText(CXCompletionString, unsigned int)
    #487 Typedef{"Pointer CXCompletionString"} clang_getCompletionChunkCompletionString(CXCompletionString, unsigned int)
    #488 UInt clang_getNumCompletionChunks(CXCompletionString)
    #489 UInt clang_getCompletionPriority(CXCompletionString)
    #490 Enum{"CXAvailabilityKind"} clang_getCompletionAvailability(CXCompletionString)
    #491 UInt clang_getCompletionNumAnnotations(CXCompletionString)
    #492 Typedef{"Record CXString"} clang_getCompletionAnnotation(CXCompletionString, unsigned int)
    #493 Typedef{"Record CXString"} clang_getCompletionParent(CXCompletionString, enum CXCursorKind *)
    #494 Typedef{"Pointer CXCompletionString"} clang_getCursorCompletionString(CXCursor)
    #499 UInt clang_defaultCodeCompleteOptions()
    #500 Pointer{"Typedef", "Record", "CXCodeCompleteResults"} clang_codeCompleteAt(CXTranslationUnit, const char *, unsigned int, unsigned int, struct CXUnsavedFile *, unsigned int, unsigned int)
    #501 Void clang_sortCodeCompletionResults(CXCompletionResult *, unsigned int)
    #502 Void clang_disposeCodeCompleteResults(CXCodeCompleteResults *)
    #503 UInt clang_codeCompleteGetNumDiagnostics(CXCodeCompleteResults *)
    #504 Typedef{"Pointer CXDiagnostic"} clang_codeCompleteGetDiagnostic(CXCodeCompleteResults *, unsigned int)
    #505 ULongLong clang_codeCompleteGetContexts(CXCodeCompleteResults *)
    #506 Enum{"CXCursorKind"} clang_codeCompleteGetContainerKind(CXCodeCompleteResults *, unsigned int *)
    #507 Typedef{"Record CXString"} clang_codeCompleteGetContainerUSR(CXCodeCompleteResults *)
    #508 Typedef{"Record CXString"} clang_codeCompleteGetObjCSelector(CXCodeCompleteResults *)

#509 Typedef{"Record CXString"} clang_getClangVersion()
@cx CXString getClangVersion {}

  # TODO: useful?
  #510 Void clang_toggleCrashRecovery(unsigned int)
  #512 Void clang_getInclusions(CXTranslationUnit, CXInclusionVisitor, CXClientData)

    #514 Typedef{"Pointer CXRemapping"} clang_getRemappings(const char *)
    #515 Typedef{"Pointer CXRemapping"} clang_getRemappingsFromFileList(const char **, unsigned int)
    #516 UInt clang_remap_getNumFiles(CXRemapping)
    #517 Void clang_remap_getFilenames(CXRemapping, unsigned int, CXString *, CXString *)
    #518 Void clang_remap_dispose(CXRemapping)

  # TODO: useful?
  #522 Void clang_findReferencesInFile(CXCursor, CXFile, CXCursorAndRangeVisitor)

    #575 Int clang_index_isEntityObjCContainerKind(CXIdxEntityKind)
    #576 Pointer{"Typedef", "Record", "CXIdxObjCContainerDeclInfo"} clang_index_getObjCContainerDeclInfo(const CXIdxDeclInfo *)
    #577 Pointer{"Typedef", "Record", "CXIdxObjCInterfaceDeclInfo"} clang_index_getObjCInterfaceDeclInfo(const CXIdxDeclInfo *)
    #578 Pointer{"Typedef", "Record", "CXIdxObjCCategoryDeclInfo"} clang_index_getObjCCategoryDeclInfo(const CXIdxDeclInfo *)
    #579 Pointer{"Typedef", "Record", "CXIdxObjCProtocolRefListInfo"} clang_index_getObjCProtocolRefListInfo(const CXIdxDeclInfo *)
    #580 Pointer{"Typedef", "Record", "CXIdxObjCPropertyDeclInfo"} clang_index_getObjCPropertyDeclInfo(const CXIdxDeclInfo *)
    #581 Pointer{"Typedef", "Record", "CXIdxIBOutletCollectionAttrInfo"} clang_index_getIBOutletCollectionAttrInfo(const CXIdxAttrInfo *)
    #582 Pointer{"Typedef", "Record", "CXIdxCXXClassDeclInfo"} clang_index_getCXXClassDeclInfo(const CXIdxDeclInfo *)
    #583 Typedef{"Pointer CXIdxClientContainer"} clang_index_getClientContainer(const CXIdxContainerInfo *)
    #584 Void clang_index_setClientContainer(const CXIdxContainerInfo *, CXIdxClientContainer)
    #585 Typedef{"Pointer CXIdxClientEntity"} clang_index_getClientEntity(const CXIdxEntityInfo *)
    #586 Void clang_index_setClientEntity(const CXIdxEntityInfo *, CXIdxClientEntity)
    #588 Typedef{"Pointer CXIndexAction"} clang_IndexAction_create(CXIndex)
    #589 Void clang_IndexAction_dispose(CXIndexAction)
    #592 Int clang_indexSourceFile(CXIndexAction, CXClientData, IndexerCallbacks *, unsigned int, unsigned int, const char *, const char *const *, int, struct CXUnsavedFile *, unsigned int, CXTranslationUnit *, unsigned int)
    #593 Int clang_indexTranslationUnit(CXIndexAction, CXClientData, IndexerCallbacks *, unsigned int, unsigned int, CXTranslationUnit)
    #594 Void clang_indexLoc_getFileLocation(CXIdxLoc, CXIdxClientFile *, CXFile *, unsigned int *, unsigned int *, unsigned int *)
    #595 Typedef{"Record CXSourceLocation"} clang_indexLoc_getCXSourceLocation(CXIdxLoc)

write_output()


text:
################################################
## Specialized routines for tridiagonal matrices
################################################

export numnegevals

"""
Computes the number of negative eigenvalues of T - σI, a.k.a. spectrum slicing

Inputs:
    T: A SymTridiagonal{<:Real} matrix
    σ: The shift parameter

Outputs:
    ν: The number of negative eigenvalues

Reference:
    B. N. Parlett, "The symmetric eigenvalue problem", Section 3.3.1, p. 52.
"""
function numnegevals(T::SymTridiagonal{S}, σ::S = zero(S)) where S
    α = T.dv
    β = T.ev
    n = length(α)
    ϵ = eps(S)
    δ = α[1] - σ
    ν = δ < 0 ? 1 : 0
    for k = 1:n - 1
        if δ == 0
        info("zero in iteration $k")
        δ = ϵ * (β[k]+ϵ) #Parlett prefers adjusting σ and starting again
    end
        δ = (α[k + 1] - σ) - β[k]*(β[k]/δ)
    ν += (δ < 0)
    end
    ν
end


text:
# wrappers of low-level functionality

function curandGetProperty(property::libraryPropertyType)
  value_ref = Ref{Cint}()
  curandGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(curandGetProperty(MAJOR_VERSION),
                          curandGetProperty(MINOR_VERSION),
                          curandGetProperty(PATCH_LEVEL))


text:
include("shortcuts.jl")
include("sequence.jl")
include("Diff.jl")
include("TrivialGate.jl")
include("Bag.jl")
include("Mod.jl")
include("ConditionBlock.jl")
include("RotBasis.jl")
include("reflect_gate.jl")
include("math_gate.jl")
include("pauli_strings.jl")
include("qft.jl")
include("FSimGate.jl")


text:
"""
    ==ₜ(x, y)

Check that _type_ and value of `x` and `y` are equal.
"""
==ₜ(_, _) = false
==ₜ(x::T, y::T) where T = x == y

"""
    ==ₛ(a::AbstractString, b::AbstractString)

Equality check ignoring white spaces
"""
==ₛ(a::AbstractString, b::AbstractString) =
    replace(a, r"\s" => "") == replace(b, r"\s" => "")

"""
    ==ᵣ(a::AbstractString, b::AbstractString)

Equality check appropriate for comparing `repr` output.
"""
==ᵣ
if VERSION >= v"1.6-"
    const ==ᵣ = ==ₛ
else
    const ==ᵣ = ==
end


text:
abstract type AbstractAlgorithm end

mutable struct Algorithm{T} <: AbstractAlgorithm
    parameters::T
    status::State
    information::Information
    options::Options
end

function Algorithm(
    parameters;
    initial_state::State = State(nothing, []),
    information::Information = Information(),
    options::Options = Options(),
)

    Algorithm(parameters, initial_state, information, options)

end

function Base.show(io::IO, alg::Algorithm)
    Base.show(io, alg.parameters)
end

mutable struct Problem <: AbstractProblem
    f::Function
    bounds::Array{Float64,2}
    f_calls::Int
end

function Problem(f::Function, bounds::Array)

    if size(bounds,1) > 2 && size(bounds,2) == 2
        bounds = Array(bounds')
    end

    Problem(f, bounds, 0)
end

function evaluate(x, problem::Problem)
    problem.f_calls += 1
    return problem.f(x)
end


text:
module Diagnostics

using ClimateMARGO.Models
using ClimateMARGO.Utils

export
    t, future_mask,
    ramp_emissions, emissions, effective_emissions,
    c, F, Tslow, Tfast, T,
    τd, B, F2x, ECS,
    discount, f, E,
    damage, cost, benefit,
    net_benefit, net_present_cost, net_present_benefit

include("carbon.jl")
include("energy_balance.jl")
include("cost_benefit.jl")
include("utils.jl")

end


text:
using Flux
using Flux: throttle, nfan, glorot_uniform, glorot_normal,
             kaiming_normal, kaiming_uniform, orthogonal,
             sparse_init, stack, unstack, Zeros, batch, unbatch
using StatsBase: var, std
using Random
using Test

@testset "Throttle" begin
  @testset "default behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=true, trailing=false)
    f()
    f()
    f()
    sleep(1.01)
    @test length(a) == 1
  end

  @testset "leading behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=true, trailing=false)
    f()
    @test length(a) == 1
    f()
    @test length(a) == 1
    sleep(1.01)
    f()
    @test length(a) == 2
  end

  @testset "trailing behaviour" begin
    a = []
    f = throttle(()->push!(a, time()), 1, leading=false, trailing=true)
    f()
    @test length(a) == 0
    f()
    @test length(a) == 0
    sleep(1.01)
    @test length(a) == 1
  end

  @testset "arguments" begin
    a = []
    f = throttle((x)->push!(a, x), 1, leading=true, trailing=true)
    f(1)
    @test a == [1]
    f(2)
    @test a == [1]
    f(3)
    @test a == [1]
    sleep(1.01)
    @test a == [1, 3]
  end
end

@testset "Initialization" begin
  # Set random seed so that these tests don't fail randomly
  Random.seed!(0)

  @testset "Fan in/out" begin
    @test nfan() == (1, 1) #For a constant
    @test nfan(100) == (1, 100) #For vector
    @test nfan(100, 200) == (200, 100) == nfan((100, 200)) #For Dense layer
    @test nfan(2, 30, 40) == (2 * 30, 2 * 40) #For 1D Conv layer
    @test nfan(2, 3, 40, 50) == (2 * 3 * 40, 2 * 3 * 50) #For 2D Conv layer
    @test nfan(2, 3, 4, 50, 60) == (2 * 3 * 4 * 50, 2 * 3 * 4 * 60) #For 3D Conv layer
  end

  @testset "glorot" begin
    # glorot_uniform and glorot_normal should both yield a kernel with
    # variance ≈ 2/(fan_in + fan_out)
    for dims ∈ [(1000,), (100, 100), (100, 400), (2, 3, 32, 64), (2, 3, 4, 32, 64)]
      for init ∈ [glorot_uniform, glorot_normal]
        v = init(dims...)
        fan_in, fan_out = nfan(dims...)
        σ2 = 2 / (fan_in + fan_out)
        @test 0.9σ2 < var(v) < 1.1σ2
        @test eltype(v) == Float32
      end
    end
  end

  @testset "kaiming" begin
    # kaiming_uniform should yield a kernel in range [-sqrt(6/n_out), sqrt(6/n_out)]
    # and kaiming_normal should yield a kernel with stddev ~= sqrt(2/n_out)
    for (n_in, n_out) in [(100, 100), (100, 400)]
      v = kaiming_uniform(n_in, n_out)
      σ2 = sqrt(6/n_out)
      @test -1σ2  < minimum(v) < -0.9σ2
      @test 0.9σ2  < maximum(v) < 1σ2
      @test eltype(v) == Float32

      v = kaiming_normal(n_in, n_out)
      σ2 = sqrt(2/n_out)
      @test 0.9σ2 < std(v) < 1.1σ2
      @test eltype(v) == Float32
    end
  end

  @testset "orthogonal" begin
    # A matrix of dim = (m,n) with m > n should produce a QR decomposition. In the other case, the transpose should be taken to compute the QR decomposition.
    for (rows,cols) in [(5,3),(3,5)]
      v = orthogonal(rows, cols)
      rows < cols ? (@test v * v' ≈ I(rows)) : (@test v' * v ≈ I(cols))
    end
    for mat in [(3,4,5),(2,2,5)]
      v = orthogonal(mat...)
      cols = mat[end]
      rows = div(prod(mat),cols)
      v = reshape(v, (rows,cols))
      rows < cols ? (@test v * v' ≈ I(rows)) : (@test v' * v ≈ I(cols))
    end
  end

  @testset "sparse_init" begin
    # sparse_init should yield an error for non 2-d dimensions
    # sparse_init should yield no zero elements if sparsity < 0
    # sparse_init should yield all zero elements if sparsity > 1
    # sparse_init should yield exactly ceil(n_in * sparsity) elements in each column for other sparsity values
    # sparse_init should yield a kernel in its non-zero elements consistent with the std parameter

    @test_throws ArgumentError sparse_init(100, 100, 100, sparsity=0.1)
    v = sparse_init(100, 100, sparsity=-0.1)
    @test sum(v .== 0) == 0
    @test eltype(v) == Float32
    v = sparse_init(100, 100, sparsity=1.1)
    @test sum(v .== 0) == length(v)
    @test eltype(v) == Float32

    for (n_in, n_out, sparsity, σ) in [(100, 100, 0.25, 0.1), (100, 400, 0.75, 0.01)]
      expected_zeros = ceil(Integer, n_in * sparsity)
      v = sparse_init(n_in, n_out, sparsity=sparsity, std=σ)
      @test all([sum(v[:,col] .== 0) == expected_zeros for col in 1:n_out])
      @test 0.9 * σ < std(v[v .!= 0]) < 1.1 * σ
      @test eltype(v) == Float32
    end
  end

  @testset "partial_application" begin
    big = 1e9

    partial_ku = kaiming_uniform(gain=big)
    @test maximum(partial_ku(8, 8)) > big / 2
    @test maximum(partial_ku(8, 8, gain=1)) < big / 2

    partial_kn = kaiming_normal(gain=big)
    @test maximum(partial_kn(8, 8)) > big / 2
    @test maximum(partial_kn(8, 8, gain=1)) < big / 2

    partial_si = sparse_init(sparsity=1)
    @test maximum(partial_si(8, 8)) == 0
    @test maximum(partial_si(8, 8, sparsity=0)) > 0
  end

  @testset "identity_init" begin
    import Flux: identity_init

    @testset "Basic" begin
      partial = identity_init(gain=3)
      @test partial(3, 3) == identity_init(3, 3; gain=3) == [3 0 0; 0 3 0; 0 0 3]
    end

    @testset "Non-identity sizes" begin
        @test identity_init(2, 3)[:, end] == zeros(Float32, 2)
        @test identity_init(3, 2; shift=1)[1, :] == zeros(Float32, 2)
        @test identity_init(1, 1, 3, 4)[:, :, :, end] == zeros(Float32, 1, 1, 3)
        @test identity_init(2, 1, 3, 3)[end, :, :, :] == zeros(Float32, 1, 3, 3)
        @test identity_init(1, 2, 3, 3)[:, end, :, :] == zeros(Float32, 1, 3, 3)
    end

    @testset "Dense ID mapping" begin
        l = Dense(3,3, init = identity_init)

        indata = reshape(collect(Float32, 1:9), 3, 3)
        @test l(indata) == indata
    end

    @testset "$layer ID mapping with kernelsize $kernelsize" for layer in (Conv, ConvTranspose, CrossCor), kernelsize in (
        (1,),
        (3,),
        (1, 3),
        (3, 5),
        (3, 5, 7))
        nch = 3
        l = layer(kernelsize, nch=>nch, init=identity_init, pad=SamePad())

        indata = randn(Float32, kernelsize..., nch, nch)
        @test l(indata) == indata
    end

    @testset "Inception identity" begin
      insize = 7
      path1 = Conv((1, 3), insize=>2; init=identity_init, pad=SamePad())
      path2 = Conv((3, 5), insize=>3; init=identity_init(shift=(0, 0, 2, 0)), pad=SamePad())
      path3 = Conv((5, 7), insize=>2; init=identity_init(shift=(0, 0, 5, 0)), pad=SamePad())
      block = Parallel((xs...) -> cat(xs...;dims=3), path1, path2, path3)

      indata = randn(Float32, 9, 9, 7, 2)
      @test block(indata) == indata
    end
  end
end

@testset "Params" begin
  m = Dense(10, 5)
  @test size.(params(m)) == [(5, 10), (5,)]
  m = RNN(10, 5)
  @test size.(params(m)) == [(5, 10), (5, 5), (5,), (5, 1)]

  # Layer duplicated in same chain, params just once pls.
  c = Chain(m, m)
  @test size.(params(c)) == [(5, 10), (5, 5), (5,), (5, 1)]

  # Self-referential array. Just want params, no stack overflow pls.
  r = Any[nothing,m]
  r[1] = r
  @test size.(params(r)) == [(5, 10), (5, 5), (5,), (5, 1)]
end

@testset "Basic Stacking" begin
  x = randn(3,3)
  stacked = stack([x, x], 2)
  @test size(stacked) == (3,2,3)
end

@testset "Precision" begin
  m = Chain(Dense(10, 5, relu), Dense(5, 2))
  x64 = rand(Float64, 10)
  x32 = rand(Float32, 10)
  @test eltype(m[1].weight) == Float32
  @test eltype(m(x32)) == Float32
  @test eltype(m(x64)) == Float64
  @test eltype(f64(m)(x32)) == Float64
  @test eltype(f64(m)(x64)) == Float64
  @test eltype(f64(m)[1].weight) == Float64
  @test eltype(f32(f64(m))[1].weight) == Float32
end

@testset "Zeros" begin
  m = Dense(3,2; bias=false)
  @test f64(m).bias === m.bias === Zeros()
  @test f32(m).bias === m.bias === Zeros()

  @testset "Gradients for broadcasted $op with sizes $s" for op in (+,-,*), s in ((1,), (2,3))
    o = ones(s)
    z = zeros(s)
    Z = Zeros()

    @testset "Explicit" begin
      gfun(args...) = gradient((x, y) -> sum(op.(x,y)), args...)
      g = gfun(o, z)
      @test gfun(o, Z) == (g[1], nothing)

      g = gfun(z, o)
      @test gfun(Z, o) == (nothing, g[2])
    end

    @testset "Implicit" begin
      gfun(args...) = gradient(() -> sum(op.(args...)), params(collect(args)))
      g = gfun(o, z)

      gres = gfun(o, Z)
      @test gres[o] == g[o]
      @test Z ∉ gres.params

      g = gfun(z, o)
      gres = gfun(Z, o)
      @test gres[o] == g[o]
      @test Z ∉ gres.params
    end
  end

  @testset "Gradients for broadcasted / with sizes $s" for s in ((1,), (2,3))
    o = ones(s)
    z = zeros(s)
    Z = Zeros() # Only defined for 0-dim

    @testset "Explicit" begin
      gfun(args...) = gradient((x, y) -> sum(x ./ y), args...)
      g = gfun(z, o)
      @test gfun(Z, o) == (nothing, g[2])
    end

    @testset "Implicit" begin
      gfun(x,y) = gradient(() -> sum(x ./ y), params([x,y]))

      g = gfun(z, o)
      gres = gfun(Z, o)
      @test gres[o] == g[o]
      @test Z ∉ gres.params
    end
  end

  @testset "Gradients for $op with sizes $s" for op in (+,-), s in (tuple(), (1,), (2,3))
    o = ones(s)
    z = zeros(s)
    Z = Zeros()

    @testset "Explicit" begin
      gfun(args...) = gradient((x, y) -> sum(op(x,y)), args...)

      g = gfun(o, z)
      @test gfun(o, Z) == (g[1], nothing)

      g = gfun(z, o)
      @test gfun(Z, o) == (nothing, g[2])
    end

    @testset "Implicit" begin
      gfun(args...) = gradient(() -> sum(op(args...)), params(collect(args)))
      g = gfun(o, z)
      gres = gfun(o, Z)
      @test gres[o] == g[o]
      @test Z ∉ gres.params

      g = gfun(z, o)
      gres = gfun(Z, o)
      @test gres[o] == g[o]
      @test Z ∉ gres.params
    end
  end
end

@testset "Stacking" begin
  stacked_array=[ 8 9 3 5; 9 6 6 9; 9 1 7 2; 7 4 10 6 ]
  unstacked_array=[[8, 9, 9, 7], [9, 6, 1, 4], [3, 6, 7, 10], [5, 9, 2, 6]]
  @test unstack(stacked_array, 2) == unstacked_array
  @test stack(unstacked_array, 2) == stacked_array
  @test stack(unstack(stacked_array, 1), 1) == stacked_array
end

@testset "Batching" begin
  stacked_array=[ 8 9 3 5
                  9 6 6 9
                  9 1 7 2
                  7 4 10 6 ]
  unstacked_array=[[8, 9, 9, 7], [9, 6, 1, 4], [3, 6, 7, 10], [5, 9, 2, 6]]
  @test unbatch(stacked_array) == unstacked_array
  @test batch(unstacked_array) == stacked_array

  # no-op for vector of non-arrays
  @test batch([1,2,3]) == [1,2,3]
  @test unbatch([1,2,3]) == [1,2,3]

  # generic iterable
  @test batch(ones(2) for i=1:3) == ones(2, 3)
  @test unbatch(ones(2, 3)) == [ones(2) for i=1:3]
end

@testset "Param remapping" begin
  ls(dims...) = reshape(collect(Float32, 1:prod(dims)), dims...) # accepts dims in reverse order to Dense
  dl(nin, nout, bias) = Dense(ls(nout, nin), bias(nout))
  dm(bias) = Chain(
    dl(3, 5, bias),
    dl(5, 4, bias),
    dl(4, 3, bias)
  )

  nobias(n) = Zeros()
  testdense(m, bt) = @testset "Check layer $i" for (i, (l1, l2)) in enumerate(zip(m, dm(bt)))
    @test l1.weight == l2.weight
    @test l1.bias == l2.bias
    @test_skip typeof(l1.bias) === typeof(l2.bias)
  end

  @testset "loadparams!" begin
    import Flux: loadparams!
    pars(w, b) = [w, b]
    import Flux: loadparams!, Zeros

    pars(w, b::Zeros) = [w, Flux.zeros32(size(w,1))]
    pars(l) = pars(l.weight, l.bias)
    pararray(m) = mapreduce(pars, vcat, m)
    weights(m) = mapreduce(l -> [l.weight], vcat, m)
    @testset "Bias type $bt" for bt in (Flux.zeros32, nobias)
      m = dm(bt)
      loadparams!(m, params(m))
      testdense(m, bt)
    end

    @testset "$b1 to $b2" for (b1, b2, be) in (
      (Flux.zeros32, Flux.ones32, Flux.ones32),   # Load ones as bias to a model with zeros as bias -> model gets ones as bias
      (Flux.ones32, nobias, Flux.zeros32), # Load Zeros as bias to a model with ones as bias-> model gets zeros as bias
      (nobias, Flux.ones32, nobias),     # Load ones as bias to a model with Zeros as bias-> model bias does not change
    )
      m1 = dm(b1)
      m2 = dm(b2)
      loadparams!(m1, b1 == nobias ? weights(m2) : pararray(m2))
      testdense(m1, be)
    end
  end

  @testset "destructure" begin
    import Flux: destructure
    @testset "Bias type $bt" for bt in (zeros, nobias)
      m = dm(bt)
      p, re = destructure(m)
      testdense(re(p), bt)
    end

    @testset "restructure in gradient" begin
      x = rand(Float32, 3, 1)
      m = dm(zeros)
      ∇m = gradient(m -> sum(m(x)), m)[1]
      p, re = destructure(m)
      ∇p = gradient(θ -> sum(re(θ)(x)), p)[1]
      @test ∇p ≈ destructure(∇m)[1]
    end
  end
end

@testset "Train and test mode" begin
  mutable struct DummyLayer
    testing::Bool
  end
  Flux.testmode!(m::DummyLayer, testing=true) = (m.testing = testing; m)

  c = Chain(DummyLayer(true))
  testmode!(c)
  @test c[1].testing
  trainmode!(c)
  @test !c[1].testing
end

@testset "modules" begin
  m1 = Conv((2,3), 4=>5; pad=6, stride=7)
  m2 = LayerNorm(8)
  m3 = m2.diag
  m4 = SkipConnection(m1, +)
  m5 =  Chain(m4, m2)
  modules = Flux.modules(m5)
  # Depth-first descent
  @test length(modules) == 5
  @test modules[1] === m5
  @test modules[2] === m4
  @test modules[3] === m1
  @test modules[4] === m2
  @test modules[5] === m3

  modules = Flux.modules(Chain(Dense(2,3), BatchNorm(3), LSTM(3,4)))
  @test length(modules) == 5

  modules = Flux.modules(Chain(SkipConnection(
                                  Conv((2,3), 4=>5; pad=6, stride=7),
                                  +),
                                LayerNorm(8)))
  @test length(modules) == 5
end

@testset "Patience triggers" begin
  @testset "patience" begin
    trigger = Flux.patience(() -> true, 3)

    @test trigger() == false
    @test trigger() == false
    @test trigger() == true

    v = [false, true, false, true, true, true]
    trigger = let v = v
      Flux.patience(i -> v[i], 3)
    end

    n_iter = 0
    for i in 1:length(v)
      trigger(i) && break
      n_iter += 1
    end

    @test n_iter == 5
  end

  @testset "early stopping" begin
    @testset "args & kwargs" begin
      es = Flux.early_stopping((x; y = 1) -> x + y, 10; min_dist=3)

      n_iter = 0
      while n_iter < 99
        es(-n_iter; y=-n_iter) && break
        n_iter += 1
      end

      @test n_iter == 9
    end

    @testset "distance" begin
      es = Flux.early_stopping(identity, 10; distance=(best_score, score) -> score - best_score)

      n_iter = 0
      while n_iter < 99
        es(n_iter) && break
        n_iter += 1
      end

      @test n_iter == 99
    end

    @testset "init_score" begin
      es = Flux.early_stopping(identity, 10; init_score=10)

      n_iter = 0
      while n_iter < 99
        es(n_iter) && break
        n_iter += 1
      end

      @test n_iter == 10
    end
  end

  @testset "plateau" begin
    f = let v = 10
      () -> v = v / abs(v) - v
    end

    trigger = Flux.plateau(f, 3, init_score=10, min_dist=18)

    n_iter = 0
    while n_iter < 99
      trigger() && break
      n_iter += 1
    end

    @test n_iter == 3
  end
end


text:
# Structured matrices
using LinearAlgebra: AbstractTriangular

# Matrix wrapper types that we know are square and are thus potentially invertible. For
# these we can use simpler definitions for `/` and `\`.
const SquareMatrix{T} = Union{Diagonal{T}, AbstractTriangular{T}}

function rrule(::typeof(/), A::AbstractMatrix{<:Real}, B::T) where T<:SquareMatrix{<:Real}
    Y = A / B
    function slash_pullback(Ȳ)
        ∂A = @thunk Ȳ / B'
        ∂B = @thunk _unionall_wrapper(T)(-Y' * (Ȳ / B'))
        return (NO_FIELDS, ∂A, ∂B)
    end
    return Y, slash_pullback
end

function rrule(::typeof(\), A::T, B::AbstractVecOrMat{<:Real}) where T<:SquareMatrix{<:Real}
    Y = A \ B
    function backslash_pullback(Ȳ)
        ∂A = @thunk _unionall_wrapper(T)(-(A' \ Ȳ) * Y')
        ∂B = @thunk A' \ Ȳ
        return NO_FIELDS, ∂A, ∂B
    end
    return Y, backslash_pullback
end

#####
##### `Diagonal`
#####

function rrule(::Type{<:Diagonal}, d::AbstractVector)
    function Diagonal_pullback(ȳ::AbstractMatrix)
        return (NO_FIELDS, diag(ȳ))
    end
    function Diagonal_pullback(ȳ::Composite)
        # TODO: Assert about the primal type in the Composite, It should be Diagonal
        # infact it should be exactly the type of `Diagonal(d)`
        # but right now Zygote loses primal type information so we can't use it.
        # See https://github.com/FluxML/Zygote.jl/issues/603
        return (NO_FIELDS, ȳ.diag)
    end
    return Diagonal(d), Diagonal_pullback
end

function rrule(::typeof(diag), A::AbstractMatrix)
    function diag_pullback(ȳ)
        return (NO_FIELDS, Diagonal(ȳ))
    end
    return diag(A), diag_pullback
end
if VERSION ≥ v"1.3"
    function rrule(::typeof(diag), A::AbstractMatrix, k::Integer)
        function diag_pullback(ȳ)
            return (NO_FIELDS, diagm(size(A)..., k => ȳ), DoesNotExist())
        end
        return diag(A, k), diag_pullback
    end

    function rrule(::typeof(diagm), m::Integer, n::Integer, kv::Pair{<:Integer,<:AbstractVector}...)
        function diagm_pullback(ȳ)
            return (NO_FIELDS, DoesNotExist(), DoesNotExist(), _diagm_back.(kv, Ref(ȳ))...)
        end
        return diagm(m, n, kv...), diagm_pullback
    end
end
function rrule(::typeof(diagm), kv::Pair{<:Integer,<:AbstractVector}...)
    function diagm_pullback(ȳ)
        return (NO_FIELDS, _diagm_back.(kv, Ref(ȳ))...)
    end
    return diagm(kv...), diagm_pullback
end

function _diagm_back(p, ȳ)
    k, v = p
    d = diag(ȳ, k)[1:length(v)] # handle if diagonal was smaller than matrix
    return Composite{typeof(p)}(second = d)
end

function rrule(::typeof(*), D::Diagonal{<:Real}, V::AbstractVector{<:Real})
    function times_pullback(Ȳ)
        return (NO_FIELDS, @thunk(Diagonal(Ȳ .* V)), @thunk(D * Ȳ))
    end
    return D * V, times_pullback
end

#####
##### `Adjoint`
#####

function rrule(::Type{<:Adjoint}, A::AbstractMatrix{<:Number})
    Adjoint_pullback(ȳ::Composite) = (NO_FIELDS, ȳ.parent)
    Adjoint_pullback(ȳ::AbstractVecOrMat) = (NO_FIELDS, adjoint(ȳ))
    return Adjoint(A), Adjoint_pullback
end

function rrule(::Type{<:Adjoint}, A::AbstractVector{<:Number})
    Adjoint_pullback(ȳ::Composite) = (NO_FIELDS, vec(ȳ.parent))
    Adjoint_pullback(ȳ::AbstractMatrix) = (NO_FIELDS, vec(adjoint(ȳ)))
    return Adjoint(A), Adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractMatrix{<:Number})
    adjoint_pullback(ȳ::Composite) = (NO_FIELDS, ȳ.parent)
    adjoint_pullback(ȳ::AbstractVecOrMat) = (NO_FIELDS, adjoint(ȳ))
    return adjoint(A), adjoint_pullback
end

function rrule(::typeof(adjoint), A::AbstractVector{<:Number})
    adjoint_pullback(ȳ::Composite) = (NO_FIELDS, vec(ȳ.parent))
    adjoint_pullback(ȳ::AbstractMatrix) = (NO_FIELDS, vec(adjoint(ȳ)))
    return adjoint(A), adjoint_pullback
end

#####
##### `Transpose`
#####

function rrule(::Type{<:Transpose}, A::AbstractMatrix{<:Number})
    Transpose_pullback(ȳ::Composite) = (NO_FIELDS, ȳ.parent)
    Transpose_pullback(ȳ::AbstractVecOrMat) = (NO_FIELDS, Transpose(ȳ))
    return Transpose(A), Transpose_pullback
end

function rrule(::Type{<:Transpose}, A::AbstractVector{<:Number})
    Transpose_pullback(ȳ::Composite) = (NO_FIELDS, vec(ȳ.parent))
    Transpose_pullback(ȳ::AbstractMatrix) = (NO_FIELDS, vec(Transpose(ȳ)))
    return Transpose(A), Transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractMatrix{<:Number})
    transpose_pullback(ȳ::Composite) = (NO_FIELDS, ȳ.parent)
    transpose_pullback(ȳ::AbstractVecOrMat) = (NO_FIELDS, transpose(ȳ))
    return transpose(A), transpose_pullback
end

function rrule(::typeof(transpose), A::AbstractVector{<:Number})
    transpose_pullback(ȳ::Composite) = (NO_FIELDS, vec(ȳ.parent))
    transpose_pullback(ȳ::AbstractMatrix) = (NO_FIELDS, vec(transpose(ȳ)))
    return transpose(A), transpose_pullback
end

#####
##### Triangular matrices
#####

function rrule(::Type{<:UpperTriangular}, A::AbstractMatrix)
    function UpperTriangular_pullback(ȳ)
        return (NO_FIELDS, Matrix(ȳ))
    end
    return UpperTriangular(A), UpperTriangular_pullback
end

function rrule(::Type{<:LowerTriangular}, A::AbstractMatrix)
    function LowerTriangular_pullback(ȳ)
        return (NO_FIELDS, Matrix(ȳ))
    end
    return LowerTriangular(A), LowerTriangular_pullback
end

function rrule(::typeof(triu), A::AbstractMatrix, k::Integer)
    function triu_pullback(ȳ)
        return (NO_FIELDS, triu(ȳ, k), DoesNotExist())
    end
    return triu(A, k), triu_pullback
end
function rrule(::typeof(triu), A::AbstractMatrix)
    function triu_pullback(ȳ)
        return (NO_FIELDS, triu(ȳ))
    end
    return triu(A), triu_pullback
end

function rrule(::typeof(tril), A::AbstractMatrix, k::Integer)
    function tril_pullback(ȳ)
        return (NO_FIELDS, tril(ȳ, k), DoesNotExist())
    end
    return tril(A, k), tril_pullback
end
function rrule(::typeof(tril), A::AbstractMatrix)
    function tril_pullback(ȳ)
        return (NO_FIELDS, tril(ȳ))
    end
    return tril(A), tril_pullback
end

_diag_view(X) = view(X, diagind(X))
_diag_view(X::Diagonal) = parent(X)  #Diagonal wraps a Vector of just Diagonal elements

function rrule(::typeof(det), X::Union{Diagonal, AbstractTriangular})
    y = det(X)
    s = conj!(y ./ _diag_view(X))
    function det_pullback(ȳ)
        return (NO_FIELDS, Diagonal(ȳ .* s))
    end
    return y, det_pullback
end

function rrule(::typeof(logdet), X::Union{Diagonal, AbstractTriangular})
    y = logdet(X)
    s = conj!(one(eltype(X)) ./ _diag_view(X))
    function logdet_pullback(ȳ)
        return (NO_FIELDS, Diagonal(ȳ .* s))
    end
    return y, logdet_pullback
end


text:
function mark_kross(r::Robot)
    a = 0
    b = 0
    A = 0
    B = 0
    k = 0
    i = o
    knA = 0
    knB = 0
    verh = 0
    niz = 0
    while !isborder(r, Nord)
        move!(r, Nord)
        a = a + 1
    end
    while !isborder(r, West)
        move!(r, West)
        b = b + 1
        while !isborder(r, Nord)
            move!(r, Nord)
            a = a + 1
        end
    end
    while !isborder(r, Sud)
        move!(r, Sud)
        A = A + 1
    end
    knA = A
    while !isborder(r, Ost)
        move!(r, Ost)
        B = B + 1
    end
    knB = B
    while B > 0
        while A > 0
            if isborder(r, Nord)==true
                break
            end
            move!(r, Nord)
        end
        if isborder(r, Nord)==true
            if (A > 0)
                niz = 1
                break
            end
        end
        A = knA
        move!(r, West)
        B = B - 1
        while A > 0
            if isborder(r, Sud)==true
                break
            end
            move!(r, Sud)
        end
        if isborder(r, Sud)==true
            if (B > 0)
                verh = 1
                break
            end
        end
        A = knA
        move!(r, West)
        B = B - 1
    end
    if (verh = 1)
        while isborder(r, Sud)==true
            move!(r, West)
        end
        move!(r, Sud)
        putmarker!(r)
        while isborder(r, Ost)==true
            move!(r, Sud)
            putmarker(r)
        end
        move!(r, Ost)
        putmarker!(r)
        while isborder(r, Nord)==true
            move!(r, Ost)
            putmarker!(r)
        end
        move!(r, Nord)
        putmarker!(r)
        while isborder(r, West)==true
            move!(r, Nord)
            putmarker(r)
        end
        move!(r, West)
        putmarker!(r)
        while isborder(r, Sud)==true
            move!(r, West)
            putmarker!(r)
        end
    end
    if (niz = 1)
        while isborder(r, Nord)==true
            move!(r, West)
        end
        move!(r, Nord)
        putmarker!(r)
        while isborder(r, Ost)==true
            move!(r, Nord)
            putmarker(r)
        end
        move!(r, Ost)
        putmarker!(r)
        while isborder(r, Sud)==true
            move!(r, Ost)
            putmarker!(r)
        end
        move!(r, Sud)
        putmarker!(r)
        while isborder(r, West)==true
            move!(r, Sud)
            putmarker(r)
        end
        move!(r, West)
        putmarker!(r)
        while isborder(r, Nord)==true
            move!(r, West)
            putmarker!(r)
        end
    end
    while !isborder(r, Nord)
        move!(r, Nord)
    end
    while !isborder(r, West)
        move!(r, West)
    end
    while a > 0
        move!(r, Sud)
        a = a -1
    end
    while b > 0
        while isborder(r, Nord)==true
            move!(r, Nord)
            a = a + 1
        end
        move!(r, Ost)
        b = b - 1
    end
    while a > 0
        move!(r, Sud)
        a = a - 1
    end
end


text:
module IncrementalPruning

using LinearAlgebra
using POMDPs, POMDPPolicies
using POMDPLinter: @POMDP_require
using JuMP, GLPK
using POMDPModelTools: ordered_states, ordered_actions, ordered_observations, StateActionReward, LazyCachedSAR

import POMDPs: Solver, Policy
import POMDPs: solve, action, value, update, initialize_belief, updater
import Base: ==, hash

export
    PruneSolver

include("vanilla.jl")

end # module


text:
# Tracker.jl Compat

# The following blocks ever constructing TrackedArrays of NamedDimArrays.
# This is not strictly forbidden (thus is commented out) but is useful for debugging things
#==
Tracker.TrackedArray(::NamedDimsArray) = error("Should not make Tracked NamedDimsArray")
Tracker.TrackedArray(::Tracker.Call, ::NamedDimsArray) = error("Should not make Tracked NamedDimsArray")
Tracker.TrackedArray(::Tracker.Call, ::NamedDimsArray, ::AbstractArray) = error("Should not make Tracked NamedDimsArray")
Tracker.TrackedArray(::Tracker.Call, ::AbstractArray, ::NamedDimsArray) = error("Should not make Tracked NamedDimsArray")
==#

function Base.BroadcastStyle(::NamedDimsStyle{A}, b::Tracker.TrackedStyle) where {A}
    return NamedDimsStyle(A(), b)
end
function Base.BroadcastStyle(a::Tracker.TrackedStyle, ::NamedDimsStyle{B}) where {B}
    return NamedDimsStyle(a, B())
end

@declare_matmul(Tracker.TrackedMatrix, Tracker.TrackedVector)

function Tracker.data(nda::NamedDimsArray{L}) where {L}
    content = Tracker.data(parent(nda))
    return NamedDimsArray{L}(content)
end

function Tracker.track(c::Tracker.Call, nda::NamedDimsArray{L}) where {L}
    content = Tracker.track(c, parent(nda))
    return NamedDimsArray{L}(content)
end

for f in (:forward, :back, :back!, :grad, :istracked, :tracker)
    @eval function Tracker.$f(nda::NamedDimsArray, args...)
        return Tracker.$f(parent(nda), args...)
    end
end


text:
"""
    TagBot(;
        file="$(contractuser(default_file("github", "workflows", "TagBot.yml")))",
        destination="TagBot.yml",
        cron="0 0 * * *",
        token=Secret("GITHUB_TOKEN"),
        ssh=nothing,
        ssh_password=nothing,
        changelog=nothing,
        changelog_ignore=nothing,
        gpg=nothing,
        gpg_password=nothing,
        registry=nothing,
        branches=nothing,
        dispatch=nothing,
        dispatch_delay=nothing,
    )

Adds GitHub release support via [TagBot](https://github.com/JuliaRegistries/TagBot).

## Keyword Arguments
- `file::AbstractString`: Template file for the workflow file.
- `destination::AbstractString`: Destination of the workflow file, relative to `.github/workflows`.
- `cron::AbstractString`: Cron expression for the schedule interval.
- `token::Secret`: Name of the token secret to use.
- `ssh::Secret`: Name of the SSH private key secret to use.
- `ssh_password::Secret`: Name of the SSH key password secret to use.
- `changelog::AbstractString`: Custom changelog template.
- `changelog_ignore::Vector{<:AbstractString}`: Issue/pull request labels to ignore in the changelog.
- `gpg::Secret`: Name of the GPG private key secret to use.
- `gpg_password::Secret`: Name of the GPG private key password secret to use.
- `registry::AbstractString`: Custom registry, in the format `owner/repo`.
- `branches::Bool`: Whether not to enable the `branches` option.
- `dispatch::Bool`: Whether or not to enable the `dispatch` option.
- `dispatch_delay::Int`: Number of minutes to delay for dispatch events.
"""
@plugin struct TagBot <: FilePlugin
    file::String = default_file("github", "workflows", "TagBot.yml")
    destination::String = "TagBot.yml"
    cron::String = "0 0 * * *"
    token::Secret = Secret("GITHUB_TOKEN")
    ssh::Union{Secret, Nothing} = nothing
    ssh_password::Union{Secret, Nothing} = nothing
    changelog::Union{String, Nothing} = nothing
    changelog_ignore::Union{Vector{String}, Nothing} = nothing
    gpg::Union{Secret, Nothing} = nothing
    gpg_password::Union{Secret, Nothing} = nothing
    registry::Union{String, Nothing} = nothing
    branches::Union{Bool, Nothing} = nothing
    dispatch::Union{Bool, Nothing} = nothing
    dispatch_delay::Union{Int, Nothing} = nothing
end

source(p::TagBot) = p.file
destination(p::TagBot) = joinpath(".github", "workflows", p.destination)

function view(p::TagBot, ::Template, ::AbstractString)
    changelog = if p.changelog === nothing
        nothing
    else
        # This magic number aligns the text block just right.
        lines = map(line -> rstrip(repeat(' ', 12) * line), split(p.changelog, "\n"))
        "|\n" * join(lines, "\n")
    end
    ignore = p.changelog_ignore === nothing ? nothing : join(p.changelog_ignore, ", ")

    return Dict(
        "BRANCHES" => p.branches === nothing ? nothing : string(p.branches),
        "CHANGELOG" => changelog,
        "CHANGELOG_IGNORE" => ignore,
        "CRON" => p.cron,
        "DISPATCH" => p.dispatch === nothing ? nothing : string(p.dispatch),
        "DISPATCH_DELAY" => p.dispatch_delay,
        "GPG" => p.gpg,
        "GPG_PASSWORD" => p.gpg_password,
        "REGISTRY" => p.registry,
        "SSH" => p.ssh,
        "SSH_PASSWORD" => p.ssh_password,
        "TOKEN" => p.token,
    )
end


text:
using BenchmarkTools
using DataFrames
using VegaLite

results = only(BenchmarkTools.load(joinpath(@__DIR__, "results.json")))

function parse_problem_label(label)
    m = match(
        r"""
        ^([^-]+)
        (?:- (?:clustersize=([^-]+)))?
        - (?:ncounters=([^-]+))
        - (?:nindices=([^-]+))
        """x,
        label,
    )
    @assert m !== nothing
    problem = Symbol(m[1])
    if m[2] === nothing
        clustersize = missing
    else
        clustersize = parse(Int, m[2])
    end
    ncounters = parse(Int, m[3])
    nindices = parse(Int, m[4])
    return (; problem, clustersize, ncounters, nindices)
end
@assert parse_problem_label("uniform-ncounters=1-nindices=2") ===
        (problem = :uniform, clustersize = missing, ncounters = 1, nindices = 2)
@assert parse_problem_label("cluster-clustersize=0-ncounters=1-nindices=2") ===
        (problem = :cluster, clustersize = 0, ncounters = 1, nindices = 2)

function parsenamed(T, prefix, str)
    @assert startswith(str, prefix)
    return parse(T, str[length(prefix)+1:end])
end

begin
    rawtable = map(leaves(results)) do (labels, trial)
        plabel = parse_problem_label(labels[1])
        ntries = nbatches = missing
        if labels[2] == "seq"
            impl = :seq
            ntasks = 1
        else
            ntasks = parsenamed(Int, "ntasks=", labels[2])
            impl = Symbol(labels[3])
            if impl === :tsx_static
                ntries = parsenamed(Int, "ntries=", labels[4])
                nbatches = parsenamed(Int, "nbatches=", labels[5])
            end
        end
        return (; plabel..., ntasks, impl, ntries, nbatches, trial)
    end

    rawdf = DataFrame(rawtable)
end

begin
    df1 = select(rawdf, Not(:trial))
    df1[!, :time] = map(trial -> minimum(trial).time ./ 1000_000, rawdf.trial)
    df1
end

function saveplot(; plt_...)
    name, plt = only(plt_)
    save(joinpath(@__DIR__, "$name.png"), plt)
    save(joinpath(@__DIR__, "$name.vegalite"), plt)
end

config = let
    fontconfig = @vlfrag(labelFontSize = 20, titleFontSize = 20)
    @vlfrag(
        axis = fontconfig,
        legend = fontconfig,
        header = fontconfig,
        title = fontconfig,
        # mark = ...
    )
end

plt_raw = @vlplot(
    mark = {type = :line, point = true},
    x = {
        field = :ncounters,
        title = "length(counters)",
        axis = {format = "e"},
        scale = {type = :log},
        type = :quantitative,
    },
    y = {field = :time, title = "time [ms]", type = :quantitative},
    color = {field = :impl, title = ""},
    column = :problem,
    row = {field = :ntasks, title = "Number of parallel tasks"},
    data = df1,
    width = 200,
    height = 100,
    config = config,
)
saveplot(; plt_raw)

function add_density!(df)
    df[!, :density_value] = df.nindices ./ df.ncounters
    df[!, :Density] = map(df.nindices ./ df.ncounters) do d
        if d >= 1
            string(round(Int, d))
        else
            "1/$(round(Int, 1/d))"
        end
    end
    return df
end

let idx_seq = df1.impl .== :seq
    cols = [:problem, :ncounters, :nindices]
    global df_seq = df1[idx_seq, [cols..., :time]]
    rename!(df_seq, :time => :time_seq)
    global df_speedup =
        leftjoin(df1[.!idx_seq, :], df_seq, on = cols, validate = (false, true))
    df_speedup[!, :speedup] = df_speedup.time_seq ./ df_speedup.time
    add_density!(df_speedup)
    df_speedup
end

plt_speedup = @vlplot(
    facet = {
        column = {field = :problem, title = "Problems"},
        row = {field = :impl, title = "Implementations"},
    },
    spec = {
        layer = [
            {mark = {type = :rule}, encoding = {y = {datum = 1}}},
            {
                mark = {type = :line, point = true},
                encoding = {
                    x = {field = :ntasks, type = :quantitative, title = "#Tasks"},
                    y = {field = :speedup, type = :quantitative},
                    color = {
                        field = :Density,
                        type = :ordinal,
                        sort = {field = :density_value},
                        scale = {scheme = :inferno}
                    },
                },
            },
        ],
        width = 120,
        height = 120,
    },
    data = df_speedup,
    config = config,
)
saveplot(; plt_speedup)

let idx_uniform = df1.problem .== :uniform
    cols = [:impl, :ncounters, :nindices, :ntasks]
    global df_uniform = df1[idx_uniform, [cols..., :time]]
    rename!(df_uniform, :time => :time_uniform)
    global df_vs_uniform =
        leftjoin(df1[.!idx_uniform, :], df_uniform, on = cols, validate = (false, true))
    df_vs_uniform[!, :relative_time] = df_vs_uniform.time ./ df_vs_uniform.time_uniform
    add_density!(df_vs_uniform)
    df_vs_uniform
end

plt_vs_uniform_on_ntasks = @vlplot(
    facet = {
        column = {field = :problem, title = "Problems"},
        row = {field = :impl, title = "Implementations"},
    },
    spec = {
        layer = [
            {mark = {type = :rule}, encoding = {y = {datum = 1}}},
            {
                mark = {type = :line, point = true},
                encoding = {
                    x = {field = :ntasks, type = :quantitative, title = "#Tasks"},
                    y = {field = :relative_time, type = :quantitative},
                    color = {
                        field = :Density,
                        type = :ordinal,
                        sort = {field = :density_value},
                    },
                },
            },
        ],
        width = 120,
        height = 120,
    },
    data = df_vs_uniform,
    config = config,
)

plt_vs_uniform = @vlplot(
    facet = {
        column = {field = :problem, title = "Problems"},
        row = {field = :impl, title = "Implementations"},
    },
    spec = {
        layer = [
            {mark = {type = :rule}, encoding = {y = {datum = 1}}},
            {
                mark = {type = :line, point = true},
                encoding = {
                    x = {
                        field = :Density,
                        type = :ordinal,
                        sort = {field = :density_value},
                    },
                    y = {field = :relative_time, type = :quantitative},
                    color = {
                        field = :ntasks,
                        type = :ordinal,
                        title = "#Tasks",
                        scale = {scheme = :viridis},
                    },
                },
            },
        ],
        width = 120,
        height = 120,
    },
    data = df_vs_uniform,
    config = config,
)
saveplot(; plt_vs_uniform)


text:
using HMMs
using Brobdingnag
using Distributions
using GCTools

filename = joinpath( dirname(dirname(pathof(HMMs))), "data", "thmm3_20040630" )
hmm = open( filename, "r" ) do file
    read( file, HMMs.HMM{HMMs.GenTDist,Brob,Float64} )
end

n = 1000
path = zeros(n)
@time Distributions.rand!( hmm, path )

GCTools.reset()
@time Distributions.rand!( hmm, path )
GCTools.print()


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
using Random;
using CPUTime;
include("../utilities/peps.jl");
include("../utilities/expfam.jl");
include("samplingrules.jl");

# Run the learning algorithm, paramterised by a sampling rule
# The stopping and recommendation rules are common
# βs must be a list of thresholds *in increasing order*

function runit(seed, sr, μs, pep, βs)
    βs = collect(βs); # mutable copy
    rng = MersenneTwister(seed);
    K = length(μs);
    N = zeros(Int64, K); # counts
    S = zeros(K);        # sum of samples
    baseline = CPUtime_us();
    # pull each arm once, and also avoid boundary ̂μ for Bernoulli
    for k in 1:K
        while (N[k] == 0) || (typeof(getexpfam(pep,k)) === Bernoulli && (S[k] == 0 || S[k] == N[k]))
            S[k] += sample(rng, getexpfam(pep, k), μs[k]); N[k] += 1;
        end
    end

    state = start(sr, N);
    R = Tuple{Int64, Array{Int64,1}, UInt64}[]; # collect return values
    while true
        t = sum(N);
        hμ = S./N; # emp. estimates
        # test stopping criterion
        Z, (_, _), (istar, ξ) = glrt(pep, N, hμ);
        while Z > βs[1](t)
            popfirst!(βs);
            push!(R, (istar, copy(N), CPUtime_us()-baseline));
            if isempty(βs)
                return R;
            end
        end
        # invoke sampling rule
        k = nextsample(state, pep, istar, ξ, N, S);
        # and actually sample
        S[k] += sample(rng, getexpfam(pep, k), μs[k]);
        N[k] += 1;
        t += 1;
    end
end