text:
# -----------------------------------------------------------------------------
# IMPORTS
# -----------------------------------------------------------------------------

# External libraries
using Logging                       # Logging is used for operation diagnostics
using Parameters                    # MetaICVIOpts are Parameters structs
using StatsBase: corspearman        # Rank correlation for cvi criterion values
using PyCall                        # PyCall object definition
using JLD                           # JLD is currently recommended for saving/loading ScikitLearn objects
using PyCallJLD                     # PyCall definition for serialization with JLD
using ScikitLearn                   # Classifiers are scikit-learn pyobjects
using ProgressMeter: @showprogress  # Data loading progress for training
using DelimitedFiles: readdlm       # Loading cvi data
# using BSON
# using ScikitLearn.Skcore: FitBit

# Custom libraries
using ClusterValidityIndices        # All Julia-implemented CVI definitions

# Local libraries
# Get the rocket kernel definitions
include("rocket.jl")
using .Rocket

# -----------------------------------------------------------------------------
# CONSTANTS
# -----------------------------------------------------------------------------

# const MetaICVIClassifier = ScikitLearn.Skcore.FitBit
const MetaICVIClassifier = PyCall.PyObject

# Top of the module for default paths
const module_dir(paths...) = joinpath(dirname(pathof(@__MODULE__)), "..", paths...)

# -----------------------------------------------------------------------------
# STRUCTURES
# -----------------------------------------------------------------------------

"""
    MetaICVIOpts()

Meta-ICVI module options.

# Examples
```julia-repl
julia> MetaICVIOpts()
```
"""
@with_kw mutable struct MetaICVIOpts @deftype Int
    # Which scikitlearn classifier to load
    classifier_selection::Symbol = :SGDClassifier
    # Scikitlearn classifier keyword arguments
    classifier_opts::NamedTuple = (loss="log", max_iter=30)
    # Size of ICVI window: [1, infty)
    icvi_window = 5; @assert icvi_window >= 1
    # Size of correlation window: [1, infty)
    correlation_window = 5; @assert correlation_window >= 1
    # Number of rocket kernels: [1, infty)
    n_rocket = 5; @assert n_rocket >= 1
    # Rocket file location
    rocket_file::String = module_dir("data", "models", "rocket.jld2")
    # Classifier file location
    classifier_file::String = module_dir("data", "models", "classifier.jld")
    # Display flag
    display::Bool = true
    # Flag to fail if any file is missing (rather than creating new objects)
    fail_on_missing::Bool = false
end # MetaICVIOpts

"""
    MetaICVIModule

Stateful information for a single MetaICVI module.

# Fields
- `opts::MetaICVIOpts`: options for construction.
- `cvis::Vector{AbstractCVI}`: list of cvis used for computing the CVIs.
- `criterion_values::RealVector`: list of outputs of the cvis used for computing correlations.
- `correlations::RealVector`: list of outputs of the rank correlations.
- `features::RealVector`: list of outputs of the rocket feature kernels.
- `rocket::RocketModule`: time-series random feature kernels module.
- `classifier::MetaICVIClassifier`: ScikitLearn classifier.
- `performance::RealFP`: final output of the most recent the Meta-ICVI step.
- `is_pretrained::Bool`: internal flag for if the classifier is trained and ready for inference.
"""
mutable struct MetaICVIModule
    opts::MetaICVIOpts
    cvis::Vector{AbstractCVI}
    criterion_values::Vector{RealVector}
    correlations::RealVector
    features::RealVector
    rocket::RocketModule
    classifier::MetaICVIClassifier
    performance::RealFP
    probabilities::RealVector
    is_pretrained::Bool
end # MetaICVIModule

"""
    MetaICVIModule(opts::MetaICVIOpts)

Instantiate a MetaICVIModule with given options.

# Arguments
- `opts::MetaICVIOpts`: options struct for the MetaICVI object.
"""
function MetaICVIModule(opts::MetaICVIOpts)
    # Load the correct classifier from the opts
    @eval @sk_import linear_model: $(opts.classifier_selection)

    # Construct the CVIs
    cvis = [
        PS(),
        GD43()
    ]

    # Initialize the empty vectors for each criterion value
    cvi_values = [Array{RealFP}(undef, 0) for _ = 1:length(cvis)]

    # Construct the rocket kernels
    if isfile(opts.rocket_file)
        # If we have a file, load the module
        @info "Loading rocket kernels: $(opts.rocket_file)"
        rocket_module = load_rocket(opts.rocket_file)
        # Correct the expected number of kernels if necessary
        if opts.n_rocket != length(rocket_module.kernels)
            @warn "Provided incorrect number of kernels with loaded file, correcting MetaICVIModule options."
            opts.n_rocket = length(rocket_module.kernels)
        end
    else
        # Otherwise, construct a module
        @warn "Missing/incorrect path for the rocket module file."
        @info "Constructing a new rocket module"
        rocket_module = RocketModule(opts.correlation_window, opts.n_rocket)
    end

    # Load the classifier
    if isfile(opts.classifier_file)
        @info "Loading classifier: $(opts.classifier_file)"
        # classifier = JLD.load(opts.classifier_file, "learner")
        classifier = load_classifier(opts.classifier_file)
        is_pretrained = true
    # If we didn't get a valid file, check for errors or construct a new object
    else
        error_msg = "Missing/incorrect path for pretrained classifier file."
        # If we want to fail on a missing file, throw an error
        if opts.fail_on_missing
            error(error_msg)
        # Otherwise, put up a warning and create a new classifier.
        else
            # If a file is even provided, construct a classifier
            if !isempty(opts.classifier_file)
                @warn error_msg
                classifier = construct_classifier(opts)
                is_pretrained = false
            # Otherwise, throw an error that no filename is provided
            else
                error("No filename provided for classifier saving/loading.")
            end
        end
    end

    # Construct and return the module
    return MetaICVIModule(
        opts,                       # opts
        cvis,                       # cvis
        cvi_values,                 # criterion_values
        Array{RealFP}(undef, 0),    # correlations
        Array{RealFP}(undef, 0),    # features
        rocket_module,              # rocket
        classifier,                 # classifier
        0.0,                        # performance
        zeros(3),                   # probabilities
        is_pretrained               # is_pretrained
    )
end # MetaICVIModule(opts::MetaICVIOpts)

"""
    MetaICVIModule()

Default constructor for the MetaICVIModule.
"""
function MetaICVIModule()
    # Create the default options
    opts = MetaICVIOpts()
    # Return the Meta-ICVI module constructed with the default options
    return MetaICVIModule(opts)
end # MetaICVIModule()

# -----------------------------------------------------------------------------
# METHODS
# -----------------------------------------------------------------------------

"""
    Base.show(io::IO, metaicvi::MetaICVIModule)

Display a metaicvi module to the command line.

# Arguments
- `io::IO`: default io stream.
- `metaicvi::MetaICVIModule`: metaicvi object about which to display info.
"""
function Base.show(io::IO, metaicvi::MetaICVIModule)
    compact = get(io, :compact, false)
    if compact
        # print(metaicvi.opts)
        print("MetaICVIModule")
    else
        # show(io, metaicvi.opts)
        pretrain_display_flag = MetaICVI.is_pretrained(metaicvi)
        opts = metaicvi.opts
        print(
            """
            MetaICVIModule:
            \tperformance: $(metaicvi.performance)
            \tpretrained: $pretrain_display_flag
            Options:
            \tclassifier: $(opts.classifier_selection)
            \tclassifier_opts: $(opts.classifier_opts)
            \ticvi_window: $(opts.icvi_window)
            \tcorrelation_window: $(opts.correlation_window)
            \tn_rocket: $(opts.n_rocket)
            \trocket_file: $(opts.rocket_file)
            \tclassifier_file: $(opts.classifier_file)
            Flags:
            \tdisplay: $(opts.display)
            \tfail_on_missing: $(opts.display)
            """
        )
    end
