text:
function run_multilevel_from_file(filename, outputdir)
    inputdict = JSON.parsefile(filename, dicttype=Dict{Symbol, Any})
    input = loadinput(MultilevelInput, inputdict[:input])
    seeds = inputdict[:seeds]
    output = inputdict[:output]
    if inputdict[:repeat] > 1
        for i in 1:inputdict[:repeat]
            rng, seed = seeds !== nothing ? (MersenneTwister(seeds[i]), seeds[i]) : (MersenneTwister(), nothing)
            population = multilevel_simulation(input, rng, Symbol(inputdict[:simtype]))
            saveoutput(population, output, outputdir, seed, i)
        end
    else
        rng = MersenneTwister(seeds)
        population = multilevel_simulation(input, rng, Symbol(inputdict[:simtype]))
        saveoutput(population, output, outputdir, seeds, 1)
    end
end

function run_multilevel_from_file(filename, outputdir, id)
    inputdict = JSON.parsefile(filename, dicttype=Dict{Symbol, Any})
    input = loadinput(MultilevelInput, inputdict[:input])
    seed = (inputdict[:seeds])[id]
    output = inputdict[:output]
    rng = MersenneTwister(seed)
    population = multilevel_simulation(input, rng, Symbol(inputdict[:simtype]))
    saveoutput(population, output, outputdir, seed, id)
end

"""
    multilevel_simulation(input::MultilevelInput, rng::AbstractRNG = MersenneTwister;
        <keyword arguments>)

    Simulate a growing population of modules, until it reaches `maxmodules`.

    Start with a single module comprised of a single cell that grows via a branching process
    to size `input.modulesize`, then switches to a Moran process.
    New modules are created with rate `input.branchrate`, by sampling cells from the parent
    module. Return output as a Population.

    - if `simtype == :normal`, simulation is
    implemented by a Gillespie algorithm and runs until the number of modules exceeds
    input.maxmodules or time exceeds input.maxtime.

    - if `simtype == :fixedtime`, a list of modules is created (initially of length 1)
    and each module is simulated independently until input.maxtime is reached. New modules
    are added to the end of the list. This implementation is faster, but only works if we
    end the simulation at a fixed time (not fixed size). If the list of modules is bigger
    than input.maxmodules and error is thrown.

"""
function multilevel_simulation(input::MultilevelInput, rng::AbstractRNG=Random.GLOBAL_RNG,
    simtype=:normal)

    #if the population dies out we start a new simulation
    while true
        populationtracker = initialize_population(
            input.modulesize,
            clonalmutations=0
        )
        if simtype == :normal || simtype == "normal"
            populationtracker =
                simulate!(
                    populationtracker,
                    input.maxtime,
                    input.maxmodules,
                    input.b,
                    input.d,
                    input.bdrate,
                    input.branchrate,
                    input.modulesize,
                    input.branchinitsize,
                    rng
                )
        elseif simtype == :fixedtime || simtype == "fixedtime"
            populationtracker =
                simulatefixedtime!(
                    populationtracker,
                    input.maxtime,
                    input.maxmodules,
                    input.b,
                    input.d,
                    input.bdrate,
                    input.branchrate,
                    input.modulesize,
                    input.branchinitsize,
                    rng
                )
        else
            error("incorrect value for simtype")
        end
        if length(populationtracker) != 0
            break
        end
    end
    populationtracker =
        processresults!(populationtracker, input.μ, input.clonalmutations, rng)
    return Population(input, populationtracker)
end

"""
    simulate!(populationtracker, maxtime, maxmodules, b, d, bdrate, branchrate,
        modulesize, branchinitsize, rng)

Run a single multilevel simulation using the Gillespie algorithm, with the
`populationtracker` giving the initial state. Simulation runs until the module population
size reaches `maxmodules` or the age of the population reaches `maxtime`.

"""
function simulate!(populationtracker, maxtime, maxmodules, b, d, bdrate, branchrate,
    modulesize, branchinitsize, rng)

    mutID = getmutID(populationtracker[1].cells) #get id of next mutation
    t = age(populationtracker)
    while t < maxtime && length(populationtracker) < maxmodules
        populationtracker, mutID, t =
            update_population!(populationtracker, b, d, bdrate, branchrate, modulesize,
                branchinitsize, mutID, t, maxtime, rng)
        #returns empty list of modules if population dies out
        if length(populationtracker) == 0
            return populationtracker
        end
    end
    return populationtracker
end

"""
    update_population!(populationtracker, b, d, bdrate, branchrate, modulesize,
        branchinitsize, mutID, t, maxtime, rng)

Perform a single Gillespie step to advance the simulation of a multilevel population, and
increase the time t → t + Δt.

Calculates the transition rates for a Moran, birth, death or module branching event; then
selects from these transitions with probability proportional to rate. Time is increased by
`Δt ~ Exp(sum(transitionrates))`. If new time exceeds `maxtime` do not perform any transition.
"""

function update_population!(populationtracker, b, d, bdrate, branchrate, modulesize,
    branchinitsize, mutID, t, maxtime, rng)

    transitionrates = get_transitionrates(populationtracker, b, d,
        bdrate, branchrate, modulesize)
    t += exptime(rng, sum(transitionrates))
    #only update the population if t < maxtime
    if t < maxtime
        #choose transition type: 1=moran, 2=birth, 3=death, 4=branch
        transitionid = sample(
            rng,
            1:4,
            ProbabilityWeights(transitionrates ./ sum(transitionrates))
        )
        populationtracker, mutID = transition!(
            populationtracker,
            transitionid,
            modulesize,
            branchinitsize,
            mutID,
            t,
            1,
            rng
        )
    end
    return populationtracker, mutID, t
end

"""
    transition!(populationtracker, transitionid, modulesize, branchinitsize, mutID, t,
        μ, rng)

Perform a single transition step on `populationtracker`, determined by `transitionid`.
"""
function transition!(populationtracker, transitionid, modulesize, branchinitsize, mutID, t,
    μ, rng)

    if transitionid == 1
        _, mutID = moranupdate!(populationtracker, modulesize, mutID, t, μ, rng)
    elseif transitionid == 2
        _, mutID = birthupdate!(populationtracker, modulesize, mutID, t, μ, rng)
    elseif transitionid == 3
        deathupdate!(populationtracker, modulesize, t, rng)
    elseif transitionid == 4
        modulebranchingupdate!(populationtracker, modulesize, branchinitsize, t, rng)
    end
    return populationtracker, mutID
end

"""
    moranupdate!(populationtracker, modulesize, mutID, t, μ, rng)

Selects a homeostatic module uniformly at random to undergo a single Moran update. From
that module one cell divides and one cell dies.
"""
function moranupdate!(populationtracker, modulesize, mutID, t, μ, rng)
    moduletracker, parentcell, deadcell =
        choose_homeostaticmodule_cells(populationtracker, modulesize, rng)
    _, mutID = celldivision!(moduletracker, parentcell, mutID, μ, rng)
    celldeath!(moduletracker, deadcell)
    updatemodulehistory!(moduletracker, 0, t)
    return populationtracker, mutID
end

"""
    birthupdate!(populationtracker, modulesize, mutID, t, μ, rng)

Selects a cell uniformly at random from all cells in non-homeostatic modules to divide.
"""
function birthupdate!(populationtracker, modulesize, mutID, t, μ, rng)
    moduletracker, parentcell =
        choose_growingmodule_cell(populationtracker, modulesize, rng)
    _, mutID = celldivision!(moduletracker, parentcell, mutID, μ, rng, fixedmu=true)
    updatemodulehistory!(moduletracker, 1, t)
    return populationtracker, mutID
end

"""
    deathupdate!(populationtracker, modulesize, t, rng)

Selects a cell uniformly at random from all cells in non-homeostatic modules to die. If
cell death results in an empty module, remove that module from the population.
"""
function deathupdate!(populationtracker, modulesize, t, rng)
    moduletracker, deadcell =
        choose_growingmodule_cell(populationtracker, modulesize, rng)
    celldeath!(moduletracker, deadcell)
    updatemodulehistory!(moduletracker, -1, t)
    if length(moduletracker) == 0
        moduledeath!(populationtracker, moduletracker)
    end
    return populationtracker
end

"""
    modulebranchingupdate!(populationtracker, modulesize, branchinitsize, t, rng)
Select a homeostatic module, uniformly at random, to undergo branching. Cells are sampled
(number of cells giving by `branchinitsize`) from the parent module to form a new module,
which is added to `populationtracker`
"""
function modulebranchingupdate!(populationtracker, modulesize, branchinitsize, t, rng)
    parentmodule = choose_homeostaticmodule(populationtracker, modulesize, rng)
    _, newmodule =
        sample_new_module!(parentmodule, length(populationtracker) + 1,
            branchinitsize, t, rng)
    push!(populationtracker, newmodule)
    return populationtracker
end

"""
    updatemodulehistory!(moduletracker, ΔN, newtime)
Adds new population size (old size + `ΔN`) and `newtime` to the end of `moduletracker.Nvec`
and `moduletracker.tvec`.

"""
function updatemodulehistory!(moduletracker, ΔN, newtime)
    push!(moduletracker.Nvec, moduletracker.Nvec[end] + ΔN)
    push!(moduletracker.tvec, newtime)
    return moduletracker
end

"""
    choose_homeostaticmodule(populationtracker, maxmodulesize, rng::AbstractRNG)
Select a homeostatic module, i.e. module of size `maxmodulesize`, uniformly at random.
"""
function choose_homeostaticmodule(populationtracker, maxmodulesize, rng::AbstractRNG)
    homeostatic_modules = filter(x -> length(x) == maxmodulesize, populationtracker)
    return rand(rng, homeostatic_modules)
end

"""
    choose_homeostaticmodule_cells(populationtracker, maxmodulesize, rng::AbstractRNG)
Select a homeostatic module and two cells from the module, uniformly at random.
"""
function choose_homeostaticmodule_cells(populationtracker, maxmodulesize, rng::AbstractRNG)
    chosenmodule = choose_homeostaticmodule(populationtracker, maxmodulesize, rng)
    return chosenmodule, rand(rng, 1:maxmodulesize), rand(rng, 1:maxmodulesize)
end

"""
    choose_growingmodule_cells(populationtracker, maxmodulesize, rng::AbstractRNG)
Select a cell unifor´mly at random from all cells in growing (non-homeostatic) modules, and
return the module and cell.
"""
function choose_growingmodule_cell(populationtracker, maxmodulesize, rng::AbstractRNG)
    modulesizes = map(length, populationtracker)
    modules = populationtracker[modulesizes .< maxmodulesize]
    modulesizes = modulesizes[modulesizes .< maxmodulesize]
    chosenmodule = sample(rng, modules, ProbabilityWeights(modulesizes ./ sum(modulesizes)))
    chosencell = rand(rng, 1:length(chosenmodule))
    return chosenmodule, chosencell
