text:
export ScikitRegressor

@with_kw mutable struct ScikitRegressor <: BaseEstimator
    pyo::PyObject = PyNULL()
    name::String = "lightgbm"
    warm_start::Bool = false
    alpha::Float32 = 1.0
    fit_intercept::Bool = true
    kernel::String = "rbf"
    gamma::Float32 = 0
    n_neighbors::Int = 10
    num_layers::Int = 1
    hidden_size::Int = 100
    lr::Float32 = 0.001
    n_estimators::Int = 20
    num_leaves::Int = 31
    colsample_bytree::Float32 = 1.0
    subsample::Float32 = 1.0
    reg_lambda::Float32 = 0.1
    min_child_samples::Int = 100
end

is_classifier(::ScikitRegressor) = false

function gridparams(m::ScikitRegressor)
    grid = [
        [
            "name" => ["lasso", "ridge", "linearsvr"],
            "alpha" => [0.001, 0.01, 0.1, 1]
        ],
        [
            "name" => ["kernelridge", "svr"],
            "alpha" => [0.001, 0.01, 0.1, 1],
            "kernel" => ["rbf", "poly", "sigmoid"],
            "gamma" => [0, 0.01, 0.1, 1, 10, 100]
        ],
        [
            "name" => ["knn"],
            "n_neighbors" => [10, 100, 1000],
        ],
        [
            "name" => ["mlp"],
            "num_layers" => [1, 2],
            "hidden_size" => [10, 100],
            "lr" => [0.001, 0.01]
        ],
        [
            "name" => ["lightgbm"],
            "num_leaves" => [15, 31],
            "min_child_samples" => [1000, 100, 5000],
            "n_estimators" => [20, 50],
            "reg_lambda" => [0.1, 1, 10],
            "subsample" => [1, 0.8],
            "colsample_bytree" => [1, 0.8],
        ]
    ]
    gridparams(grid)
end

function fit!(m::ScikitRegressor, x, y, w = nothing; columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    x, y, w = pymat(x), vec(y), vec(w)
    @unpack name, warm_start, alpha, fit_intercept, kernel, gamma, n_neighbors = m
    @unpack num_layers, hidden_size, lr, n_estimators, num_leaves = m
    @unpack colsample_bytree, subsample, reg_lambda, min_child_samples = m
    if name == "ridge"
        @from sklearn.linear_model imports Ridge
        pyo = Ridge(alpha = alpha)
    elseif name == "lasso"
        @from sklearn.linear_model imports Lasso
        pyo = Lasso(alpha = alpha)
    elseif name == "ridgecv"
        @from sklearn.linear_model imports RidgeCV
        pyo = RidgeCV(fit_intercept = fit_intercept)
    elseif name == "lassocv"
        @from sklearn.linear_model imports LassoCV
        pyo = LassoCV(fit_intercept = fit_intercept)
    elseif name == "kernelridge"
        @from sklearn.kernel_ridge imports KernelRidge
        pyo = KernelRidge(alpha = alpha, kernel = kernel, gamma = gamma == 0 ? nothing : gamma)
    elseif name == "knn"
        @from sklearn.neighbors imports KNeighborsRegressor
        pyo = KNeighborsRegressor(n_neighbors = n_neighbors)
    elseif name == "linearsvr"
        @from sklearn.svm imports LinearSVR
        pyo = LinearSVR(dual = false, C = 0.5 / alpha, loss = "squared_epsilon_insensitive")
    elseif name == "svr"
        use_thunder = true
        try
            @from thundersvm imports SVR
            pyo = SVR(kernel = kernel, C = 0.5 / alpha, gamma = gamma == 0 ? "auto" : gamma, n_jobs = 20)
        catch e
            use_thunder = false
            @from sklearn.svm imports SVR
            pyo = SVR(kernel = kernel, C = 0.5 / alpha, gamma = gamma == 0 ? "auto" : gamma)
        end
    elseif name == "mlp"
        @from sklearn.neural_network imports MLPRegressor
        pyo = MLPRegressor(
            hidden_layer_sizes = ntuple(i -> hidden_size, num_layers),
            alpha = alpha,
            learning_rate_init = lr,
            verbose = true,
            warm_start = true,
        )
        if warm_start && occursin("MLP", string(m.pyo))
            @unpack pyo = m
        end
    elseif name == "lightgbm"
        @from lightgbm imports LGBMRegressor
        if warm_start
            @imports lightgbm
            @from unidecode imports unidecode
            dset = lightgbm.Dataset(data = x, label = y, weight = w)
            columns′ = [filter(!isspace, unidecode(c)) for c in columns]
            dset.set_feature_name(columns′)
            params = Dict(pairs((n_jobs = 20, max_bin = 15, num_leaves = num_leaves,
                        subsample = subsample, colsample_bytree = colsample_bytree,
                        reg_lambda = reg_lambda, min_child_samples = min_child_samples)))
            init_model = isfile("lightgbm") ? "lightgbm" : nothing
            pyo = lightgbm.train(params, dset, num_boost_round = n_estimators, init_model = init_model)
            pyo.save_model("lightgbm")
        else
            pyo = LGBMRegressor(n_jobs = 20, max_bin = 15, num_leaves = num_leaves,
                        subsample = subsample, colsample_bytree = colsample_bytree,
                        reg_lambda = reg_lambda, min_child_samples = min_child_samples)
            pyo.fit(x, y, sample_weight = w)
        end
    elseif name == "autosk"
        @from autosklearn.regression imports AutoSklearnRegressor
        AutoSklearnRegressor(resampling_strategy = "cv", resampling_strategy_arguments = Dict("folds" => 5))
    end
    if name == "svr" && use_thunder || occursin(r"mlp|lasso|knn", name)
        pyo.fit(x, y)
    elseif name != "lightgbm"
        pyo.fit(x, y, sample_weight = w)
    end
    @pack! m = pyo
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

predict(m::ScikitRegressor, x) = m.pyo.predict(pymat(x))

istree(m::ScikitRegressor) = occursin(r"tree|gbm"i, m.name)

reset!(m::ScikitRegressor) = rm("lightgbm", force = true)

modelhash(m::ScikitRegressor) = hash(m.pyo)

function visualize(m::ScikitRegressor, columns)
    @unpack name, pyo = m
    if name == "lightgbm"
        gbm = GbmModel(pyo = pyo)
        visualize(gbm, columns)
    elseif occursin(r"^ridge|lasso", name)
        write_feaimpt(pyo.coef_, columns)
    end
end


text:
include("utilities/misc.jl")
include("utilities/group_actions.jl")
include("utilities/sorting_points.jl")
include("utilities/linear_algebra.jl")
include("utilities/inner.jl")
include("utilities/polynomials.jl")

ext: jl


text:
# ---
# title: 852. Peak Index in a Mountain Array
# id: problem852
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Binary Search
# link: <https://leetcode.com/problems/peak-index-in-a-mountain-array/description/>
# hidden: true
# ---
#
# Let's call an array `arr` a **mountain**  if the following properties hold:
#
#   * `arr.length >= 3`
#   * There exists some `i` with `0 < i < arr.length - 1` such that:
#     * `arr[0] < arr[1] < ... arr[i-1] < arr[i] `
#     * `arr[i] > arr[i+1] > ... > arr[arr.length - 1]`
#
# Given an integer array arr that is **guaranteed** to be a mountain, return any
# `i` such that `arr[0] < arr[1] < ... arr[i - 1] < arr[i] > arr[i + 1] > ... >
# arr[arr.length - 1]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [0,1,0]
#     Output: 1
#
#
# **Example 2:**
#
#
#
#     Input: arr = [0,2,1,0]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: arr = [0,10,5,2]
#     Output: 1
#
#
# **Example 4:**
#
#
#
#     Input: arr = [3,4,5,1]
#     Output: 2
#
#
# **Example 5:**
#
#
#
#     Input: arr = [24,69,100,99,79,78,67,36,26,19]
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `3 <= arr.length <= 104`
#   * `0 <= arr[i] <= 106`
#   * `arr` is **guaranteed** to be a mountain array.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# Remember to install HDF5 and JLD packages
using HDF5, JLD, Nervoso

# Initializing our examples and network
@load "data.jld"
net = FFNNet(400, 25, 10)

# Just checking our in-sample error now to compare it later
println("In-Sample Mean Error before training: ", meanerror(net, ins, outs))
println("In-Sample Classification Error before training: ", classerror(net, ins, outs))

# Let's train out network 100 times using our examples
train!(net, ins, outs, α=0.1, η=0., epochs=100, batchsize=1)

# Now we compare our new in-sample error
println("In-Sample Mean Error after training: ", meanerror(net, ins, outs))
println("In-Sample Classification Error after training: ", classerror(net, ins, outs))


text:
module RecSysChallenge

export FileSpec, read_input, Inputs, prepdata, ratings, train

end # module


text:
using FastaIO
using DataFrames
using Lumberjack
using GZip

include(Pkg.dir("PseudoGenomes", "src","pseudogenome-misc.jl"))
#
# This is a file formatted in the way you get it
# from WTSI see here: ftp://ftp-mouse.sanger.ac.uk/current_svs/
#
function read_SV_file(filename::ASCIIString;gzip=true)
    Lumberjack.info("Start Reading $filename")
    num_lines=0
    heading = Regex("^#CHROM")
    line_itr = nothing
    if gzip
       #line_itr = eachline(GZip.open(filename) ) - GZip is broken for some versions of gzip - gives an extraline - see bug report
       line_itr  = eachline(`zcat $filename`)
    else
       line_itr = eachline(open(filename) )
    end

    # check for comment lines and header
    df = DataFrame()
    header=[]
    for line in line_itr
       line=chomp(line)
       if ismatch(heading, line)
           Lumberjack.info("DETECT HEADING: $line")
           header = split(line,'\t')
           species_names=header[5:end]
           Lumberjack.info("The species names are $species_names")
           break
       end
       # anything above the line heading is a comment
       Lumberjack.info("DETECT COMMENT: $line")
    end

    Lumberjack.info("Load first line into dataframe")
    first_data_line = first(line_itr)
    first_data_line=chomp(first_data_line)
    first_data_fields = split(first_data_line,'\t')
    df[:seq_id]     = first_data_fields[1]
    df[:start]      = first_data_fields[2]
    df[:stop]       = first_data_fields[3]
    df[:formats]    = first_data_fields[4]
    for field_idx   = 5:length(header)
        df[symbol(header[field_idx])] = first_data_fields[field_idx]
    end
    Lumberjack.info("First dataline: $df")

    df_num_cols = ncol(df)
    for line in line_itr
         line=chomp(line)
         fields = split(line,'\t')
         num_fields = length(fields)
         if num_fields == df_num_cols
             push!(df,fields)
         else
             error("Number of fields parsed inconsistent $fields from $line.\nNumber fields parsed: $num_fields Num Cols Expected $df_num_cols\nLine number: $(num_lines+1)")
         end
         num_lines += 1
    end
    #DELS = SVs[ SVs[:SV_type] .== "DEL", :]
    Lumberjack.info( "Total number of data lines: $(num_lines+1)")
    return df
end

function get_strain_SVs(SVs::DataFrame, strain_id::ASCIIString,sv_type::Regex; seq_id_format="ucsc")

    # get the strain
    strain_column = SVs[symbol(strain_id)]

    chrs=ASCIIString[]
    starts=Int64[]
    stops=Int64[]

    # and parse the first field to get the actual coordinates for the structural vairant
    for row in strain_column
        # The fields are POS::CLASSIFICATION::BREAKPPOINT::REPEAT_TYPE
        # it seems there used to be a field called PEM_pattern that crops up.
        # The fields are specified  in the format field of the tab file for each row (POS:CL:BP:TY)
        # so this prob should be take this into account in case they change the fields in the future.

        fields = split(row, ';')

        if length(fields) == 1
            continue
        elseif length(fields) < 4
            Lumberjack.warn("Less than expected number of fields")
        end

        if !ismatch(sv_type,fields[2])
            continue
        end
        coordinates=fields[1]
        (chr, position) = split(coordinates,':')
        (start,stop)=split(position,'-')
        start_int = parseint(start)
        stop_int  = parseint(stop)
        if seq_id_format == "ucsc"
            push!(chrs, ens2ucsc(chr) )
        else
            push!(chrs, chr)
        end
        push!(starts,start_int)
        push!(stops,stop_int)
    end
    return DataFrame(seq_id=chrs,start=starts,stop=stops)
end

function mask_region!(sequence,SV_deletions_df)
   num_masked=0
   for row=1:nrow(SV_deletions_df)
        startpos = SV_deletions_df[row,:start]
        stoppos  = SV_deletions_df[row,:stop]
        for i=startpos:stoppos
            sequence[i]='N'
            num_masked +=1
        end
    end
    Lumberjack.info("Done masking")
    return num_masked
end

function mask_fasta_file_with_strain_SV_deletions(filename,outfilename,SV_deletions::DataFrame)
    # roll through the fastafile and then roll through the fastafile switch to N.

    fr = FastaReader{Vector{Char}}(filename)
    fw = FastaWriter(outfilename)
    Lumberjack.info("Start masking fasta file")
    for (desc, seq) in fr
         seqlength = length(seq)
         Lumberjack.info("processing $desc - length: $seqlength")
         seqid_df = SV_deletions[ SV_deletions[:seq_id] .== desc,:]
         if nrow(seqid_df) > 0
             num_masked = mask_region!(seq, seqid_df )
             Lumberjack.info("Masked $num_masked nucleotides")
         else
             Lumberjack.info("No deletion coordinates for $desc")
         end

         writeentry(fw, desc, seq)
    end
   Lumberjack.info("read $(fr.num_parsed) entries")
   close(fw)
end

function mask_fasta_file_with_strain_SV_deletions(fastafile::ASCIIString,outfilename::ASCIIString, SVfile::ASCIIString,strain::ASCIIString)
    SV_df = read_SV_file(SVfile)
    strain_SVs_df = get_strain_SVs(SV_df, strain,Regex("DEL"))
    mask_fasta_file_with_strain_SV_deletions(fastafile,outfilename,strain_SVs_df)
end


text:
# NOTE: Currently, not used.
module TestGibbsPT

import Pkg; Pkg.activate("../..")  # CytofRepFAM
import CytofRepFAM.MCMC, Random
using Distributions, Distributed

mutable struct Param
  mu
  sig2
  w
  lam
end

struct Priors
  mu
  sig2
  w
end

Priors(K) = Priors(Normal(0, 1), InverseGamma(3, 2), Dirichlet(K, 1/K))

Param(N, priors) = let
  K = length(priors.w.alpha)
  Param(rand(priors.mu, K), rand(priors.sig2),
        rand(priors.w), rand(1:K, N))
end

# Set up number of cores to use
maxcores = 20

# if length(workers()) != maxcores
#   rmprocs(filter(w -> w > 1, workers()))
#   addprocs(maxcores)
#   @everywhere begin
#     import Pkg; Pkg.activate("../")
#     using CytofRepFAM
#   end
# end

function update_mu(k, theta, priors, y)
  idx_k = findall(theta.lam .== k)
  nk = length(idx_k)
  yk = y[idx_k]
  newPrec = 1/var(priors.mu) + nk/theta.sig2
  newVar = 1/newPrec
  newMean = newVar * (mean(priors.mu)/var(priors.mu) + sum(yk)/theta.sig2)
  theta.mu[k] = rand(Normal(newMean, sqrt(newVar)))
end

function update_mu(theta, priors, y)
  for k in 1:length(theta.mu)
    update_mu(k, theta, priors, y)
  end
end

function update_sig2(theta, priors, y)
  N = length(y)
  ss = sum((y - theta.mu[theta.lam]) .^ 2)
  newShape = N / 2 + shape(priors.sig2)
  newScale = ss / 2 + scale(priors.sig2)
  theta.sig2 = rand(InverseGamma(newShape, newScale))
end

function update_w(theta, priors, y)
  N = length(y)
  K = length(theta.w)
  counts = zeros(Int, K)
  for i in 1:N
    counts[theta.lam[i]] += 1
  end
  newAlpha = priors.w.alpha + counts
  theta.w = rand(Dirichlet(newAlpha))
end

function update_lam(theta, priors, y)
  N = length(y)
  for i in 1:N
    update_lam(i, theta, priors, y)
  end
end

function update_lam(i, theta, priors, y)
  lls = MCMC.lpdf_normal.(y[i], theta.mu, sqrt(theta.sig2))
  lp = log.(theta.w)
  logprob = lls + lp
  theta.lam[i] = MCMC.wsample_logprob(logprob)
end

function update(theta, priors, y)
  update_mu(theta, priors, y)
  update_sig2(theta, priors, y)
  update_lam(theta, priors, y)
  update_w(theta, priors, y)
end

begin
  mu_true = [-2, 3, 1]
  sig2_true = 0.01
  w_true = [.3, .5, .2]
  K_true = length(mu_true)
  nobs = 100
  mus = wsample(mu_true, w_true, nobs)
  y = rand.(Normal.(mus, sqrt(sig2_true)))

  K = 5
  priors = Priors(5)
  theta = Param(nobs, priors)

  nburn = 100
  nmcmc = 1000
  states = Param[]

  for i in 1:nburn
    update(theta, priors, y)
  end

  for i in 1:nmcmc
    update(theta, priors, y)
    append!(states, [deepcopy(theta)])
  end

  extract(sym) = [getfield(s, sym) for s in states]

  mus = extract(:mu)
  sig2s = extract(:sig2)
  ws = extract(:w)

  println("mu: ", mean(mus))
  println("sig2: ", mean(sig2s))
  println("w: ", mean(ws))
end

end


text:
# NOTE(tim):
# value V will be placed in bin B if V ∈ [Bₗ Bᵣ)
# (or V ∈ [Bₗ Bᵣ] if B is the rightmost bin)

# TODO(tim):
# handle NA
# handle Nullable
# handle NaN
# ability to specify decoding behavior

const DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST = true

immutable LinearDiscretizer{N<:Real, D<:Integer} <: AbstractDiscretizer{N,D}
    binedges :: Vector{N}   # list of bin edges, sorted smallest to largest
    nbins    :: Int
    i2d      :: Dict{Int,D} # maps bin index to discrete label
    d2i      :: Dict{D,Int} # maps discrete label to bin index

    force_outliers_to_closest :: Bool # if true, real values outside of the bin ranges will be forced to the nearest bin, otherwise will throw an error
end

function LinearDiscretizer{N<:Real, D<:Integer}( binedges::Vector{N}, i2d::Dict{Int,D};
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST )

    length(binedges) > 1 || error("bin edges must contain at least 2 values")

    if N <: AbstractFloat
        findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)) == 0 ||
            error("Bin edges must be sorted in increasing order")
    else # for integers, bins of unit width require repeated values
        (findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)-1) == 0 &&
         findfirst(i->binedges[i-1] > binedges[i], 2:length(binedges)) == 0 ) ||
            error("Bin edges must be sorted in increasing order")
    end

    for i = 1 : length(binedges)-1
        haskey(i2d, i) || error("i2d must contain all necessary keys")
    end

    d2i = Dict{D,Int}()
    for (k,v) in i2d
        d2i[v] = k
    end

    LinearDiscretizer{N,D}(binedges, length(binedges)-1, i2d, d2i, force_outliers_to_closest)