end # Base.show(io::IO, metaicvi::MetaICVIModule)

"""
    construct_classifier(opts::MetaICVIOpts)

Construct a new classifier for the MetaICVI module with metaprogramming.

# Arguments
- `opts::MetaICVIOpts`: options containing the classifier type and options for instantiation.
"""
function construct_classifier(opts::MetaICVIOpts)
    @info "Constructing a new classifier"
    @eval classifier = $(opts.classifier_selection)(;$(opts.classifier_opts)...)
    return MetaICVIClassifier(classifier)
end # construct_classifier(opts::MetaICVIOpts)

"""
    safe_save_classifier(metaicvi::MetaICVIModule)

Error handle saving of the metaicvi classifier.

# Arguments
- `metaicvi::MetaICVIModule`: metaicvi module containing the classifier and path for saving.
"""
function safe_save_classifier(metaicvi::MetaICVIModule)
    # If we specified a file but none was there, then save to that file
    if !isempty(metaicvi.opts.classifier_file)
        @info "Saving trained classifier"
        save_classifier(metaicvi.classifier, metaicvi.opts.classifier_file)
    else
        error("No filename provided for classifier file saving/loading.")
    end
end # safe_save_classifier(metaicvi::MetaICVIModule)

"""
    safe_save_rocket(metaicvi::MetaICVIModule)

Error handle saving of the metaicvi rocket kernels.

# Arguments
- `metaicvi::MetaICVIModule`: metaicvi module containing the classifier and path for saving.
"""
function safe_save_rocket(metaicvi::MetaICVIModule)
    # If we specified a file but none was there, then save to that file
    if !isempty(metaicvi.opts.rocket_file)
        @info "Saving current rocket kernels"
        save_rocket(metaicvi.rocket, metaicvi.opts.rocket_file)
    else
        error("No filename provided for rocket file saving/loading.")
    end
end # safe_save_rocket(metaicvi::MetaICVIModule)

"""
    load_classifier(filepath::String)

Load the classifier at the filepath.

# Arguments
- `filepath::String`: location of the classifier .jld file.
"""
function load_classifier(filepath::String)
    return MetaICVIClassifier(JLD.load(filepath, "classifier"))
    # return BSON.load(filepath, @__MODULE__)["classifier"]
end # load_classifier(filepath::String)

"""
    save_classifier(classifier::MetaICVIClassifier, filepath::String)

Save the classifier at the filepath.

# Arguments
- `classifier::MetaICVIClassifier`: classifier object to save.
- `filepath::String`: name/path to save the classifier .jld file.
"""
function save_classifier(classifier::MetaICVIClassifier, filepath::String)
    JLD.save(filepath, "classifier", classifier)
    # bson(filepath, Dict("classifier" => classifier))
end # save_classifier(classifier::MetaICVIClassifier, filepath::String)

"""
    train_and_save(metaicvi::MetaICVIModule, x::RealMatrix, y::IntegerVector)

Train the classifier on x/y and save the kernels and classifier.

# Arguments
- `metaicvi::MetaICVIModule`: metaicvi module to save with.
- `x::RealMatrix`: features to train on.
- `y::IntegerVector`: correct/over/under partition targets.
"""
function train_and_save(metaicvi::MetaICVIModule, x::RealMatrix, y::IntegerVector)
    # Create a new classifier
    # TODO: option to train/retrain loaded/serialized models
    classifier = construct_classifier(metaicvi.opts)
    metaicvi.classifier = classifier

    # Train the classifier
    @info "Training classifier"
    fit!(metaicvi.classifier, x, y)
    metaicvi.is_pretrained = true

    # Save the rocket kernels used
    safe_save_rocket(metaicvi)
    # Save the classifier used
    safe_save_classifier(metaicvi)
end # train_and_save(metaicvi::MetaICVIModule, x::RealMatrix, y::IntegerVector)

"""
    get_icvis(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

Compute and store the icvi criterion values.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
- `sample::RealVector`: the sample used for clustering.
- `label::Integer`: the label prescribed to the sample by the clustering algorithm.
"""
function get_icvis(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)
    # Update all of the cvis incrementally
    for ix = 1:length(metaicvi.cvis)
        # Compute and push the criterion value
        value = get_icvi!(metaicvi.cvis[ix], sample, label)
        push!(metaicvi.criterion_values[ix], value)
        # FIFO the list to size
        while length(metaicvi.criterion_values[ix]) > metaicvi.opts.icvi_window
            popfirst!(metaicvi.criterion_values[ix])
        end
    end
end # get_icvis(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

"""
    get_correlations(metaicvi::MetaICVIModule)

Compute and store the rank correlations from the cvi values.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
"""
function get_correlations(metaicvi::MetaICVIModule)
    # If the cvi window is big enough, compute the correlations
    if length(metaicvi.criterion_values[1]) >= metaicvi.opts.icvi_window
        # Get the spearman correlation
        correlation = corspearman(metaicvi.criterion_values[1], metaicvi.criterion_values[2])
        # Sanitize a potential NaN response
        # metaicvi.performance = isequal(performance, NaN) ? 0 : performance
        push!(metaicvi.correlations, correlation)
        # FIFO the list to size
        while length(metaicvi.correlations) > metaicvi.opts.correlation_window
            popfirst!(metaicvi.correlations)
        end
    end
end # get_correlations(metaicvi::MetaICVIModule)

"""
    get_rocket_features(metaicvi::MetaICVIModule)

Compute and store the rocket features.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
"""
function get_rocket_features(metaicvi::MetaICVIModule)
    # If there are enough correlations, compute compute the meta-icvi value
    if length(metaicvi.correlations) >= metaicvi.opts.correlation_window
        metaicvi.features = apply_kernels(metaicvi.rocket, metaicvi.correlations)[:, 1]
    else
        metaicvi.features = zeros(metaicvi.opts.n_rocket)
    end
end # get_rocket_features(metaicvi::MetaICVIModule)

"""
    is_pretrained(metaicvi::MetaICVIModule)

Checks if the classifier is pretrained to permit inference.

# Arguments
- `metaicvi::MetaICVIModule`: metaicvi module containing the classifier to check.
"""
function is_pretrained(metaicvi::MetaICVIModule)
    # Check if the model is pretrained with the isfit function
    # return ScikitLearn.Utils.isfit(metaicvi.classifier)
    return metaicvi.is_pretrained
end # is_pretrained(metaicvi::MetaICVIModule)

"""
    get_probability(metaicvi::MetaICVIModule)

Compute and store the metaicvi value from the classifier.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
"""
function get_probability(metaicvi::MetaICVIModule)
    # If we have previously computed features
    if !isempty(metaicvi.features) && is_pretrained(metaicvi)
        # Compute the class 'probabilities'
        probs = predict_proba(metaicvi.classifier, transpose([metaicvi.features]))
        metaicvi.probabilities = probs[:]

        # Store only the probability of correct partitioning
        metaicvi.performance = probs[1]
    else
        metaicvi.performance = 0.0
    end
end # get_probability(metaicvi::MetaICVIModule)