end

"""
    get_transitionrates(populationtracker, b, d, bdrate, branchrate, modulesize)
Compute the rates for moran update, cell division, death and module branching and return as
a Vector.
"""
function get_transitionrates(populationtracker, b, d, bdrate, branchrate, modulesize)
    rates = zeros(Float64, 4)
    number_homeostatic_modules = 0
    for (i, moduletracker) in enumerate(populationtracker)
        N = length(moduletracker)
        if N < modulesize
            rates[2] += N * b
            rates[3] += N * d
        elseif N == modulesize
            number_homeostatic_modules += 1
        else
            error("module size exceeds homeostatic size")
        end
    end
    rates[1] = number_homeostatic_modules * bdrate * modulesize
    rates[4] = number_homeostatic_modules * branchrate
    return rates
end

"""
    moduledeath!(populationtracker, moduletracker)
Remove the given `moduletracker` from `populationtracker`.
"""
function moduledeath!(populationtracker, moduletracker)
    return deleteat!(populationtracker, findall(x->x==moduletracker, populationtracker))
end

"""
    simulatefixedtime!(populationtracker, maxtime, maxmodules, b, d, bdrate, branchrate,
        modulesize, branchinitsize, rng)

Run a single multilevel simulation until `maxtime`, with the `populationtracker` giving the
initial state. Faster than `simulate!`, but cannot run to fixed population size.
"""
function simulatefixedtime!(populationtracker, maxtime, maxmodules, b, d, bdrate, branchrate,
    modulesize, branchinitsize, rng)

    for moduletracker in populationtracker
        check_module_number(populationtracker, maxmodules)
        while true
            moduletracker, newmoduletracker =
                module_simulate_to_branching!(
                    moduletracker,
                    maxtime,
                    b,
                    d,
                    bdrate,
                    branchrate,
                    modulesize,
                    branchinitsize,
                    length(populationtracker) + 1,
                    rng
                )

            if newmoduletracker !== nothing
                push!(populationtracker, newmoduletracker)
            elseif length(moduletracker) == 0
                moduledeath!(populationtracker, moduletracker)
                break
            else
                break
            end
        end
    end
    return populationtracker
end

"""
    check_module_number(populationtracker, maxmodules)
Throws an error if the number of modules in `populationtracker` exceeds `maxmodules`.
"""
function check_module_number(populationtracker, maxmodules)
    numbermodules = length(populationtracker)
    if numbermodules >= maxmodules
        error("population size exceeds maxmodules: ", numbermodules, " > ", maxmodules)
    end
end

"""
    module_simulate_to_branching(moduletracker, maxtime, b, d, bdrate, branchrate,
        modulesize, branchinitsize, newmoduleid, rng)

    Simulates moduletracker dynamics until a module branching event occurs, or maxtime is
    reached.

    If moduletracker is smaller than `modulesize` it grows according to a branching
    process, with birthrate `b` and deathrate `d`, until it reaaches `modukesize`. Once it
    reaches `modulesize` the dynamic switches to a Moran process with event rate `bdrate`.

"""

function module_simulate_to_branching!(moduletracker, maxtime, b, d, bdrate,
    branchrate, modulesize, branchinitsize, newmoduleid, rng)

    #branching process until module reaches homeostatic size
    if length(moduletracker) < modulesize
        moduletracker =
            branchingprocess!(moduletracker, b, d, modulesize, 1,
                rng, numclones=0, fixedmu=true, maxclonesize=Inf, tmax=maxtime)
        if length(moduletracker) == 0
            return moduletracker, nothing
        end
    end

    #get time of next branching event
    branchtime = 1 / branchrate .* exptime(rng) + moduletracker.tvec[end]

    #if module branching will occur before maxime is reached, simulate moran process until
    #branchtime and sample new module
    if branchtime < maxtime
        moduletracker =
            moranprocess!(moduletracker, bdrate, branchtime, 1, rng, numclones=0,
                fixedmu=true)

        moduletracker, newmoduletracker =
            sample_new_module!(moduletracker, newmoduleid, branchinitsize, branchtime, rng)

        return moduletracker, newmoduletracker

    #if module branching will occur after maxime is reached, simulate moran process until
    #maxtime
    else
        moduletracker =
            moranprocess!(moduletracker, bdrate, maxtime, 1, rng, numclones=0,
                fixedmu=true)
        return moduletracker, nothing
    end
end

"""
    sample_new_module(moduletracker, newmoduleid, branchinitsize, branchtime,
        rng::AbstractRNG)

Sample cells, uniformly at random, from `moduletracker` to form `newmoduletracker`. The
number of cells is given by `branchinitsize` and the newmodule is given initial time
`branchtime`.

"""
function sample_new_module!(moduletracker, newmoduleid, branchinitsize, branchtime,
    rng::AbstractRNG)

    sampleids = sample(rng, 1:length(moduletracker.cells), branchinitsize, replace=false)
    newmoduletracker =
        initializesim_from_cells(moduletracker.cells[sampleids], moduletracker.subclones,
            newmoduleid, moduletracker.id, inittime=branchtime)
    moduletracker = celldeath!(moduletracker, sampleids)
    push!(moduletracker.Nvec, moduletracker.Nvec[end] - branchinitsize)
    push!(moduletracker.tvec, branchtime)

    return moduletracker, newmoduletracker
end

"""
    initialize_population([modulesize=nothing]; clonalmutations=0)

Create a Vector{ModuleTracker} of length 1 that contains a single module comprising a single
cell at time t=0.
"""
function initialize_population(modulesize=nothing; clonalmutations=0)

    #population is defined by ModuleTracker vector with each entry coprresponding to
    #an individual module
    population = ModuleTracker[]

    #initialise population with a single module conisting of a single cell
    initialmodule = initializesim_branching(
        modulesize,
        clonalmutations=clonalmutations
    )
    push!(population, initialmodule)

    return population
end


text:
## stress test clenshaw evaluation against naive evaluate of Chebyshev polynomials.

@everywhere using Distributed
@everywhere using SharedArrays

@everywhere import Random
@everywhere import LinearAlgebra
@everywhere import Printf

@everywhere include("../src/minimax/interpolators.jl")

@everywhere include("../src/misc/utilities.jl")

Random.seed!(25)

fig_num = 1

N = 1000 # number of positions used in each test.
N_tests = 1 # number of tests.

function stresstestget𝑤scaledfast(N_tests, N)
    discrepancy = SharedArray{Float64}(N_tests)

    @sync @distributed for n = 1:N_tests
        X = collect( convert(BigFloat,rand()*2-1) for i = 1:N )

        w = get𝑤scaled(X)
        w2 = get𝑤scaledfast(X)

        discrepancy[n] = LinearAlgebra.norm(w-w2)
    end

    return convert(Vector{Float64},discrepancy)
end

discrepancy = stresstestget𝑤scaledfast(N_tests, N)
println("norm(discrepancy) = ", LinearAlgebra.norm(discrepancy))
println()

#### speed comparison.
X = collect( convert(BigFloat,rand()*2-1) for i = 1:N )

w = get𝑤scaled(X)
w2 = get𝑤scaledfast(X)

d = LinearAlgebra.norm(w-w2)

println("discrepancy between naive and parallel: ", d)

ext: jl

max_stars_repo_name: RoyCCWang/FIRRemez

max_issues_repo_name: RoyCCWang/FIRRemez

max_forks_repo_name: RoyCCWang/FIRRemez


text:
struct Rastrigin{B<:Vector{Float64}, P<:Real} <: AbstractDensity
    bounds::B
    a::P
    b::P
    dist
end

function Rastrigin(;bounds=[-3.,3.],a=30,b=10)
    integral_val = -a*bounds[1] + b*bounds[1] + (bounds[1]^3)/3 + a*bounds[2] - b*bounds[2] - (bounds[2]^3)/3 + (b*(-sin(2*bounds[1]*pi) + sin(2*bounds[2]*pi)))/(2*pi)
    dist = x-> (1/integral_val) * (a - b - x^2 + b*cos(2*pi*x))
    Rastrigin(bounds,a,b,dist)
end

BAT.eval_logval_unchecked(density::Rastrigin,v::Any) = log(broadcast(density.dist,v))

"""
    Rastrigin <: AbstractDensity
"""

struct SineSquared{B<:Vector{Float64}} <: AbstractDensity
    bounds::B
    dist
end

function SineSquared(;bounds=[-0.,25.])
    integral_val = 1/40*(-4*(bounds[1])^5 + 10*(2*(bounds[1])^2 - 3)*(bounds[1])*cos(2*(bounds[1])) + 5*(2*(bounds[1])^4 - 6*(bounds[1])^2 + 3)*sin(2*(bounds[1])) + 4*(bounds[2])^5 + 10*(bounds[2])*(3 - 2*(bounds[2])^2)*cos(2*(bounds[2])) - 5*(2*(bounds[2])^4 - 6*(bounds[2])^2 + 3)*sin(2*(bounds[2])))
    dist = x-> (1/integral_val) * (x^4*sin(x)^2)
    SineSquared(bounds,dist)
end

BAT.eval_logval_unchecked(density::SineSquared,v::Any) = log(broadcast(density.dist,v))

"""
    SineSquared <: AbstractDensity
"""

struct HoelderTable{B<:Vector{Float64}} <: AbstractDensity
    bounds::B
    dist
end
export HoelderTable

function HoelderTable()
    bounds = [-10.,10.]
    integral_val = (-18*exp(1)^2*π^2 - 12*exp(1)^3*π^2 + exp(1)*(500 + 488*π^2) - 6*exp(1)^(10/π)*π*(π*cos(10) - sin(10)))/(3*exp(1)*(1 + π^2))
    dist = x-> (1/integral_val) * (10-((x^2)/20) - abs( sin(x)*exp(abs(1-((sqrt(x^2))/(pi))))) )
    HoelderTable(bounds,dist)
end

BAT.eval_logval_unchecked(density::HoelderTable,v::Any) = log(broadcast(density.dist,v))

"""
    HoelderTable <: AbstractDensity
    Fixed bounds, integral not easily calculable
"""


text:
# Copyright 2018-2019 Eric S. Tellez <eric.tellez@infotec.mx>
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
export KernelClassifier, KConfigurationSpace, predict, predict_one

using KernelMethods
import SimilaritySearch: l2_distance, angle_distance