end
function LinearDiscretizer{N<:Real, D<:Integer}( binedges::Vector{N}, ::Type{D} = Int;
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST )

    length(binedges) > 1 || error("bin edges must contain at least 2 values")

    if N <: AbstractFloat
        findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)) == 0 ||
            error("Bin edges must be sorted in increasing order")
    else # for integers, bins of unit width require repeated values
        (findfirst(i->binedges[i-1] ≥ binedges[i], 2:length(binedges)-1) == 0 &&
         findfirst(i->binedges[i-1] > binedges[i], 2:length(binedges)) == 0 ) ||
            error("Bin edges must be sorted in increasing order")
    end

    i2d = Dict{Int,D}()
    for i in 1 : length(binedges)-1
        i2d[i] = convert(D, i)
    end

    LinearDiscretizer(binedges,i2d,force_outliers_to_closest=force_outliers_to_closest)
end
function LinearDiscretizer{N<:Real, D<:Integer}(arr::AbstractVector{N}, ::Type{D}=Int;
    force_outliers_to_closest::Bool = DEFAULT_LIN_DISC_FORCE_OUTLIERS_TO_CLOSEST)
    LinearDiscretizer(convert(Vector{N}, arr), D, force_outliers_to_closest=force_outliers_to_closest)
end

function supports_encoding{N<:Real,D<:Integer}(ld::LinearDiscretizer{N,D}, x::N)
    if ld.force_outliers_to_closest
        return true
    else
        return ld.binedges[1] ≤ x ≤ ld.binedges[end]
    end
end
supports_decoding{N<:Real,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D) = 1 ≤ d ≤ ld.nbins

function encode{N,D<:Integer}(ld::LinearDiscretizer{N,D}, x::N)
    if isnan(x)
        error("cannot encode NaN values")
    end
    if x < ld.binedges[1]
        return ld.force_outliers_to_closest ? ld.i2d[1] : throw(BoundsError())
    elseif x > ld.binedges[end]
        return ld.force_outliers_to_closest ? ld.i2d[ld.nbins] : throw(BoundsError())
    end

    # run bisection search
    binedges = ld.binedges
    a, b = 1, length(binedges)
    va, vb = binedges[a], binedges[b]
    if isapprox(vb, x)
        return ld.i2d[b-1]
    end
    while b-a > 1
        c = div(a+b, 2)
        vc = binedges[c]
        if isapprox(x, vc)
            return ld.i2d[c]
        end
        if x < vc
            b, vb = c, vc
        else
            a, va = c, vc
        end
    end
    return ld.i2d[a]
end
encode{N,D}(ld::LinearDiscretizer{N,D}, x) = encode(ld, convert(N, x))::D
function encode{N,D<:Integer}(ld::LinearDiscretizer{N,D}, data::AbstractArray)
    arr = Array(D, length(data))
    for (i,x) in enumerate(data)
        arr[i] = encode(ld, x)
    end
    reshape(arr, size(data))
end

function decode{N<:AbstractFloat,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniform)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    convert(N, lo + rand()*(hi-lo))
end
function decode{N<:AbstractFloat,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleBinCenter)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    convert(N, (hi + lo)/2)
end
function decode{N<:AbstractFloat,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniformZeroBin)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]

    if lo ≤ 0.0 ≤ hi
        return zero(N)
    end
    return convert(N, lo + rand()*(hi-lo))
end
decode{N<:AbstractFloat,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D) = decode(ld, d, SAMPLE_UNIFORM)

function decode{N<:Integer,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniform)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    if hi != ld.binedges[end]
        retval = rand(lo:hi-1)::Int
    else
        retval = rand(lo:hi)::Int
    end
    convert(N, retval)
end
function decode{N<:Integer,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleBinCenter)
    ind = ld.d2i[d]
    lo = ld.binedges[ind]
    hi = ld.binedges[ind+1]
    convert(N, div(lo+hi,2))
end
function decode{N<:Integer,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D, ::SampleUniformZeroBin)
    ind = ld.d2i[d]
    lo = round(Int, ld.binedges[ind])
    hi = round(Int, ld.binedges[ind+1])

    if lo ≤ 0 ≤ hi
        retval = 0
    elseif hi != ld.binedges[end]
        retval = rand(lo:hi-1)::Int
    else
        retval = rand(lo:hi)::Int
    end

    convert(N, retval)
end
decode{N<:Integer,D<:Integer}(ld::LinearDiscretizer{N,D}, d::D) = decode(ld, d, SAMPLE_UNIFORM)

decode{N<:Real,D<:Integer,I<:Integer}(ld::LinearDiscretizer{N,D}, d::I, method::AbstractSampleMethod=SAMPLE_UNIFORM) =
    decode(ld, convert(D,d), method)

function decode{N,D<:Integer}(ld::LinearDiscretizer{N,D}, data::AbstractArray{D}, ::AbstractSampleMethod=SAMPLE_UNIFORM)
    arr = Array(N, length(data))
    for (i,d) in enumerate(data)
        arr[i] = decode(ld, d)
    end
    reshape(arr, size(data))
end

Base.max(ld::LinearDiscretizer) = ld.binedges[end]
Base.min(ld::LinearDiscretizer) = ld.binedges[1]
function extrema{N,D}(ld::LinearDiscretizer{N,D})
    lo  = ld.binedges[1]
    hi  = ld.binedges[end]
    (lo, hi)
end
function extrema{N<:AbstractFloat,D}(ld::LinearDiscretizer{N,D}, d::D)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    (lo, hi)
end
function extrema{N<:Integer,D}(ld::LinearDiscretizer{N,D}, d::D)
    ind = ld.d2i[d]
    lo  = ld.binedges[ind]
    hi  = ld.binedges[ind+1]
    if hi  == ld.binedges[end]
        (lo, hi)
    else
        (lo, hi-1)
    end
end
function totalwidth(ld::LinearDiscretizer)
    ex = extrema(ld)
    ex[2] - ex[1]
end

nlabels(ld::LinearDiscretizer) = ld.nbins
binedges(ld::LinearDiscretizer) = ld.binedges
bincenters{N<:AbstractFloat,D}(ld::LinearDiscretizer{N,D}) = (0.5*(ld.binedges[1:ld.nbins] + ld.binedges[2:end]))::Vector{Float64}
function bincenters{N<:Integer,D}(ld::LinearDiscretizer{N,D})
    retval = Array(Float64, ld.nbins)
    for i = 1 : length(retval)-1
        retval[i] = 0.5(ld.binedges[i+1]-1 + ld.binedges[i])
    end
    retval[end] = 0.5(ld.binedges[end] + ld.binedges[end-1])
    retval
end
binwidth{N<:AbstractFloat,D}(ld::LinearDiscretizer{N,D}, d::D) = ld.binedges[d+1] - ld.binedges[d]
binwidths{N<:AbstractFloat,D}(ld::LinearDiscretizer{N,D}) = ld.binedges[2:end] - ld.binedges[1:end-1]


text:
using Mimi

include("analysis/allscc/main_model_annual.jl")

m = getpage()
run(m)
explore(m)

## to run PAGE2020 on an annual timestep
# using Mimi
include("../src/main_model_annual.jl")
# include("../src/main_model_annual.jl")

# m = getpage()
# run(m)
# explore(m)


text:
abstract type FragMatrix end

struct FragMatrixPair{T} <: FragMatrix

    name::String

    numregions::Int
    totalfrags::Int
    numfields::Int
    fields::Vector{String}
    totalecfrags::Vector{Int}

    fragstats::Dict{Symbol, Int}

    index::Dict{String, UnitRange{Int}}
    chromindex::Dict{Int, String}

    FM::T
end

struct FragMatrixSingle{T} <: FragMatrix

    name::String
    readlength::Int
    numregions::Int
    totalfrags::Int
    numfields::Int
    fields::Vector{String}
    totalecfrags::Vector{Int}

    fragstats::Dict{Symbol, Int}

    index::Dict{String, UnitRange{Int}}
    chromindex::Dict{Int, String}

    FM::T
end

function Base.show(io::IO, fm::FragMatrixPair)
    if get(io, :compact, false)
        ### compact
        print(io, "FM Paired Indexed: ", fm.name)
        print(io, "    ", fm.numregions)
        print(io, "    ", fm.totalfrags)
    else
        println(io, "Index FM      :   ", fm.name)
        #println(io, "paired        :   ", fm.paired)
        println(io, "numregions    :   ", fm.numregions)
        println(io, "totalfrags    :   ", fm.totalfrags)
        println(io, "numfields     :   ", fm.numfields)
        println(io, "fields        :   ", fm.fields)
        println(io, "totalecfrags  :   ", fm.totalecfrags)
        println(io, "index         :   ", fm.index)
        println(io, "chromindex    :   ", fm.chromindex)
        println(io, "FM            :   ", summary(fm.FM))
    end