"""
    get_features(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

Compute only the features on the sample and label without classifier inference.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
- `sample::RealVector`: the sample used for clustering.
- `label::Integer`: the label prescribed to the sample by the clustering algorithm.
"""
function get_features(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)
    # Compute the icvi values
    get_icvis(metaicvi, sample, label)

    # Compute the rank correlations
    get_correlations(metaicvi)

    # Compute the rocket features
    get_rocket_features(metaicvi)
end # get_features(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

"""
    get_metaicvi(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

Compute and return the meta-icvi value.

# Arguments
- `metaicvi::MetaICVIModule`: the Meta-ICVI module.
- `sample::RealVector`: the sample used for clustering.
- `label::Integer`: the label prescribed to the sample by the clustering algorithm.
"""
function get_metaicvi(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)
    # If the sample was not misclassified
    if label > 0
        # Compute the rocket features
        get_features(metaicvi, sample, label)

        # Compute and store the
        get_probability(metaicvi)
    else
        # Default to 0
        metaicvi.performance = 0.0
    end

    return metaicvi.performance
end # get_metaicvi(metaicvi::MetaICVIModule, sample::RealVector, label::Integer)

"""
    get_cvi_data(data_file::String)

Get the cvi data specified by the data_file path.

# Arguments
- `data_file::String`: file containing clustered data for cvi processing.
"""
function get_cvi_data(data_file::String)
    # Parse the data
    data = readdlm(data_file, ',')
    data = permutedims(data)
    train_x = convert(Matrix{Float64}, data[1:2, :])
    train_y = convert(Vector{Integer}, data[3, :])

    return train_x, train_y
end # get_cvi_data(data_file::String)

"""
    get_training_features(metaicvi::MetaICVIModule, data_path::String)

# Arguments

"""
function get_training_features(metaicvi::MetaICVIModule, data_path::String)
    # Point to data
    data_dir(args...) = joinpath(data_path, args...)

    # Load the data
    correct_x, correct_y = get_cvi_data(data_dir("correct_partition.csv"))
    under_x, under_y = get_cvi_data(data_dir("under_partition.csv"))
    over_x, over_y = get_cvi_data(data_dir("over_partition.csv"))

    # Package the data conveniently
    data = Dict(
        "correct" => Dict(
            "x" => correct_x,
            "y" => correct_y
        ),
        "under" => Dict(
            "x" => under_x,
            "y" => under_y
        ),
        "over" => Dict(
            "x" => over_x,
            "y" => over_y
        )
    )

    # Create the target containers
    data_lengths = [length(correct_y), length(under_y), length(over_y)]
    offset_lengths = [0, length(correct_y), length(under_y)]
    data_length = sum(data_lengths)
    features_data = zeros(metaicvi.opts.n_rocket, data_length)
    features_targets = zeros(Int, data_length)

    # Mapping of datatype to numeral target for classification
    type_to_num = Dict(
        "correct" => 1,
        "under" => 2,
        "over" => 3
    )

    # Itereate over all data to get features
    for (type, subdata) in data
        # Get the offset directly from the type mapping and explicit definition
        data_offset = offset_lengths[type_to_num[type]]

        # Iterate over all entries of the type of partitioning (i.e., correct, under, and over)
        @showprogress for i = 1:length(subdata["y"])
            # Extract the sample and label
            sample = subdata["x"][:, i]
            label = subdata["y"][i]

            # Compute and retrieve the features
            get_features(metaicvi, sample, label)
            features = metaicvi.features

            # Save the results
            features_data[:, i + data_offset] = features
            features_targets[i + data_offset] = type_to_num[type]
        end
    end

    return features_data, features_targets
end # get_training_features(metaicvi::MetaICVIModule, data_path::String)


text:
using Grid
using Cubature
using Distributions
using Iterators

iters = 0
function f(x)
    global iters
    iters += 1
    prod(pdf(Poisson(3.0), floor(x)))
end

function gen_grid(mins, maxs, N=30)
    D = length(mins)
    stepsizes = ceil(float(maxs .- mins) ./N)
    ranges = [mins[i]:stepsizes[i]:maxs[i] for i=1:D]
    grid = hcat(map(collect, collect(product(ranges...)))...)'
    return map(length,ranges), stepsizes, grid
end

const tol = 0.001
const upp = 100

println("####################")
@time val = hquadrature(f, [0], [upp], abstol=tol)
println(val)
println("1D used $iters iters")

println("####################")
iters = 0
@time val = hcubature(f, zeros(2), zeros(2).+upp, abstol=tol)
println(val)
println("2D used $iters iters")

println("####################")
iters = 0
@time val = hcubature(f, zeros(3), zeros(3).+upp, abstol=tol)
println(val)
println("3D used $iters iters")

println("####################")
println("Now using Grid...")
println("####################")

valsi = nothing
iiters = 0
function fi(x)
    global iters, iiters
    v,g,h = valgradhess(valsi,x...)
    if abs(prod(g)) > 1e-15 || abs(det(h)) > 1e-15
        f(x)
    else
        iiters += 1
        valsi[x...]
    end
end

N = 20

println("####################")
_,_,pts = gen_grid([0],[upp], N)
@time vals = [f(pt) for pt in pts]
valsi = InterpGrid(vals, BCnearest, InterpQuadratic)
iters = 0
iiters = 0
@time val = hquadrature(fi, [0], [upp], abstol=tol)
println(val)
println("1D used $iters actual iters, $iiters interp iters")

println("####################")
lens,_,pts = gen_grid(zeros(2), zeros(2).+upp, N)
@time vals = [f(pts[i,:]) for i=1:size(pts,1)]
vals = reshape(vals, lens...)
valsi = InterpGrid(vals, BCnearest, InterpQuadratic)
iters = 0
iiters = 0
@time val = hcubature(fi, zeros(2), zeros(2).+upp, abstol=tol)
println(val)
println("2D used $iters actual iters, $iiters interp iters")

#println("####################")
#lens,_,pts = gen_grid(zeros(3), zeros(3).+upp, N)
#@time vals = [f(pts[i,:]) for i=1:size(pts,1)]
#vals = reshape(vals, lens...)
#valsi = InterpGrid(vals, BCnearest, InterpQuadratic)
#iters = 0
#iiters = 0
#@time val = hcubature(fi, zeros(3), zeros(3).+upp, abstol=tol)
#println(val)
#println("3D used $iters actual iters, $iiters interp iters")

hess = reshape([det(valgradhess(valsi, pts[i,:]...)[3]) for i=1:size(pts,1)], lens...);
grad = reshape([prod(valgrad(valsi, pts[i,:]...)[2]) for i=1:size(pts,1)], lens...);

ext: jl


text:
"""
Created on Wed Apr 22 2020

@author: Yoann Pradat

Parameters for nmf algorithms in the subfolder.
"""