struct KConfigurationSpace
    normalize
    distances
    kdistances
    sampling
    kernels
    reftypes
    classifiers

    function KConfigurationSpace(;
        normalize=[true, false],
        distances=[l2_distance, cosine_distance],
        kdistances=[l2_distance, cosine_distance],
        sampling=vcat(
            [(method=fftraversal, stop=x) for x in (sqrt_criterion, log_criterion, change_criterion)],
            [(method=dnet, kfun=x) for x in (log2, x -> min(x, log2(x)^2))]
        ),
        #kernels=[linear_kernel, gaussian_kernel, sigmoid_kernel, cauchy_kernel, tanh_kernel],
        kernels=[linear_kernel, gaussian_kernel, tanh_kernel],
        #reftypes=[:centroids, :centers],
        reftypes=[:centroids, :centers],
        #classifiers=[NearNeighborClassifier, NaiveBayesClassifier]
        classifiers=[NearNeighborClassifier]
        )
        new(normalize, distances, kdistances, sampling, kernels, reftypes, classifiers)
    end

end

struct KConfiguration
    normalize
    dist
    kdist
    kernel
    net
    reftype
    classifier
end

function randconf(space::KConfigurationSpace)
    kdist = rand(space.kdistances)
    normalize = (kdist in (cosine_distance, angle_distance)) || rand(space.normalize)

    c = KConfiguration(
        normalize,
        rand(space.distances),
        kdist,
        rand(space.kernels),
        rand(space.sampling),
        rand(space.reftypes),
        rand(space.classifiers)
    )

    c
end

function randconf(space::KConfigurationSpace, num::Integer)
    [randconf(space) for i in 1:num]
end

struct KernelClassifierType{ItemType}
    kernel
    refs::Vector{ItemType}
    classifier
    conf::KConfiguration
end

"""
Searches for a competitive configuration in a parameter space using random search
    """
    function KernelClassifier(X, y;
        folds=3,
        score=recall,
        size=32,
        ensemble_size=3,
        space=KConfigurationSpace()
        )

        bestlist = []
        tabu = Set()
        dtype = typeof(X[1])

        for conf in randconf(space, size)
            if conf in tabu
                continue
            end

            @info "testing configuration $(conf), data-type $(typeof(X))"
            push!(tabu, conf)
            dist = conf.dist
            kdist = conf.kdist
            refs = Vector{dtype}()
            dmax = 0.0

            if conf.kernel in (cauchy_kernel, gaussian_kernel, sigmoid_kernel, tanh_kernel)
                kernel = conf.kernel(dist, dmax/2)
            elseif conf.kernel == linear_kernel
                kernel = conf.kernel(dist)
            else
                kernel = conf.kernel
            end

            if conf.net.method == fftraversal
                function pushcenter1(c, _dmax)
                    push!(refs, X[c])
                    dmax = _dmax
                end

                fftraversal(pushcenter1, dist, X, conf.net.stop())
                # after fftraversal refs is populated
                R = fit(Sequential, refs)
                @info "computing kmap, conf: $conf"
                if conf.reftype == :centroids
                    a = [centroid!(X[plist]) for plist in invindex(dist, X, R) if length(plist) > 0]
                    M = kmap(X, kernel, a)
                else
                    M = kmap(X, kernel, refs)
                end
            elseif conf.net.method == dnet
                function pushcenter2(c, dmaxlist)
                    if conf.reftype == :centroids
                        a = vcat([X[c]], X[[p.objID for p in dmaxlist]]) |> centroid!
                        push!(refs, a)
                    else
                        push!(refs, X[c])
                    end

                    dmax += last(dmaxlist).dist
                end

                k = conf.net.kfun(length(X)) |> ceil |> Int
                dnet(pushcenter2, dist, X, k)
                if k == length(refs)
                    @info "$k != $(length(refs))"
                    @info refs
                    error("incorrect number of references, $conf")
                end
                @info "computing kmap, conf: $conf"
                M = kmap(X, kernel, refs)
                dmax /= length(refs)
            end

            if conf.normalize
                for m in M
                    normalize!(m)
                end
            end

            @info "creating and optimizing classifier, conf: $conf"
            if conf.classifier == NearNeighborClassifier
                classifier = NearNeighborClassifier(kdist, M, y)
                best = optimize!(classifier, score, folds=folds)[1]
            else
                classifier = NaiveBayesClassifier(M, y)
                best = optimize!(classifier, M, y, score, folds=folds)[1]
            end

            model = KernelClassifierType(kernel, refs, classifier, conf)
            push!(bestlist, (best[1], model))
            @info "score: $(best[1]), conf: $conf"
            sort!(bestlist, by=x->-x[1])

            if length(bestlist) > ensemble_size
                bestlist = bestlist[1:ensemble_size]
            end
        end

        @info "final scores: ", [b[1] for b in bestlist]
        # @show [b[1] for b in bestlist]
        [b[2] for b in bestlist]
    end

    function predict(kmodel::AbstractVector{KernelClassifierType{ItemType}}, vector) where {ItemType}
        [predict_one(kmodel, x) for x in vector]
    end

    function predict_one(kmodel::AbstractVector{KernelClassifierType{ItemType}}, x) where {ItemType}
        C = Dict()
        for m in kmodel
            label = predict_one(m, x)
            C[label] = get(C, label, 0) + 1
        end

        counter = [(c, label) for (label, c) in C]
        sort!(counter, by=x->-x[1])

        # @show counter
        counter[1][end]
    end

    function predict_one(kmodel::KernelClassifierType{ItemType}, x) where {ItemType}
        kernel = kmodel.kernel
        refs = kmodel.refs

        vec = Vector{Float64}(undef, length(refs))
        for i in 1:length(refs)
            vec[i] = kernel(x, refs[i])
        end

        if kmodel.conf.normalize
            normalize!(vec)
        end

        predict_one(kmodel.classifier, vec)
    end


text:
"""
    parallactic_angles(filelist; latitude=19.82636)

Calculate the parallactic angles (in degrees) for the given FITS files in the middle of the exposure. Uses the given
`latitude` which defaults to the Keck II site.
"""
parallactic_angles(filelist; kwargs...) = map(f -> parallactic_angle(read_header(f); kwargs...), filelist)

parallactic_angles(filelist::DataFrame; kwargs...) = map(r -> parallactic_angle(r; kwargs...), eachrow(filelist))

function parallactic_angle(row; latitude=19.82636)
    # compute the hour angle at the middle of the frame
    exp_start = Time(row["EXPSTART"])
    exp_stop = Time(row["EXPSTOP"])
    exp_mid = Dates.value((exp_stop - exp_start) / 2)
    # convert exp_mid from nanoseconds of time to degrees
    Δha = exp_mid / 1e9 / 24 / 10
    ha_mid = row["HA"] + Δha

    # precess the star coordinates to the appropriate epoch
    ra = deg2rad(row["RA"])
    dec = deg2rad(row["DEC"])
    coord = FK5Coords{2000}(ra, dec)
    obs_epoch = Date(row["DATE-OBS"])
    yr = year(obs_epoch) + (dayofyear(obs_epoch) - 1) / daysinyear(obs_epoch)
    coord_curr = convert(FK5Coords{yr}, coord)

    # derive true PA at the middle of the frame
    ra = coord_curr.ra
    dec = coord_curr.dec
    PA = -rad2deg(atan(-sind(ha_mid), cos(dec) * tand(latitude) - sin(dec) * cosd(ha_mid)))
    # apply instrumental effects
    PA += row["ROTPOSN"] - row["INSTANGL"] + (row["PARANG"] - row["PARANTEL"])

    return PA
end

max_stars_repo_name: mileslucas/sirius-b

max_issues_repo_name: mileslucas/sirius-b

max_forks_repo_name: mileslucas/sirius-b


text:
# Conjuagtes for
#
#	Gamma - Exponential
#

function posterior_canon(prior::Gamma, ss::ExponentialStats)
	return Gamma(shape(prior) + ss.sw, 1.0 / (rate(prior) + ss.sx))
end

complete(G::Type{Exponential}, pri::Gamma, θ::Float64) = Exponential(1.0 / θ)


text:
# This is a convenience function for libraries like MaBLAS and PaddedMatrices to use explicitly, and for others to force more type inference in precompilation.
let GEMMLOOPSET = LoopVectorization.LoopSet(
    :(for m ∈ 1:size(A,1), n ∈ 1:size(B,2)
          Cₘₙ = zero(eltype(C))
          for k ∈ 1:size(A,2)
              Cₘₙ += A[m,k] * B[k,n]
          end
          C[m,n] += Cₘₙ
      end)
    );
    order = LoopVectorization.choose_order(GEMMLOOPSET)
    mr = order[5]
    nr = last(order)
    @eval const mᵣ = $mr
    @eval const nᵣ = $nr
end


text:
using SparseArrays: sparse
using LinearAlgebra: cholesky, ldiv!
# import Convex
# import ECOS

"""Fit Ax = b using Huber loss in ADMM (alternatinve directions method of multiplers)
Adapted from https://web.stanford.edu/~boyd/papers/admm/huber/huber_fit.html
"""
# TODO: differentiate this from the regularization alpha
function huber_fit(
    A,
    b,
    rho = 1.0,
    alpha = 1.0;
    lu_tuple = nothing,
    quiet = true,
    max_iter = 1000,
    abstol = 1e-4,
    reltol = 1e-2,
)
    m, n = size(A)
    Atb = A' * b

    x = zeros(n)
    z = zeros(m)
    zold = zeros(m)
    u = zeros(m)
    # If you prefactor the A matrix
    if !isnothing(lu_tuple)
        L, U = lu_tuple
    else
        L, U = factor(A)
    end

    idx = 0

    kappa = 1 + 1 / rho  # For shrinkage
    r_factor = rho / (1 + rho)
    s_factor = 1 / (1 + rho)
    while idx < max_iter
        # x update
        zold .= z
        q = Atb .+ A' * (z - u)
        x .= U \ (L \ q)

        # x update with relaxation
        Ax_hat = alpha .* A * x + (1 - alpha) * (z .+ b)
        tmp = Ax_hat - b + u
        z .= (r_factor .* tmp) + (s_factor .* shrinkage(tmp, kappa))
        u += (Ax_hat - z - b)

        # Stopping check
        r_norm = norm(A * x - z - b)

        s_norm = norm(BLAS.gemv('T', -rho, A, (z - zold)))
        # ^^ equivalent to this below,  but faster
        # s_norm  = norm(-rho * A' * (z - zold))

        eps_pri = sqrt(n) * abstol + reltol * maximum([norm(A * x), norm(-z), norm(b)])
        eps_dual = sqrt(n) * abstol + reltol * norm(rho * u)

        if !quiet
            @show r_norm, eps_pri, s_norm, eps_dual, objective(z)
        end

        if r_norm < eps_pri && s_norm < eps_dual
            if !quiet
                println("Number of iterations to converge: $idx")
            end
            return x
        else
            idx += 1
        end

    end
    println("Caution: problem did not converge within $max_iter iterations to within $abstol, $reltol (alpha: $alpha, rho: $rho, sizes: $(size(A)), $(size(b)))")
    return x
end

"""Evaluate the huber loss with a bandwidth of `M`"""
function huber(x; M = 1)
    y = max.(x, 0)
    z = min.(y, M)
    return z .* (2 .* y .- z)