end

function Base.show(io::IO, fm::FragMatrixSingle)
    if get(io, :compact, false)
        ### compact
        print(io, "FM Paired Indexed: ", fm.name)
        print(io, "    ", fm.numregions)
        print(io, "    ", fm.totalfrags)
    else
        println(io, "Index FM      :   ", fm.name)
        println(io, "Readlength    :   ", fm.readlength)
        #println(io, "paired        :   ", fm.paired)
        println(io, "numregions    :   ", fm.numregions)
        println(io, "totalfrags    :   ", fm.totalfrags)
        println(io, "numfields     :   ", fm.numfields)
        println(io, "fields        :   ", fm.fields)
        println(io, "totalecfrags  :   ", fm.totalecfrags)
        println(io, "index         :   ", fm.index)
        println(io, "chromindex    :   ", fm.chromindex)
        println(io, "FM            :   ", summary(fm.FM))
    end
end

##### Find indexed frags
##### Return empty view if chrom is missing
function get_frags(chrom, location::UnitRange{Int}, fm::FragMatrix)

    !haskey(fm.index, chrom) && return view(fm.FM, :, 1:0)

    ind = fm.index[chrom]
    starts = view(fm.FM, 1, ind)
    stops  = view(fm.FM, 2, ind)

    st = min(searchsortedfirst(starts, location.start), searchsortedfirst(stops, location.start))
    (st > size(starts, 1)) && return view(fm.FM, :, 1:0)

    en = min(max(searchsortedlast(starts, location.stop), searchsortedlast(stops, location.stop)), length(stops))

    view(fm.FM, :, ind.start .+ (st:en) .- 1)
end

function get_frags(chrom, fm::FragMatrix)
    !haskey(fm.index, chrom) && return view(fm.FM, :, 1:0)
    ind = fm.index[chrom]
    view(fm.FM, :, ind)
end

#### sorts chroms in numerical order then X, Y, M
function chrom_lt(chrA, chrB)
    chrA[4] == 'M' && return false
    chrB[4] == 'M' && return true

    if isnumber(chrA[4])
        !isnumber(chrB[4]) && return true
        return isless(parse(Int, chrA[4:end]), parse(Int, chrB[4:end]))
    else
        isnumber(chrB[4]) && return false
        return isless(chrA[4:end], chrB[4:end])
    end
end

function load_frag_matrix(file, paired=true)

     io = open(file)

    paired, numregions, totalfrags, numfields, fields, totalecfrags, datatype, index, chromindex, fragstats = read_header_footer(io, paired)

    #  @show paired, numregions, totalfrags, numfields, fields, totalecfrags, datatype, index, chromindex, fragstats

    a = Mmap.mmap(io, Vector{UInt8}, sizeof(datatype)*numfields*totalfrags, position(io))
    # Create an array of the desired eltype and size:
    ua = UnalignedVector{datatype}(a)

    F = reshape(ua, (numfields, totalfrags))

    ft = typeof(F)
    if paired
        FragMatrixPair{ft}(file, numregions, totalfrags, numfields, fields, totalecfrags, fragstats, index, chromindex, F)
    else
        ### estimate readlength
        readlength = estreadlength(F)
        FragMatrixSingle{ft}(file, readlength, numregions, totalfrags, numfields, fields, totalecfrags, fragstats, index, chromindex, F)
    end
end

### for single reads
function estreadlength(F, nr=min(500, size(F, 2)))
    readlengths = counter(Int)
    for i = 1:nr
        readlength = F[2, i] - F[1, i] + 1
        push!(readlengths, readlength)
    end
    rls = collect(keys(readlengths))
    c = [readlengths[k] for k in rls]
    si = sortperm(c, rev=true)

    if rls[si[1]] != maximum(rls)
        println("Warning: maxfraglength ≠ most common frag length: ", rls[si[1]], " ≠ ", maximum(rls))
    end
    return rls[si[1]]
end

function read_header_footer(io, paired=true)

    numregions = -1
    numfields = -1
    totalfrags = -1
    fmfields = String[]
    totalecfrags = Int[]
    dataencoding = Int32
    index = Dict{String, UnitRange{Int}}()
    chromindex = Dict{Int,  String}()

    fragstats = Dict{Symbol, Int}()

    bottom_header = false
    first_line = readline(io)

    if chomp(first_line) == "#bottom_header"
        headerpos = read(io, Int64)
        bottom_header = true
        seek(io, headerpos)
    else
        seekstart(io)
    end

    ### assumes index is final field
    while true
         fields = split(readline(io), '\t')
        (fields[1] == "#paired")     && (paired     = parse(Bool, chomp(fields[2]))::Bool)
        (fields[1] == "#numregions")   && (numregions   = parse(Int,  fields[2])::Int)
        (fields[1] == "#totalfrags") && (totalfrags = parse(Int,  fields[2])::Int)
        (fields[1] == "#numfields")  && (numfields  = parse(Int,  fields[2])::Int)
        (fields[1] == "#fields")     && (fmfields     = split(chomp(fields[2]), [',', ' ', '[', ']'], keepempty=false))
        (fields[1] == "#totalecfrags") && (totalecfrags = parse.([Int], split(chomp(fields[2]), [',', ' ', '[', ']'], keepempty=false))::Vector{Int})
        (fields[1] == "#dataencoding")   && (dataencoding = get_datatype(chomp(fields[2])))

        if (fields[1] == "#numregions") || (fields[1] == "#totalfrags") || (fields[1] == "#totaldupfrags")
            fragstats[Symbol(fields[1][2:end])] = parse(Int,  fields[2])::Int
        end

        if fields[1] == "#index"
            index_str = split(replace(chomp(fields[2]), r"\(|\"" => ""), r"\)[,]*", keepempty=false)
            index_vec = parse_chrom_range.(index_str)

            index = Dict(index_vec)
            chromindex = Dict(i => c[1] for (i, c) in enumerate(index_vec))
            break
        end
    end

    if isempty(totalecfrags)
        totalecfrags = [totalfrags]
    end

    if bottom_header
        seekstart(io)
        readline(io)
        read(io, Int64)
    end

    paired, numregions, totalfrags, numfields, fmfields, totalecfrags, dataencoding, index, chromindex, fragstats

end

function parse_chrom_range(s)
    chrom, rind = split(s, ",")
    sstart, sstop = split(rind, ":")
    chrom, parse(Int, sstart):parse(Int, sstop)
end

function get_datatype(s)
    (s == "Int32")  && return Int32
    (s == "UInt32") && return UInt32
    (s == "Int64")  && return Int64
    (s == "UInt64") && return UInt64

    error("Unsupported data type: $s")
end

####### Encode Decode of strand_chrom
set_strand_chrom_enc(strand, chromindex, ::Type{UInt32}, topbit=0x80000000) = (strand == STRAND_NEG) ?  UInt32(chromindex) | topbit : UInt32(chromindex)
get_strand_chrom_enc(enc::UInt32, topbit=0x80000000) = (enc & topbit) > 0 ? (STRAND_NEG, enc $ topbit) : (STRAND_POS, enc)
set_strand_chrom_enc(strand, chromindex, ::Type{UInt64}, topbit=0x8000000000000000) = (strand == STRAND_NEG) ?  UInt64(chromindex) | topbit : UInt64(chromindex)
get_strand_chrom_enc(enc::UInt64, topbit=0x8000000000000000) = (enc & topbit) > 0 ? (STRAND_NEG, enc $ topbit) : (STRAND_POS, enc)
set_strand_chrom_enc(strand, chromindex, ::Type{Int32}) = (strand == STRAND_NEG) ? -Int32(chromindex) : Int32(chromindex)
get_strand_chrom_enc(enc::Int32) = (enc < 0) ? (STRAND_NEG, -enc) : (STRAND_POS, enc)
set_strand_chrom_enc(strand, chromindex, ::Type{Int64}) = (strand == STRAND_NEG) ? -Int64(chromindex) : Int64(chromindex)
get_strand_chrom_enc(enc::Int64) = (enc < 0) ? (STRAND_NEG, -enc) : (STRAND_POS, enc)

isneg(enc::UInt32, topbit=0x80000000)         = (enc & topbit) > 0 ? true : false
isneg(enc::UInt64, topbit=0x8000000000000000) = (enc & topbit) > 0 ? true : false
isneg(enc::Int32)                             =  enc           < 0 ? true : false
isneg(enc::Int64)                             =  enc           < 0 ? true : false

# chromsome, strand = GenomeFrags2.get_strand_chrom_enc(v[3, i])
# STRAND_POS
# STRAND_NEG

#### Accessor
totalfrags(FM) = FM.totalfrags
totalfrags(FM, s::Symbol) = totalfrags(FM, string(s))
function totalfrags(FM, s::String)
    ind = findfirst(f -> f == s, FM.fields)
    (ind == 0) && error("Field: $s not found")
    FM.totalecfrags[s - 3]
end


text:
using RxNav
using Test

@test rcui("ibuprofen") == "5640"

@test contains(first(drugs("naproxen")), "Oral")

@test contains(last(interaction("61148"; ONCHigh = false)).description, "creased")

@test first(interaction("sumatriptan")).severity == "high"

@test contains(last(interaction_within_list(["207106", "656659"])).description, "metabolism")

@test last(interaction_within_list(["divalproex", "lamotrigine"])).severity == "N/A"

prescribable(true)
RxNav.getSpellingSuggestions("nortriptelene") == ["nortriptyline"]
prescribable(false)
RxNav.getSpellingSuggestions("nortriptelene") == ["nortriptyline", "Nortriptylina"]


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBase.jl/blob/master/LICENSE

using FEMBase
using FEMBase: get_local_coordinates, inside
using Base.Test

@testset "inverse isoparametric mapping" begin
    el = Element(Quad4, [1, 2, 3, 4])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    update!(el, "geometry", X)
    time = 0.0
    X1 = el("geometry", [0.1, 0.2], time)
    xi = get_local_coordinates(el, X1, time)
    X2 = el("geometry", xi, time)
    info("X1 = $X1, X2 = $X2")
    @test isapprox(X1, X2)
end

@testset "inside of linear element" begin
    el = Element(Quad4, [1, 2, 3, 4])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [1.0, 1.0],
        4 => [0.0, 1.0])
    update!(el, "geometry", X)
    time = 0.0
    @test inside(el, [0.5, 0.5], time) == true
    @test inside(el, [1.0, 0.5], time) == true
    @test inside(el, [1.0, 1.0], time) == true
    @test inside(el, [1.01, 1.0], time) == false
    @test inside(el, [1.0, 1.01], time) == false
end

@testset "inside of quadratic element" begin
    el = Element(Tri6, [1, 2, 3, 4, 5, 6])
    X = Dict{Int64, Vector{Float64}}(
        1 => [0.0, 0.0],
        2 => [1.0, 0.0],
        3 => [0.0, 1.0],
        4 => [0.5, 0.2],
        5 => [0.8, 0.6],
        6 => [-0.2, 0.5])
    update!(el, "geometry", X)
    p = [0.94, 0.3] # visually checked to be inside
    @test inside(el, p, 0.0) == true
    p = [-0.2, 0.8] # visually checked to be outside
    @test inside(el, p, 0.0) == false
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Algorithm test for SGP4 orbit propagator. All tests are based on [1].
#
#   Notice that only the tests related to the SGP4 in [1] are considered,
#   because the SDP4 (orbit propagator for deep space) are not available yet in
#   SatelliteToolbox.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A., Crawford, P., Hujsak, R., Kelso, T. S (2006). Revisiting
#       Spacetrack Report #3: Rev1. AIAA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                   Test 01
################################################################################

@testset "Tests from the paper AIAA 2006-6753" begin
    # Read all TLEs that will be used to test.
    tles = read_tle("./sgp4_tests/sgp4_tests.tle")

    # Array with the satellite numbers that must be tested.
    # sat_num_array = [5; 6251; 28057; 28350; 28872; 29141; 29238; 88888]

    # Find the TLEs with the respective satellite numbers.
    # tles_test_01 = tles[find( (tle->tle.sat_num in sat_num_array), tles)]

    # Test all the TLEs.
    tles_test_01 = tles

    st_sgp4_result = []
    SGP4_results = []

    for tle in tles_test_01
        filename = @sprintf("./sgp4_tests/aiaa-2006-6753/sgp4_tle_%d_result.txt",
                            tle.sat_num)
        SGP4_results = readdlm(filename; comments=true)

        # Initialize the orbit propagator.
        orbp = init_orbit_propagator(Val(:sgp4), tle; sgp4_gc = sgp4_gc_wgs72)

        # Propagate the orbit.
        t = SGP4_results[:,1]*60
        r_TEME, v_TEME = propagate!(orbp, t)

        # Compare the results.
        @inbounds for k = 1:length(t)
            # Assemble the result vector.
            st_sgp4_result = vcat(t[k]/60, r_TEME[k]/1000, v_TEME[k]/1000)

            # Compare the values.
            @test st_sgp4_result[1] == SGP4_results[k,1]
            @test st_sgp4_result[2] ≈  SGP4_results[k,2] atol=1e-8
            @test st_sgp4_result[3] ≈  SGP4_results[k,3] atol=1e-8
            @test st_sgp4_result[4] ≈  SGP4_results[k,4] atol=1e-8
            @test st_sgp4_result[5] ≈  SGP4_results[k,5] atol=1e-9
            @test st_sgp4_result[6] ≈  SGP4_results[k,6] atol=1e-9
            @test st_sgp4_result[7] ≈  SGP4_results[k,7] atol=1e-9
        end
    end

    # Test using the function `propagate_to_epoch!`
    # ==========================================================================
    #
    # In this case, testing the function using only one case is enough.

    tle = read_tle_from_string(
        "1 23599U 95029B   06171.76535463  .00085586  12891-6  12956-2 0  2905",
        "2 23599   6.9327   0.2849 5782022 274.4436  25.2425  4.47796565123555")

    filename     = "./sgp4_tests/aiaa-2006-6753/sgp4_tle_23599_result.txt"
    SGP4_results = readdlm(filename; comments=true)

    # Initialize the orbit propagator.
    orbp = init_orbit_propagator(Val(:sgp4), tle[1]; sgp4_gc = sgp4_gc_wgs72)

    # Propagate the orbit.
    t = SGP4_results[:,1]   # [min.]
    r_TEME, v_TEME = propagate_to_epoch!(orbp, tle[1].epoch .+ t/1440)

    # Compare the results.
    @inbounds for k = 1:length(t)
        # Assemble the result vector.
        st_sgp4_result = vcat(t[k], r_TEME[k]/1000, v_TEME[k]/1000)

        # Compare the values.
        #
        # It was necessary to relax the tolerance due to the precision in the
        # conversion from Julian Day to seconds.
        @test st_sgp4_result[1] == SGP4_results[k,1]
        @test st_sgp4_result[2] ≈  SGP4_results[k,2] atol=5e-4
        @test st_sgp4_result[3] ≈  SGP4_results[k,3] atol=5e-4
        @test st_sgp4_result[4] ≈  SGP4_results[k,4] atol=5e-4
        @test st_sgp4_result[5] ≈  SGP4_results[k,5] atol=1e-6
        @test st_sgp4_result[6] ≈  SGP4_results[k,6] atol=1e-6
        @test st_sgp4_result[7] ≈  SGP4_results[k,7] atol=1e-6
    end