"""
    CParams([,α, max_iter, stopping_crit, stopping_tol, α_H, α_W, l₁ratio_H, l₁ratio_W, μ_H, μ_W, α_sparse_H,
    α_sparse_W, ω, verbose])

Parameters passed to `nmf_α`.

# Arguments
- `α::Real`               : parameter of the α-divergence
- `max_iter::Integer`     : maximum number of iterations (in main procedure)
- `stopping_crit::Symbol` : criterion for stopping iterations. One of :none or :cost
- `stopping_tol::Real`    : tolerance for the stopping crit
- `scale_W::Bool`         : scale the components of W to have unit l₁ norm
- `α_H::Real`             : constant that multiplies the regularization term on H
- `α_W::Real`             : constant that multiplies the regularization term on W
- `l₁ratio_H::Real`       : l₁/l₂ regularization balance for H (in [0; 1]). 1 is l₁ only, 0 is l₂ only
- `l₁ratio_W::Real`       : l₁/l₂ regularization balance for W (in [0; 1]). 1 is l₁ only, 0 is l₂ only
- `μ_H::Real`             : constant that multiplies the orthogonalisation term on H
- `μ_W::Real`             : constant that multiplies the orthogonalisation term on W
- `α_sparse_H::Real`      : (1+α_sparse_H) is used as global exponent in the update of H
- `α_sparse_W::Real`      : (1+α_sparse_W) is used as global exponent in the update of W
- `ω::Real`               : overrelaxation parameter, typically in (0,2)
- `verbose::Bool`         : whether to show intermediate information
"""
mutable struct CParams
    α::Real
    max_iter::Integer
    stopping_crit::Symbol
    stopping_tol::Real
    α_H::Real
    α_W::Real
    l₁ratio_H::Real
    l₁ratio_W::Real
    μ_H::Real
    μ_W::Real
    α_sparse_H::Real
    α_sparse_W::Real
    ω::Real
    verbose::Bool

    function CParams(;
        α::Real               = 2,
        max_iter::Integer     = 1000,
        stopping_crit::Symbol = :none,
        stopping_tol::Real    = 1e-5,
        α_H::Real             = 0,
        α_W::Real             = 0,
        l₁ratio_H::Real       = 0,
        l₁ratio_W::Real       = 0,
        μ_H::Real             = 0,
        μ_W::Real             = 0,
        α_sparse_H::Real      = 0,
        α_sparse_W::Real      = 0,
        ω::Real               = 1,
        verbose::Bool         = false)

        stopping_criteria = [:none, :cost]

        valid_iter = max_iter > 1
        valid_stopping_crit = stopping_crit in stopping_criteria
        valid_stopping_tol = stopping_tol > 0
        valid_α_H = α_H >= 0
        valid_α_W = α_W >= 0
        valid_l₁ratio_H = l₁ratio_H >= 0 && l₁ratio_H <= 1
        valid_l₁ratio_W = l₁ratio_W >= 0 && l₁ratio_W <= 1
        valid_μ_H = μ_H >= 0
        valid_μ_W = μ_W >= 0
        valid_α_sparse_H = α_sparse_H >= 0
        valid_α_sparse_W = α_sparse_W >= 0
        valid_ω = ω > 0

        valid_iter  || throw(ArgumentError("max_iter must be at least 2"))
        valid_stopping_crit || throw(ArgumentError(string("stopping_crit should be one of ", stopping_criteria)))
        valid_stopping_tol || throw(ArgumentError("stopping_tol must be > 0"))
        valid_α_H || throw(ArgumentError("α_H must be > 0"))
        valid_α_W || throw(ArgumentError("α_W must be > 0"))
        valid_l₁ratio_H || throw(ArgumentError("l₁ratio_H must be in [0,1]"))
        valid_l₁ratio_W || throw(ArgumentError("l₁ratio_W must be in [0,1]"))
        valid_μ_H || throw(ArgumentError("μ_H must be > 0"))
        valid_μ_W || throw(ArgumentError("μ_W must be > 0"))
        valid_α_sparse_H || throw(ArgumentError("α_sparse_H must be > 0"))
        valid_α_sparse_W || throw(ArgumentError("α_sparse_W must be > 0"))
        valid_ω || throw(ArgumentError("ω must be > 0"))

        new(α, max_iter, stopping_crit, stopping_tol, α_H, α_W, l₁ratio_H, l₁ratio_W, μ_H, μ_W, α_sparse_H, α_sparse_W, ω, verbose)
    end
end


text:
#
# Data archiving for state estimators.
#
using Plots

"""
    EstimatorHistory(target_estimator::Estimator)

Estimation history archival type.
"""
immutable EstimatorHistory{T<:AbstractFloat, S<:AbstractUncertainState{T}}
    states::Vector{S}
    residuals::Vector{S}

    function EstimatorHistory(target_estimator::Estimator{T,S}) where {T,S}
        new{T,S}([], [])
    end
end

"""
    vectorize_state_history(state_history::Vector)
"""
function vectorize_state_history(state_history::Vector)
    state_length = length(state_history[1].x)
    num_states = length(state_history)

    times = fill(0*state_history[1].t, num_states)
    states = fill(0*state_history[1].x[1], (num_states, state_length))
    sigmas = fill(0*state_history[1].x[1], (num_states, state_length))

    for i = 1:num_states
        times[i] .= state_history[i].t
        states[i,:] .= state_history[i].x
        sigmas[i,:] .= sqrt.(diag(state_history[i].P))
    end

    return times, states, sigmas
end

"""
    plot_state_history(archive::EstimatorHistory)

Plot the state history in an estimator history archive.
"""
function plot_state_history(archive::EstimatorHistory)
    t, data, cov = vectorize_state_history(archive.states)
    scatter(t, data, markersize=1, markercolor="black",
            layout=(length(archive.states[1].x), 1), legend=false)
    plot!(t, data+3*cov, linestyle=:dot, color="gray")
    plot!(t, data-3*cov, linestyle=:dot, color="gray")
end

"""
    plot_residuals(archive::EstimatorHistory)

Plot the residuals in an estimator history archive.
"""
function plot_residuals(archive::EstimatorHistory)
    t, data, cov = vectorize_state_history(archive.residuals)
    scatter(t, data, markersize=1, markercolor="black",
            layout=(length(archive.states[1].x), 1), legend=false)
    plot!(t, data+3*cov, linestyle=:dot, color="gray")
    plot!(t, data-3*cov, linestyle=:dot, color="gray")
end


text:
struct OpenGLGUIBackend <: AbstractGUIBackend
    function OpenGLGUIBackend()
        GLFW.Init()
        new()
    end
end

mutable struct OpenGLWindow <: AbstractGUIWindow
    screen::GLWindow.Screen
    events::Queue{Any}
    mouse_position::Tuple{Float64,Float64}

    function OpenGLWindow(w,h,name)
        screen =  GLWindow.Screen(name, resolution=(w,h))
        GLVisualize.add_screen(screen)
        GLWindow.add_complex_signals!(screen) #add the drag events and such
        GLFW.MakeContextCurrent(GLWindow.nativewindow(screen))

        events = Queue(Any)
        win = new(screen,events,(0.,0.))

        GLFW.SetKeyCallback(screen.glcontext.window, (_, key, scancode, action, mods) -> begin
            enqueue!(events,("Key", key, scancode, action, mods))
        end)
        GLFW.SetMouseButtonCallback(screen.glcontext.window, (_, button, action, mods) -> begin
            enqueue!(events,("Button", button, action, mods))
        end)
        GLFW.SetCursorPosCallback(screen.glcontext.window, (_, x, y) -> begin
            win.mouse_position = (x,y)
            enqueue!(events,("Cursor", x, y))
        end)

        win
    end
end

width(win::OpenGLWindow)  = GLWindow.widths(win.screen)[1]
height(win::OpenGLWindow) = GLWindow.widths(win.screen)[2]

window(b::OpenGLGUIBackend,w,h,name) = OpenGLWindow(w,h,name)
window(b::OpenGLGUIBackend,w,h) =  window(b,w,h,"")

destroy(win::OpenGLWindow) = GLFW.DestroyWindow(win.screen.glcontext.window)

function update!(win::OpenGLWindow)

    GLWindow.render_frame(win.screen)
    GLWindow.swapbuffers(win.screen)
    GLWindow.poll_glfw()
    yield()