end

objective(z) = sum(huber(z)) / 2

# Sparse factoring
function factor(A)
    ch = cholesky(A' * A)
    return sparse(ch.L), sparse(ch.U)
end

pos(x) = max.(x, 0)

# Soft threshold operator (section 4.4.3)
shrinkage(x, kappa) = pos(1 .- kappa ./ abs.(x)) .* x

# ### Alternate: Solve the L1 fitting problem using Convex
# # Note: this had a memory leak earlier when run on many pixels,
# # possibly from the ECOS library
# function l1_fit(A, b, x::Union{Convex.Variable, Nothing}=nothing, quiet=true)
#     solver = ECOS.ECOSSolver(verbose=0)
#     if isnothing(x)
#         x = Convex.Variable(size(A, 2))
#     end
#
#     problem = Convex.minimize(norm(A*x - b, 1))
#     Convex.solve!(problem, solver)
#     return x.value
# end

#### IRLS algorithm functions: ####
l1_objective(A, x, b) = sum(abs.(A * x - b))

"""Iteratively reweighted least squares (IRLS), used to solve L1 minimization
Source: https://en.wikipedia.org/wiki/Iteratively_reweighted_least_squares"""
function irls(
    A::AbstractArray{<:AbstractFloat,2},
    b::AbstractArray{<:AbstractFloat,1};
    p::Int = 1,
    iters = 50,
)
    M, N = size(A)

    # Use Float64 to avoid roundoff NaNs
    x = Array{Float64,1}(undef, N)
    W = zeros(Float64, (size(A, 1), size(A, 1)))

    ep = sqrt(eps(eltype(x)))
    W .= diagm(0 => (abs.(b - A * x) .+ ep).^(p - 2))

    for ii = 1:iters
        x .= (A' * W * A) \ (A' * W * b)
        W .= diagm(0 => (abs.(b - A * x) .+ ep).^(p - 2))
    end
    # println("objective: ", l1_objective(A, x, b))
    return Float32.(x)
end

# # Function compatible for run_sbas
# function invert_pixel(pixel::AbstractArray{T, 1}, B::AbstractArray{T,2}; iters=50, p=1) where {T<:AbstractFloat}
#     return irls(B, pixel, iters=iters, p=p)
# end


text:
# Autogenerated wrapper script for Xyce_jll for i686-linux-musl-cxx11
export Xyce, libxyce

using Trilinos_jll
using SuiteSparse_jll
using OpenBLAS32_jll
using FFTW_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Xyce")
JLLWrappers.@declare_executable_product(Xyce)
JLLWrappers.@declare_library_product(libxyce, "libxyce.so.0")
function __init__()
    JLLWrappers.@generate_init_header(Trilinos_jll, SuiteSparse_jll, OpenBLAS32_jll, FFTW_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_executable_product(
        Xyce,
        "bin/Xyce",
    )

    JLLWrappers.@init_library_product(
        libxyce,
        "lib/libxyce.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

plot(x=rand(100), y=rand(100),
     color=rand([colorant"red",colorant"blue",colorant"green"], 100),
     shape=rand([Shape.circle,Shape.square,Shape.utriangle], 100),
     size=rand([1mm, 2mm, 3mm], 100),
     Geom.point)


text:
function tri_zero_spurious_modes!(A::AbstractMatrix)
    M, N = size(A)
    @inbounds for j = 2:N
        @simd for i = M-j+2:M
            A[i,j] = 0
        end
    end
    A
end

function trirand(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = rand(T)
        end
    end
    A
end

function trirandn(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = randn(T)
        end
    end
    A
end

function triones(::Type{T}, m::Int, n::Int) where T
    A = zeros(T, m, n)
    for j = 1:n
        for i = 1:m+1-j
            A[i,j] = one(T)
        end
    end
    A
end

trizeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, n)

"""
Pointwise evaluation of triangular harmonic:

```math
\\tilde{P}_{\\ell,m}^{(\\alpha,\\beta,\\gamma)}(x,y).
```
"""
trievaluate(x, y, L, M, α, β, γ) = trievaluate(x, L, M, α, β, γ)*trievaluate(x, y, M, β, γ)

function trievaluate(x::Number, L::Integer, M::Integer, α::Number, β::Number, γ::Number)

end

function trievaluate(x::Number, y::Number, M::Integer, β::Number, γ::Number)

end


text:
@defcomp CO2Cycle begin
    # Permafrost outputs
    permte0_permafrostemissions0 = Parameter(unit="Mtonne")
    permte_permafrostemissions = Parameter(index=[time], unit="Mtonne")

    e_globalCO2emissions = Parameter(index=[time], unit="Mtonne/year")
    e0_globalCO2emissions = Parameter(unit="Mtonne/year", default=41223.85968577856)

    pic_preindustconcCO2 = Parameter(unit="ppbv", default=278000.)
    exc_excessconcCO2 = Variable(unit="ppbv")
    c0_CO2concbaseyr = Parameter(unit="ppbv", default=400859.5833333334)
    re_remainCO2 = Variable(index=[time], unit="Mtonne")
    re_remainCO2base = Variable(unit="Mtonne")
    renoccf_remainCO2wocc = Variable(index=[time], unit="Mtonne")
    air_CO2fractioninatm = Parameter(unit="%", default=62.00)
    stay_fractionCO2emissionsinatm = Parameter(default=0.2341802168612297)# percent co2 stay but not in percent
    ccf_CO2feedback = Parameter(unit="%/degreeC", default=0.0)
    ccfmax_maxCO2feedback = Parameter(unit="%", default=20.0)
    ce_0_basecumCO2emissions = Parameter(unit="Mtonne", default=2040000.)
    conoccf_concentrationCO2wocc = Variable(index=[time], unit="ppbv")
    y_year = Parameter(index=[time], unit="year")
    y_year_0 = Parameter(unit="year")
    res_CO2atmlifetime = Parameter(unit="year", default=73.3333333333333)
    den_CO2density = Parameter(unit="Mtonne/ppbv", default=7.8)
    rt_g0_baseglobaltemp = Parameter(unit="degreeC", default=0.9461666666666667)
    rt_g_globaltemperature = Parameter(index=[time], unit="degreeC")

    # Parameters for components of CO2
    thist_timescaleco2hist = Parameter(unit="year", default=49.59054463320648)
    corrf_correctionfactorco2_0 = Parameter(default=0.8177307839027955)
    a1_percentco2oceanlong = Parameter(unit="%", default=22.97348182941324)
    a2_percentco2oceanshort = Parameter(unit="%", default=26.64454562239661)
    a3_percentco2land = Parameter(unit="%", default=26.96395086206718)
    t1_timeco2oceanlong = Parameter(unit="year", default=312.54206687556456)
    t2_timeco2oceanshort = Parameter(unit="year", default=34.873449489671636)
    t3_timeco2land = Parameter(unit="year", default=4.257701532922862)

    # Projections for components of CO2
    asymptote_co2_hist = Variable(index=[time], unit="Mtonne")
    ocean_long_uptake_component_hist = Variable(index=[time], unit="Mtonne")
    ocean_short_uptake_component_hist = Variable(index=[time], unit="Mtonne")
    land_uptake_co2hist = Variable(index=[time], unit="Mtonne")

    asymptote_co2_proj = Variable(index=[time], unit="Mtonne")
    ocean_long_uptake_component_proj = Variable(index=[time], unit="Mtonne")
    ocean_short_uptake_component_proj = Variable(index=[time], unit="Mtonne")
    land_uptake_co2_proj = Variable(index=[time], unit="Mtonne")

    te_totalemissions = Variable(index=[time], unit="Mtonne/year")
    c_CO2concentration = Variable(index=[time], unit="ppbv")

    function run_timestep(p, v, d, t)

        te0_totalemissions0 = p.e0_globalCO2emissions + p.permte0_permafrostemissions0
        if is_first(t)
            # adapted from eq.1 in Hope(2006) - calculate excess concentration in base year
            v.exc_excessconcCO2 = p.c0_CO2concbaseyr - p.pic_preindustconcCO2
            # Eq. 2 from Hope (2006) - base-year remaining emissions
            v.re_remainCO2base = v.exc_excessconcCO2 * p.den_CO2density

            v.asymptote_co2_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * p.thist_timescaleco2hist * p.stay_fractionCO2emissionsinatm
            v.ocean_long_uptake_component_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t1_timeco2oceanlong)) * (p.a1_percentco2oceanlong / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t1_timeco2oceanlong)
            v.ocean_short_uptake_component_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t2_timeco2oceanshort)) * (p.a2_percentco2oceanshort / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t2_timeco2oceanshort)
            v.land_uptake_co2hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t3_timeco2land)) * (p.a3_percentco2land / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t3_timeco2land)

            v.te_totalemissions[t] = p.e_globalCO2emissions[t] + p.permte_permafrostemissions[t]
            v.asymptote_co2_proj[t] = 0.5 * (te0_totalemissions0 + v.te_totalemissions[t]) * (p.y_year[t] - p.y_year_0) * p.stay_fractionCO2emissionsinatm
            v.ocean_long_uptake_component_proj[t] = 0.5 * (te0_totalemissions0 + v.te_totalemissions[t]) * p.t1_timeco2oceanlong * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.t1_timeco2oceanlong)) * (p.a1_percentco2oceanlong / 100)
            v.ocean_short_uptake_component_proj[t] = 0.5 * (te0_totalemissions0 + v.te_totalemissions[t]) * p.t2_timeco2oceanshort * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.t2_timeco2oceanshort)) * (p.a2_percentco2oceanshort / 100)
            v.land_uptake_co2_proj[t] = 0.5 * (te0_totalemissions0 + v.te_totalemissions[t]) * p.t3_timeco2land * (1 - exp(-(p.y_year[t] - p.y_year_0) / p.t3_timeco2land)) * p.a3_percentco2land / 100
        else
            v.te_totalemissions[t] = p.e_globalCO2emissions[t] + p.permte_permafrostemissions[t]

            v.asymptote_co2_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * p.thist_timescaleco2hist * p.stay_fractionCO2emissionsinatm
            v.ocean_long_uptake_component_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t1_timeco2oceanlong)) * (p.a1_percentco2oceanlong / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t1_timeco2oceanlong)
            v.ocean_short_uptake_component_hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t2_timeco2oceanshort)) * (p.a2_percentco2oceanshort / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t2_timeco2oceanshort)
            v.land_uptake_co2hist[t] = p.corrf_correctionfactorco2_0 * te0_totalemissions0 * (p.thist_timescaleco2hist / (1 + p.thist_timescaleco2hist / p.t3_timeco2land)) * (p.a3_percentco2land / 100) * exp(-(p.y_year[t] - p.y_year_0) / p.t3_timeco2land)

            v.asymptote_co2_proj[t] = v.asymptote_co2_proj[t - 1] + 0.5 * (v.te_totalemissions[t - 1] + v.te_totalemissions[t]) * (p.y_year[t] - p.y_year[t - 1]) * p.stay_fractionCO2emissionsinatm
            v.ocean_long_uptake_component_proj[t] = v.ocean_long_uptake_component_proj[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t1_timeco2oceanlong) + 0.5 * (v.te_totalemissions[t - 1] + v.te_totalemissions[t]) * p.t1_timeco2oceanlong * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t1_timeco2oceanlong)) * (p.a1_percentco2oceanlong / 100)
            v.ocean_short_uptake_component_proj[t] = v.ocean_short_uptake_component_proj[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t2_timeco2oceanshort) + 0.5 * (v.te_totalemissions[t - 1] + v.te_totalemissions[t]) * p.t2_timeco2oceanshort * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t2_timeco2oceanshort)) * (p.a2_percentco2oceanshort / 100)
            v.land_uptake_co2_proj[t] = v.land_uptake_co2_proj[t - 1] * exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t3_timeco2land) + 0.5 * (v.te_totalemissions[t - 1] + v.te_totalemissions[t]) * p.t3_timeco2land * (1 - exp(-(p.y_year[t] - p.y_year[t - 1]) / p.t3_timeco2land)) * p.a3_percentco2land / 100
        end

        # remaining emmissions C02 before ccf
        v.renoccf_remainCO2wocc[t] = v.asymptote_co2_hist[t] + v.ocean_long_uptake_component_hist[t] + v.ocean_short_uptake_component_hist[t] + v.land_uptake_co2hist[t] + v.asymptote_co2_proj[t] + v.ocean_long_uptake_component_proj[t] + v.ocean_short_uptake_component_proj[t] + v.land_uptake_co2_proj[t]
        # CO2 concentration CO2 before CCF
        v.conoccf_concentrationCO2wocc[t] = p.pic_preindustconcCO2 + v.exc_excessconcCO2 * (v.renoccf_remainCO2wocc[t] * v.re_remainCO2base)
        v.re_remainCO2[t] = v.renoccf_remainCO2wocc[t]

        # eq.11 from Hope(2006) - CO2 concentration
        v.c_CO2concentration[t] = p.pic_preindustconcCO2 + v.exc_excessconcCO2 * v.re_remainCO2[t] / v.re_remainCO2base
    end