end


text:
module Hubbard

using LinearMaps: FunctionMap

"""
    hops(n)

Computes neighbor bit masks for the `n` site square lattice with
periodic boundary conditions.

# Examples
```jldoctest
julia> [digits(h, base = 2, pad = 8) for h ∈ hops(8)]
16-element Array{Array{Int64,1},1}:
 [1, 0, 0, 0, 1, 0, 0, 0]
 [0, 1, 1, 0, 0, 0, 0, 0]
 [0, 0, 1, 1, 0, 0, 0, 0]
 [0, 0, 0, 1, 0, 0, 0, 1]
 [0, 0, 0, 0, 1, 1, 0, 0]
 [0, 0, 0, 0, 0, 1, 1, 0]
 [1, 0, 0, 0, 0, 0, 1, 0]
 [0, 1, 0, 0, 0, 0, 0, 1]
 [1, 0, 1, 0, 0, 0, 0, 0]
 [0, 1, 0, 0, 1, 0, 0, 0]
 [0, 0, 1, 0, 0, 1, 0, 0]
 [0, 0, 0, 1, 0, 0, 1, 0]
 [0, 0, 0, 1, 1, 0, 0, 0]
 [0, 0, 0, 0, 0, 1, 0, 1]
 [0, 1, 0, 0, 0, 0, 1, 0]
 [1, 0, 0, 0, 0, 0, 0, 1]

```
"""
function hops(n::Integer)::Vector
    a = first([u, v] for u=isqrt(n):-1:1 for v=0:u if sum(abs2, [u, v]) == n)
    b = [0 -1; 1 0]a  # orthogonal basis vector
    # apply periodic boundary conditions
    fold(x) = isequal(reduce((x, y) -> x - y * fld(x'y, n), [a, b], init = x))  # '
    inside(x) = fold(x)(x)
    sites = [[x, y] for y=0:sum(a) for x=b[1]:b[2] if inside([x, y])]
    masks(y) = [(1<<(i-1))|(1<<(findfirst(fold(x + y), sites)-1)) for (i, x) in enumerate(sites)]
    return [masks([1, 0]); masks([0, 1])]
end

"""
    index(m, c)

Compute the index of the combination of `count_ones(mask)` items
in lexicographic order

# Examples
```jldoctest
julia> c = [binomial(i, j) for i=0:3, j=1:2];

julia> [index(parse(Int, b, base = 2), c) for b ∈ ["0011", "0101", "0110", "1001", "1010", "1100"]]
6-element Array{Int64,1}:
 1
 2
 3
 4
 5
 6

```
"""
function index(m::Integer, c::Matrix)::Integer
    i = 1
    @inbounds for k ∈ 1:size(c, 2)
        i += c[trailing_zeros(m) + 1, k]
        m &= m - 1
    end
    return i
end

"""
    mask(i, c)

Compute the binary mask of the combination of `k` items out of `n`
at position `index` in lexicographic order

# Examples
```jldoctest
julia> c = [binomial(i, j) for i=0:3, j=1:2];

julia> [digits(mask(i, c), base = 2, pad = 4) for i=1:6]
6-element Array{Array{Int64,1},1}:
 [1, 1, 0, 0]
 [1, 0, 1, 0]
 [0, 1, 1, 0]
 [1, 0, 0, 1]
 [0, 1, 0, 1]
 [0, 0, 1, 1]

```
"""
function mask(i::Integer, c::Matrix)::Integer
    k = size(c, 2)
    m = 0
    @inbounds for n ∈ size(c, 1):-1:1
        j = i - c[n, k]
        if j ≥ 1
            m |= 1
            k == 1 && return m << (n - 1)
            i = j
            k -= 1
        end
        m <<= 1
    end
end

"""
    hamiltonian(n, k, t, U)

Implements the [Hubbard hamiltonian](https://en.wikipedia.org/wiki/Hubbard_model)
as a [mutating linear map](https://github.com/Jutho/LinearMaps.jl).

# Arguments
- `n::Integer`: sites in the square lattice with periodic boundary conditions
- `k::Integer`: spin up and spin down fermions
- `t::Real`: hopping integral
- `U::Real`: strength of on-site interaction

# Examples
```jldoctest
julia> h = hamiltonian(8, 4, 1., 4.)
LinearMaps.FunctionMap{Float64}(multiply!, 4900, 4900; ismutating=true, issymmetric=true, ishermitian=true, isposdef=false)

julia> using Arpack: eigs

julia> λ, ϕ = eigs(h, nev=1, which=:SR);

julia> λ
1-element Array{Float64,1}:
 -11.775702792136244

julia> size(ϕ)
(4900, 1)

```
"""
function hamiltonian(n::Integer, k::Integer, t::T, U::T) where T <: Real
    c = [binomial(i, j) for i=0:(n-1), j=1:k]
    nCk = binomial(n, k)
    hop = hops(n)
    function mult!(y, x)
        fill!(y, 0.0)
        xm = reshape(x, nCk, :)
        ym = reshape(y, nCk, :)
        @inbounds for i ∈ CartesianIndices(xm)
            w = xm[i]
            w == 0 && continue

            l, h = i.I
            u, d = mask(h, c), mask(l, c)
            ym[i] += U * w * count_ones(u & d)

            tw = t * w
            for p ∈ hop
                a = u ⊻ p
                if a ≠ u && count_ones(a) == k
                    ym[l, index(a, c)] -= tw
                end
                a = d ⊻ p
                if a ≠ d && count_ones(a) == k
                    ym[index(a, c), h] -= tw
                end
            end
        end
    end
    return FunctionMap{T}(mult!, abs2(nCk), ismutating=true, issymmetric=true)
end

end

max_stars_repo_name: jburgy/hubbard

max_issues_repo_name: jburgy/hubbard

max_forks_repo_name: jburgy/hubbard


text:
using Tabben: num_examples, num_inputs, num_outputs, num_classes, task,
has_extras, has_extra, license, bibtex

@testset "classification dataset tests" begin
    ns = [5, 10, 17, 100, 1000]
    ms = [5, 10, 17, 100, 1000, 5000]
    ks = [1, 2, 3, 5, 10]
    cs = [2, 3, 5, 10, 100]

    for n in ns, m in ms, k in ks, c in cs
        ds = random_classification_dataset(n, m, k, c)

        @test num_examples(ds) == n == length(ds) == lastindex(ds)
        @test num_inputs(ds) == m
        @test num_outputs(ds) == k
        @test num_classes(ds) == c

        @test task(ds) == "classification"

        @test has_extras(ds)

        for extra_name in ("license", "bibtex", "column-names-attributes", "column-names-target")
            @test has_extra(ds, extra_name)
        end
        @test !has_extra(ds, "other")
        @test !has_extra(ds, "blah")

        @test license(ds) == "test license"
        @test bibtex(ds) == "test bibtex"

        count = 0
        for i in 1:length(ds)
            single_input, single_output = ds[i]

            @test length(single_input) == m
            @test length(single_output) == k
            @test all(val in 0:c-1 for val in single_output)

            count += 1
        end
        @test count == n

        count = 0
        for (single_input, single_output) in ds
            @test length(single_input) == m
            @test length(single_output) == k
            @test all(val in 0:c-1 for val in single_output)

            count += 1
        end
        @test count == n
    end

end


text:
module Segmenter
    export Segment

    struct Segment
        type::Symbol
        x_start::Int
        x_end::Int
        value::Any
    end

    Base.length(segment::Segment) = segment.x_end - segment.x_start + 1

    function segment(type, y; segment_by_value = false)
        x_start = 1
        segments = []
        while true
            idx = findfirst(x->x !== nothing, y[x_start:end])
            idx !== nothing || break
            x_start += idx - 1

            x_end = findfirst(x->(segment_by_value ? x != y[x_start] : x === nothing), y[x_start:end])
            if x_end === nothing
                x_end = length(y) - x_start + 1
            else
                x_end -= 1
            end
            x_end += x_start - 1

            value = (segment_by_value ? y[x_start] : y[x_start:x_end])

            push!(segments, Segment(type, x_start, x_end, value))

            x_end != length(y) || break
            x_start = x_end + 1
        end
        return segments
    end

    function consolidate_segments(segments)
        idx = 1
        consolidated = []
        segments_to_merge = [idx]
        while idx <= length(segments)
            while true
                idx += 1
                idx <= length(segments) && segments[idx].x_start - segments[idx-1].x_end == 1 || break # is contiguous
                push!(segments_to_merge, idx)
            end

            if length(segments_to_merge) > 1
                x_start = segments[segments_to_merge[1]].x_start
                x_end = segments[segments_to_merge[end]].x_end
                value = zeros(length(x_start:x_end))
                for j in segments_to_merge
                    copyto!(value, segments[j].x_start - x_start + 1, [segments[j].value for i in 1:length(segments[j])], 1)
                end
                push!(consolidated, Segment(segments[segments_to_merge[1]].type, x_start, x_end, value))
            else
                push!(consolidated, segments[idx-1])
            end

            segments_to_merge = [idx]
        end
        return consolidated
    end
end


text:
immutable AirspeedType{str} end
kcas = AirspeedType{:kcas}()
keas = AirspeedType{:keas}()
ktas = AirspeedType{:ktas}()
mach = AirspeedType{:mach}()

type Airspeed{T}
    val::Float64
    Airspeed(v) = typeof(T)<:AirspeedType ? new{T}(v) : error("Invalid Airspeed Type")
end

show{T}(io::IO,::AirspeedType{T}) = print(io,T)
show{T}(io::IO,a::Airspeed{T}) = print(io,a.val,T)
*{T}(x::Real,t::AirspeedType{T}) = Airspeed{t}(x)

immutable Env
    atm::Atm
    calibrated_airspeed::@quantity(Float64,m/s)
    equivalent_airspeed::@quantity(Float64,m/s)
    true_airspeed::@quantity(Float64,m/s)
    mach_number::Float64
    q::@quantity(Float64,N/m^2)
end

function calc_calibrated(mach_number::Float64,pressure_ratio::Float64)
     sqrt(5*((((0.2*mach_number^2+1)^3.5 - 1)*pressure_ratio+1)^(1/3.5) - 1))*SeaLevelSpeedOfSound
end

function calc_mach(calibrated_airspeed::@quantity(Float64,m/s),pressure_ratio::Float64)
    sqrt(5*((((0.2*(calibrated_airspeed/SeaLevelSpeedOfSound)^2+1)^3.5 - 1)/pressure_ratio+1)^(1/3.5) - 1))
end

function Env(airspeed::Airspeed{kcas},altitude::@quantity(Float64,m),delta_temperature::@quantity(Float64,K) = 0.0K)
    atm = Atm(altitude,delta_temperature)
    calibrated_airspeed = convert(SIQuantity,airspeed.val*knots)
    mach_number = calc_mach(calibrated_airspeed,atm.pressure_ratio)
    equivalent_airspeed = atm.speed_of_sound*mach_number
    true_airspeed = equivalent_airspeed/sqrt(atm.density_ratio)
    q = 0.5*atm.density*true_airspeed^2
    Env(atm,calibrated_airspeed,equivalent_airspeed,true_airspeed,mach_number,q)
end

function Env(airspeed::Airspeed{mach},altitude::@quantity(Float64,m),delta_temperature::@quantity(Float64,K) = 0.0K)
    atm = Atm(altitude,delta_temperature)
    mach_number = airspeed.val
    calibrated_airspeed = calc_calibrated(mach_number,atm.pressure_ratio)
    equivalent_airspeed = atm.speed_of_sound*mach_number
    true_airspeed = equivalent_airspeed/sqrt(atm.density_ratio)
    q = 0.5*atm.density*true_airspeed^2
    Env(atm,calibrated_airspeed,equivalent_airspeed,true_airspeed,mach_number,q)
end

function Env(airspeed::Airspeed{keas},altitude::@quantity(Float64,m),delta_temperature::@quantity(Float64,K) = 0.0K)
    atm = Atm(altitude,delta_temperature)
    equivalent_Airspeed = convert(SIQuantity,airspeed.val*knots)
    mach_number = equivalent_airspeed/atm.speed_of_sound
    calibrated_airspeed = calc_calibrated(mach_number,atm.pressure_ratio)
    true_airspeed = equivalent_airspeed/sqrt(atm.density_ratio)
    q = 0.5*atm.density*true_airspeed^2
    Env(atm,calibrated_airspeed,equivalent_airspeed,true_airspeed,mach_number,q)
end

# function Env{N1<:Real,N2<:Real,U1<:LengthUnit,U2<:TemperatureUnit}(airspeed::Airspeed,altitude::Quantity{N1,U1},delta_temperature::Quantity{N2,U2} = 0.0K)
#     Env(airspeed,convert(SIQuantity{Float64},altitude),convert(SIQuantity{Float64},delta_temperature))
# end


text:
""" Functions that were used with NEID data products prior to shipping to KPNO """
module PreShip
using DataFrames, FITSIO

""" Read NEID (non-solar) data from FITS file, and return in a Spectra2DBasic object."""
function read_data   end

function read_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    @assert read_key(f[1],"SKY-OBJ")[1] != "Solar"
    λ, flux, var  = FITSIO.read(f["SCIWAVE"]), FITSIO.read(f["Sci Flux"]), FITSIO.read(f["Sci Variance"])
    metadata[:normalization] = :raw
    Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
end

function read_data(fn::String)
    f = FITS(fn)
    @assert read_key(f[1],"SKY-OBJ")[1] != "Solar"
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    λ, flux, var  = FITSIO.read(f["SCIWAVE"]), read(f["Sci Flux"]), read(f["Sci Variance"])
    metadata[:normalization] = :raw
    Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
end

function read_data(dfr::DataFrameRow{DataFrame,DataFrames.Index})
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_data(fn,metadata)
end

""" Read NEID Solar data from FITS file, and return in a Spectra2DBasic object."""
function read_solar_data
end

function read_solar_data(fn::String, metadata::Dict{Symbol,Any} )
    f = FITS(fn)
    @assert read_key(f[1],"SKY-OBJ")[1] == "Solar"
    λ, flux, var  = read(f["SKYWAVE"]), read(f["Sky Flux"]), read(f["Sky Variance"])
    Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
end

function read_solar_data(fn::String)
    f = FITS(fn)
    @assert read_key(f[1],"SKY-OBJ")[1] == "Solar"
    hdr = FITSIO.read_header(f[1])
    metadata = Dict(zip(map(k->Symbol(k),hdr.keys),hdr.values))
    λ, flux, var  = FITSIO.read(f["SKYWAVE"]), FITSIO.read(f["Sky Flux"]), FITSIO.ead(f["Sky Variance"])
    Spectra2DBasic(λ, flux, var, NEID2D(), metadata=metadata)
end

function read_solar_data(dfr::DataFrameRow{DataFrame,DataFrames.Index})
    fn = dfr.Filename
    metadata = Dict(zip(keys(dfr),values(dfr)))
    read_solar_data(fn,metadata)
end

""" Read CSV of NEID drift corrections, interpolate to bjd's in df and insert into df[:,drift]. """
function read_drift_corrections!(fn::String, df::DataFrame, df_time_col::Symbol = :bjd)
    drift_corrections = CSV.read(fn, DataFrame, header=["bjd", "sci_drift", "cal_drift"]);
    @assert any(isequal(:bjd),propertynames(drift_corrections))
    @assert any(isequal(:cal_drift),propertynames(drift_corrections))
    drift_interp = LinearInterpolation(drift_corrections[!,:bjd],drift_corrections[!,:cal_drift])
    df[!,:drift] = -drift_interp.(df[!,df_time_col])
    return df
end

""" Read CSV of NEID barycentric corrections, interpolate to bjd's in df and insert into df[:,ssb_rv]. """
function read_barycentric_corrections!(fn::String, df::DataFrame, df_time_col::Symbol = :bjd)
    ssb_corrections = CSV.read(fn, DataFrame, header=["bjd","rv_ssb"], select=[1,2], types=types=[Float64,Float64], datarow=2, silencewarnings=true);
    @assert any(isequal(:bjd),propertynames(ssb_corrections))
    @assert any(isequal(:rv_ssb),propertynames(ssb_corrections))
    ssb_interp = LinearInterpolation(ssb_corrections.bjd, ssb_corrections.rv_ssb)
    df[!,:ssb_rv] = ssb_interp.(df[!,df_time_col])
    return df
end

end  # module PreShip

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using jInv.Mesh
using jInv.Utils
using jInv.InverseSolve

using EikonalInv
using MAT

#############################################################################################################
plotting = false;
if plotting
	using jInvVis
	using  PyPlot
	close("all");
end
#############################################################################################################
dataDir = pwd();
resultsDir = pwd();

###############################################################################################################
# this filename can be used to save the model during the iterations. If set to "", no images or files will be saved.
# resultsFilenamePrefix = "";
resultsFilenamePrefix = string(resultsDir,"/shortTravelTimeInvSEG");
################################################################################################################

modelDir = pwd();
include("../drivers/prepareTravelTimeDataFiles.jl");
include("../drivers/readModelAndGenerateMeshMref.jl");
include("../drivers/setupTravelTimeTomography.jl");

########################################################################################################
######################################## for SEG 256X512 ###############################################
#######################################################################################################
dim     = 2;
pad     = 0;
jump    = 5;
newSize = [256,128];
offset  = ceil(Int64,(newSize[1]*(6/13.5)));;
(m,Minv,mref,boundsHigh,boundsLow) = readModelAndGenerateMeshMref(modelDir,"SEGmodel2Dsalt.dat",dim,pad,[0.0,13.5,0.0,4.2],newSize,1.752,2.9);
dataFilenamePrefix = string(dataDir,"/shortDATA_SEG",tuple((Minv.n+1)...));

limits = (minimum(m),maximum(m));

if plotting
	figure()
	plotModel(m,true,Minv,pad,limits);
	figure()
	plotModel(mref,true,Minv,pad,limits);
end

prepareTravelTimeDataFiles(m,Minv,mref,boundsHigh,boundsLow,dataFilenamePrefix,pad,jump,offset);

(Q,P,pMisRFs,SourcesSubInd,contDiv,Iact,sback,mref,boundsHigh,boundsLow,resultsFilename,Dobs) = setupTravelTimeTomography(m,dataFilenamePrefix, resultsFilenamePrefix);

########################################################################################################
##### Set up Inversion #################################################################################
########################################################################################################
############### Inversion for the velocity:
# modfun = velocityToSlowSquared;
################################################################################################################
############### USE A BOUND MODEL For velocity #################################################################
################################################################################################################
# maxStep=0.1*maximum(boundsHigh);
# a = minimum(boundsLow);
# b = maximum(boundsHigh);
# function modfun(m)
	# bm,dbm = getBoundModel(m,a,b);
	# bs,dbs = velocityToSlowSquared(bm);
	# dsdm = dbm*dbs;
	# return bs,dsdm;
# end
# mref = getBoundModelInv(mref,a,b);
# boundsHigh = boundsHigh*10000000.0;
# boundsLow = -boundsLow*100000000.0;
################################################################################################################
################ Inversion for the squared slowness: ###########################################################
################################################################################################################

mref 		= velocityToSlowSquared(mref)[1];
t    		= copy(boundsLow*0.9);
boundsLow 	= velocityToSlowSquared(boundsHigh*1.1)[1];
boundsHigh 	= velocityToSlowSquared(t)[1]; t = 0;
###################################################
a = minimum(boundsLow);
b = maximum(boundsHigh);
maxStep=0.1*maximum(boundsHigh);

############### Standard bounded optimization:
modfun 		= identityMod;

############### USE A BOUND MODEL For slow squared #############################################################

# modfun(m) = getBoundModel(m,a,b);
# mref = getBoundModelInv(mref,a,b);
# boundsHigh = boundsHigh*10000000.0;
# boundsLow = -boundsLow*100000000.0

################################################################################################################
################################################################################################################
################################################################################################################

cgit  = 8;
maxit = 10;
alpha = 5e-1;
pcgTol = 1e-1;

HesPrec=getSSORCGRegularizationPreconditioner(1.0,1e-5,100)

################################################# GIT VERSION OF JINV #################################################

regparams = [1.0,1.0,1.0,1e-6];
#### Use smoothness regularization
regfun(m, mref, M) = wdiffusionRegNodal(m, mref, M, Iact=Iact, C = regparams);
#### Use TV regularization
# regfun(m, mref, M) = wTVRegNodal(m, mref, M, Iact=Iact, C = regparams);

function dump(mc,Dc,iter,pInv,pMis)
	if resultsFilename!=""
		ntup = tuple((pInv.MInv.n+1)...);
		fullMc = slowSquaredToVelocity(reshape(Iact*pInv.modelfun(mc)[1] + sback,ntup))[1];
		Temp = splitext(resultsFilename);
		Temp = string(Temp[1],"_GN",iter,Temp[2]);
		writedlm(Temp,convert(Array{Float16},fullMc));
		if plotting
			close(888);
			figure(888);
			# plotModel(fullMc,false,[],0,limits,splitdir(Temp)[2]); # no information on axes
			plotModel(fullMc,true,Minv,pad,limits,splitdir(Temp)[2]); # with domain information on axes
		end
	end
end

pInv = getInverseParam(Minv,modfun,regfun,alpha,mref[:],boundsLow,boundsHigh,
                     maxStep=maxStep,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate=1e-3, maxIter = maxit,HesPrec=HesPrec);

D0RFs, = computeMisfit(pInv.modelfun(mref[:])[1],pMisRFs,false);

#### Projected Gauss Newton
mc,Dc,flag,his = projGNCG(copy(mref[:]),pInv,pMisRFs,indCredit = [],dumpResults = dump);
#### Barrier Gauss Newton
# mc,Dc,flag,his = barrierGNCG(copy(mref[:]),pInv,pMisRFs,indCredit = [],dumpResults = dump);

if plotting
	figure()
	Fvals = his.F;
	Jvals = his.Jc;
	numiter = length(Fvals);
	semilogy(0:numiter-1,Fvals,"-*");
	xlabel("Gauss Newton iterations")
	ylabel("Misfit value")
end

D0 = Array(Array{Float64,2},length(pMisRFs))
for k = 1:length(pMisRFs)
	D0[k] = fetch(D0RFs[k]);
end
D0 = arrangeRemoteCallDataIntoLocalData(D0); D0RFs = 0;

Dpred = Array(Array{Float64,2},length(pMisRFs))
for k = 1:length(pMisRFs)
	Dpred[k] = fetch(Dc[k]);
end
Dpred = arrangeRemoteCallDataIntoLocalData(Dpred); Dc = 0;

if plotting
	figure()
	imshow(Dpred');
	colorbar();
	xlabel("Sources")
	ylabel("Receivers")
	figure()
	imshow(D0');
	colorbar();
	xlabel("Sources")
	ylabel("Receivers")
	figure()
	imshow(abs(Dpred - Dobs)');
	colorbar();
	xlabel("Sources")
	ylabel("Receivers")
	figure()
	imshow(abs(D0 - Dobs)');
	colorbar();
	xlabel("Sources")
	ylabel("Receivers")
end

if resultsFilename!=""
	Temp = splitext(resultsFilename);
	writedlm(string(Temp[1],"_predictedData",Temp[2]),convert(Array{Float16},Dpred));
	writedlm(string(Temp[1],"_recoveredModel",Temp[2]),convert(Array{Float16},reshape(Iact*modfun(mc)[1] + sback,tuple((pInv.MInv.n+1)...))));
end

#############################################################################################


text:
using BinaryBuilder

# Collection of sources required to build libiconv
sources = [
    "https://ftp.gnu.org/pub/gnu/libiconv/libiconv-1.15.tar.gz" =>
    "ccf536620a45458d26ba83887a983b96827001e92a13847b45e4925cc8913178",
]

script = raw"""
cd $WORKSPACE/srcdir
cd libiconv-1.15/
./configure --prefix=$prefix --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless
# further platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = prefix -> [
    LibraryProduct(prefix,"libcharset",:libcharset),
    LibraryProduct(prefix,"libiconv",:libiconv),
    ExecutableProduct(prefix,"iconv",:iconv)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

build_tarballs(ARGS, "IConv", sources, script, platforms, products, dependencies)


text:
function neighpairs(x, q, pairlist, pairnum, atomtype, ilist, alist)

dim = size(x,1);
ncq = size(q,1);
inum = length(ilist);
ijnum = sum(pairnum[ilist.+1]-pairnum[ilist]);

ai = zeros(eltype(ilist),ijnum);
aj = zeros(eltype(ilist),ijnum);
ti = zeros(eltype(ilist),ijnum);
tj = zeros(eltype(ilist),ijnum);
xij = zeros(dim,ijnum);
qi = zeros(ncq,ijnum);
qj = zeros(ncq,ijnum);

count = 0;
for ii=1:inum
    i = ilist[ii]; # atom i
    n1 = pairnum[i];
    n2 = pairnum[i+1];
    n = n2-n1;
    g = pairlist[(n1+1):n2]; # a list of neighbors around atom i
    j = alist[g];
    ind = (count+1):(count+n);
    count = count + n;

    ai[ind] .= i;
    aj[ind] .= j;
    ti[ind] .= atomtype[i];
    tj[ind] .= atomtype[j];
    xij[:,ind] = x[:,g] .- x[:,i];

    if ncq>0
        qi[ind] .= q[i];
        qj[ind] .= q[j];
    end
end

return xij, ai, aj, ti, tj, qi, qj

end

function neighpairs(x, pairlist, pairnum, atomtype, ilist, alist)

    dim = size(x,1);
    inum = length(ilist);
    ijnum = sum(pairnum[ilist.+1]-pairnum[ilist]);

    ai = zeros(eltype(ilist),ijnum);
    aj = zeros(eltype(ilist),ijnum);
    ti = zeros(eltype(ilist),ijnum);
    tj = zeros(eltype(ilist),ijnum);
    xij = zeros(dim,ijnum);

    count = 0;
    for ii=1:inum
        i = ilist[ii]; # atom i
        n1 = pairnum[i];
        n2 = pairnum[i+1];
        n = n2-n1;
        g = pairlist[(n1+1):n2]; # a list of neighbors around atom i
        j = alist[g];
        ind = (count+1):(count+n);
        count = count + n;

        ai[ind] .= i;
        aj[ind] .= j;
        ti[ind] .= atomtype[i];
        tj[ind] .= atomtype[j];
        xij[:,ind] = x[:,g] .- x[:,i];

    end

    return xij, ai, aj, ti, tj

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
@hierarchical KernelMatrix BarycentricMatrix2D Matrix

import LinearAlgebra: matprod
function (*)(H::AbstractKernelMatrix{T}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(zeros(TS, size(H, 1)), H, x)
end
function (*)(H::AbstractKernelMatrix{T}, x::AbstractMatrix{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    mul!(zeros(TS, size(H, 1), size(x, 2)), H, x)
end

LinearAlgebra.mul!(u::Vector, H::AbstractKernelMatrix, v::AbstractVector) = mul!(u, H, v, 1, 1)

@generated function mul!(u::Vector{S}, H::KernelMatrix{S}, v::AbstractVector{S}, istart::Int, jstart::Int) where S
    L = length(fieldnames(H))-1
    T = fieldname(H, 1)
    str = "
    begin
        M, N = blocksize(H)
        p = 0
        for m = 1:M
            q = 0
            for n = 1:N
                Hmn = H.assigned[m,n]
                if Hmn == 1
                    mul!(u, getindex(H.$T, m, n), v, istart + p, jstart + q)"
    for l in 2:L
        T = fieldname(H, l)
        str *= "
                elseif Hmn == $l
                    mul!(u, getindex(H.$T, m, n), v, istart + p, jstart + q)"
    end
    str *= "
                end
                q += blocksize(H, 1, n, 2)
            end
            p += blocksize(H, m, N, 1)
        end
        return u
    end"
    return Meta.parse(str)
end

KernelMatrix(f::Function, x::Vector{T}, y::Vector{T}, a::T, b::T, c::T, d::T) where T = KernelMatrix(f, x, y, 1:length(x), 1:length(y), a, b, c, d)

function KernelMatrix(f::Function, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}, a::T, b::T, c::T, d::T) where T
    ir1, ir2 = indsplit(x, ir, a, b)
    jr1, jr2 = indsplit(y, jr, c, d)
    ab2 = half(T)*(a+b)
    cd2 = half(T)*(c+d)

    if length(ir1) < BLOCKSIZE(T) && length(ir2) < BLOCKSIZE(T) && length(jr1) < BLOCKSIZE(T) && length(jr2) < BLOCKSIZE(T)
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = T[f(x[i], y[j]) for i in ir1, j in jr1]
        H[Block(1), Block(2)] = T[f(x[i], y[j]) for i in ir1, j in jr2]
        H[Block(2), Block(1)] = T[f(x[i], y[j]) for i in ir2, j in jr1]
        H[Block(2), Block(2)] = T[f(x[i], y[j]) for i in ir2, j in jr2]
        H
    else
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = KernelMatrix(f, x, y, ir1, jr1, a, ab2, c, cd2)
        H[Block(1), Block(2)] = KernelMatrix1(f, x, y, ir1, jr2, a, ab2, cd2, d)
        H[Block(2), Block(1)] = KernelMatrix2(f, x, y, ir2, jr1, ab2, b, c, cd2)
        H[Block(2), Block(2)] = KernelMatrix(f, x, y, ir2, jr2, ab2, b, cd2, d)
        H
    end
end

function KernelMatrix1(f::Function, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}, a::T, b::T, c::T, d::T) where T
    ir1, ir2 = indsplit(x, ir, a, b)
    jr1, jr2 = indsplit(y, jr, c, d)
    ab2 = half(T)*(a+b)
    cd2 = half(T)*(c+d)

    if length(ir1) < BLOCKSIZE(T) && length(ir2) < BLOCKSIZE(T) && length(jr1) < BLOCKSIZE(T) && length(jr2) < BLOCKSIZE(T)
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = BarycentricMatrix2D(T, f, a, ab2, c, cd2, x, y, ir1, jr1)
        H[Block(1), Block(2)] = BarycentricMatrix2D(T, f, a, ab2, cd2, d, x, y, ir1, jr2)
        H[Block(2), Block(1)] = T[f(x[i], y[j]) for i in ir2, j in jr1]
        H[Block(2), Block(2)] = BarycentricMatrix2D(T, f, ab2, b, cd2, d, x, y, ir2, jr2)
        H
    else
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = BarycentricMatrix2D(T, f, a, ab2, c, cd2, x, y, ir1, jr1)
        H[Block(1), Block(2)] = BarycentricMatrix2D(T, f, a, ab2, cd2, d, x, y, ir1, jr2)
        H[Block(2), Block(1)] = KernelMatrix1(f, x, y, ir2, jr1, ab2, b, c, cd2)
        H[Block(2), Block(2)] = BarycentricMatrix2D(T, f, ab2, b, cd2, d, x, y, ir2, jr2)
        H
    end
end

function KernelMatrix2(f::Function, x::Vector{T}, y::Vector{T}, ir::UnitRange{Int}, jr::UnitRange{Int}, a::T, b::T, c::T, d::T) where T
    ir1, ir2 = indsplit(x, ir, a, b)
    jr1, jr2 = indsplit(y, jr, c, d)
    ab2 = half(T)*(a+b)
    cd2 = half(T)*(c+d)

    if length(ir1) < BLOCKSIZE(T) && length(ir2) < BLOCKSIZE(T) && length(jr1) < BLOCKSIZE(T) && length(jr2) < BLOCKSIZE(T)
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = BarycentricMatrix2D(T, f, a, ab2, c, cd2, x, y, ir1, jr1)
        H[Block(1), Block(2)] = T[f(x[i], y[j]) for i in ir1, j in jr2]
        H[Block(2), Block(1)] = BarycentricMatrix2D(T, f, ab2, b, c, cd2, x, y, ir2, jr1)
        H[Block(2), Block(2)] = BarycentricMatrix2D(T, f, ab2, b, cd2, d, x, y, ir2, jr2)
        H
    else
        H = KernelMatrix(T, 2, 2)
        H[Block(1), Block(1)] = BarycentricMatrix2D(T, f, a, ab2, c, cd2, x, y, ir1, jr1)
        H[Block(1), Block(2)] = KernelMatrix2(f, x, y, ir1, jr2, a, ab2, cd2, d)
        H[Block(2), Block(1)] = BarycentricMatrix2D(T, f, ab2, b, c, cd2, x, y, ir2, jr1)
        H[Block(2), Block(2)] = BarycentricMatrix2D(T, f, ab2, b, cd2, d, x, y, ir2, jr2)
        H
    end
end


text:
@testset "coercion between cyclotomic fields" begin
  F2, z2 = CyclotomicField(2)
  @test Hecke.force_coerce_cyclo(F2, z2) == z2

  F1, z1 = CyclotomicField(1)
  up = Hecke.force_coerce_cyclo(F2, z1)
  @test Hecke.force_coerce_cyclo(F1, up) == z1

  choices = fmpz[collect(-5:5)...]

  # coerce first up and then down
  for n in 1:15
    Fn, zn = CyclotomicField(n)
    for m in 1:15
      nm = n*m
      Fnm, znm = CyclotomicField(nm)
      x = rand(Fn, choices)
      x_up = Hecke.force_coerce_cyclo(Fnm, x)
      x_down = Hecke.force_coerce_cyclo(Fn, x_up)
      @test x_down == x
    end
  end

  # coerce first down and then up
  for n in 1:15
    Fn, zn = CyclotomicField(n)
    for g in divisors(n)
      Fg, zg = CyclotomicField(g)
      for m in 1:15
        gm = g*m
        Fgm, zgm = CyclotomicField(gm)
        x = rand(Fg, choices)
        x_up = Hecke.force_coerce_cyclo(Fgm, x)
        x_n = Hecke.force_coerce_cyclo(Fn, x_up)
        @test x_n == Hecke.force_coerce_cyclo(Fn, x)
      end
    end
  end

  # impossible coercions
  for n in 1:45
    Fn, zn = CyclotomicField(n)
    for m in 1:45
      if n % m != 0 && ! (isodd(n) && (2*n) % m == 0)
        Fm, zm = CyclotomicField(m)
        @test_throws ErrorException Hecke.force_coerce_cyclo(Fn, zm)
        @test Hecke.force_coerce_cyclo(Fn, zm, Val{false}) == nothing
      end
    end
  end

end


text:
module Overlap

using LinearAlgebra: dot
using StaticArrays: SVector
using OffsetArrays: OffsetArray
using Common: Float, Coordinate
using BasisFunctions: Gaussian

export buildCartesianOverlaps

"""
    gaussianOverlap(Rμ::Coordinate, ζμ::Float, Rν::Coordinate, ζν::Float)

Calculates the overlap integral ∫ μ⋅ν dr between the gaussian functions μ and ν having centers Rμ, Rν and exponents ζμ, ζν.
"""
function gaussianOverlap(Rμ::SVector{3, Float}, ζμ::Float, Rν::SVector{3, Float}, ζν::Float)
    dR = Rμ - Rν
    p = ζμ + ζν
    return sqrt(π / p) * exp(- ζμ * ζν / p * dot(dR, dR))
end

"""
    buildCartesianOverlaps(Rμ::SVector{3}, ζμ::Float, Rν::SVector{3}, ζν::Float, lμ::Int, lν::Int)

Calculates all of the primitive cartesian overlap integrals up to angular momentums lμ, lν between gaussians at centers Rμ, Rν with exponents ζμ, ζν.
"""
function buildCartesianOverlaps(Rμ::SVector{3, Float}, ζμ::Float, lμ::Int, Rν::SVector{3, Float}, ζν::Float, lν::Int)::AbstractArray
    p = ζμ + ζν
    twop = 0.5 / p
    S = OffsetArray(Array(zeros(Float, lμ+1, lν+1, 3)), -1, -1, 0)
    S[0,0,:] .= gaussianOverlap(Rμ, ζμ, Rν, ζν)
    for n=1:3 #  Over cartesian directions
        P = (ζμ * Rμ[n] + ζν * Rν[n]) / p
        XPμ = P - Rμ[n]
        XPν = P - Rν[n]
        if lμ > 0
            S[1,0,n] = XPμ * S[0,0,n]
        end
        if lν > 0
            S[0,1,n] = XPν * S[0,0,n]
        end
        if lμ > 0 && lν > 0
            S[1,1,n] = XPμ * S[0,1,n] + twop * S[0,0,n]
        end
        for i=1:lμ-1
            S[i+1,0,n] = XPμ * S[i,0,n] + twop * i * S[i-1,0,n]
        end
        for j=1:lν-1
            S[0,j+1,n] = XPν * S[0,j,n] + twop * j * S[0,j-1,n]
        end
        if lμ > 0
            for j=1:lν-1
                S[1,j+1,n] = XPν * S[1,j,n] + twop * (S[0,j,n] + j * S[1,j-1,n])
            end
        end
        for i =1:lμ-1
            for j=1:lν
                S[i+1,j,n] = XPμ * S[i,j,n] + twop * (i * S[i-1,j,n] + j * S[i,j-1,n])
            end
        end
    end
    return S
end

end


text:
using BadukGoWeiqiTools
using Test

@testset "BadukGoWeiqiTools.jl" begin

end


text:
include("VelocityScale.jl")

VelocityScale(125,750,-17.07,0.223,300.9113701541558,0.04974920123242903)


text:
module ImageFiltering

using FFTW
using ImageCore, MappedArrays, FFTViews, OffsetArrays, StaticArrays, ComputationalResources, TiledIteration
using Statistics, LinearAlgebra
using ColorVectorSpace  # for filtering RGB arrays
using Base: Indices, tail, fill_to_length, @pure, depwarn, @propagate_inbounds
using OffsetArrays: IdentityUnitRange   # using the one in OffsetArrays makes this work with multiple Julia versions
using Requires

export Kernel, KernelFactors,
    Pad, Fill, Inner, NA, NoPad,
    BorderArray,
    Algorithm,
    imfilter, imfilter!,
    mapwindow, mapwindow!,
    imgradients, padarray, centered, kernelfactors, reflect,
    freqkernel, spacekernel

FixedColorant{T<:Normed} = Colorant{T}
StaticOffsetArray{T,N,A<:StaticArray} = OffsetArray{T,N,A}
OffsetVector{T} = OffsetArray{T,1}

# Add a fix that should have been included in julia-1.0.3
if isdefined(Broadcast, :_sametype) && !isdefined(Broadcast, :axistype)
    axistype(a::T, b::T) where T = a
    axistype(a, b) = UnitRange{Int}(a)
    Broadcast._bcs1(a, b) = Broadcast._bcsm(b, a) ? axistype(b, a) : (Broadcast._bcsm(a, b) ? axistype(a, b) : throw(DimensionMismatch("arrays could not be broadcast to a common size")))
end

# Needed for type-stability
function Base.transpose(A::StaticOffsetArray{T,2}) where T
    inds1, inds2 = axes(A)
    OffsetArray(transpose(parent(A)), inds2, inds1)
end

module Algorithm
    # deliberately don't export these, but it's expected that they
    # will be used as Algorithm.FFT(), etc.
    abstract type Alg end
    "Filter using the Fast Fourier Transform" struct FFT <: Alg end
    "Filter using a direct algorithm" struct FIR <: Alg end
    "Cache-efficient filtering using tiles" struct FIRTiled{N} <: Alg
        tilesize::Dims{N}
    end
    "Filter with an Infinite Impulse Response filter" struct IIR <: Alg end
    "Filter with a cascade of mixed types (IIR, FIR)" struct Mixed <: Alg end

    FIRTiled() = FIRTiled(())
end
using .Algorithm: Alg, FFT, FIR, FIRTiled, IIR, Mixed

Alg(r::AbstractResource{A}) where {A<:Alg} = r.settings

include("utils.jl")
include("kernelfactors.jl")
using .KernelFactors: TriggsSdika, IIRFilter, ReshapedOneD, iterdims, kernelfactors

ReshapedVector{T,N,Npre,V<:AbstractVector} = ReshapedOneD{T,N,Npre,V}
ArrayType{T} = Union{AbstractArray{T}, ReshapedVector{T}}
ReshapedIIR{T,N,Npre,V<:IIRFilter} = ReshapedOneD{T,N,Npre,V}
AnyIIR{T} = Union{IIRFilter{T}, ReshapedIIR{T}}
ArrayLike{T} = Union{ArrayType{T}, AnyIIR{T}}

include("kernel.jl")
using .Kernel
using .Kernel: Laplacian, reflect, ando3, ando4, ando5, scharr, bickley, prewitt, sobel, gabor

NDimKernel{N,K} = Union{AbstractArray{K,N},ReshapedOneD{K,N},Laplacian{N}}

include("border.jl")
include("borderarray.jl")

BorderSpec{T} = Union{Pad{0}, Fill{T,0}, Inner{0}}
BorderSpecNoNa{T} = Union{Pad{0}, Fill{T,0}, Inner{0}}
const BorderSpecAny = Union{BorderSpec,NA,NoPad}

const ProcessedKernel = Tuple

include("imfilter.jl")
include("specialty.jl")

include("mapwindow.jl")
using .MapWindow

function __init__()
    # See ComputationalResources README for explanation
    push!(LOAD_PATH, dirname(@__FILE__))
    # if haveresource(ArrayFireLibs)
    #     @eval using DummyAF
    # end
    pop!(LOAD_PATH)
    @require AxisArrays="39de3d68-74b9-583c-8d2d-e117c070f3a9" begin
        centered(ax::AxisArrays.Axis{name}) where name = AxisArrays.Axis{name}(centered(ax.val))
        centered(a::AxisArrays.AxisArray) = AxisArrays.AxisArray(centered(a.data), centered.(AxisArrays.axes(a)))
    end
    @require ImageMetadata="bc367c6b-8a6b-528e-b4bd-a4b897500b49" begin
        centered(a::ImageMetadata.ImageMeta) = ImageMetadata.ImageMeta(centered(ImageMetadata.arraydata(a)), ImageMetadata.properties(a))
    end
end

end # module


text:
module MLJInterface

import MLJModelInterface
import LightGBM

const LGBM_METRICS = (
    "None", "l1", "l2", "rmse", "quantile", "mape", "huber", "fair", "poisson", "gamma", "gamma_deviance",
    "tweedie", "ndcg", "lambdarank", "map", "mean_average_precision", "auc", "average_precision", "binary_logloss",
    "binary_error", "auc_mu", "multi_logloss", "multi_error", "cross_entropy", "xentropy", "multi_logloss","cross_entropy_lambda",
    "xentlambda", "kullback_leibler", "kldiv",
)

const CLASSIFICATION_OBJECTIVES = (
    "binary", "multiclass", "softmax",
)

const REGRESSION_OBJECTIVES = (
    "regression", "regression_l2", "l2", "mean_squared_error", "mse", "l2_root", "root_mean_squared_error", "rmse",
    "regression_l1", "l1", "mean_absolute_error", "mae", "huber", "fair", "poisson", "quantile", "mape",
    "mean_absolute_percentage_error", "gamma", "tweedie",
)

MLJModelInterface.@mlj_model mutable struct LGBMRegressor <: MLJModelInterface.Deterministic

    # Hyperparameters, see https://lightgbm.readthedocs.io/en/latest/Parameters.html for defaults
    boosting::String = "gbdt"::(_ in ("gbdt", "goss", "rf", "dart"))
    num_iterations::Int = 10::(_ >= 0)
    learning_rate::Float64 = 0.1::(_ > 0.)
    num_leaves::Int = 31::(1 < _ <= 131072)
    max_depth::Int = -1;#::(_ != 0);
    tree_learner::String = "serial"::(_ in ("serial", "feature", "data", "voting"))
    histogram_pool_size::Float64 = -1.0;#::(_ != 0.0);
    min_data_in_leaf::Int = 20::(_ >= 0)
    min_sum_hessian_in_leaf::Float64 = 1e-3::(_ >= 0.0)
    max_delta_step::Float64 = 0.0
    lambda_l1::Float64 = 0.0::(_ >= 0.0)
    lambda_l2::Float64 = 0.0::(_ >= 0.0)
    min_gain_to_split::Float64 = 0.0::(_ >= 0.0)
    feature_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    feature_fraction_bynode::Float64 = 1.0::(0.0 < _ <= 1.0)
    feature_fraction_seed::Int = 2
    bagging_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    bagging_freq::Int = 0::(_ >= 0)
    bagging_seed::Int = 3
    early_stopping_round::Int = 0
    extra_trees::Bool = false
    extra_seed::Int = 6
    max_bin::Int = 255::(_ > 1)
    bin_construct_sample_cnt = 200000::(_ > 0)
    init_score::String = ""
    drop_rate::Float64 = 0.1::(0.0 <= _ <= 1.0)
    max_drop::Int = 50
    skip_drop:: Float64 = 0.5::(0.0 <= _ <= 1)
    xgboost_dart_mode::Bool
    uniform_drop::Bool
    drop_seed::Int = 4
    top_rate::Float64 = 0.2::(0.0 <= _ <= 1.0)
    other_rate::Float64 = 0.1::(0.0 <= _ <= 1.0)
    min_data_per_group::Int = 100::(_ > 0)
    max_cat_threshold::Int = 32::(_ > 0)
    cat_l2::Float64 = 10.0::(_ >= 0)
    cat_smooth::Float64 = 10.0::(_ >= 0)

    # Model properties
    objective::String = "regression"::(_ in REGRESSION_OBJECTIVES)
    categorical_feature::Vector{Int} = Vector{Int}()
    data_random_seed::Int = 1
    is_sparse::Bool = true
    is_unbalance::Bool = false
    boost_from_average::Bool = true
    use_missing::Bool = true

    # Metrics
    metric::Vector{String} = ["l2"]::(all(in.(_, (LGBM_METRICS, ))))
    metric_freq::Int = 1::(_ > 0)
    is_training_metric::Bool = false
    ndcg_at::Vector{Int} = Vector{Int}([1, 2, 3, 4, 5])::(all(_ .> 0))

    # Implementation parameters
    num_machines::Int = 1::(_ > 0)
    num_threads::Int  = 0::(_ >= 0)
    local_listen_port::Int = 12400::(_ > 0)
    time_out::Int = 120::(_ > 0)
    machine_list_file::String = ""
    save_binary::Bool = false
    device_type::String = "cpu"::(_ in ("cpu", "gpu"))
    force_col_wise::Bool = false
    force_row_wise::Bool = false

end

MLJModelInterface.@mlj_model mutable struct LGBMClassifier <: MLJModelInterface.Probabilistic

    # Hyperparameters, see https://lightgbm.readthedocs.io/en/latest/Parameters.html for defaults
    boosting::String = "gbdt"::(_ in ("gbdt", "goss", "rf", "dart"))
    num_iterations::Int = 10::(_ >= 0)
    learning_rate::Float64 = 0.1::(_ > 0.)
    num_leaves::Int = 31::(1 < _ <= 131072)
    max_depth::Int = -1;#::(_ != 0);
    tree_learner::String = "serial"::(_ in ("serial", "feature", "data", "voting"))
    histogram_pool_size::Float64 = -1.0;#::(_ != 0.0);
    min_data_in_leaf::Int = 20::(_ >= 0)
    min_sum_hessian_in_leaf::Float64 = 1e-3::(_ >= 0.0)
    max_delta_step::Float64 = 0.0
    lambda_l1::Float64 = 0.0::(_ >= 0.0)
    lambda_l2::Float64 = 0.0::(_ >= 0.0)
    min_gain_to_split::Float64 = 0.0::(_ >= 0.0)
    feature_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    feature_fraction_bynode::Float64 = 1.0::(0.0 < _ <= 1.0)
    feature_fraction_seed::Int = 2
    bagging_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    pos_bagging_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    neg_bagging_fraction::Float64 = 1.0::(0.0 < _ <= 1.0)
    bagging_freq::Int = 0::(_ >= 0)
    bagging_seed::Int = 3
    early_stopping_round::Int = 0
    extra_trees::Bool = false
    extra_seed::Int = 6
    max_bin::Int = 255::(_ > 1)
    bin_construct_sample_cnt = 200000::(_ > 0)
    init_score::String = ""
    drop_rate::Float64 = 0.1::(0.0 <= _ <= 1.0)
    max_drop::Int = 50
    skip_drop:: Float64 = 0.5::(0.0 <= _ <= 1)
    xgboost_dart_mode::Bool
    uniform_drop::Bool
    drop_seed::Int = 4
    top_rate::Float64 = 0.2::(0.0 <= _ <= 1.0)
    other_rate::Float64 = 0.1::(0.0 <= _ <= 1.0)
    min_data_per_group::Int = 100::(_ > 0)
    max_cat_threshold::Int = 32::(_ > 0)
    cat_l2::Float64 = 10.0::(_ >= 0)
    cat_smooth::Float64 = 10.0::(_ >= 0)

    # For documentation purposes: A calibration scaling factor for the output probabilities for binary and multiclass OVA
    # Not included above because this is only present for the binary model in the FFI wrapper, hence commented out
    # sigmoid::Float64 = 1.0::(_ > 0.0 )

    # Model properties
    objective::String = "multiclass"::(_ in CLASSIFICATION_OBJECTIVES)
    categorical_feature::Vector{Int} = Vector{Int}();
    data_random_seed::Int = 1
    is_sparse::Bool = true
    is_unbalance::Bool = false
    boost_from_average::Bool = true
    scale_pos_weight = 1.0
    use_missing::Bool = true

    # Metrics
    metric::Vector{String} = ["None"]::(all(in.(_, (LGBM_METRICS, ))))
    metric_freq::Int = 1::(_ > 0)
    is_training_metric::Bool = false
    ndcg_at::Vector{Int} = Vector{Int}([1, 2, 3, 4, 5])::(all(_ .> 0))

    # Implementation parameters
    num_machines::Int = 1::(_ > 0)
    num_threads::Int  = 0::(_ >= 0)
    local_listen_port::Int = 12400::(_ > 0)
    time_out::Int = 120::(_ > 0)
    machine_list_file::String = ""
    save_binary::Bool = false
    device_type::String = "cpu"::(_ in ("cpu", "gpu"))
    force_col_wise::Bool = false
    force_row_wise::Bool = false

end

MODELS = Union{LGBMClassifier, LGBMRegressor}

function mlj_to_kwargs(model::LGBMRegressor)

    return Dict{Symbol, Any}(
        name => getfield(model, name)
        for name in fieldnames(typeof(model))
    )

end

function mlj_to_kwargs(model::LGBMClassifier, classes)

    num_class = length(classes)
    if model.objective == "binary"
        if num_class != 2
            throw(ArgumentError("binary objective number of classes $num_class greater than 2"))
        end
        # munge num_class for LightGBM
        num_class = 1
    end

    retval = Dict{Symbol, Any}(
        name => getfield(model, name)
        for name in fieldnames(typeof(model))
    )

    retval[:num_class] = num_class
    return retval

end

# X and y and w must be untyped per MLJ docs
function fit(mlj_model::MODELS, verbosity::Int, X, y, w=AbstractFloat[])

    # MLJ docs are clear that 0 means silent. but 0 in LightGBM world means "warnings"
    # and < 0 means fatal logs only, so we put intended silence to -1 (which is probably the closest we get)
    verbosity = if verbosity == 0; -1 else verbosity end

    y_lgbm, classes = prepare_targets(y, mlj_model)
    model = model_init(mlj_model, classes)
    X = MLJModelInterface.matrix(X)
    # The FFI wrapper wants Float32 for these
    w = Float32.(w)
    train_metrics = LightGBM.fit!(model, X, y_lgbm; verbosity=verbosity, weights=w)

    fitresult = (model, classes, deepcopy(mlj_model))
    # because update needs access to the older version of training metrics we keep them in the cache
    # so the update can merge old and additional metrics as necessary.
    cache = (num_boostings_done=[LightGBM.get_iter_number(model)], training_metrics=train_metrics)
    report = user_fitreport(model, train_metrics)

    return (fitresult, cache, report)

end

function update(mlj_model::MLJInterface.MODELS, verbosity::Int, fitresult, cache, X, y, w=AbstractFloat[])

    old_lgbm_model, old_classes, old_mlj_model = fitresult

    # we can continue boosting if and only if num_iterations has changed
    if !MLJModelInterface.is_same_except(old_mlj_model, mlj_model, :num_iterations)
        # if we get here it's just means we need to call fit directly
        return MLJInterface.fit(mlj_model, verbosity, X, y, w)
    end

    additional_iterations = mlj_model.num_iterations - old_mlj_model.num_iterations

    if additional_iterations < 0
        # less iterations isn't very valid so re-fit from scratch
        # TODO: I think theres a LightGBM API where you can prune
        # the boosting, so we would just do that instead of wasting time re-fitting
        # although we'd still have to consider fit-report etc
        return MLJInterface.fit(mlj_model, verbosity, X, y, w)
    end

    if verbosity >= 1
        @info("Not refitting from scratch", additional_iterations)
    end

    # splice the data into the estimator -- we need to update num_iterations,
    # taking into account it may have stopped early previously
    # It might well stop early again too, but give it a chance
    # Also ideally early stopping would be implemented via this mechanism and MLJ anyway
    num_iterations = sum(cache.num_boostings_done)
    old_lgbm_model.num_iterations = num_iterations + additional_iterations

    # eh this is ugly, possibly prompts a need for some refactoring
    train_metrics = LightGBM.train!(old_lgbm_model, additional_iterations, String[], verbosity, LightGBM.Dates.now())
    fitresult = (old_lgbm_model, old_classes, deepcopy(mlj_model))

    final_num_iter = LightGBM.get_iter_number(old_lgbm_model)
    iteration_history = deepcopy(cache.num_boostings_done)
    push!(iteration_history, final_num_iter - num_iterations)
    report = user_fitreport(old_lgbm_model, cache.training_metrics, train_metrics)
    newcache = (num_boostings_done=iteration_history, training_metrics=report.training_metrics)

    return (fitresult, newcache, report)

end

# This does prep for classification tasks
@inline function prepare_targets(targets, model::LGBMClassifier)

    classes = MLJModelInterface.classes(first(targets))
    # -1 because these will be 1,2 and LGBM uses the 0/1 boundary
    # -This also works for multiclass because the classes are 0 indexed
    targets = Float64.(MLJModelInterface.int(targets) .- 1)
    return targets, classes

end
# This does prep for Regression, which is basically a no-op (or rather, just creation of an empty placeholder classes object
@inline prepare_targets(targets::AbstractVector, model::LGBMRegressor) = targets, []

function predict_classifier((fitted_model, classes, _), Xnew)

    Xnew = MLJModelInterface.matrix(Xnew)
    predicted = LightGBM.predict(fitted_model, Xnew)
    # when the objective == binary, lightgbm internally has classes = 1 and spits out only probability of positive class
    if size(predicted, 2) == 1
        predicted = hcat(1. .- predicted, predicted)
    end

    return MLJModelInterface.UnivariateFinite(classes, predicted)

end

function predict_regression((fitted_model, classes, _), Xnew)

    Xnew = MLJModelInterface.matrix(Xnew)
    return dropdims(LightGBM.predict(fitted_model, Xnew), dims=2)

end

# This function returns a user accessible report and therefore needs to not be a source of breaking changes
# Keep all of the functinality required to generate this into a single function so that it can be
# kept under control via tests and so on. It's not a lot for now but is likely to grow.
function user_fitreport(estimator::LightGBM.LGBMEstimator, eval_metrics::Dict)

    importance = (gain = LightGBM.gain_importance(estimator), split = LightGBM.split_importance(estimator))
    return (training_metrics = deepcopy(eval_metrics), importance = importance)

end
function user_fitreport(estimator::LightGBM.LGBMEstimator, eval_metrics::Dict, new_eval_metrics::Dict)

    metrics = LightGBM.merge_scores(eval_metrics, new_eval_metrics)
    return user_fitreport(estimator, metrics)

end

# multiple dispatch the various signatures for each model and args combo
MLJModelInterface.fit(model::MLJInterface.MODELS, verbosity::Int, X, y) = MLJInterface.fit(model, verbosity, X, y)
MLJModelInterface.fit(model::MLJInterface.MODELS, verbosity::Int, X, y, w::Nothing) = MLJInterface.fit(model, verbosity, X, y)
MLJModelInterface.fit(model::MLJInterface.MODELS, verbosity::Int, X, y, w) = MLJInterface.fit(model, verbosity, X, y, w)

# Add `update` method
MLJModelInterface.update(model::MLJInterface.MODELS, verbosity::Int, fitresult, cache, X, y) = MLJInterface.update(model, verbosity, fitresult, cache, X, y)
MLJModelInterface.update(model::MLJInterface.MODELS, verbosity::Int, fitresult, cache, X, y, w::Nothing) = MLJInterface.update(model, verbosity, fitresult, cache, X, y)
MLJModelInterface.update(model::MLJInterface.MODELS, verbosity::Int, fitresult, cache, X, y, w) = MLJInterface.update(model, verbosity, fitresult, cache, X, y, w)

MLJModelInterface.predict(model::MLJInterface.LGBMClassifier, fitresult, Xnew) = MLJInterface.predict_classifier(fitresult, Xnew)
MLJModelInterface.predict(model::MLJInterface.LGBMRegressor, fitresult, Xnew) = MLJInterface.predict_regression(fitresult, Xnew)

# multiple dispatch the model initialiser functions
model_init(mlj_model::MLJInterface.LGBMClassifier, classes) = LightGBM.LGBMClassification(; mlj_to_kwargs(mlj_model, classes)...)
model_init(mlj_model::MLJInterface.LGBMRegressor, targets) = LightGBM.LGBMRegression(; mlj_to_kwargs(mlj_model)...)

# metadata
MLJModelInterface.metadata_pkg.(
    (LGBMClassifier, LGBMRegressor); # end positional args
    name="LightGBM",
    uuid="7acf609c-83a4-11e9-1ffb-b912bcd3b04a",
    url="https://github.com/IQVIA-ML/LightGBM.jl",
    julia=false,
    license="MIT Expat",
    is_wrapper=false,
)

MLJModelInterface.metadata_model(
    LGBMClassifier; # end positional args
    path="LightGBM.MLJInterface.LGBMClassifier",
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{<:MLJModelInterface.Finite},
    weights=true,
    descr="Microsoft LightGBM FFI wrapper: Classifier",
)

MLJModelInterface.metadata_model(
    LGBMRegressor; # end positional args
    path="LightGBM.MLJInterface.LGBMRegressor",
    input=MLJModelInterface.Table(MLJModelInterface.Continuous),
    target=AbstractVector{MLJModelInterface.Continuous},
    weights=true,
    descr="Microsoft LightGBM FFI wrapper: Regressor",
)

end # module


text:
function alloc_check(
    title,
    dependencies,
    targets,
    pre_alloc,
    alloc;
    path=pwd(),
    labeller=:version,
    threads=nothing,
)
    @info "Tracking allocations: $title"

    # cd to path if valid
    isdir(path) && cd(path)

    # add a proc (id == p) that track allocations
    p = first(
        if isnothing(threads)
            addprocs(1; exeflags=["--track-allocation=user", "--project=$path"])
        else
            addprocs(
                1; exeflags=["--track-allocation=user", "--project=$path", "-t $threads"]
            )
        end,
    )

    @eval @everywhere $p using Pkg
    @eval @everywhere $p Pkg.instantiate()
    @eval @everywhere $p using Profile

    # @info read("Project.toml", String)
    # @warn read("Manifest.toml", String)

    for d in dependencies
        @eval @everywhere $p using $(Symbol(d))
    end

    @eval @everywhere $p $pre_alloc()
    @eval @everywhere $p Profile.clear_malloc_data()
    @eval @everywhere $p $alloc()

    rmprocs(p)

    myallocs = Coverage.analyze_malloc(map(dirname ∘ pathof ∘ eval, targets))

    for t in dependencies, d in walkdir(dirname(pathof(t))), f in d[end]
        splitext(f)[2] == ".mem" && rm(joinpath(d[1], f))
    end
    for d in walkdir(path), f in d[end]
        splitext(f)[2] == ".mem" && rm(joinpath(d[1], f))
    end

    if isempty(myallocs)
        @warn "No allocations was found in " targets
        return nothing
    end

    # Smart paths
    common, specifics = smart_paths(map(a -> a.filename, Iterators.reverse(myallocs)))
    # @info "sizes" map(a -> a.filename, Iterators.reverse(myallocs)) specifics
    slash = Sys.iswindows() ? "\\" : "/"

    # Make the allocations data readable through a dataframe
    df = DataFrame()
    df.bytes = map(a -> a.bytes, Iterators.reverse(myallocs))
    df[!, "ratio (%)"] = round.(df.bytes / sum(df.bytes) * 100; digits=2)
    df[!, "filename: [$common$slash"] = map(first ∘ splitext ∘ first ∘ splitext, specifics)
    df.linenumber = map(a -> a.linenumber, Iterators.reverse(myallocs))

    # Save it as a CSV file
    label = ""
    if labeller == :oid
        label = oid2string(map(p -> joinpath(dirname(pathof(p)), ".."), targets))
    elseif labeller == :version
        label = version2string(map(p -> joinpath(dirname(pathof(p)), ".."), targets))
    end
    mkpath("mallocs")
    CSV.write(joinpath(path, "mallocs/mallocs$label.csv"), df)

    # Visualize a pretty table
    pretty_table(df)

    return nothing
end

function pie_filter(df, threshold=5.0)
    i = findfirst(x -> x < threshold, df[!, 2])
    X = df[1:(i - 1), 3] .* " line " .* map(string, df[1:(i - 1), 4])
    push!(X, "others .< $threshold")
    Y = df[1:(i - 1), 2]
    push!(Y, sum(df[i:end, 2]))
    return X, Y
end

function alloc_plot(
    targets; formats=["pdf", "svg", "png"], backend=Plots.GRBackend, seriestype=:step
)
    backend()
    for target in targets
        path = normpath(joinpath(dirname(pathof(target)), "..", "perf", "mallocs"))
        versions = Vector{VersionNumber}()
        for f in readdir(path; join=true)
            st = splitext(basename(f))
            last(st) == ".csv" || continue
            push!(versions, VersionNumber(first(st)[9:end]))
        end
        sort!(versions)

        bytes = Dict{String,Vector{Int}}()
        for (i, version) in enumerate(versions)
            csv_path = joinpath(path, "mallocs-$(string(version)).csv")
            df = DataFrame(CSV.File(csv_path))
            v = get!(bytes, "Total", zeros(Int, length(versions)))
            v[i] = sum(df.bytes)
            for (b, f) in zip(df[:, 1], df[:, 3])
                w = get!(bytes, "$f", zeros(Int, length(versions)))
                w[i] += b
            end

            X, Y = pie_filter(df)
            pie(X, Y; title="Mallocs for $target.jl@v$(string(version))", l=0.5)
            for format in formats
                savefig(joinpath(path, "mallocs-$version.$format"))
            end
        end

        X = map(string, versions)
        Y = reshape(
            collect(Iterators.flatten(values(bytes))), length(versions), length(bytes)
        )
        L = reshape(collect(keys(bytes)), 1, length(bytes))
        plot(
            X,
            Y;
            xlabel="version",
            ylabel="bytes",
            markershape=:circle,
            # seriestype,
            title="Mallocs evolution in\n$target.jl",
            l=(0.5, 2),
            label=L,
            # yaxis=:log,
        )
        for format in formats
            savefig(joinpath(path, "mallocs-evolutions.$format"))
        end
    end
end


text:
using Test
using ExproniconLite
using ComoniconTypes
using ComoniconCast

args = [
    # name, type, require, vararg, default
    JLArgument(:arg1, Any,     true, false, nothing),
    JLArgument(:arg2, Int,     false, false, "::Int64"),
    JLArgument(:arg3, String,  false, false, "::String"),
]

flags = [
    JLFlag(:flag1),
    JLFlag(:flag2),
]

options = [
    # name, type, hint
    JLOption(:option1, Any, "::Nothing"),
    JLOption(:option2, Int, "::Int64"),
    JLOption(:option3, String, "::String"),
]

"""
    foo(arg1, arg2::Int, arg3::String; kw...)

a test function.

# Args

- `arg1`: test argument.
- `arg2`: test argument.
- `arg3`: test argument.

# Options

- `--option1, -o`: test option.
- `--option2 <int>`: test option.
- `--option3=<name>`: test option.

# Flags
- `--flag1, -f`: test flag.
- `--flag2`: test flag.
"""
function foo(arg1, arg2::Int, arg3::String;
        option1=nothing, option2::Int=1, option3::String="abc",
        flag1::Bool=false, flag2::Bool=false)
end

@testset "cast(::Function, ...)" begin
    cmd = cast(foo, "foo", args, options, flags)

    @test cmd.fn === foo
    @test cmd.name == "foo"
    @test cmd.args[1].name == "arg1"
    @test cmd.args[1].vararg == false
    @test cmd.args[1].require == true
    @test cmd.args[1].default === nothing
    @test cmd.args[1].description.brief == "test argument."
    @test cmd.vararg === nothing
    @test cmd.flags["flag1"].name == "flag1"
    @test cmd.flags["f"].name == "flag1"
    @test cmd.flags["f"].short == true
    @test cmd.options["option1"].name == "option1"
    @test cmd.options["option2"].name == "option2"
    @test cmd.options["o"].short == true
    @test cmd.options["o"].hint === nothing

    @test cmd.options["option2"].short == false
    @test cmd.options["option2"].hint == "int"
    @test cmd.options["option2"].type === Int
    @test cmd.description == Description("a test function.")
    @test cmd.line == LineNumberNode(0)
end

@testset "split_leaf_command" begin
    def = @expr JLFunction function foo(arg1, arg2::Int=1, arg3::String="abc";
        option1=nothing, option2::Int=1, option3::String="abc",
        flag1::Bool=false, flag2::Bool=false)
    end

    args′, options′, flags′ = split_leaf_command(def)
    args′, options′, flags′ = eval(args′), eval(options′), eval(flags′)

    @test args′ == args
    @test options == options′
    @test flags == flags′

    def = @expr JLFunction function foo(;option1::Bool=true)
    end
    @test_throws ErrorException split_leaf_command(def)

    def = @expr JLFunction function foo(;option1::Bool)
    end
    @test_throws ErrorException split_leaf_command(def)

    def = @expr JLFunction function foo(;option1)
    end
    @test_throws ErrorException split_leaf_command(def)
end

@test_throws ErrorException eval(:(
    module TestA
    using ComoniconCast
    @cast module Foo end
    end
))

module TestB
using ComoniconCast

@cast module Foo
    using ComoniconCast
    @cast foo(x) = 1
end
end

@testset "cast module" begin
    @test TestB.CASTED_COMMANDS["foo"].name == "foo"
    @test TestB.CASTED_COMMANDS["foo"] isa NodeCommand
    @test TestB.Foo.CASTED_COMMANDS["foo"].name == "foo"
    @test TestB.Foo.CASTED_COMMANDS["foo"] isa LeafCommand
end

module TestC
    using ComoniconCast
    @cast function command_a(a, b::String, c::Int; option_a="abc")
    end
end

@testset "replace _ => -" begin
    @test haskey(TestC.CASTED_COMMANDS, "command-a")
    cmd = TestC.CASTED_COMMANDS["command-a"]
    @test cmd isa LeafCommand
    option = cmd.options["option-a"]
    @test option.name == "option-a"
    @test option.hint == "::String"
end