end

function rectangle(win::OpenGLWindow,x,y,w,h,color=RGB(0.9,0.9,0.9))
    _view( visualize(GLVisualize.SimpleRectangle{Float32}(x,height(win)-y-h,w,h), color=color), win.screen )
end

function clear!(win::OpenGLWindow,color=RGB(1.0,1.0,1.0))

    ModernGL.glClearColor(red(color), green(color), blue(color), 1)
    GLWindow.clear_all!(win.screen)
    GLWindow.empty!(win.screen)
end

function text(win::OpenGLWindow,x,y,str::String,color=RGB(0.0,0.0,0.0))
    _view( visualize(
        str, color=color,
        model = GLAbstraction.translationmatrix(GeometryTypes.Vec3f0(x,height(win)-y,0)),
    ), win.screen )
end

mouse_position(win::OpenGLWindow) = win.mouse_position

function events(win::OpenGLWindow)
    out = []
    ev = win.events
    while !isempty(ev)
        push!(out,dequeue!(ev))
    end
    out
end

##


text:
"""
```
binarize(Intermodes(), img)
```

Under the assumption that the image histogram is bimodal the histogram is
smoothed using a length-3 mean filter until two modes remain. The binarization
threshold is then set to the average value of the two modes.

# Output

Returns the binarized image as an `Array{Gray{Bool},2}`.

# Arguments

The function argument is described in more detail below.

##  `img`

An `AbstractArray` representing an image. The image is automatically converted
to `Gray` in order to construct the requisite graylevel histogram.

# Example

Binarize the "cameraman" image in the `TestImages` package.

```julia
using TestImages, ImageBinarization

img = testimage("cameraman")
img_binary = binarize(Intermodes(), img)
```

## Reference

C. A. Glasbey, “An Analysis of Histogram-Based Thresholding Algorithms,” CVGIP: Graphical Models and Image Processing, vol. 55, no. 6, pp. 532–537, Nov. 1993. [doi:10.1006/cgip.1993.1040](https://doi.org/10.1006/cgip.1993.1040)
"""
function binarize(algorithm::Intermodes,  img::AbstractArray{T,2}) where T <: Colorant
  img₀₁ = zeros(Gray{Bool}, axes(img))
  edges, counts = build_histogram(img,  256)
  t = find_threshold(HistogramThresholding.Intermodes(), counts[1:end], edges)
  for i in CartesianIndices(img)
    img₀₁[i] = img[i] < t ? 0 : 1
  end
  img₀₁
end


text:
using NESSie.Format
using NESSie.Format: readmsms_nodes, readmsms_elements

testfiles = ((mktemp()..., mktemp()..., (0,0)), # empty file
             (mktemp()..., mktemp()..., (3,2))) # dummy file

write(testfiles[2][2], """
# MSMS solvent excluded surface vertices
#vertex #sphere density probe_r
  x  x  x  x
  0.0 0.0 0.0 x x x x x x
  1.0 2.0 3.0 x x x x x x
  -2.00001 1.337 42.0 x x x x x x
""")

write(testfiles[2][4], """
# MSMS solvent excluded surface faces
#faces  #sphere density probe_r
  x  x  x  x
  1 2 3 x x
  2 3 1 x x
""")

@testset "readmsms" begin
    try
        for (fnamev, fv, fnamef, ff, len) in testfiles, fn in ("single", "bulk"), T in testtypes
            seekstart(fv)
            seekstart(ff)
            nodes = (); elements = ()
            if fn == "single"
                # check single methods
                nodes = readmsms_nodes(fv, T)
                elements = readmsms_elements(ff, nodes)
            elseif fn == "bulk"
                # check bulk method (implicit single methods with atstart=true)
                model = readmsms(fv, ff, T)
                nodes, elements = (model.nodes, model.elements)
            end
            # check return types
            @test typeof(nodes) == Vector{Vector{T}}
            @test typeof(elements) == Vector{Triangle{T}}
            # check lengths
            @test length(nodes) == len[1]
            @test length(elements) == len[2]

            # check content
            if fnamev == testfiles[2][1]
                @test nodes[1] == T[0, 0, 0]
                @test nodes[2] == T[1, 2, 3]
                @test nodes[3] == T[-2.00001, 1.337, 42]
                @test elements[1].v1 === nodes[1]
                @test elements[1].v2 === nodes[2]
                @test elements[1].v3 === nodes[3]
                @test elements[2].v1 === nodes[2]
                @test elements[2].v2 === nodes[3]
                @test elements[2].v3 === nodes[1]
            end
        end
    finally
        [begin close(fv); rm(fnamev); close(ff); rm(fnamef) end for (fnamev, fv, fnamef, ff, _) in testfiles]
    end
end


text:
"""
    maximum_error(mw::MWState)

Compute maximum error of synthetic data on the query set. Result is not
differentially private.
"""
function maximum_error(mw::MWState)
    norm(evaluate(mw.queries, mw.synthetic) - mw.real_answers, Inf)
end

"""
    mean_squared_error(mw::MWState)

Compute mean_squared_error of synthetic data on the query set. Result is not
differentially private.
"""
function mean_squared_error(mw::MWState)

    errors = evaluate(mw.queries, mw.synthetic) - mw.real_answers
    (norm(errors)^2)/length(errors)
end

function total_squared_error(mw::MWState)

    errors = evaluate(mw.queries, mw.synthetic) - mw.real_answers
    (norm(errors)^2)/length(errors)
end

"""
   kl_divergence(a, b)

Compute KL-divergence between two normalized histograms a and b.
"""
function kl_divergence(a::Array{Float64, 1}, b::Array{Float64, 1})
    r = 0.0
    for i = 1 : length(a)
        @inbounds ai = a[i]
        @inbounds bi = b[i]
        if ai > 0 && bi > 0
            r += ai * log(2, ai / bi)
        end
    end
    r
end

function kl_divergence_error(mw::MWState)

    #error = kl_divergence(mw.synthetic.weights, mw.real.weights)
    error = kl_divergence(mw.real.weights, mw.synthetic.weights)
    error
end


text:
using Dates
include("util.jl")

lines = getlines("inputs/4.txt")

function get_minute(line)
    m = match(r"^\[\d+-\d+-\d+ \d+:(\d+)\]", line)
    return parse(Int64, m[1])
end

function parse_guard_records(lines)
    m = match(r"^.*Guard #(\d+) begins shift$", lines[1])
    id = m[1]
    asleep::Array{Bool} = []
    curr_minute = 0
    for line in lines[2:end]
        if occursin("falls asleep", line)
            m = get_minute(line)
            for _ in curr_minute:m-1
                push!(asleep, false)
            end
            curr_minute = m
        elseif occursin("wakes up", line)
            m = get_minute(line)
            for _ in curr_minute:m-1
                push!(asleep, true)
            end
            curr_minute = m
        end
    end
    for _ in curr_minute:59
        push!(asleep, false)
    end
    return (parse(Int64, id), asleep)
end

function parse_lines(lines)
    lines = order_lines(lines)
    guard_info = [lines[1]]
    records = []
    for line in lines[2:end]
        if occursin("begins shift", line)
            push!(records, parse_guard_records(guard_info))
            guard_info = [line]
        else
            push!(guard_info, line)
        end
    end
    push!(records, parse_guard_records(guard_info))
    return records
end

function order_lines(lines)
    lines_dict::Dict{DateTime, String} = Dict()
    for line in lines
        pattern = "[yyyy-mm-dd HH:MM"
        datepart = split(line, "] ")[1]
        date = Dates.DateTime(datepart, pattern)
        lines_dict[date] = line
    end
    return [i[2] for i in sort(collect(lines_dict), by=x -> x[1])]