end

function addco2cycle(model::Model, use_permafrost::Bool)
    co2cycle = add_comp!(model, CO2Cycle)

    if use_permafrost
        co2cycle[:permte0_permafrostemissions0] = 692.2476420621228
        co2cycle[:corrf_correctionfactorco2_0] = 0.8177307839027955
        co2cycle[:thist_timescaleco2hist] = 49.59054463320648
    else
        co2cycle[:permte0_permafrostemissions0] = 0
        co2cycle[:permte_permafrostemissions] = zeros(10)
        co2cycle[:corrf_correctionfactorco2_0] = 0.834514918731254
        co2cycle[:thist_timescaleco2hist] = 49.36461591688456
    end

    co2cycle
end


text:
module TRLWE
    using ..MulFFT
    using ..Utils
    function trlweSymEncrypt(p::Array{Float64,1},α::Float64,key::Array{UInt32,1},plan::MulFFT.TwistFFTPlan)
        a::Array{UInt32,1} = rand(0:1,length(key))
        b::Array{UInt32,1} = Utils.gaussian32(p,α) + MulFFT.PolyMul(signed.(a),signed.(key),plan)
        return [a,b]
    end
    function trlweSymDecrypt(c::Array{Array{UInt32,1},1},key::Array{UInt32,1},plan::MulFFT.TwistFFTPlan)
        return -((c[2] - MulFFT.PolyMul(signed.(c[1]),signed.(key),plan)).>UInt32(2)^31) .+ 1
    end

    function SampleExtractIndex(trlwe::Array{Array{UInt32,1},1},index::Int32)
        N = length(trlwe[1])
        return vcat([trlwe[1][index - i] for i in 0:index-1],[-trlwe[1][N+index-i] for i in index:N-1],trlwe[2][index])
    end
end


text:
"""
    iscallalist(calls, i)

    true if the call starting at `i` is a simple call without childs.
    In the example below, `add_float` is a call list, but `+` is not
    becaus it is `+`.
```
+: 8.73396258e8
 add_float: 8.73396817e8
 add_float: 8.73396957e8
+: 8.73396999e8
```
"""
function iscallalist(calls, i)
    calls.startstop[i] == :stop && return(false)
    calls.i[]  calls.stamps[i],
        :stop => calls.stamps[i+1],
        :name => calls.event[i],
        )
    new_node, i+2
end

"""
    function tape2structure(calls)

    convert flat list of logs of calls and convert it to
    structured format
"""
function tape2structure(calls::Events)
    if calls.i[] >= length(calls.stamps)
        @warn "The recording buffer was too small, consider increasing it"
    end
    root = Dict{Symbol,Any}(:children => [])
    stack = Stack{Dict}()
    node = root;
    i = 1
    while i <= min(calls.i[], length(calls.stamps))
        if iscallalist(calls, i)
            new_node, i = parselist(calls, i)
            push!(node[:children], new_node)
            continue
        end
        if calls.startstop[i] == :start
            push!(stack, node)
            new_node = Dict(
                :start => calls.stamps[i],
                :name => calls.event[i],
                :children => []
                )
            # println(i, " new node ", calls.event[i])
            push!(node[:children], new_node)
            node = new_node
            i += 1
            continue
        end

        if calls.startstop[i] == :stop
            # println(i, " end node ", calls.event[i], " expecting: ", node[:name])
            @assert node[:name] == calls.event[i]
            node[:stop] = calls.stamps[i]
            i += 1
            node = pop!(stack)
            continue
        end
    end
    root[:children]
end
tape2structure() = [tape2structure(t) for t in to];

_visualize(filename;scaling = 1) = _visualize(filename, tape2structure(); scaling)
function _visualize(filename::String, events; basetime= events[1][:start], scaling = 1)
    open(filename, "w") do fio
        header(fio)
        _visualize(fio, events, Dict(); basetime, scaling)
        footer(fio)
    end
end

function _visualize(io::IO, events, mapping ; basetime, scaling)
    for event in events
        name = event[:name]
        start = scaling * event[:start]
        stop = scaling * event[:stop]
        expandable = haskey(event, :children) ? "expandable" : ""
        duration = round((event[:stop] - event[:start]) / 10^9, digits=2)

        if name ∉ keys(mapping)
            mapping[name] = Dict(:id => rand(Int), :color => random_color(), :counter => 0, :children => Dict())
        end
        m = mapping[name]

        println(io, """

                $(!haskey(event, :children) ? "" : """ """)$(name)

                    $(duration)s)

        """)
        m[:counter] += 1

        if haskey(event, :children)
            _visualize(io, event[:children], m[:children] ; basetime=event[:start], scaling=scaling)
        end

        println(io, """

""")
    end
end

function random_color()
    c = Colors.color_names
    return("#"*hex(RGB((c[rand(keys(c))]./255)...)))
end

function header(io::IO)
    print(io, """

  <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.15.2/css/all.min.css" />
  <link rel="stylesheet" href="https://unpkg.com/purecss@2.0.5/build/pure-min.css" integrity="sha384-LTIDeidl25h2dPxrB2Ekgc9c7sEC3CWGM6HeFmuDNUjX76Ert4Z4IY714dhZHPLd" crossorigin="anonymous">

""")
end

function footer(io::IO)
    print(io, """

          This call
          Overall

          Time
          Memory
          GC Time
          CPU util
          Time
          Memory
          GC Time
          CPU util

""")
end

max_stars_repo_name: pevnak/LoggingProfiler

max_issues_repo_name: pevnak/LoggingProfiler

max_forks_repo_name: pevnak/LoggingProfiler


text:
using ArgParse
using MLDataUtils
using ForwardDiff
using Distances
using LinearAlgebra
using ReverseDiff
using AutoGrad
using Random
#using CuArrays
using Flux
using Flux.Tracker
using Flux: @epochs
using Base.Iterators: partition
using DelimitedFiles

arguments = ArgParseSettings()
@add_arg_table arguments begin
    "--input", "-i"
    help = "input file containing normalized OWL EL axioms; usually the output of Normalize.groovy"
    required = true
    "--output", "-o"
    help = "output file containing class, relation, and instance coordinates"
    required = true
    "--epochs", "-e"
    help = "number of epochs to train"
    arg_type = Int
    default = 1000
    "--lr", "-l"
    help = "learning rate"
    arg_type = Float64
    default = 0.01
    "--dim", "-d"
    help = "input dimensions"
    arg_type = Int
    default = 20
    "--gpu", "-g"
    action = :store_true
    help = "use GPU accelleration"
    "--margin", "-m"
    arg_type = Float64
    default = 0.1
    help = "margin parameter"
end
args = parse_args(ARGS, arguments)
if args["gpu"] == true
    using CuArrays
end

# Model parameters

LR = args["lr"] # learning rate
EPOCHS = args["epochs"] # number of epochs
DIM = args["dim"] # dimensionality of embeddings

# DL Constants; not valid indices of an array
BOTTOM = -1
TOP = -2

# Setting up the data

classes = Dict()
relations = Dict()
ccount = 1
rcount = 1

classes["owl:Nothing"] = BOTTOM
classes["<http://www.w3.org/2002/07/owl#Nothing>"] = BOTTOM
classes["owl:Thing"] = TOP
classes["<http://www.w3.org/2002/07/owl#Thing>"] = TOP

function filldict(x)
    global ccount
    if haskey(classes, x)
        return classes[x]
    else
        classes[x] = ccount
        ccount += 1
        return ccount - 1
    end
end

function rfilldict(x)
    global rcount
    if haskey(relations, x)
        return relations[x]
    else
        relations[x] = rcount
        rcount += 1
        return rcount - 1
    end
end

nf1 = Set()
nf2 = Set()
nf3 = Set()
nf4 = Set()

# Reading the file of normalized axioms (generated by Normalize.groovy)
open(args["input"]) do file
    for line in eachline(file)
        if (startswith(line, "SubClassOf")) # ignore subproperty axioms
            if occursin(r"ObjectIntersectionOf", line) # normal form 2
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: C and D SubClassOf: E
                if m != nothing
                    c, d, e = map(filldict, m.captures)
                    push!(nf2, (c,d,e))
                end
            elseif occursin(r"SubClassOf.<[^\s]*>[\s]*<[^\s]*>.$", line) # normal form 1
                m = match(r"(<.*>).*(<.*>)", line) # 2 captures: C SubClassOf: D
                c,d = map(filldict, m.captures)
                push!(nf1, (c,d))
            elseif occursin(r"SubClassOf.ObjectSomeValuesFrom", line)  # normal form 4
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: R some C SubClassOf: D
                r,c,d = m.captures
                c = filldict(c)
                d = filldict(d)
                r = rfilldict(r)
                push!(nf4, (c,d,r))
            elseif occursin(r"SubClassOf.<[^\s]*>[\s]*ObjectSomeValuesFrom", line)  # normal form 3
                m = match(r"(<.*>).*(<.*>).*(<.*>)", line) # 3 captures: C SubClassOf: R some D
                c,r,d = m.captures
                c = filldict(c)
                d = filldict(d)
                r = rfilldict(r)
                push!(nf3, (c,d,r))
            end
        end
    end
end

# there are now rcount-1 relations and ccount-1 classes; generate random embeddings; offload to GPU; make them trainable
cvec = rand(ccount-1, DIM + 1)
rvec = rand(rcount-1, DIM)
#cvec = rand(ccount-1, DIM + 1)
#rvec = rand(rcount-1, DIM)

@inline function radius(x::Array{AbstractFloat,1})
    return x[end]
end

@inline function centerpoint(x::Array{AbstractFloat,1})
    return x[1:end-1]
end

@inline function radius(x)
    return x[end]
end

@inline function centerpoint(x)
    return x[1:end-1]
end

# Radius: 0 if BOTTOM; infinity if TOP; otherwise to be optimized
@inline function radius(x::Int)
    if x == BOTTOM
        return 0
    elseif x == TOP
        return Inf
    else
        return radius(cvec[x,:])
    end
end

# Centerpoint: (0,...,0) if BOTTOM or TOP; otherwise to be optimized
@inline function centerpoint(x::Int)
    if x == BOTTOM
        return zeros(DIM)
    elseif x == TOP
        return zeros(DIM)
    else
        return centerpoint(cvec[x,:])
    end
end

@inline function v(x::Int)
    return rvec[x,:]
end

# loss functions that work with class indices
function loss1(c::Int, d::Int)
    return max(0, euclidean(centerpoint(c), centerpoint(d)) + abs(radius(c)) - abs(radius(d)) + args["margin"])
end
function loss1(x::Tuple{Int64,Int64})
    return loss1(x[1],x[2])
end
function loss1(x::Array{Tuple{Int64,Int64}})
    return loss1.(x)
end

function loss2(c1::Int, c2::Int, d::Int) # loss for c1 and c2 SubClassOf: d
    dist = euclidean(centerpoint(c1), centerpoint(c2))
    if dist > abs(radius(c1)) + abs(radius(c2)) # no solution, circles are separate
        rad = 0 # radius is 0
        cent = centerpoint(c1) + (centerpoint(c2) - centerpoint(c1)) / (dist - radius(c2)) # half-way between outer edges of the class
        # loss1 between the class with radius 0 and centerpoint half-way between c1 and c2, plus a penalty for the distance between c1 and c2
        return max(0, euclidean(cent, centerpoint(d)) + abs(rad) - abs(radius(d)) + (dist - abs(radius(c1)) - abs(radius(c2)) + args["margin"]))
    elseif dist < abs(radius(c1)-radius(c2)) # no solution, one circle contained in other
        if abs(radius(c1)) < abs(radius(c2))
            return loss1(c1, d) # return loss1 of smaller class (c1) and d
        else
            return loss1(c2, d) # return loss1 of smaller class (c2) and d
        end
    elseif dist == 0 && radius(c1) == radius(c2) ## circles are coincident
        return loss1(c1, d) # choose one
    else
        a = (radius(c1)^2 - radius(c2)^2 + dist^2)/(2 * dist)
        rad = sqrt(radius(c1)^2 - a^2)
        cent = centerpoint(c1) + a * (centerpoint(c2) - centerpoint(c1)) / dist
        return max(0, euclidean(cent, centerpoint(d)) + abs(rad) - abs(radius(d)) + args["margin"])

    end
end

function loss2(x::Tuple{Int64,Int64,Int64})
    return loss2(x[1], x[2], x[3])
end

function loss3(c::Int, d::Int, r::Int) # normal form 3
    return max(0, euclidean(centerpoint(d) - v(r), centerpoint(c)) + abs(radius(c)) - abs(radius(d)) + args["margin"])
end

function loss3(x::Tuple{Int64,Int64,Int64})
    return loss3(x[1], x[2], x[3])
end

function loss4(c::Int, d::Int, r::Int) # normal form 4
    return max(0, euclidean(centerpoint(c) + v(r), centerpoint(d)) + abs(radius(c)) - abs(radius(d)) + args["margin"])
end

function loss4(x::Tuple{Int64,Int64,Int64})
    return loss4(x[1], x[2], x[3])
end

# generate arrays for our training data, easier to move to GPU
nf1arr =  shuffle(collect(nf1))
nf2arr =  shuffle(collect(nf2))
nf3arr =  shuffle(collect(nf3))
nf4arr =  shuffle(collect(nf4))

#BATCHSIZE = 10
#nf1arr = Flux.chunk(nf1arr, ceil(length(nf1arr)/BATCHSIZE))
#nf2arr = Flux.chunk(nf2arr, ceil(length(nf2arr)/BATCHSIZE))
#nf3arr = Flux.chunk(nf3arr, ceil(length(nf3arr)/BATCHSIZE))
#nf4arr = Flux.chunk(nf4arr, ceil(length(nf4arr)/BATCHSIZE))

# Training now...

if args["gpu"] == true
    cvec = param(cvec) |> gpu
    rvec = param(rvec) |> gpu
else
    cvec = param(cvec)
    rvec = param(rvec)
end

opt = SGD([cvec, rvec], LR)
if (length(nf1) < 10) || (length(nf2) < 10) || (length(nf3) < 10) || (length(nf4) < 10)
    for i in 1:EPOCHS
        Flux.train!(loss1, nf1arr, opt)
        Flux.train!(loss2, nf2arr, opt)
        Flux.train!(loss3, nf3arr, opt)
        Flux.train!(loss4, nf4arr, opt)
    end
else
    TESTSIZE = 20
    test1 = collect(RandomBatches(nf1arr, TESTSIZE, 1))
    test2 = collect(RandomBatches(nf2arr, TESTSIZE, 1))
    test3 = collect(RandomBatches(nf3arr, TESTSIZE, 1))
    test4 = collect(RandomBatches(nf4arr, TESTSIZE, 1))

    #evalcb1 = throttle(() -> @show(loss1([1])), 5)
    evalcb1 = Flux.throttle(() -> @show(sum(loss1.(test1[1]))), 20)
    evalcb2 = Flux.throttle(() -> @show(sum(loss2.(test2[1]))), 20)
    evalcb3 = Flux.throttle(() -> @show(sum(loss3.(test3[1]))), 20)
    evalcb4 = Flux.throttle(() -> @show(sum(loss4.(test4[1]))), 20)

    for i in 1:EPOCHS
        Flux.train!(loss1, nf1arr, opt, cb=evalcb1)
        Flux.train!(loss2, nf2arr, opt, cb=evalcb2)
        Flux.train!(loss3, nf3arr, opt, cb=evalcb3)
        Flux.train!(loss4, nf4arr, opt, cb=evalcb4)
    end
end

#println("Loss is ", sum(loss1.(nf1arr)) + sum(loss2.(nf2arr)) + sum(loss3.(nf3arr)) + loss4.(nf4arr))

cvec = Flux.Tracker.data(cvec)
rvec = Flux.Tracker.data(rvec)
open(args["output"], "w") do file
    write(file, "IRI")
    for i in 1:DIM
        write(file, "\tV$i")
    end
    write(file, "\tr\n")
    for (c,ind) in classes
        if ind > 0
            write(file, c)
            vec = cvec[ind,:]
            for val in vec
                write(file,"\t$val")
            end
            write(file, "\n")
        end
    end
    for (r,ind) in relations
        if ind > 0
            write(file, r)
            vec = rvec[ind,:]
            for val in vec
                write(file,"\t$val")
            end
            write(file, "\n")
        end
    end
end

max_stars_repo_name: leechuck/julia-stuff

max_issues_repo_name: leechuck/julia-stuff

max_forks_repo_name: leechuck/julia-stuff


text:
using GR
using Base.Test

@test tick(1.2,3.14) == 0.5

tests = [ "ex", "griddata" ]

for t in tests
    tp = joinpath(Pkg.dir("GR"), "test", "$(t).jl")
    println("running $(tp) ...")
    inline("pdf")
    include(tp)
    @test isfile("gks.pdf")
    rm("gks.pdf")
end


text:
# Numerical solutions of planetary structure models

# Solver types

abstract type Integrator end
abstract type FixedStepIntegrator <: Integrator end

const FloatOrVecFloat = Union{Float64, Vector{Float64}}

"RK4 solver"
type RK4{ResultType<:FloatOrVecFloat} <: FixedStepIntegrator
    func
    y0::ResultType
    tgrid::AbstractVector
    ks::Vector{ResultType}
    y::ResultType
end
RK4(f, y0::FloatOrVecFloat, tgrid) = RK4(f, y0, tgrid, [zero(y0) for _=1:4], copy(y0))
RK4(f, y0::Real, tgrid) = RK4(f, Float64(y0), tgrid)
RK4{R<:Real}(f, y0::Vector{R}, tgrid) = RK4(f, collect(Float64, y0), tgrid)

const ScalarRK4 = RK4{Float64}
const VectorRK4 = RK4{Vector{Float64}}

# Working with integrator and state types

function reset!(integrator::ScalarRK4)
    integrator.y = integrator.y0
    integrator.ks = [zero(integrator.y0) for _=1:4]
    return nothing
end
function reset!(integrator::VectorRK4)
    integrator.y[:] = integrator.y0
    integrator.ks = [zero(integrator.y0) for _=1:4]
    return nothing
end
nvars(integrator) = length(integrator.y0)
ntimes(integrator) = length(integrator.tgrid)
Base.length(integrator::RK4) = ntimes(integrator)
time(integrator, n) = integrator.tgrid[n]
nexttime(integrator, n) = integrator.tgrid[n + 1]
y(integrator::Integrator) = integrator.y
y0(integrator) = integrator.y0
atfinaltime(integrator, n) = (n == ntimes(integrator))
pastfinaltime(integrator, n) = (n > ntimes(integrator))
Base.eltype(::Type{ScalarRK4}) = Float64
Base.eltype(::Type{VectorRK4}) = Vector{Float64}

# Numerical steps

const rk4_dy_coeffs = [1/3, 1/6, 1/6, 1/3]

"Do a RK4 step of function `F`, initial conds `x0` and `t0`, step size `dt`"
function ode4_step!(i::ScalarRK4, t0::Float64, t1::Float64)
    let ks = i.ks, F = i.func, order = 4, dt = t1 - t0

        ks[1] = dt*F(t0,        i.y)
        ks[2] = dt*F(t0 + dt/2, i.y + ks[1]/2)
        ks[3] = dt*F(t0 + dt/2, i.y + ks[2]/2)
        ks[4] = dt*F(t0 + dt,   i.y + ks[3])

        for d = 1:order
            i.y += ks[d] * rk4_dy_coeffs[d]
        end

        nothing
    end