end

function time_asleep_fn(arr::Array{Bool})::Int64
    x = map(x -> begin
        if x
            return 1
        else
            return 0
        end
    end, arr)
    return sum(x)
end

function guard_asleep_at_minute_n(guard_id::Int64, results)
    guard_records = filter(x -> x[1]==guard_id, results)
    days_asleep_at_minute = zeros(60)
    for record in guard_records
        for minute in 0:59
            if record[2][minute+1]
                days_asleep_at_minute[minute+1] += 1
            end
        end
    end
    return (maximum(days_asleep_at_minute), guard_id, argmax(days_asleep_at_minute)-1)
end

function part_1(results)
    time_asleep::Dict{Int64,Int64} = Dict()
    for guard in results
        time_asleep[guard[1]] = 0
    end

    for result in results
        time_asleep[result[1]] += time_asleep_fn(result[2])
    end

    guard_asleep_most = sort(collect(time_asleep), by=x->-x[2])[1][1]

    guard_records = filter(x -> x[1]==guard_asleep_most, results)
    days_asleep_at_minute = zeros(60)
    for record in guard_records
        for minute in 0:59
            if record[2][minute+1]
                days_asleep_at_minute[minute+1] += 1
            end
        end
    end
    return guard_asleep_most*(argmax(days_asleep_at_minute)-1)
end

function part_2(results)
    part_2_results = []
    for guard_id in Set([r[1] for r in results])
        push!(part_2_results, guard_asleep_at_minute_n(guard_id, results))
    end
    a = sort(collect(part_2_results), by=x->-x[1])
    return a[1][2]*a[1][3]
end

results = parse_lines(lines)
time_asleep_list = part_1(results)
part_2(results)


text:
"""
    labeled_accuracy(prediction, gold)

Accuracy score for dependency arcs, including the labels.
"""
function labeled_accuracy end

"""
    unlabeled_accuracy(prediction, gold)

Accuracy score for dependency arcs, not including the labels.
"""
function unlabeled_accuracy end

const labelled_accuracy = labeled_accuracy
const unlabelled_accuracy = unlabeled_accuracy

function labeled_accuracy(tree::DependencyTree, gold::DependencyTree)
    correct, total = 0, 0
    for (t, g) in zip(tree, gold)
        t.head == g.head && t.label == g.label && (correct += 1)
        total += 1
    end
    return correct / total
end

function labeled_accuracy(trees, goldtrees)
    correct, total = 0, 0
    for (tree, gold) in zip(trees, goldtrees), (t, g) in zip(tree, gold)
        t.head == g.head && t.label == g.label && (correct += 1)
        total += 1
    end
    return correct / total
end

function unlabeled_accuracy(tree::DependencyTree, gold::DependencyTree)
    correct, total = 0, 0
    for (t, g) in zip(tree, gold)
        t.head == g.head && (correct += 1)
        total += 1
    end
    return correct / total
end

function unlabeled_accuracy(trees, goldtrees)
    correct, total = 0, 0
    for (tree, gold) in zip(trees, goldtrees), (t, g) in zip(tree, gold)
        t.head == g.head && (correct += 1)
        total += 1
    end
    return correct / total
end


text:
using StokesParticles, Test
import StokesParticles: is_inside_rect, Particle2D, Point2D

# test isInsideRect
sampleRect    = [0.0, 0.6, -0.2, 0.4]
pointIn       = Point2D(0.1,0.3)
pointOut      = Point2D(1.1,0.3)
pointOnEdge   = Point2D(0.6,0.3)
pointOnCorner = Point2D(0.6,0.4)
@test is_inside_rect(sampleRect,pointIn) == true
@test is_inside_rect(sampleRect,pointOut) == false
@test is_inside_rect(sampleRect,pointOnEdge) == true
@test is_inside_rect(sampleRect,pointOnCorner) == true

# test generateCellList
n_particles = 3
particleList1 = [Particle2D(Point2D(rand(),rand()),rand(),1) for i=1:n_particles]
TotalBounds = [0.0,1.0,0.0,1.0]
L = 0.1
cl = generate_cell_list(particleList1,TotalBounds,L)
sum = 0
for i=1:length(cl.cells)
  global sum += length(cl.cells[i].particleIDList)
end
@test n_particles == sum

# test updateCellList!
n_particles = 1_000
particleList2 = [Particle2D(Point2D(rand(),rand()),rand(),1) for i=1:n_particles]
TotalBounds = [0.0,1.0,0.0,1.0]
L = 0.1
cl = generate_cell_list(particleList2,TotalBounds,L)
particleList3 = [Particle2D(Point2D(rand(),rand()),rand(),1) for i=1:n_particles]
update_cell_list!(cl,particleList3)
sum = 0
for i=1:length(cl.cells)
  global sum += length(cl.cells[i].particleIDList)
end
@test n_particles == sum


text:
"""
    simulate a single path of drilling for section `i` given
        `uv` and `prob_uv`
"""
function simulate_lease!(simprim::SimulationPrimitives,
        lease::DrillLease, sim::SimulationDraw, weight::Number=1)

    simtmp = SimulationTmp(simprim)
    m = _model(DataDrill(lease))
    wp = statespace(m)
    θt = theta_drill(simprim)

    isfinite(_ψ2(sim)) || throw(DomainError("_ψ2(sim) = $(_ψ2(sim)) is not finite"))
    isfinite(_ψ1(sim)) || throw(DomainError("_ψ1(sim) = $(_ψ1(sim)) is not finite"))
    isfinite(weight)   || throw(DomainError("  weight = $(weight) is not finite"))

    # reset sa and sb, which hold the current state
    # we start out in period 1 at state0
    reset!(simtmp, first(_x(lease)))

    # iterate over calendar periods from the lease signing to the
    #    end of the price vector
    for (obs,zt) in LeaseCounterfactual(lease)

        # Define temporary scalars
        @declarezero Float64 d0 d1 epsdeq1 epsdgt1 Prdeq1 Prdgt1
        @declarezero Float64 _profit surp cost rev extension

        update_sparse_state_transition!(simtmp, obs, sim, θt)
        eur, ext_cost = update!(simtmp, obs, sim, θt)

        tday, tmrw = today_tomorrow(simtmp, zt)
        _Eeps = dot(tday, Eeps(simtmp))
        @assert isfinite(_Eeps) || throw(error("_Eeps = $_Eeps is not finite!"))

        for si in _x(obs):length(wp) #_x(obs) is state if never drilled
            PrTday = tday[si]
            @assert 0 <= PrTday <= 1 || PrTday ≈ 1

            if PrTday > 0
                action_probs = uview_col(Pprime(simtmp), si)
                @assert sum(action_probs) ≈ 1
                actions = actionspace(wp, si)
                dmx = _dmax(wp,si)
                numactions = dmx+1

                # add to extension if we are expiring
                if si == end_ex1(wp)
                    extension += ext_cost * first(action_probs) * PrTday
                end

                # ϵ cost shocks so we can add them to observed drilling cost later
                if dmx >= 1
                    Pr1 = action_probs[2]
                    epsdeq1 += PrTday * xlogx(Pr1)
                    Prdeq1  += PrTday * Pr1
                    if dmx >= 2
                        Pr2 = view(action_probs, 3:numactions)
                        twodmx = 2:dmx
                        epsdgt1 += PrTday * sum(xlogx.(Pr2) ./ twodmx )
                        Prdgt1  += PrTday * sum(Pr2)
                    end
                end

                # expected drilling for D=0 vs D=1
                expected_d = dot(actions, action_probs) * PrTday
                if si <= end_ex0(wp) # before final expiration?
                    d0 += expected_d
                else
                    d1 += expected_d
                end

                # Expected payoffs
                dp1space_t = OneTo(numactions) # 0:dmax
                _profit += dot(action_probs, view(profit(simtmp),    dp1space_t)) * PrTday
                surp    += dot(action_probs, view(surplus(simtmp),   dp1space_t)) * PrTday
                cost    += dot(action_probs, view(drillcost(simtmp), dp1space_t)) * PrTday
                rev     += dot(action_probs, view(revenue(simtmp),   dp1space_t)) * PrTday
            end # PrTday > 0
        end # states

        # update our tracking vectors
        ds = SharedSimulations(simprim)

        update!(ds, zt, uniti(lease), weight,
            d0, d1, _ψ2(sim), eur,
            epsdeq1, epsdgt1, Prdeq1, Prdgt1,
            _Eeps, _profit, surp, rev, cost, extension
        )

        # Full transition
        mul!(tmrw, Pprime(simtmp), tday)

        # Distribution Pr(Dₜ = D|T) for D ∈ 1:Dmax after actions in T (eg sₜ₊₁)
        if zt == Tstop(simprim)-1
            i = uniti(lease)
            @inbounds for D in 1:_Dmax(wp)
                ds.D_at_T[D,i]  += weight * sum( view(tmrw, s_of_D(wp, D)) )
            end
        end

    end # loop over time