end
function ode4_step!(i::VectorRK4, t0::Float64, t1::Float64)
    let ks = i.ks, y = i.y, F = i.func, n = nvars(i), order = 4, dt = t1 - t0

        ks[1] = dt*F(t0,        y)
        ks[2] = dt*F(t0 + dt/2, y + ks[1]/2)
        ks[3] = dt*F(t0 + dt/2, y + ks[2]/2)
        ks[4] = dt*F(t0 + dt,   y + ks[3])

        for d = 1:order, j=1:n
            y[j] += ks[d][j] * rk4_dy_coeffs[d]
        end

        nothing
    end
end
function ode4_step!(i::Integrator, n)
    let t = time(i, n), tn = nexttime(i, n)
        ode4_step!(i, t, tn)
    end

    nothing
end
function nexty!(i::Integrator, n)
    if atfinaltime(i, n)
        # can't integrate past the end point
        nothing
    else
        # do an integration
        ode4_step!(i, n)
    end
end

# iterator setup
Base.start(rk4::RK4) = (reset!(rk4); return 1)
# stepper
function Base.next(rk4::RK4, n)
    yout = copy(y(rk4))
    nexty!(rk4, n)
    return yout, n + 1
end
# termination condition
Base.done(fs::FixedStepIntegrator, n) = pastfinaltime(fs, n)

# for planets, if a value becomes unphysical (such as r<0 or P<0) then the
# structure equations will yield zeroes for each remaining step. This allows us
# to step "over" the core and leaves us with a slightly negative value at the
# centre. This negative sign signals the radius search function that we should
# increase our radius.

# Dense solvers

""" Solve the ODE dx/dt = `F` with initial conds `x0` on time grid `tgrid`.

    Returns dense output (an array of solutions at each point in `tgrid`) """
function ode4_dense(F, y0::Real, tgrid)
    integrator = RK4(F, y0, tgrid)

    # return a 1D array
    [y for y in integrator]
end
# version for multi-valued functions
function ode4_dense(F, x0::Vector, tgrid)
    integrator = RK4(F, x0, tgrid)

    # return a 2D array
    [y[i] for y in integrator, i in 1:nvars(integrator)]
end

# Solutions

"Solve a `PlanetSystem`, writing the solution to a given `PlanetStructure`."
function solve!(sys::PlanetSystem, soln::PlanetStructure)
    # initial values
    y0 = nonmass(sys.boundary_values)
    tgrid = sys.solution_grid

    # solve and write to solution array
    solver = RK4(sys.structure_equations, y0, tgrid)
    for (i, y) in enumerate(solver)
        setmass!(soln, i, tgrid[i])
        setnonmass!(soln, i, y)
    end

    return soln
end

"Generate a new `PlanetStructure` object and solve a `PlanetSystem`"
function solve(sys::PlanetSystem)
    soln = blank_structure(sys)
    solve!(sys, soln)

    return soln
end

# Radius searching

# helper functions
"Has the radius passed the central point?"
hit_the_centre(R::Radius) = R < 0
hit_the_centre(ps::PlanetStructure) = hit_the_centre(radius(centre(ps)))
"Is the radius not yet close enough to the centre?"
not_far_enough(R::Radius) = R > 1000
not_far_enough(ps::PlanetStructure) = not_far_enough(radius(centre(ps)))
unacceptable(ps) = hit_the_centre(ps) || not_far_enough(ps)
"Is the structural solution acceptable?"
acceptable(system, structure) = !unacceptable(structure)

"Update the radius search bracket using a bisection search"
function refine_r_bracket!(system::PlanetSystem, previousresult::PlanetStructure)
    if hit_the_centre(previousresult)
        minR = currentradiusguess(system)
        maxR = maxradius(system)
        refine_r_bracket!(system, minR, maxR)
    elseif not_far_enough(previousresult)
        minR = minradius(system)
        maxR = currentradiusguess(system)
        refine_r_bracket!(system, minR, maxR)
    end
end

"Update the planet's boundary conditions using the results of a test solution"
function refine_boundary_conditions!(system::PlanetSystem, previousresult::PlanetStructure)
    refine_r_bracket!(system, previousresult)
    refine_boundary_r!(system)
    refine_surface!(system)
end

""" Repeatedly solve a `PlanetSystem` until its radius is suitable

    Recursively search the radius search bracket of a `PlanetSystem` to
    get a result that has an acceptable error at the centre. Returns a
    `PlanetSystem` with the appropriate radius. """
function converge!(system)
    # first make sure the system has the correct boundary conditions to start with
    refine_boundary_r!(system)
    refine_surface!(system)

    # then get the structure for this radius guess
    result = solve(system)

    # and recursively iterate
    converge!(system, result)
    return system
end
function converge!(system, result)
    # we have already done an iteration and have a result to work with
    if !acceptable(system, result)
        ΔR = system.radius_search_bracket[2] - system.radius_search_bracket[1]
        R0 = radius(centre(result))
        if ΔR < 0.1
            info("Incomplete centre convergence ($R0 m)")
            return system
        end
        # need to refine and try again
        refine_boundary_conditions!(system, result)
        solve!(system, result)
        converge!(system, result)
    else
        # we're done
        return system
    end
end

# basic functions for geting radii and structures
"Get the correct internal structure for a `PlanetSystem`"
find_structure!(system) = solve(converge!(system))
"Find the radius of a `PlanetSystem` by solving the system"
find_radius!(system) = radius(surface(find_structure!(system)))

"Solve a `PlanetSystem` and give both the total radius and the structure"
function find_structure_and_radius!(system::PlanetSystem)
    structure = find_structure!(system)
    r = radius(surface(structure))

    structure, r
end

# higher level functions for doing MR diagrams
"Create and solve for a `PlanetSystem`'s radius'"
function R(M::Mass, structure_equations::EquationSet, Psurf::Pressure,
           solution_grid, R_bracket)

    R_guess = mean(R_bracket)
    system = PlanetSystem(M, R_guess, Psurf, structure,
                          solution_grid, R_bracket)

    radius = find_radius!(system)
end

"Find the radius of a solid sphere of mass `M` using an given `EOS`."
function R(M, eos; in_earth_units=false)
    Mscale = in_earth_units ? M_earth : 1.
    Rscale = in_earth_units ? 1/R_earth : 1.

    sys = DefaultPlanetSystem(M * Mscale, eos)
    radius = find_radius!(sys) * Rscale
end

# temperature-dependent version
function R(M, eos, Cₚ; in_earth_units=false)
    Mscale = in_earth_units ? M_earth : 1.
    Rscale = in_earth_units ? 1/R_earth : 1.

    sys = DefaultPlanetSystem(M * Mscale, eos, Cₚ)
    radius = find_radius!(sys) * Rscale
end

# vectorized form of the above
function R{T<:Real}(ms::AbstractVector{T}, args...; in_earth_units=false)
    R_withEOS(M) = R(M, args...; in_earth_units=in_earth_units)
    radii = map(R_withEOS, ms)
end


text:
module TestML

using FactCheck
importall GradientBoost.ML
importall GradientBoost.LossFunctions

instances = [
  1.0 1.0;
  1.0 8.0;
  1.0 10.0
]
labels = [
  0.0;
  1.0;
  1.0;
]

facts("Machine Learning API") do
  context("not implemented functions throw an error") do
    gbl = GBLearner(GBDT(), :regression)
    instances = 1
    labels = 1

    @fact_throws fit!(gbl, instances, labels)
    @fact_throws predict!(gbl, instances)
  end

  context("fit! on Float64 arrays works") do
    gbl = GBLearner(GBDT(), :regression)
    @fact gbl.model => nothing
    fit!(gbl, instances, labels)
    @fact gbl.model => not(nothing)
  end

  context("predict! on Float64 arrays works") do
    gbl = GBLearner(GBDT(;loss_function=BinomialDeviance()), :class)
    fit!(gbl, instances, labels)
    predictions = predict!(gbl, instances)
    @fact eltype(predictions) => Float64
  end

  context("postprocess_pred works") do
    predictions = [-Inf, 0.0, Inf]
    expected = [0.0, 1.0, 1.0]
    actual = ML.postprocess_pred(:class, BinomialDeviance(), predictions)
    @fact actual => expected

    predictions = [-Inf, 0.0, Inf]
    actual = ML.postprocess_pred(:class_prob, BinomialDeviance(), predictions)
    @fact all(i -> (0 <= i <= 1), actual) => true

    predictions = [-Inf, 0.0, Inf]
    expected = predictions
    actual = ML.postprocess_pred(:regression, LeastSquares(), predictions)
    @fact actual => expected

    predictions = [-Inf, 0.0, Inf]
    @fact_throws ML.postprocess_pred(:class, LeastSquares(), predictions)
  end

  context("logistic works") do
    x = [-Inf, -1, 0, 1, Inf]
    expected = [0.0, 0.0, 1.0, 1.0, 1.0]

    actual = round(ML.logistic(x))
    @fact actual => expected
  end
end

end # module


text:
# Useful for CDD, LRS and matrix representations

# HyperPlanes and HalfSpaces mixed
abstract type MixedHRep{T} <: HRepresentation{T} end
# SymPoints and Points mixed and/or Lines and Rays mixed
abstract type MixedVRep{T} <: VRepresentation{T} end

const MixedRep{T} = Union{MixedHRep{T}, MixedVRep{T}}

# This could be defined for all types Indices but it is not recommended for reps other than Mixed to use this implementation of length as it is inefficient
# a MethodError is more helpful than a hidden inefficiency
function mixedlength(idxs::Indices)
    count = 0
    for idx in idxs
        count += 1
    end
    count
end
Base.length(idxs::Indices{T, ElemT, <:MixedRep{T}}) where {T, ElemT} = mixedlength(idxs)

function mixednext(rep::MixedRep{T}, idx::IdxT) where {T, ElemT, IdxT<:Index{T, ElemT}}
    idx = IdxT(idx.value+1)
    while !done(Indices{T, ElemT}(rep), idx) && !isvalid(rep, idx)
        idx = IdxT(idx.value+1)
    end
    if done(Indices{T, ElemT}(rep), idx)
        return nothing
    else
        return idx
    end
end
startindex(idx::Indices{T, ElemT, <:MixedRep{T}}) where {T, ElemT} = mixednext(idx.rep, Index{T, ElemT}(0))
nextindex(rep::MixedRep{T}, idx::Index{T}) where {T} = mixednext(rep, idx)


text:
function parse_registry(
    registry_dir = joinpath(homedir(), ".julia", "registries", "General")
)
    registry_dict = TOML.parsefile(joinpath(registry_dir, "Registry.toml"))
    pkgs = Dict{String, Package}()
    for (uuid_str, pkg_lookup) in registry_dict["packages"]
        pkg_name = pkg_lookup["name"]::String
        rel_path = pkg_lookup["path"]::String
        abs_path = joinpath(registry_dir, rel_path)
        pkg = parse_registry_package_dir(abs_path)
        pkgs[pkg_name] = pkg
    end
    pkgs
end

function parse_registry_package_dir(path::AbstractString)
    info = TOML.parsefile(joinpath(path, "Package.toml"))
    deps_data = Operations.load_package_data_raw(
        UUID, joinpath(path, "Deps.toml")
    )
    compat_data = Operations.load_package_data_raw(
        VersionSpec, joinpath(path, "Compat.toml")
    )
    toml_versions = Operations.load_versions(path)
    versions = Dict{VersionNumber, Version}()
    for v in keys(toml_versions)
        reqdict = Dict{String, Require}()
        ver = Version(nothing, reqdict)
        for (vr, dd) in compat_data
            if v in vr
                for (reqname, vs) in dd
                    reqdict[reqname] = Require(convert(VersionSet, vs), Symbol[])
                end
            end
        end
        for (vr, dd) in deps_data
            if v in vr
                for (reqname, uuid) in dd
                    if ! haskey(reqdict, reqname)
                        reqdict[reqname] = Require(VersionSet(), Symbol[])
                    end
                end
            end
        end
        versions[v] = ver
    end
    return Package(
        UUID(info["uuid"]),
        info["repo"],
        versions
    )
end

convert(::Type{VersionNumber}, v::VersionBound) = VersionNumber(v.t...)
function convert(::Type{VersionSet}, vs::VersionSpec)
    VersionSet(
        convert(VersionNumber, vs.ranges[1].lower),
        next_semver(convert(VersionNumber, vs.ranges[1].upper))
    )
end

function next_semver(v::VersionNumber)
    if v.major == 0
        return VersionNumber(0, v.minor + 1, 0)
    else
        return VersionNumber(v.major + 1, 0, 0)
    end
end


text:
module Atoms

using MultiIndices
using AtomicLevels
using JSON

import Base: show

const elements = JSON.parse(open(joinpath(Pkg.dir("Atoms"), "periodic-table", "data.json")))
const periodic_table = Dict(Symbol(e["symbol"]) => e
                            for e in elements)

abstract type Atom end

space(a::Atom) = a.𝔓

function show(io::IO, a::A) where A<:Atom
    write(io, "$(a.name) ($(typeof(a))")
    space_descr(io, a)
    write(io, ")\n")
    write(io, "N = $(a.N)")
    if(a.Z ≠ a.N)
        write(io, ", Z = $(a.Z)")
    end
    write(io, ", ground state: $(a.gst_config)")
end

struct BreitPauliAtom <: Atom
    𝔓::Space
    name::String
    gst_config::Config
    N::Integer # Number of electrons
    Z::Integer # Nuclear charge
end

"""
    BreitPauliAtom(N, ℓmax; Z=N, mmax=ℓmax)

Construct an atom in the Breit–Pauli approximation (non-relativistic,
but with spin) with `N` electrons, optionally specifying the atomic
number `Z` (for ions).
"""
function BreitPauliAtom(N::Integer, ℓmax::Integer;
                        Z::Integer=N,
                        mmax::Integer = ℓmax)
    element = elements[N]
    gst_config = ref_set_list(element["electronicConfiguration"])
    name = element["name"]
    if N ≠ Z
        name = "$(name)-like $(elements[Z]["name"])"
    end

    𝕷 = RotationSpace(0)
    for ℓ = 1:ℓmax
        𝕷 = 𝕷 ⊕ RotationSpace(ℓ,mmax)
    end
    𝔖 = Spin{1//2}()
    𝕵 = 𝕷 ⊗ 𝔖
    𝔓 = 𝕵 ⊗ N
    BreitPauliAtom(𝔓, name, gst_config, N, Z)
end
BreitPauliAtom(N::Symbol, args...; Z::Symbol=N, kwargs...) =
    BreitPauliAtom(periodic_table[N]["atomicNumber"], args...;
                   Z = periodic_table[Z]["atomicNumber"],
                   kwargs...)

function space_descr(io::IO, a::BreitPauliAtom)
    𝕵 = space(a).𝔖
    𝕷 = 𝕵.𝔄
    𝔖 = 𝕵.𝔅
    nopws = size(𝕷)
    ℓmax = length(𝕷.𝔖)
    write(io, "; ℓmax = $(ℓmax)")
    if nopws ≠ ℓmax
        write(io, ", $(nopws) partial waves")
    else
        write(io, ", 2d")
    end
    write(io, " ⊗ 2 spins")
end

struct SAEAtom <: Atom
    𝔓::Space
    name::String
    gst_config::Config
    N::Integer # Number of electrons
    Z::Integer # Nuclear charge
end

"""
    SAEAtom(N, ℓmax; Z=N, mmax=ℓmax)

Construct an atom in the single-active electron approximation
(non-relativistic, no spin) with `N` electrons (of which, only one
valence electron is active), optionally specifying the atomic number
`Z` (for ions).
"""
function SAEAtom(N::Integer, ℓmax::Integer;
                 Z::Integer=N,
                 mmax::Integer = ℓmax)
    element = elements[N]
    gst_config = ref_set_list(element["electronicConfiguration"])
    name = element["name"]
    if N ≠ Z
        name = "$(name)-like $(elements[Z]["name"])"
    end

    𝕷 = RotationSpace(0)
    for ℓ = 1:ℓmax
        𝕷 = 𝕷 ⊕ RotationSpace(ℓ,mmax)
    end
    SAEAtom(𝕷, name, gst_config, N, Z)
end
SAEAtom(N::Symbol, args...; Z::Symbol=N, kwargs...) =
    SAEAtom(periodic_table[N]["atomicNumber"], args...;
            Z = periodic_table[Z]["atomicNumber"],
            kwargs...)

SAEAtom2D(args...; kwargs...) = SAEAtom(args...; mmax=0, kwargs...)

function space_descr(io::IO, a::SAEAtom)
    𝕷 = space(a)
    nopws = size(𝕷)
    ℓmax = length(𝕷.𝔖)
    write(io, "; ℓmax = $(ℓmax)")
    if nopws ≠ ℓmax
        write(io, ", $(nopws) partial waves")
    else
        write(io, ", 2d")
    end
end

export Atom, space, BreitPauliAtom, SAEAtom, SAEAtom2D

end # module


text:
using DashBootstrapComponents
using HTTP

LOREM = rstrip(
    String(
        HTTP.get(
            "https://raw.githubusercontent.com/facultyai/dash-bootstrap-components/main/docs/components_page/components/modal/lorem.txt",
        ).body,
    ),
    ['\n'],
)

modal = html_div([
    dbc_button("Scrolling modal", id = "open-scroll", className = "me-1", n_clicks = 0),
    dbc_button("Modal with scrollable body", id = "open-body-scroll", n_clicks = 0),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Scrolling modal")),
            dbc_modalbody(LOREM),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id = "close-scroll",
                    className = "ms-auto",
                    n_clicks = 0,
                ),
            ),
        ],
        id = "modal-scroll",
        is_open = false,
    ),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Modal with scrollable body")),
            dbc_modalbody(LOREM),
            dbc_modalfooter(
                dbc_button(
                    "Close",
                    id = "close-body-scroll",
                    className = "ms-auto",
                    n_clicks = 0,
                ),
            ),
        ],
        id = "modal-body-scroll",
        scrollable = true,
        is_open = false,
    ),
]);

function toggle_modal(n1, n2, is_open)
    return n1 > 0 || n2 > 0 ? is_open == 0 : is_open
end;

callback!(
    app,
    Output("modal-scroll", "is_open"),
    Input("open-scroll", "n_clicks"),
    Input("close-scroll", "n_clicks"),
    State("modal-scroll", "is_open"),
) do n1, n2, is_open
    return toggle_modal(n1, n2, is_open)
end;

callback!(
    app,
    Output("modal-body-scroll", "is_open"),
    Input("open-body-scroll", "n_clicks"),
    Input("close-body-scroll", "n_clicks"),
    State("modal-body-scroll", "is_open"),
) do n1, n2, is_open
    return toggle_modal(n1, n2, is_open)
end;

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
module Statsd

using Sockets

mutable struct Client
    host::IPv4
    port::Integer
    sock::IO
    prefix::String
    function Client(host::String="localhost", port::Integer=8125, prefix="")
        host = getaddrinfo(host)
        sock = UDPSocket()
        new(host, port, sock, prefix)
    end
end

# Send the metrics
function sc_send(sc::Client, data::String)
    send(sc.sock, sc.host, sc.port, data)
end

# Counter Functions
incr(sc::Client, metric, rate=nothing) = sc_metric(sc, "c", metric, 1, rate)
decr(sc::Client, metric, rate=nothing) = sc_metric(sc, "c", metric, -1, rate)

# Gauge, Timers, and Set
gauge(sc::Client, metric, value, rate=nothing) = sc_metric(sc, "g", metric, value, rate)
timing(sc::Client, metric, value, rate=nothing) = sc_metric(sc, "ms", metric, value, rate)
set(sc::Client, metric, value, rate=nothing) = sc_metric(sc, "s", metric, value, rate)

# Generate the metric and call send func
function sc_metric(sc::Client, type, metric, value, rate)
    if rate == nothing
        sc_send(sc, (isempty(sc.prefix) ? "" : "$(sc.prefix).") * "$metric:$value|$type")
    else
        sc_send(sc, (isempty(sc.prefix) ? "" : "$(sc.prefix).") * "$metric:$value|$type@$rate")
    end
end

end # module


text:
const intptr_t = UInt
const time_t = Int

const AV_NOPTS_VALUE = 0x8000000000000000 % Int64
AV_VERSION_INT(a, b, c) = (a << 16 | b << 8) | c
AV_VERSION(a, b, c) = nothing
const AVPROBE_SCORE_MAX = 100
MKTAG(a, b, c, d) = (UInt32(a) | (UInt32(b) << 8) | (UInt32(c) << 16) | UInt32(d) << 24)
FFERRTAG(a, b, c, d) = -MKTAG(a, b, c, d)
macro AV_PIX_FMT_NE(be, le)
    return Symbol("AV_PIX_FMT_" * string(le))
end


text:
export Path, InteractionSequence, InteractionSequenceSample

const Path{T} = Vector{T} where {T<:Union{Int, String}}
const InteractionSequence{T} = Vector{Vector{T}} where {T<:Union{Int, String}}
const InteractionSequenceSample{T} = Vector{Vector{Vector{T}}} where {T<:Union{Int, String}}