end # function

function update!(ds::SharedSimulations, t::Int, i::Int, weight,
        d0, d1, ψ2, eur,
        epsdeq1, epsdgt1, Prdeq1, Prdgt1,
        Eeps, profit, surplus, revenue, drillcost, extension
    )

    eursq = eur^2
    eurcub = eursq*eur

    # drilling, eur, ψ conditional on D==0
    d0 *= weight

    ds.d0[t,i]       += d0
    ds.d0psi[t,i]    += d0 * ψ2
    ds.d0eur[t,i]    += d0 * eur
    ds.d0eursq[t,i]  += d0 * eursq
    ds.d0eurcub[t,i] += d0 * eurcub

    # drilling, eur, ψ conditional on D>0
    d1 *= weight

    ds.d1[t,i]       += d1
    ds.d1psi[t,i]    += d1 * ψ2
    ds.d1eur[t,i]    += d1 * eur
    ds.d1eursq[t,i]  += d1 * eursq
    ds.d1eurcub[t,i] += d1 * eurcub

    # expected average cost shocks
    ds.epsdeq1[t,i]  += epsdeq1 * weight
    ds.epsdgt1[t,i]  += epsdgt1 * weight
    ds.Prdeq1[t,i]   += Prdeq1  * weight
    ds.Prdgt1[t,i]   += Prdgt1  * weight

    # profitability stuff
    ds.Eeps[t,i]      += Eeps      * weight
    ds.profit[t,i]    += profit    * weight
    ds.surplus[t,i]   += surplus   * weight
    ds.revenue[t,i]   += revenue   * weight
    ds.drillcost[t,i] += drillcost * weight
    ds.extension[t,i] += extension * weight

end


text:
println("checking for libgtksourceview...")
@static if is_linux()
    try
        strip(readall(pipeline(`ldconfig -p`, `grep libgtksourceview-3`, `cut -d'>' -f2`)))
    catch
        run(`sudo apt-get install libgtksourceview-3.0-1`)
    end
end
@static if is_apple()
    if !isfile( Pkg.dir() * "/Homebrew/deps/usr/lib/libgtksourceview-3.0.dylib" )
        using Homebrew
        Homebrew.add("gtksourceview3")
    end
end


text:
using IRTools, Test
using IRTools: IR

function f(x)
    N = 4
    for i1 in 1:3    # single loop works without issue
        for i2 in 1:N   # needs to be a variable `1:4` works fine
        end
    end
    0.0 # same error with `x`
end

@test @code_ir(f(1)) isa IR

function f(a, b)
    u = 1
    while true
        if true
        end
    end
    f(u)
end

@test @code_ir(f(1, 2)) isa IR


text:
using Printf

function degree2compasspoint(d::Float64)
    majors = ("north", "east", "south", "west", "north", "east", "south", "west")
    quarter1 = ("N", "N by E", "N-NE", "NE by N", "NE", "NE by E", "E-NE", "E by N")
    quarter2 = map(p -> replace(p, "NE" => "EN"), quarter1)

    d = d % 360 + 360 / 64
    imajor, minor = divrem(d, 90)
    iminor = div(minor * 4, 45)
    imajor += 1
    iminor += 1
    p1, p2 = majors[imajor:imajor+1]
    q = p1 in ("north", "south") ? quarter1 : quarter2
    titlecase(replace(replace(q[iminor], 'N' => p1), 'E' => p2))
end

for i in 0:32
    d = i * 11.25
    i % 3 == 1 && (d += 5.62)
    i % 3 == 2 && (d -= 5.62)
    @printf("%2i %-17s %10.2f°\n", i % 32 + 1, degree2compasspoint(d), d)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
"""
Four times the size of the outermost cache (rule of thumb "laid down by Dr. Bandwidth").
"""
default_vector_length() = Int(4 * last(cachesize()) / sizeof(Float64))

_nthreads_string() = avxt() ? "@avxt" : string(nthreads())

function _run_kernels(
    copy,
    scale,
    add,
    triad;
    verbose=true,
    N=default_vector_length(),
    evals_per_sample=5,
    write_allocate=true,
)
    # initialize
    A, B, C, D, s = zeros(N), zeros(N), zeros(N), zeros(N), rand()

    α = write_allocate ? 24 : 16
    β = write_allocate ? 32 : 24

    f = t -> N * α * 1e-6 / t
    g = t -> N * β * 1e-6 / t

    # COPY
    t_copy = @belapsed copy($C, $A) samples = 10 evals = evals_per_sample
    bw_copy = f(t_copy)
    verbose && println("╟─ COPY:  ", round(bw_copy; digits=1), " MB/s")

    # SCALE
    t_scale = @belapsed scale($B, $C, $s) samples = 10 evals = evals_per_sample
    bw_scale = f(t_scale)
    verbose && println("╟─ SCALE: ", round(bw_scale; digits=1), " MB/s")

    # ADD
    t_add = @belapsed add($C, $A, $B) samples = 10 evals = evals_per_sample
    bw_add = g(t_add)
    verbose && println("╟─ ADD:   ", round(bw_add; digits=1), " MB/s")

    # TRIAD
    t_triad = @belapsed triad($A, $B, $C, $s) samples = 10 evals = evals_per_sample
    bw_triad = g(t_triad)
    verbose && println("╟─ TRIAD: ", round(bw_triad; digits=1), " MB/s")

    # statistics
    values = [bw_copy, bw_scale, bw_add, bw_triad]
    calc = f -> round(f(values); digits=1)

    return (median=calc(median), minimum=calc(minimum), maximum=calc(maximum))
end

"""
    benchmark(; kwargs...)

Runs a comprehensive STREAM benchmark. Returns the median, minimum, and maximum
of single- and multi-threaded measurements of the memory_bandwidth (in MB/s).
"""
function benchmark(; kwargs...)
    println("╔══╡ Single-threaded:")
    nt_single = _run_kernels(copy, scale, add, triad; kwargs...)
    println("╟─────────────────────")
    println("║ Median: ", nt_single.median, " MB/s")
    println("╚═════════════════════")
    println()
    println("╔══╡ Multi-threaded:")
    println("╠══╡ ($(_nthreads_string()) threads)")
    nt_multi = _run_kernels(
        copy_threaded, scale_threaded, add_threaded, triad_threaded; kwargs...
    )
    println("╟─────────────────────")
    println("║ Median: ", nt_multi.median, " MB/s")
    println("╚═════════════════════")
    println()
    return (single=nt_single, multi=nt_multi)
end

"""
    memory_bandwidth(; verbose=false, multithreading=nthreads() > 1, kwargs...)

Measure the memory bandwidth in megabytes per second (MB/s). Returns a named tuple
indicating the median, minimum, and maximum of the measurements in this order.
"""
function memory_bandwidth(; verbose=false, multithreading=nthreads() > 1, kwargs...)
    c = multithreading ? copy_threaded : copy
    s = multithreading ? scale_threaded : scale
    a = multithreading ? add_threaded : add
    t = multithreading ? triad_threaded : triad
    if verbose
        multithreading ? println("╔══╡ Multi-threaded:") : println("╔══╡ Single-threaded:")
        println("╠══╡ ($(_nthreads_string()) threads)")
    end
    nt = _run_kernels(c, s, a, t; verbose, kwargs...)
    if verbose
        println("╟─────────────────────")
        println("║ Median: ", nt.median, " MB/s")
        println("╚═════════════════════")
    end
    return nt
end

"""
    vector_length_dependence(; n=4, evals_per_sample=1) -> Dict

Measure the memory bandwidth for multiple vector lengths corresponding to
factors of the size of the outermost cache.
"""
function vector_length_dependence(; n=4, evals_per_sample=1)
    outer_cache_size = CpuId.cachesize()[end] / sizeof(Float64)
    Ns = floor.(Int, range(1, 4; length=n) .* outer_cache_size)
    membws = Dict{Int,Float64}()
    for (i, N) in pairs(Ns)
        m, _, _ = memory_bandwidth(; N=N, evals_per_sample=evals_per_sample)
        membws[N] = m
        println(i, ": ", N => m)
    end
    return membws
end


text:
using Flux

struct Embedding
    table
end

@Flux.treelike Embedding

(e::Embedding)(x) = e.table[x, :]

"""
    Embedding(num_vocab, embed_dim)

Construct an embedding layer for Flux models. This
is a mapping from the dimension of the model's
`num_vocab`, the number of unique vocabulary
terms in a corpus, to some embedding dimension, `embed_dim`.

# Examples
```julia-repl
julia> m = Embeddig(10, 5)

```
"""
Embedding(num_vocab, embed_dim) = Embedding(Flux.glorot_normal(num_vocab, embed_dim))

"""
Forward pass of the Embedding layer.

Given an array of size (input_dim, batch_size), this will
embed the information to (input_dim, embed_dim, batch_size).

# Examples
```julia-repl
julia> input = rand(1:100, (10, 32))
julia> m = Pilot.Embedding(100, 16)
julia> out = m(input)
julia> size(out)
(10, 16, 32)
```
"""
function (e::Embedding)(x::AbstractArray{T, 2}) where {T}
    out = e.table[x, :]
    permutedims(out, (1, 3, 2))
end


text:
module _Float

include("pi_iterative_gauss-legendre_float_simple.jl")
include("pi_iterative_gauss-legendre_float_inplace.jl")

end


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

using Markdown
using InteractiveUtils

# ╔═╡ c5b16d82-6281-11ec-250b-8be67db3a2e7
using Flux

# ╔═╡ b75114fe-f046-4681-9348-139318cf8ba8
xs = rand(Float32, 28, 28, 1, 50);

# ╔═╡ f90bf01b-06be-4685-9ad9-6dca0ac87b88
m = Chain(
	Conv((5,5), 1 => 7),
	MeanPool((5,5), pad=SamePad())
)

# ╔═╡ 73ce15dc-d40d-4df8-abe7-4e7b3019f8bc
# Chain also supports indexing and slicing, e.g. m[2] or m[1:end-1]. m[1:3](x) will calculate the output of the first three layers.

# ╔═╡ 596bc7ab-b3bc-4d02-b37c-04d34c3be108
m[1](xs) |> size

# ╔═╡ 2d0d9a40-6b92-49b2-ba07-bee22f8d3d0d
m[3](xs) |> size

# ╔═╡ fc8ec7eb-63ea-4f5c-9600-eba205eb3f70
m(xs) |> size

# ╔═╡ 0eb0d2b6-6c67-479b-a6f5-1a4fa7aa9e8e
begin
imgsize=(28,28,1)
nclasses=8
out_conv_size = (imgsize[1]÷4 - 3, imgsize[2]÷4 - 3, 16)
end

# ╔═╡ a79c9b40-56db-4056-bbf6-4d385f254139
model = Chain(
            Flux.normalise,
	Conv((3, 3), imgsize[end]=>32, relu),
            Conv((3, 3), 32=>64, relu),
            MaxPool((2, 2)),
            Dropout(0.25),
            flatten,
            Dense(9216, 120, relu),
            Dropout(0.5),
            Dense(120, nclasses)
          )

# ╔═╡ 25c1d1ae-3950-47d7-8060-38b3b13c68c0
Flux.outdims(Conv((3, 3), 1 => 10), (28, 28, 1, 10))

# ╔═╡ 62b4d2c2-a30f-49f6-ab0a-3e61e2a3ac45
prod(out_conv_size)

# ╔═╡ a74b1a5a-af15-4c99-ba47-e4f4dc9f8396
Flux.outdims(Chain(
		Conv((3, 3), imgsize[end]=>32, relu),
		Conv((3, 3), 32=>64, relu),
		MaxPool((2,2)),Dropout(0.25), flatten),

	(28, 28, 1, 10))

# ╔═╡ 46faca74-b9c6-470f-b916-329f566e17a9
imgsize[end]

# ╔═╡ b1b604ef-a722-441c-b58a-709602dd3e71
model[1:end](xs) |> size

# ╔═╡ 86ec66b0-bb94-4df7-ba27-91c6d89fac05
m[3](xs) |> size

# ╔═╡ Cell order:
# ╠═c5b16d82-6281-11ec-250b-8be67db3a2e7
# ╠═b75114fe-f046-4681-9348-139318cf8ba8
# ╠═f90bf01b-06be-4685-9ad9-6dca0ac87b88
# ╠═73ce15dc-d40d-4df8-abe7-4e7b3019f8bc
# ╠═596bc7ab-b3bc-4d02-b37c-04d34c3be108
# ╠═2d0d9a40-6b92-49b2-ba07-bee22f8d3d0d
# ╠═fc8ec7eb-63ea-4f5c-9600-eba205eb3f70
# ╠═0eb0d2b6-6c67-479b-a6f5-1a4fa7aa9e8e
# ╠═a79c9b40-56db-4056-bbf6-4d385f254139
# ╠═25c1d1ae-3950-47d7-8060-38b3b13c68c0
# ╠═62b4d2c2-a30f-49f6-ab0a-3e61e2a3ac45
# ╠═a74b1a5a-af15-4c99-ba47-e4f4dc9f8396
# ╠═46faca74-b9c6-470f-b916-329f566e17a9
# ╠═b1b604ef-a722-441c-b58a-709602dd3e71
# ╠═86ec66b0-bb94-4df7-ba27-91c6d89fac05

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil