text:
using KernelAbstractions.Extras.LoopInfo: @unroll

##### Kernels that ensure 'no-flux' from second- and fourth-order diffusion operators.
##### Kernel functions that ensure 'no-flux' from second- and fourth-order diffusion operators.
##### Note that flux divergence associated with a flux boundary condition is added
##### in a separate step.
#####
##### We implement two sets of kernel functions: one for filling one boundary at a time, and
##### a second that fills both boundaries at the same as a performance optimization.
#####

#####
##### Low-level functions that set data
#####

@inline _fill_flux_west_halo!(i, j, k, grid, c) = @inbounds c[1-i, j, k] = c[i, j, k]
@inline _fill_flux_east_halo!(i, j, k, grid, c) = @inbounds c[grid.Nx+i, j, k] = c[grid.Nx+1-i, j, k]

@inline _fill_flux_south_halo!(i, j, k, grid, c) = @inbounds c[i, 1-j, k] = c[i, j, k]
@inline _fill_flux_north_halo!(i, j, k, grid, c) = @inbounds c[i, grid.Ny+j, k] = c[i, grid.Ny+1-j, k]

@inline _fill_flux_bottom_halo!(i, j, k, grid, c) = @inbounds c[i, j, 1-k] = c[i, j, k]
@inline _fill_flux_top_halo!(i, j, k, grid, c)    = @inbounds c[i, j, grid.Nz+k] = c[i, j, grid.Nz+1-k]

#####
#####
##### Combined halo filling functions
#####

@inline _fill_west_halo!(j, k, grid, c, ::FBC, args...)   =   _fill_flux_west_halo!(1, j, k, grid, c)
@inline _fill_east_halo!(j, k, grid, c, ::FBC, args...)   =   _fill_flux_east_halo!(1, j, k, grid, c)
@inline _fill_south_halo!(i, k, grid, c, ::FBC, args...)  =  _fill_flux_south_halo!(i, 1, k, grid, c)
@inline _fill_north_halo!(i, k, grid, c, ::FBC, args...)  =  _fill_flux_north_halo!(i, 1, k, grid, c)
@inline _fill_bottom_halo!(i, j, grid, c, ::FBC, args...) = _fill_flux_bottom_halo!(i, j, 1, grid, c)
@inline _fill_top_halo!(i, j, grid, c, ::FBC, args...)    =    _fill_flux_top_halo!(i, j, 1, grid, c)

#####
##### Single halo filling kernels
#####

@kernel function fill_flux_west_halo!(c, grid)
    j, k = @index(Global, NTuple)
    _fill_flux_west_halo!(1, j, k, grid, c)
end

@kernel function fill_flux_south_halo!(c, grid)
    i, k = @index(Global, NTuple)
    _fill_flux_south_halo!(i, 1, k, grid, c)
end

@kernel function fill_flux_bottom_halo!(c, grid)
    i, j = @index(Global, NTuple)
    _fill_flux_bottom_halo!(i, j, 1, grid, c)
end

@kernel function fill_flux_east_halo!(c, grid)
    j, k = @index(Global, NTuple)
    _fill_flux_east_halo!(1, j, k, grid, c)
end

@kernel function fill_flux_north_halo!(c, grid)
    i, k = @index(Global, NTuple)
    _fill_flux_north_halo!(i, 1, k, grid, c)
end

@kernel function fill_flux_top_halo!(c, grid)
    i, j = @index(Global, NTuple)
    _fill_flux_top_halo!(i, j, 1, grid, c)
end

#####
##### Kernel launchers for flux boundary conditions
#####

fill_west_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :yz, fill_flux_west_halo!, c, grid; dependencies=dep, kwargs...)
fill_east_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :yz, fill_flux_east_halo!, c, grid; dependencies=dep, kwargs...)
fill_south_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :xz, fill_flux_south_halo!, c, grid; dependencies=dep, kwargs...)
fill_north_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :xz, fill_flux_north_halo!, c, grid; dependencies=dep, kwargs...)
fill_bottom_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :xy, fill_flux_bottom_halo!, c, grid; dependencies=dep, kwargs...)
fill_top_halo!(c, bc::FBC, arch, dep, grid, args...; kwargs...) =
            launch!(arch, grid, :xy, fill_flux_top_halo!, c, grid; dependencies=dep, kwargs...)


text:
using Documenter, WorkshopWizard

makedocs(
    modules = [WorkshopWizard],
    doctest = true,
    sitename = "WorkshopWizard.jl",
    pages = [
        "Introduction" => "index.md",
        "Usage" => "usage.md",
        "Troubleshooting" => "troubleshooting.md",
        "Functions" => "functions.md",
    ],
    # assets = ["assets/custom.css", "assets/custom.js"]
)

deploydocs(repo = "github.com/crstnbr/WorkshopWizard.jl.git", push_preview = true)

ext: jl


text:
"""
```
struct InpStiffness{dim, N, TF, TI, TBool, Tch <: ConstraintHandler, GO, TInds <: AbstractVector{TI}, TMeta<:Metadata} <: StiffnessTopOptProblem{dim, TF}
    inp_content::InpContent{dim, TF, N, TI}
    geom_order::Type{Val{GO}}
    ch::Tch
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end
```

- `dim`: dimension of the problem
- `TF`: number type for computations and coordinates
- `N`: number of nodes in a cell of the grid
- `inp_content`: an instance of [`InpContent`](@ref) which stores all the information from the ``.inp` file.
- `geom_order`: a field equal to `Val{GO}` where `GO` is an integer representing the order of the finite elements. Linear elements have a `geom_order` of `Val{1}` and quadratic elements have a `geom_order` of `Val{2}`.
- `metadata`: Metadata having various cell-node-dof relationships
- `black`: a `BitVector` of length equal to the number of elements where `black[e]` is 1 iff the `e`^th element must be part of the final design
- `white`:  a `BitVector` of length equal to the number of elements where `white[e]` is 1 iff the `e`^th element must not be part of the final design
- `varind`: an `AbstractVector{Int}` of length equal to the number of elements where `varind[e]` gives the index of the decision variable corresponding to element `e`. Because some elements can be fixed to be black or white, not every element has a decision variable associated.
"""
struct InpStiffness{
    dim,N,TF,TI,TBool,Tch<:ConstraintHandler,GO,TInds<:AbstractVector{TI},TMeta<:Metadata
} <: StiffnessTopOptProblem{dim,TF}
    inp_content::InpContent{dim,TF,N,TI}
    geom_order::Type{Val{GO}}
    ch::Tch
    black::TBool
    white::TBool
    varind::TInds
    metadata::TMeta
end

"""
    InpStiffness(filepath::AbstractString; keep_load_cells = false)

Imports stiffness problem from a .inp file, e.g. `InpStiffness("example.inp")`. The `keep_load_cells` keyword argument will enforce that any cell with a load applied on it is to be part of the final optimal design generated by topology optimization algorithms.
"""
function InpStiffness(filepath_with_ext::AbstractString; keep_load_cells=false)
    problem = Parser.extract_inp(filepath_with_ext)
    return InpStiffness(problem; keep_load_cells=keep_load_cells)
end
function InpStiffness(problem::Parser.InpContent; keep_load_cells=false)
    ch = Parser.inp_to_ferrite(problem)
    black, white = find_black_and_white(ch.dh)
    metadata = Metadata(ch.dh)
    geom_order = Ferrite.getorder(ch.dh.field_interpolations[1])
    if keep_load_cells
        for k in keys(problem.cloads)
            for (c, f) in metadata.node_cells[k]
                black[c] = 1
            end
        end
    end
    varind = find_varind(black, white)
    return InpStiffness(problem, Val{geom_order}, ch, black, white, varind, metadata)
end

getE(p::InpStiffness) = p.inp_content.E
getν(p::InpStiffness) = p.inp_content.ν
nnodespercell(::InpStiffness{dim,N}) where {dim,N} = N
getgeomorder(p::InpStiffness{<:Any,<:Any,<:Any,<:Any,<:Any,<:Any,GO}) where {GO} = GO
getdensity(p::InpStiffness) = p.inp_content.density
getpressuredict(p::InpStiffness) = p.inp_content.dloads
getcloaddict(p::InpStiffness) = p.inp_content.cloads
getfacesets(p::InpStiffness) = p.inp_content.facesets


text:
module Extension
using MLStyle.Err
export use, @use, used

@nospecialize
macro _depwarn(sym)
    :($Base.depwarn("No need to use this function anymore since MLStyle v0.4", $sym))
end

function used(ext::Symbol, mod::Module)::Bool
    @_depwarn :used
    false
end

function use(ext::Symbol, mod::Module)
    @_depwarn :use
    nothing
end

macro use(exts...)
    ln = __source__
    @warn "No need to use this function anymore since MLStyle v0.4 at $(ln.file):$(ln.line)"
    use(:_, @__MODULE__)
end
@specialize

end


text:
function countdict(xs)
    counts = Dict()
    for x in xs
        if haskey(counts, x)
            counts[x] += 1
        else
            counts[x] = 1
        end
    end
    return counts
end

get_words(sentence) = map(lowercase, split(sentence, " "))
build_ngram(words, n) = map(i -> join(words[i:i+n-1], " "), 1:length(words)-n+1)
ngram_counts(sentence, n) = countdict(build_ngram(get_words(sentence), n))


text:
# # [A classification task: the prediction of  plant species from floreal measures (the iris dataset)](@id clustering_tutorial)
# The task is to estimate the species of a plant given some floreal measurements. It use the classical "Iris" dataset.
# Note that in this example we are using clustering approaches, so we try to understand the "structure" of our data, without relying to actually knowing the true labels ("classes" or "factors"). However we have chosen a dataset for which the true labels are actually known, so to compare the accuracy of the algorithms we use, but these labels will not be used during the algorithms training.

#
# Data origin:
# - dataset description: [https://en.wikipedia.org/wiki/Iris_flower_data_set](https://en.wikipedia.org/wiki/Iris_flower_data_set)
# - data source we use here: [https://github.com/JuliaStats/RDatasets.jl](https://github.com/JuliaStats/RDatasets.jl)

# ## Library and data loading
# We load the Beta Machine Learning Toolkit as well as some other packages that we use in this tutorial
using BetaML
using Random, Statistics, Logging, BenchmarkTools, RDatasets, Plots, DataFrames

# We are also going to compare our results with two other leading packages in Julia for clustering analysis, [`Clustering.jl`](https://github.com/JuliaStats/Clustering.jl) that provides (inter alia) kmeans and kmedoids algorithms and [`GaussianMixtures.jl`](https://github.com/davidavdav/GaussianMixtures.jl) that provides, as the name says, Gaussian Mixture Models. So we import them (we "import" them, rather than "use", not to bound their full names into namespace as some would collide with BetaML).
import Clustering, GaussianMixtures
using  Test     #src

# We do a few tweeks for the Clustering and GaussianMixtures packages. Note that in BetaML we can also control both the random seed and the verbosity in the algorithm call, not only globally
Random.seed!(123)
#logger  = Logging.SimpleLogger(stdout, Logging.Error); global_logger(logger); ## For suppressing GaussianMixtures output

# Differently from the [regression tutorial](@ref regression_tutorial), we load the data here from [`RDatasets`](https://github.com/JuliaStats/RDatasets.jl](https://github.com/JuliaStats/RDatasets.jl), a package providing standard datasets.
iris = dataset("datasets", "iris")
describe(iris)

# The iris dataset  provides floreal measures in columns 1 to 4 and the assigned species name in column 5. There are no missing values

# ## Data preparation
# The first step is to prepare the data for the analysis. We collect the first 4 columns as our _feature_ `x` matrix and the last one as our `y` label vector.
# As we are using clustering algorithms, we are not actually using the labels to train the algorithms, we'll behave like we do not know them, we'll just let the algorithm "learn" fro mthe structure of the data itself. We'll however use it to judge the accuracy that they did reach.
x       = Matrix{Float64}(iris[:,1:4]);
yLabels = unique(iris[:,5]);
# As the labels are expressed as strings, the first thing we do is encode them as integers for our analysis using the function [`integerEncoder`](@ref).
y       = integerEncoder(iris[:,5],factors=yLabels);

# The dataset from RDatasets is ordered by species, so we need to shuffle it to avoid biases.
# Shuffling happens by default in crossValidation, but we are keeping here a copy of the shuffled version for later.
# Note that the version of [`shuffle`](@ref) that is included in BetaML accepts several n-dimensional arrays and shuffle them (by default on rows, by we can specify the dimension) keeping the association  between the various arrays in the shuffled output.
(xs,ys) = shuffle([x,y]);

# ## Main analysis

# We will try 3 BetaML models ([`kmeans`](@ref), [`kmedoids`](@ref) and [`gmm`](@ref)) and we compare them with `kmeans` from Clusterings.jl and `GMM` from GaussianMixtures.jl
# `Kmeans` and `kmedoids` works by first initialising the centers of the k-clusters (the "representative" (step a ) . For `kmeans` they must be selected within one of the data, for kmeans they are the geometrical center) n a nutshell. Then ( b ) iterate for each point to assign the point to the cluster of the closest representative (according with a user defined distance metric, default to Euclidean), and ( c ) move each representative at the center of its newly acquired cluster (where "center" depends again from the metric). Steps ( b ) and ( c ) are reiterated until the algorithm converge, i.e. the tentative k representative points (and their relative clusters) don't move any more. The result (output of the algorithm) is that each point is assigned to one of the clusters (classes).
# The `gmm` algorithm is similar in that it employs an iterative approach (the Expectation_Minimisation algorithm, "em") but here we make the hipothesis that the data points are the observed outcomes of some _mixture_ probabilistic models where we have first a k-categorical variables whose outcomes are the (unobservble) parameters of a probabilistic distribution from which the data is finally drawn. Because the parameters of each of the k-possible distributions is unobservable this is also called a model with latent variables.
# Most `gmm` models use the Gaussain distribution as the family of the mixture components, so we can tought the `gmm` acronym to indicate _Gaussian Mixture Model_. In BetaML we do implemented only Gaussain components, but any distribution could be used by just subclassing `AbstractMixture` and implementing a couple of methids (you are invited to contribute or just ask for a distribution family you are interested), so I prefer to think "gmm" as an acronym for _Generative Mixture Model_.
# The algorithm try to find the mixture that maximises the likelihood that the data has been generated indeed from such mixture, where the "E" step refers to computing the probability that each point belongs to each of the k-composants (somehow similar to the step _b_ in the kmeans/kmedoids algorithm), and the "M" step estimates, giving the association probabilities in step "M", the parameters of the mixture and of the individual components (similar to step _c_).
# The result here is that each point has a categorical distribution (PMF) representing the probabilities that it belongs to any of the k-components (our classes or clusters). This is interesting, as `gmm` can be used for many other things that clustering. It forms the backbone of the [`predictMissing`](@ref) function to impute missing values (on some or all dimensions) based to how close the record seems to its pears. For the same reasons, `predictMissing` can also be used to predict user's behaviours (or users' appreciation) according to the behaviour/ranking made by pears ("collaborative filtering").
# While the result of `gmm` is a vector of PMFs (one for each record), error measures and reports with the true values (if known) can be directly applied, as in BetaML they internally call `mode()` to retrieve the class with the highest probability for each record.

# As we are here, we also try different versions of the BetaML models, even if the default "versions" should be fine. For `kmeans` and `kmedoids` we will try different initialisation strategies ("gird", the default one, "random" and "shuffle"), while for the `gmm` model we'll choose different distributions of the Gaussain family (`SphericalGaussian` - where the variance is a scalar, `DiagonalGaussian` - with a vector variance, and `FullGaussian`, where the covariance is a matrix).

# As the result would depend on stochasticity both in the data selected and in the random initialisation, we use a cross-validation approach to run our models several times (with different data) and then we average their results.
# Cross-Validation in BetaML is very flexible and it is done using the [`crossValidation`](@ref) function.
# crossValidation works by calling the function `f`, defined by the user, passing to it the tuple `trainData`, `valData` and `rng` and collecting the result of the function f. The specific method for which `trainData`, and `valData` are selected at each iteration depends on the specific `sampler`.
# We start by selectign a k-fold sampler that split our data in 5 different parts, it uses 4 for training and 1 part (not used here) for validation. We run the simulations twice and, to be sure to have replicable results, we fix the random seed (at the whole crossValidaiton level, not on each iteration).
sampler = KFold(nSplits=5,nRepeats=3,shuffle=true, rng=copy(FIXEDRNG))

# We can now run the cross-validation with our models. Note that instead of defining the function `f` and then calling `crossValidation[f(trainData,testData,rng),[x,y],...)` we use the Julia `do` block syntax and we write directly the content of the `f` function in the `do` block.
# Also, by default crossValidation already returns the mean and the standard deviation of the output of the user-provided `f` function (or the `do` block). However this requires that the `f` function return a single scalar. Here we are returning a vector of the accuracies of the different models (so we can run the cross-validation only once), and hence we indicate with `returnStatistics=false` to crossValidation not to attempt to generate statistics but rather report the whole output.
# We'll compute the statistics ex-post.

# Inside the `do` block we do 4 things:
# - we recover from `trainData` (a tuple, as we passed a tuple to `crossValidation` too) the `xtrain` features and `ytrain` labels;
# - we run the various clustering algorithms
# - we use the real labels to compute the model accuracy. Note that the clustering algorithm know nothing about the specific label name or even their order. This is why [`accuracy`](@ref) has the parameter `ignoreLabels` to compute the accuracy oven any possible permutation of the classes found.
# - we return the various models' accuracies
cOut = crossValidation([x,y],sampler,returnStatistics=false) do trainData,testData,rng
          ## For unsupervised learning we use only the train data.
          ## Also, we use the associated labels only to measure the performances
         (xtrain,ytrain)  = trainData;
         ## We run the clustering algorithm...
         clusteringOut     = kmeans(xtrain,3,rng=rng) ## init is grid by default
         ## ... and we compute the accuracy using the real labels
         kMeansAccuracy    = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = kmeans(xtrain,3,rng=rng,initStrategy="random")
         kMeansRAccuracy   = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = kmeans(xtrain,3,rng=rng,initStrategy="shuffle")
         kMeansSAccuracy   = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = kmedoids(xtrain,3,rng=rng)   ## init is grid by default
         kMedoidsAccuracy  = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = kmedoids(xtrain,3,rng=rng,initStrategy="random")
         kMedoidsRAccuracy = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = kmedoids(xtrain,3,rng=rng,initStrategy="shuffle")
         kMedoidsSAccuracy = accuracy(clusteringOut[1],ytrain,ignoreLabels=true)
         clusteringOut     = gmm(xtrain,3,mixtures=[SphericalGaussian() for i in 1:3], verbosity=NONE, rng=rng)
         gmmSpherAccuracy  = accuracy(clusteringOut.pₙₖ,ytrain,ignoreLabels=true, rng=rng)
         clusteringOut     = gmm(xtrain,3,mixtures=[DiagonalGaussian() for i in 1:3], verbosity=NONE, rng=rng)
         gmmDiagAccuracy   = accuracy(clusteringOut.pₙₖ,ytrain,ignoreLabels=true, rng=rng)
         clusteringOut     = gmm(xtrain,3,mixtures=[FullGaussian() for i in 1:3], verbosity=NONE, rng=rng)
         gmmFullAccuracy   = accuracy(clusteringOut.pₙₖ,ytrain,ignoreLabels=true, rng=rng)
         ## For comparision with Clustering.jl
         clusteringOut     = Clustering.kmeans(xtrain', 3)
         kMeans2Accuracy   = accuracy(clusteringOut.assignments,ytrain,ignoreLabels=true)
         ## For comparision with GaussianMistures.jl - sometimes GaussianMistures.jl em! fails with a PosDefException
         dGMM              = GaussianMixtures.GMM(3, xtrain; method=:kmeans, kind=:diag)
         GaussianMixtures.em!(dGMM, xtrain)
         gmmDiag2Accuracy  = accuracy(GaussianMixtures.gmmposterior(dGMM, xtrain)[1],ytrain,ignoreLabels=true)
         fGMM              = GaussianMixtures.GMM(3, xtrain; method=:kmeans, kind=:full)
         GaussianMixtures.em!(fGMM, xtrain)
         gmmFull2Accuracy  = accuracy(GaussianMixtures.gmmposterior(fGMM, xtrain)[1],ytrain,ignoreLabels=true)
         ## Returning the accuracies
         return kMeansAccuracy,kMeansRAccuracy,kMeansSAccuracy,kMedoidsAccuracy,kMedoidsRAccuracy,kMedoidsSAccuracy,gmmSpherAccuracy,gmmDiagAccuracy,gmmFullAccuracy,kMeans2Accuracy,gmmDiag2Accuracy,gmmFull2Accuracy
 end

## We transform the output in matrix for easier analysis
accuracies = fill(0.0,(length(cOut),length(cOut[1])))
[accuracies[r,c] = cOut[r][c] for r in 1:length(cOut),c in 1:length(cOut[1])]
μs = mean(accuracies,dims=1)
σs = std(accuracies,dims=1)

@test all(μs .> 0.7) #src

@test μs[1] > 0.89 &&  μs[4] > 0.89 &&  μs[9] > 0.96 #src
modelLabels=["kMeansG","kMeansR","kMeansS","kMedoidsG","kMedoidsR","kMedoidsS","gmmSpher","gmmDiag","gmmFull","kMeans (Clustering.jl)","gmmDiag (GaussianMixtures.jl)","gmmFull (GaussianMixtures.jl)"]
report = DataFrame(mName = modelLabels, avgAccuracy = dropdims(round.(μs',digits=3),dims=2), stdAccuracy = dropdims(round.(σs',digits=3),dims=2))

#src plot(modelLabels,μs',seriestype=:scatter)
#src yerror=collect(zip(rand(12), rand(12))

# ### BetaML model accuracies

# From the output We see that the gmm models perform for this dataset generally better than kmeans or kmedoids algorithms, also with very low variances.
# In detail, it is the (default) `grid` initialisation that leads to the better results for `kmeans` and `kmedoids`, while for the `gmm` models it is the `FullGaussian` to perform better.

# ### Comparisions with `Clustering.jl` and `GaussianMixtures.jl`
# For this specific case, both `Clustering.jl` and `GaussianMixtures.jl` report substantially worst accuracies, and with very high variances. But we maintain the ranking that Full Gaussian gmm > Diagonal Gaussian > Kmeans accuracy.
# I suspect the reason that BetaML gmm works so weel is in relation to the usage of kmeans algorithm with itself the grid initialisation.
# The grid initialisation "guarantee" indeed that the initial means of the mixture components are well spread across the multidimensional space defined by the data, and it helps avoiding the EM algoritm to converge to a bad local optimus.

# ## Working without the labels

# Up to now we used the real labels to compare the model accuracies. But in real clustering examples we don't have the true classes, or we wouln't need to do clustering in the first instance, so we don't know the number of classes to use.
# There are several methods to judge clusters algorithms goodness, perhaps the simplest one, at least for the expectation-maximisation algorithm employed in `gmm` to fit the data to the unknown mixture, is to use a information criteria that trade the goodness of the lickelyhood with the parameters used to do the fit.
# BetaML provide by default in the gmm clustering outputs both the _Bayesian information criterion_  ([`BIC`](@ref bic)) and the _Akaike information criterion_  ([`AIC`](@ref aic)), where for both a lower value is better.

# We can then run the model with different number of classes and see which one leads to the lower BIC or AIC.
# We run hence `crossValidation` again with the `FullGaussian` gmm model
# Note that we use the BIC/AIC criteria here for establishing the "best" number of classes but we could have used it also to select the kind of Gaussain distribution to use. This is one example of hyper-parameter tuning that we developed more in detail (but without using cross-validation) in the [regression tutorial](@ref regression_tutorial).

# Let's try up to 8 possible classes:

K = 8
sampler = KFold(nSplits=5,nRepeats=2,shuffle=true, rng=copy(FIXEDRNG))
cOut = crossValidation([x,y],sampler,returnStatistics=false) do trainData,testData,rng
    (xtrain,ytrain)  = trainData;
    clusteringOut  = [gmm(xtrain,k,mixtures=[FullGaussian() for i in 1:k], verbosity=NONE, rng=rng) for k in 1:K]
    BICS           = [clusteringOut[i].BIC for i in 1:K]
    AICS           = [clusteringOut[i].AIC for i in 1:K]
    return (BICS,AICS)
end

## Transforming the output in matrices for easier analysis
Nit = length(cOut)

BICS = fill(0.0,(Nit,K))
AICS = fill(0.0,(Nit,K))
[BICS[r,c] = cOut[r][1][c] for r in 1:Nit,c in 1:K]
[AICS[r,c] = cOut[r][2][c] for r in 1:Nit,c in 1:K]

μsBICS = mean(BICS,dims=1)
#-
σsBICS = std(BICS,dims=1)
#-
μsAICS = mean(AICS,dims=1)
#-
σsAICS = std(AICS,dims=1)
#-
plot(1:K,[μsBICS' μsAICS'], labels=["BIC" "AIC"], title="Information criteria by number of classes", xlabel="number of classes", ylabel="lower is better")

# We see that following the "lowest AIC" rule we would indeed choose three classes, while following the "best AIC" criteria we would have choosen only two classes. This means that there is two classes that, concerning the floreal measures used in the database, are very similar, and opur models are unsure about them. Perhaps the biologists will end up one day with the conclusion that it is indeed only one specie :-).

# We could study this issue more in detail by analysing the [`ConfusionMatrix`](@ref), but the one used in BetaML does not account for the ignoreLabels option (yet).

# ## Benchmarking computational efficiency

# We now benchmark the time and memory required by the various models by using the `@btime` macro of the `BenchmarkTools` package:

# ```
# @btime kmeans($xs,3);
# # 261.540 μs (3777 allocations: 442.53 KiB)
# @btime kmedoids($xs,3);
# 4.576 ms (97356 allocations: 10.42 MiB)
# @btime gmm($xs,3,mixtures=[SphericalGaussian() for i in 1:3], verbosity=NONE);
# # 5.498 ms (133365 allocations: 8.42 MiB)
# @btime gmm($xs,3,mixtures=[DiagonalGaussian() for i in 1:3], verbosity=NONE);
# # 18.901 ms (404333 allocations: 25.65 MiB)
# @btime gmm($xs,3,mixtures=[FullGaussian() for i in 1:3], verbosity=NONE);
# # 49.257 ms (351500 allocations: 61.95 MiB)
# @btime Clustering.kmeans($xs', 3);
# # 17.071 μs (23 allocations: 14.31 KiB)
# @btime begin dGMM = GaussianMixtures.GMM(3, $xs; method=:kmeans, kind=:diag); GaussianMixtures.em!(dGMM, $xs) end;
# # 530.528 μs (2088 allocations: 488.05 KiB)
# @btime begin fGMM = GaussianMixtures.GMM(3, $xs; method=:kmeans, kind=:full); GaussianMixtures.em!(fGMM, $xs) end;
# # 4.166 ms (58910 allocations: 3.59 MiB)
# ```
# (_note: the values reported here are of a local pc, not of the GitHub CI server, as sometimes - depending on data and random initialisation - `GaussainMixtures.em!`` fails with a `PosDefException`. This in turln would lead the whole documentation to fail to compile_)

# Like for supervised models, dedicated models are much better optimized than BetaML models, and are order of magnitude more efficient. However even the slowest BetaML clusering model (gmm using full gaussians) is realtively fast and can handle mid-size datasets (tens to hundreds of thousand records) without significant slow downs.

# ## Conclusions

# We have shown in this tutorial how we can easily run clustering almgorithms in BetaML with just one line of code `choosenModel(x,k)`, but also how can we use cross-validation in order to help the model or parameter selection, with or whithout knowing the real classes.
# We retrieve here what we observed with supervised models. Globally the accuracy of BetaML models are comparable to those of leading specialised packages (in this case they are even better), but there is a significant gap in computational efficiency that restricts the pratical usage of BetaML to mid-size datasets. However we trade this relative inefficiency with very flexible model definition and utility functions (for example the BetaML gmm works with missing data, allowing it to be used as the backbone of the [`predictMissing`](@ref) missing imputation function, or for collaborative reccomendation systems).

max_issues_repo_name: sylvaticus/IntroSPMLJuliaCourse

max_forks_repo_name: sylvaticus/IntroSPMLJuliaCourse


text:
"joincmds(cmds): stitch together an array of commands"
joincmds(cmds) = foldl((cmd1,cmd2)->`$cmd1 $cmd2`, ``, cmds)

# Filesystem stuff
"mkdir_ifnot(dir): Create directoy if it does not exist"
mkdir_ifnot(dir) = if (!isdir(dir)) mkdir(dir) end

"Decimal numbers from danish software often use , rather than . as separator"
parse_danish_float(x) = parse(Float64, replace(x, ",", "."))

"Count the number of lines in a file"
count_lines(file) = length(split(readstring(file),'\n'))

"Filters array to include only variables that have been defined (vectorized operation)"
vals_if_defined(x) = map(x->@eval($x), filter(isdefined,x))


text:
module PointerBasedTrees
using ClusterTrees

import ClusterTrees.ChildIterator
import ClusterTrees: start, next, done

mutable struct Node{T}
    data::T
    num_children::Int
    next_sibling::Int
    parent::Int
    first_child::Int
end

data(n::Node) = n.data

abstract type APBTree end
struct PointerBasedTree{N<:Node} <: APBTree
    nodes::Vector{N}
    root::Int
end

getnode(tree::PointerBasedTree, node_idx) = tree.nodes[node_idx]
nextsibling(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).next_sibling
ClusterTrees.parent(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).parent
firstchild(tree::PointerBasedTree, node_idx) = getnode(tree, node_idx).first_child

start(itr::ChildIterator{<:APBTree}) = (0, firstchild(itr.tree, itr.node))
function done(itr::ChildIterator{<:APBTree}, state)
    prev, this = state
    this < 1 && return true
    sibling_par = ClusterTrees.parent(itr.tree, this)
    sibling_par != itr.node && return true
    return false
end

function next(itr::ChildIterator{<:APBTree}, state)
    prev, this = state
    nxt = nextsibling(itr.tree, this)
    return (this, (this, nxt))
end

Base.iterate(itr::ChildIterator{<:APBTree}, st = start(itr)) = done(itr,st) ? nothing : next(itr,st)

# Base.IteratorSize(cv::ChildIterator) = Base.SizeUnknown()

ClusterTrees.root(tree::PointerBasedTree) = tree.root
ClusterTrees.children(tree::APBTree, node=ClusterTrees.root(tree)) = ChildIterator(tree, node)
ClusterTrees.haschildren(tree::APBTree, node) = (firstchild(tree,node) >= 1)
ClusterTrees.data(tree::PointerBasedTree, node=ClusterTrees.root(tree)) = data(getnode(tree, node))

# """
#     insert!(tree, parent, data)
#
# Insert a node carrying `data` as the first child of `parent`
# """
# function ClusterTrees.insert!(tree::PointerBasedTree, parent_idx, data)
#     next = firstchild(tree, parent_idx)
#     ClusterTrees.insert!(tree, data, parent=parent_idx, next=next, prev=0)
# end

setfirstchild!(node::Node, child) = Node(node.data, node.num_children, node.next_sibling, node.parent, child)
setfirstchild!(tree::PointerBasedTree, node, child) = tree.nodes[node] = setfirstchild!(getnode(tree, node), child)

setnextsibling!(node::Node, next) = Node(node.data, node.num_children, next, node.parent, node.first_child)
setnextsibling!(tree::PointerBasedTree, node, next) = tree.nodes[node] = setnextsibling!(getnode(tree, node), next)

"""
    insert!(child_iterator, item, position)

Insert item at 'position' in 'child_iterator'. Here, 'position' is a state generated
using 'start/next/done' acting on 'child_iterator', which in turn is created by
'children(tree, node)'.
"""
function Base.insert!(chd_itr::ChildIterator{<:APBTree}, item, state)
    prev, next = state
    parent = chd_itr.node

    tree = chd_itr.tree
    push!(tree.nodes, Node(item, 0, next, parent, 0))
    this = lastindex(tree.nodes)
    if prev < 1
        setfirstchild!(tree, parent, this)
    else
        setnextsibling!(tree, prev, this)
    end
    return this
end

end


text:
"""
predict(𝒱::FluxData, model)

#Description
Returns a tuple of (1) the model predictions for the variable associated with the object 𝒱
and (2) the truth data for the same variable.

#Arguments
- 𝒱: (FluxData) object containing the training data for the associates variable
- model: model returned by gp_model or nn_model function
"""
function predict(𝒱, model; subsampled_only=false, scaled=false)

    if scaled
        if subsampled_only
            predictions = (model(𝒱.training_data[i][1]) for i in 1:length(𝒱.training_data))
            targets = (𝒱.training_data[i][2] for i in 1:length(𝒱.training_data))
            return (cat(predictions...,dims=2), cat(targets...,dims=2))
        end

        predictions = (model(𝒱.uvT_scaled[:,t]) for t in 1:size(𝒱.uvT_scaled,2))
        return (cat(predictions...,dims=2),  𝒱.scaled)

    else
    if subsampled_only
        predictions = (𝒱.unscale_fn(model(𝒱.training_data[i][1])) for i in 1:length(𝒱.training_data))
        targets = (𝒱.unscale_fn(𝒱.training_data[i][2]) for i in 1:length(𝒱.training_data))
        return (cat(predictions...,dims=2), cat(targets...,dims=2))
    end

    predictions = (𝒱.unscale_fn(model(𝒱.uvT_scaled[:,t])) for t in 1:size(𝒱.uvT_scaled,2))
    return (cat(predictions...,dims=2),  𝒱.coarse)
end
end


text:
using ParticleAccelerations
using Test, SafeTestsets

@testset "ParticleAccelerations.jl" begin
    @safetestset "Molly integration tests" begin
        include("molly.jl")
    end
end


text:
##### Beginning of file

module Utils # Begin submodule MirrorUpdater.Utils

__precompile__(true)

include("git.jl")

include("misc.jl")

end # End submodule MirrorUpdater.Utils

##### End of file


text:
struct Empty end

# TODO: Rename to `Moving`?
"""
    Moved{Value}

A "tag" used for making that this slot is being moved to a new `slots` vector.
"""
struct Moved{Value}
    value::Value
end

const MOVED_EMPTY = Moved{Empty}(Empty())

# TODO: Support "UInt62"? (i.e., `UInt64` but two bits used as tag)
# TODO: Get rid of Moved{Value}
const ValueUnion{Value} = Union{Value,Empty,Moved{Value},Moved{Empty}}

# Since Julia does not have atomics for arrays, `LinearProbingDict` uses an
# array filled with mutable struct `DictSlot` instead. This has the disadvantage
# that it creates one more indirection. But a somewhat nice effect is that it
# makes copying a key-value pair cheap and easy to do with non-atomic write.
# Re-using small mutable objects like this can be bad for memory locality but it
# seems that the effect is neglegible. The performance is much better than
# creating small objects and copy the key and value separately.
#
# Due to the type tag, storing two UInt64 takes space as large as storing three
# UInt64. So, we can store hash for "free".
# TODO: Check if it is actually better to store the hash in an `Vector{UInt}`.
# Even though it will waste some memory, it will reduce indirection and should
# make probing much more efficient.
mutable struct DictSlot{
    Key,
    Value,
    KeyStorage>:Union{Key,Empty},
    ValueStorage>:ValueUnion{Value},
}
    @atomic hash::UInt
    @atomic key::KeyStorage
    @atomic value::ValueStorage

    # Defining constructor manually to workaround:
    # https://github.com/JuliaLang/julia/issues/42353
    DictSlot{Key,Value,KeyStorage,ValueStorage}(
        hash,
        key,
        value,
    ) where {Key,Value,KeyStorage,ValueStorage} =
        new{Key,Value,KeyStorage,ValueStorage}(hash, key, value)
end

DictSlot{Key,Value,KeyStorage,ValueStorage}() where {Key,Value,KeyStorage,ValueStorage} =
    DictSlot{Key,Value,KeyStorage,ValueStorage}(zero(UInt), Empty(), Empty())

mutable struct LinearProbingDict{Key,Value,KeyStorage,ValueStorage} <:
               ConcurrentDict{Key,Value}
    @atomic slots::Union{Vector{DictSlot{Key,Value,KeyStorage,ValueStorage}},Nothing}
    migration::ReentrantLock
    # TODO: per-thread non-atomic counter for approximating deleted elements
    nadded::Threads.Atomic{Int}
    ndeleted::Threads.Atomic{Int}
end
# TODO: inline the counters (and pad them)?
# TODO: relaxed counter

slots_type(
    ::LinearProbingDict{Key,Value,KeyStorage,ValueStorage},
) where {Key,Value,KeyStorage,ValueStorage} =
    Vector{DictSlot{Key,Value,KeyStorage,ValueStorage}}

ConcurrentCollections.length_upper_bound(dict::LinearProbingDict) =
    dict.nadded[] - dict.ndeleted[]

ConcurrentCollections.length_lower_bound(dict::LinearProbingDict) =
    dict.nadded[] - dict.ndeleted[]

ConcurrentCollections.ConcurrentDict{Key,Value}() where {Key,Value} =
    LinearProbingDict{Key,Value}()

# Extracted as a function so that `T` is captured as a type parameter in the
# closure for the generator.
fillwith(::Type{T}, n) where {T} = [T() for _ in 1:n]

function LinearProbingDict{Key,Value}() where {Key,Value}
    capacity = 4  # TODO: customizable init size?
    if sizeof(Some{Union{Key,Empty}}) > sizeof(UInt)
        KeyStorage = Any
    else
        KeyStorage = Union{Key,Empty}
    end
    if sizeof(Some{ValueUnion{Value}}) > sizeof(UInt)
        ValueStorage = Any
    else
        ValueStorage = ValueUnion{Value}
    end
    slots = fillwith(DictSlot{Key,Value,KeyStorage,ValueStorage}, capacity)
    return LinearProbingDict{Key,Value,KeyStorage,ValueStorage}(
        slots,
        ReentrantLock(),
        Threads.Atomic{Int}(0),
        Threads.Atomic{Int}(0),
    )
end

struct ValueRef{Value}
    value::Value
end
# Note on `modify!` design: It looks like (even relaxed) atomic load is not
# eliminated when the value is not used (<https://godbolt.org/z/abE7sGfjG>).
# So, let's pass a `Ref`-like object to `modify!` and so that load is not issued
# when the user does not look at the value.

@inline Base.getindex(ref::ValueRef) = ref.value

function Base.getindex(d::LinearProbingDict{Key}, key) where {Key}
    y = maybeget(d, key)
    if y === nothing
        throw(KeyError(key))
    else
        return something(y)
    end
end

function Base.haskey(d::LinearProbingDict, key)
    @inline f(::Nothing) = nothing
    @inline f(_) = Keep(nothing)
    y = modify!(f, d, key)
    return y !== nothing
end

Base.get(d::LinearProbingDict, key, default) =
    something(ConcurrentCollections.maybeget(d, key), default)

function Base.get!(d::LinearProbingDict{<:Any,V}, key, default) where {V}
    @inline f(::Nothing) = Some{V}(default)
    @inline f(x) = Keep(x[])
    y = modify!(f, d, key)
    if y isa Some
        return something(y)::V
    else
        return y.value::V
    end
end

function Base.get!(
    f::F,
    d::LinearProbingDict{<:Any,V},
    key,
) where {F<:Union{Function,Type},V}
    value = Ref{V}()
    isset = Ref(false)
    @inline function onvalue(::Nothing)
        if isset[]
            Some{V}(value[])
        else
            local y = f()
            local v = convert(V, y)
            value[] = v
            isset[] = true
            Some{V}(v)
        end
    end
    @inline onvalue(x) = Keep(x[])
    y = modify!(onvalue, d, key)
    if y isa Some
        return something(y)::V
    else
        return y.value::V
    end
end

function ConcurrentCollections.maybeget(d::LinearProbingDict{<:Any,V}, key) where {V}
    @inline f(::Nothing) = nothing
    @inline f(x) = Keep(x[])
    y = modify!(f, d, key)
    if y === nothing
        return nothing
    else
        return Some{V}(y.value::V)
    end
end

function Base.setindex!(d::LinearProbingDict{Key,Value}, v, k) where {Key,Value}
    v′ = convert(Value, v)
    @inline f(_) = Some(v′)
    modify!(f, d, k)
    return d
end

Base.pop!(d::LinearProbingDict, key, default) = something(maybepop!(d, key), default)
function Base.pop!(d::LinearProbingDict, key)
    value = maybepop!(d, key)
    value === nothing && throw(KeyError(key))
    return something(value)
end

function ConcurrentCollections.maybepop!(d::LinearProbingDict, key)
    @inline f(::Nothing) = nothing
    @inline f(ref) = Delete(ref[])
    y = modify!(f, d, key)
    if y === nothing
        return nothing
    else
        return Some(y.value)
    end
end

function record_added!(dict)
    Threads.atomic_add!(dict.nadded, 1) + 1
end

function record_deleted!(dict)
    Threads.atomic_add!(dict.ndeleted, 1) + 1
end

# The main linearization points are the change in the value, not the key.  The
# CAS on key only determines the internal "physical" location of the slot.
function ConcurrentCollections.modify!(
    f,
    dict::LinearProbingDict{Key,Value},
    key,
) where {Key,Value}
    key = convert(Key, key)
    slots = load_slots(dict)
    h = hash(key)

    if 1.5 * length_upper_bound(dict) > length(slots)
        slots = expand!(dict, slots)
    end

    while true
        @label restart
        c = length(slots)
        offset = reinterpret(Int, h) & (c - 1)  # h % c
        nprobes = 0
        while true
            index = offset + 1
            offset = (offset + 1) & (c - 1)  # (offset + 1) % c

            ref = slots[index]
            hstored = @atomic :monotonic ref.hash
            if (!iszero(hstored)) & (h != hstored)
                nprobes += 1
                if nprobes > c ÷ 4
                    @goto too_many_probes
                end
                continue  # already occupied; continue probing
            end
            # It's possible to reach here even when the key matches (i.e., the
            # maching key has been stored but the hash has not). This is fine
            # since the key is compared directly with `isequal` below anyway.

            kstored = (@atomic ref.key)::Union{Empty,Key}
            if kstored isa Empty
                reply = f(nothing)::Union{Nothing,Some}
                reply === nothing && return nothing
                oldkey, ok = @atomicreplace ref.key Empty() => key
                if ok
                    @atomic :monotonic ref.hash = h
                else
                    hstored = @atomic :monotonic ref.hash
                    if (!iszero(hstored) & (h != hstored)) || !isequal(oldkey, key)
                        continue  # failed to obtain this slot; continue probing
                    end
                    # Reaching here means two tasks simultaneously inserted the
                    # same key.  The value must be inserted to this slot.
                end
                v = convert(Value, something(reply))
                vstored, ok = @atomicreplace ref.value Empty() => v
                if ok
                    record_added!(dict)
                    return reply
                end
                vstored = vstored::Union{Value,Moved{Value},Moved{Empty}}
            elseif isequal(kstored, key)
                vstored = (@atomic ref.value)::ValueUnion{Value}
            else
                continue  # already occupied; continue probing
            end
            # Here, we know that this is the correct slot but the value is not
            # stored yet.
            while true
                if vstored isa Empty
                    reply = f(nothing)::Union{Nothing,Some}
                    reply === nothing && return nothing
                    v = convert(Value, something(reply))
                    vstored, ok = @atomicreplace ref.value Empty() => v
                    if ok
                        record_added!(dict)
                        return reply
                    end
                    vstored = vstored::Union{Value,Moved{Value},Moved{Empty}}
                end
                if vstored isa Moved
                    slots = finishmove!(dict, slots)
                    @goto restart
                end
                vref = ValueRef{Value}(vstored)
                reply = f(vref)::Union{Keep,Nothing,Delete,Some}
                if reply isa Keep
                    # @show nprobes
                    return reply
                end
                if reply isa Union{Nothing,Delete}
                    # Deletion:
                    vstored, ok = @atomicreplace ref.value vstored => Empty()
                    if ok
                        ndeleted = record_deleted!(dict)
                        approx_len = dict.nadded[] - ndeleted
                        half_len = length(slots) ÷ 2
                        if length(slots) > 8 && approx_len < half_len
                            shrink!(dict, slots)
                        end
                        return reply
                    end
                else
                    # Update:
                    new = convert(Value, something(reply::Some))
                    vstored, ok = @atomicreplace ref.value vstored => new
                    ok && return reply
                end
                vstored = vstored::ValueUnion{Value}
            end

            nprobes += 1
            if nprobes > c ÷ 4
                @label too_many_probes
                let newslots = @atomic dict.slots
                    if (slots === newslots) || newslots === nothing
                        # @info "expand: length(slots) ≈ 2^$(floor(Int, log2(length(slots))))"
                        # global DICT = dict
                        # TODO: Check the approximated table size here. It's
                        # possible that the hash table needs cleanup but not
                        # resize (i.e., too many deleted slots).
                        slots = expand!(dict, slots)
                    else
                        slots = load_slots(dict)
                    end
                end
                @goto restart
            end
        end  # ...of the inner "CAS loop"
    end  # ...of the restart loop
end

expand!(dict, oldslots) = migrate!(dict, oldslots, true)
shrink!(dict, oldslots) = migrate!(dict, oldslots, false)

function migrate!(dict::LinearProbingDict, expand::Bool; basesize = nothing)
    slots = load_slots(dict)
    return migrate!(dict, slots, expand; basesize)
end

function load_slots(dict)
    while true
        slots = @atomic dict.slots
        if slots === nothing
            return finishmove!(dict, slots)
        else
            return slots
        end
    end
end

function migrate!(dict::LinearProbingDict, oldslots, expand; basesize = nothing)
    # Since the migration is parallelized, workers running tasks blocked by the
    # lock actually will contribute to the forward progress of the entire
    # system. (Though the OS may suspend this worker thread before the tasks are
    # spawned.)
    lock(dict.migration) do
        slots = (@atomic dict.slots)::slots_type(dict)
        if slots !== oldslots
            return slots
        end
        @atomic dict.slots = nothing
        newslots = similar(slots, expand ? length(slots) * 2 : length(slots) ÷ 2)
        if expand
            expand_parallel!(newslots, slots, basesize)
        else
            if migrate_serial!(newslots, slots, Val(true)) === nothing
                # Shrinking the slots have failed and rolled back. Keep using
                # the old slots:
                @atomic dict.slots = slots
                return slots
            end
        end
        # TODO: parallelize `shrink!`

        # This is the atomic "publishing" operation that makes the `newslots`
        # accessible to any tasks (including the ones that are/were not trying
        # to acquire the `migration` lock).
        @atomic dict.slots = newslots

        return newslots
    end
end

function finishmove!(dict::LinearProbingDict, oldslots)
    lock(dict.migration) do
        slots = (@atomic dict.slots)::slots_type(dict)
        # TODO: `oldslots` was used for a sanity check that `slots` is updated.
        # However, now that shrink can rollback migration, it's not correct any
        # more. Still keeping this code since it might make sense to actually
        # not support shrinking.
        #=
        @assert oldslots !== slots
        =#
        return slots
    end
end

struct Stopped
    i::Int
    nadded::Int
end

"""
    expand_parallel_basecase!(...) -> (nadded, seen_empty)

Process all clusters started within `start0:(start0 + basesize)` (mod `length(slots)`).

That is to say, _try_ to process `start0:(start0 + basesize)` but make sure to
avoid stepping into other base cases by using an empty slot to delimite the base
case boundaries. This uses the property of linear probing dictionary that
non-overlapping clusters (i.e., continuous regions of slots that are non-empty)
are mapped to non-overlapping regions when the `slots` array is doubled in size.

More precisely:

1. Process all clusters started within `start0:(start0 + basesize - 1)`.
2. If at least one cluster is processed, process the cluster that contains the
   start position of the next chunk `start0 + basesize` (mod `length(slots)`).

Unlike the original trick mentioned in Maier et al. (2019), there is a
difference due to that `DictSlot` is reused in the `newslots`.  Since the slot
is re-used, we can't use it as the reliable marker for delimiting the clusters.
Other tasks can empty a slot at any moment.  So, instead, each basecase sets its
own promise `firstempties[ichunk]::Promise{Int}` to the index of the first empty
slot that it sees and successfully CAS'ed to `Moved{Empty}`.  If no empty slot
is observed, the promise is set to -1.  When the basecase reaches the end of
chunk (`start0 + basesize`), it confirms the end of the cluster it owns by
finding out the first valid `firstempties[mod1((ichunk+i),end)]` with the
smallest `i > 0`. If consecutive basecase tasks are started at the same time, by
the time a basecase task needs the promise, it is likely already is ready since
that's the first thing that the other task does.

---

This trick is mentioned in:

> Maier, Tobias, Peter Sanders, and Roman Dementiev. “Concurrent Hash Tables:
> Fast and General(?)!” ACM Transactions on Parallel Computing 5, no. 4
> (February 22, 2019): 16:1–16:32. https://doi.org/10.1145/3309206.
"""
function expand_parallel_basecase!(newslots, slots, basesize, start0, ichunk, firstempties)
    c = length(slots)
    stop0 = min(start0 + basesize - 1, c)
    stpd = migrate_serial_nofill!(nothing, slots, start0, stop0, Val(true))
    if stpd isa Int
        @assert stpd == 0
        # This chunk does not own any clusters.
        put!(firstempties[ichunk], -1)
        return (0, false)
    end

    migrate_between(start, stop) = migrate_serial_nofill!(newslots, slots, start, stop)::Int

    # An empty slot is observed. There is at least one cluster started within
    # this chunk.
    stpd::Stopped
    @assert stpd.nadded == 0
    put!(firstempties[ichunk], stpd.i)
    nadded = migrate_between(stpd.i + 1, stop0)

    # Process the cluster that includes `start0 + basesize` (if any).
    chunk_starts = (
        start0+basesize:basesize:c,
        1:basesize:start0-1,  # wrap around
    )
    chunk_indices = (
        ichunk+1:length(firstempties),
        1:ichunk-1,  # wrap around
    )
    for half in 1:2
        for (jchunk, start) in zip_strict(chunk_indices[half], chunk_starts[half])
            nextempty = fetch(firstempties[jchunk])
            if nextempty == -1  # next chunk is owned by this task
                stop = min(start + basesize - 1, c)
                nadded += migrate_between(start, stop)
            else
                nadded += migrate_between(start, nextempty - 1)
                return (nadded, true)
            end
        end
    end
    # Reaching here means there was only one empty slot.
    nadded += migrate_between(start0, stpd.i - 1)
    return (nadded, true)
end

plus_or((a, b), (c, d)) = (a + c, b | d)

# See`BenchDictMigration` for benchmarking this:
const LINEAR_PROBING_DICT_EXPAND_BASESIZE = Ref(2^13)

function expand_parallel!(newslots, slots, basesize)
    # TODO: Make the default `basesize` configurable?
    basesize = something(basesize, LINEAR_PROBING_DICT_EXPAND_BASESIZE[])
    @assert length(newslots) > length(slots)

    if length(slots) <= basesize || Threads.nthreads() == 1
        return migrate_serial!(newslots, slots)
    end
    basesize = min(basesize, cld(length(slots), Threads.nthreads()))

    # TODO: since `nadded` from this is not used any more, `|` can be used
    # instead of `plus_or`.
    c = length(slots)
    chunkstarts = 1:basesize:c
    firstempties = [Promise{Int}() for _ in chunkstarts]
    nadded, seen_empty = threaded_typed_mapreduce(
        Tuple{Int,Bool},  # return type of the do block and the doamin of `plus_or`
        plus_or,
        eachindex(chunkstarts),
    ) do ichunk
        start0 = chunkstarts[ichunk]
        return expand_parallel_basecase!(
            newslots,
            slots,
            basesize,
            start0,
            ichunk,
            firstempties,
        )
    end
    if seen_empty
        # fill_undef!(newslots)  # TODO: parallel?
        Threads.@threads for i in eachindex(newslots)
            @inbounds fill_undef_at!(newslots, i)
        end
        return nadded
    else
        # The `slots` are all non-empty. Fallback to serial migration:
        return migrate_serial!(newslots, slots)
    end
end

function migrate_serial!(
    newslots,
    slots,
    rollback_on_error::Union{Val{false},Val{true}} = Val(false),
)
    nadded = migrate_serial_nofill!(
        newslots,
        slots,
        1,
        length(slots),
        Val(false),
        rollback_on_error,
    )
    if nadded === nothing
        rollback_on_error::Val{true}
        return nothing
    end
    nadded = nadded::Int
    fill_undef!(newslots)
    return nadded
end

Base.@propagate_inbounds function fill_undef_at!(newslots, i)
    if !isassigned(newslots, i)
        newslots[i] = eltype(newslots)()
    end
end

function fill_undef!(newslots)
    for i in eachindex(newslots)
        @inbounds fill_undef_at!(newslots, i)
    end
    return newslots
end

# TODO: Create "defrag" version that actually copies keys/values into pre-filled
# (continuously allocated) slots?
"""
Migrate `DictSlots` entries to `slots` to `newslots`.
"""
function migrate_serial_nofill!(
    newslots::Union{AbstractVector{Slot},Nothing},
    slots::AbstractVector{Slot},
    start::Int,
    stop::Int,
    # TODO: use custom singletong rather than Val:
    stop_on_empty::Union{Val{false},Val{true}} = Val(false),
    rollback_on_error::Union{Val{false},Val{true}} = Val(false),
) where {K,V,Slot<:DictSlot{K,V}}
    nadded = 0
    for i in start:stop
        ref = slots[i]
        while true
            local value
            value = @atomic ref.value
            if value isa Empty
                _, ok = @atomicreplace ref.value Empty() => MOVED_EMPTY
                ok || continue
                stop_on_empty == Val(true) && return Stopped(i, nadded)
                break  # next index
            elseif value isa Moved{Empty}
                @static_error("unexpected Moved{Empty}")
            end
            @goto move
        end
        continue
        @label move

        newslots === nothing && continue

        h = @atomic :monotonic ref.hash
        if iszero(h)
            key = @atomic ref.key
            key = key::K
            h = hash(key)
        end

        # Insertion to `newslots` does not have to use atomics since
        # it's protected by the `.migration` lock.
        c = length(newslots)
        offset = reinterpret(Int, h) & (c - 1)  # h % c
        nprobes = 0
        while true
            index = offset + 1
            offset = (offset + 1) & (c - 1)  # (offset + 1) % c

            if !isassigned(newslots, index)
                newslots[index] = ref
                nadded += 1
                break
            end

            nprobes += 1
            if nprobes > c
                if rollback_on_error === Val(true)
                    rollback_migration!(slots, start, i)
                    return nothing
                end
                @static_error "unreachable: too many probings during migration"
            end
        end
    end
    return nadded
end

function rollback_migration!(slots, start, stop)
    for i in start:stop
        ref = slots[i]
        value = @atomic ref.value
        while true
            if value isa Moved{Empty}
                value, ok = @atomicreplace ref.value value => Empty()
                ok && break
                # TODO: maybe simple store is OK, since no one should be touching it?
            elseif value isa Moved
                @static_error("unexpected Moved{Value}")
            else
                break
            end
        end
    end
end

Base.IteratorSize(::Type{<:LinearProbingDict}) = Base.SizeUnknown()
Base.IteratorSize(::Type{<:Base.KeySet{<:Any,<:LinearProbingDict}}) = Base.SizeUnknown()
Base.IteratorSize(::Type{<:Base.ValueIterator{<:LinearProbingDict}}) = Base.SizeUnknown()

function Base.iterate(dict::LinearProbingDict)
    slots = load_slots(dict)
    return iterate(dict, (slots, 1))
end

function Base.iterate(::LinearProbingDict, (slots, index))
    index < 1 && return nothing
    while true
        index > length(slots) && return nothing
        ref = slots[index]
        value = @atomic ref.value
        if value isa Empty
        elseif value isa Moved{Empty}
        else
            if value isa Moved
                value = value.value
            end
            key = @atomic ref.key
            return (key => value), (slots, index + 1)
        end
        index += 1
    end
end

kvtype(::Type{Pair{Key,Value}}) where {Key,Value} = (Key, Value)
kvtype(::Type{Tuple{Key,Value}}) where {Key,Value} = (Key, Value)
kvtype(::Type) = nothing

# TODO: batch and/or racy construction
function LinearProbingDict{Key,Value}(kvs) where {Key,Value}
    dict = LinearProbingDict{Key,Value}()
    for (k, v) in kvs
        dict[k] = v
    end
    return dict
end

function ConcurrentCollections.ConcurrentDict(kvs)
    if IteratorEltype(kvs) isa HasEltype
        typ = kvtype(eltype(kvs))
        if typ !== nothing
            return LinearProbingDict{first(typ),last(typ)}(kvs)
        end
    end
    kvs = [k => v for (k, v) in kvs]
    typ = kvtype(eltype(kvs))
    return LinearProbingDict{first(typ),last(typ)}(kvs)
end

ConcurrentCollections.ConcurrentDict{Key,Value}(kvs) where {Key,Value} =
    LinearProbingDict{Key,Value}(kvs)

function ConcurrentCollections.ConcurrentDict(kv::Pair, kvs::Pair...)
    ks = promote(first(kv), map(first, kvs)...)
    vs = promote(last(kv), map(last, kvs)...)
    allkvs = map(=>, ks, vs)
    Key = typeof(first(ks))
    Value = typeof(first(vs))
    return LinearProbingDict{Key,Value}(allkvs)
end

function ConcurrentCollections.ConcurrentDict{Key,Value}(
    kv::Pair,
    kvs::Pair...,
) where {Key,Value}
    allkvs = map((kv, kvs...)) do (k, v)
        convert(Key, k) => convert(Value, v)
    end
    return LinearProbingDict{Key,Value}(allkvs)
end

function Base.summary(io::IO, dict::LinearProbingDict{Key,Value}) where {Key,Value}
    @nospecialize dict
    print(io, ConcurrentCollections.ConcurrentDict)
    print(io, '{')
    print(io, Key)
    print(io, ", ")
    print(io, Value)
    print(io, '}')
end

function Base.show(io::IO, ::MIME"text/plain", dict::LinearProbingDict)
    @nospecialize dict
    summary(io, dict)
    n = 0
    for kv in dict
        n += 1
        println(io)
        if n > 3
            print(io, "  ⋮")
            break
        end
        print(io, "  ", kv)
    end
end

"""
    clusters(dict) -> ranges
    clusters(slots) -> ranges

Compute clusters in the slots. Used for performance debugging.

```julia
using ConcurrentCollections
using ConcurrentCollections.Implementations: clusters
d = ConcurrentDict(1:1000 .=> 0)
cs = clusters(d.slots)

using UnicodePlots
histogram(map(length, cs))

using StatsBase
describe(map(length, cs))
```
"""
clusters(d::LinearProbingDict) = clusters(d.slots)
function clusters(slots::AbstractVector)
    cs = typeof(1:2)[]
    i = 1
    while true
        while true
            i > length(slots) && return cs
            ref = slots[i]
            key = @atomic ref.key
            i += 1
            if key isa Union{Empty,Moved{Empty}}
                break
            end
        end
        start = i - 1
        while true
            if i > length(slots)
                push!(cs, start:i-1)
                return cs
            end
            ref = slots[i]
            key = @atomic ref.key
            i += 1
            if key isa Union{Empty,Moved{Empty}}
                break
            end
        end
        push!(cs, start:i-1)
    end
end


text:
# build a factor graph from PS3 eye acoustic data

# Load required libraries
using Caesar, RoMEPlotting
# using Distributions
# using ApproxManifoldProducts

## Environment setup -- in Juno you can run cells with Alt+Enter

# true pose position (assume hidden)  - [X Y θ]
p0 = [0.0;0.0;0.0]

# noise level for Gaussian assumption -- see https://github.com/JuliaStats/Distributions.jl/issues/584
pcov = Matrix(Diagonal([0.01; 0.01; 0.05].^2))
lcov = Matrix(Diagonal([0.01; 0.01].^2))

# true landmark positions
l0 = [0.4841;0.559/2]
l1 = [0.4841;-0.559/2]

# d0 = l0-p0[1:2]
# m0 = atan(d0[2],d0[1])
#
# d1 = l1-p0[1:2]
# m1 = atan(d1[2],d1[1])

# mstd = 0.1

## build the factor graph

fg = initfg()

# add two beacons
addVariable!(fg, :l0, Point2, tags=[:BEACON;])
# addFactor!(fg, [:l0], Prior(MvNormal(l0, lcov)))

addVariable!(fg, :l1, Point2, tags=[:BEACON;])
# addFactor!(fg, [:l1], Prior(MvNormal(l1, lcov)))

addFactor!(fg, [:l0; :l1], Point2Point2(MvNormal(l1-l0, [0.01 0; 0 0.01].^2)))

# add unknown pose location
addVariable!(fg, :x0, Pose2, tags=[:POSE;])
addFactor!(fg, [:x0], PriorPose2(MvNormal(p0, pcov)))

# add two bearing only measurements
addFactor!(fg, [:x0; :l0], Pose2Point2BearingRange(pcLeft, Rayleigh(2.0)))
addFactor!(fg, [:x0; :l1], Pose2Point2BearingRange(pcRight, Rayleigh(2.0)))

writeGraphPdf(fg, show=true)

##

solveTree!(fg)

##

drawPosesLandms(fg, spscale=0.5)

getKDEMean(getKDE(fg, :l0))
getKDEMax(getKDE(fg, :l0))

getKDEMean(getKDE(fg, :l1))
getKDEMax(getKDE(fg, :l1))

## OPTION 2

fg = initfg()

# add two beacons
addVariable!(fg, :l0, Point2, tags=[:BEACON;])
addFactor!(fg, [:l0], Prior(MvNormal(l0, lcov)))

addVariable!(fg, :l1, Point2, tags=[:BEACON;])
addFactor!(fg, [:l1], Prior(MvNormal(l1, lcov)))

# add unknown pose location
addVariable!(fg, :x0, Pose2, tags=[:POSE;])
addFactor!(fg, [:x0], PartialPrior(Normal(0.0, 0.01), (2,)) )
# addFactor!(fg, [:x0], PartialPrior(Uniform(-1.0, 0.5), (1,)) )
# addFactor!(fg, [:x0], PriorPose2(MvNormal(p0, pcov)))

# add two bearing only measurements
addFactor!(fg, [:x0; :l0], Pose2Point2BearingRange(pcLeft, Rayleigh(1.0)))
addFactor!(fg, [:x0; :l1], Pose2Point2BearingRange(pcRight, Rayleigh(1.0)))

writeGraphPdf(fg, show=true)

##

getSolveParams(fg).N = 100
solveTree!(fg)

## look at one mode only

pts = getVal(fg, :x0)
pts0 = pts[:,pts[1,:] .< 1.0]
X00 = kde!(pts0[1,:])
@show getKDEMax(X00)[1]
plotKDE(X00)

##

pl = plotPose(Pose2(), [getKDE(fg, :x0);])

drawPosesLandms(fg, spscale=0.5, meanmax=:mean, xmin=-3.0, xmax=3.0)

getKDEMean(getKDE(fg, :x0))
getKDEMax(getKDE(fg, :x0))

getKDEMean(getKDE(fg, :l0))
getKDEMax(getKDE(fg, :l0))

getKDEMean(getKDE(fg, :l1))
getKDEMax(getKDE(fg, :l1))

##

pl |> PDF("/tmp/test.pdf",30cm,20cm)
pl |> SVG("/tmp/test.svg",15cm,6cm)

##

@async run(`evince /tmp/test.pdf`)

##


text:
function col_iter(x)
s=zero(eltype(x))
for i = 1:size(x, 2)
for j = 1:size(x, 1)
s = s + x[j, i] ^ 2
x[j, i] = s
end
end
end
function row_iter(x)
s=zero(eltype(x))
for i = 1:size(x, 1)
for j = 1:size(x, 2)
s = s + x[i, j] ^ 2
x[i, j] = s
end
end
end

max_stars_repo_path: Module 2/Chapter05/Timestep.jl

max_issues_repo_path: Module 2/Chapter06/col_iter.jl

max_forks_repo_path: Module 2/Chapter06/col_iter.jl


text:
"""
    $(TYPEDSIGNATURES)

Finds all leaves of MIOTree.
"""
find_leaves(mt::MIOTree) = [nd for nd in allnodes(mt) if is_leaf(nd)]

""" Checks if a MIOTree is trained. """
function check_if_trained(mt::MIOTree)
    for nd in allnodes(mt)
        (is_leaf(nd) && isnothing(nd.label)) && return false
        (!is_leaf(nd) && (isnothing(nd.a) || isnothing(nd.b))) && return false
    end
    return true
end

""" Prints problems with MIOTree for debugging. """
function debug_if_trained(mt::MIOTree)
    for nd in allnodes(mt)
        (is_leaf(nd) && isnothing(nd.label)) && println("Leaf $(nd.idx) has no label.")
        (!is_leaf(nd) && (isnothing(nd.a) || isnothing(nd.b))) && println("Split $(nd.idx) has no hyperplane data.")
    end
end

"""
    $(TYPEDSIGNATURES)

Returns the hyperplane data from an MIOTree,
in the format Dict[leaf_number] containing [B0, B].
"""
function trust_region_data(mt::MIOTree)
    all_leaves = find_leaves(mt)
    upperDict = Dict()
    lowerDict = Dict()
    for leaf in all_leaves
        line = [leaf, lineage(leaf)...]
        ups = [line[i] for i = 2:length(line) if line[i-1] == line[i].left]
        lows = [line[i] for i = 2:length(line) if line[i-1] == line[i].right]
        upperDict[leaf.idx] = [[nd.b, nd.a] for nd in ups]
        lowerDict[leaf.idx] = [[nd.b, nd.a] for nd in lows]
    end
    return upperDict, lowerDict
end

"""
    $(TYPEDSIGNATURES)

Returns the regression weights from an OptimalTreeLearner by leaf,
i.e. Dict[leaf_number] containing [B0, B].
"""
function pwl_constraint_data(lnr::MIOTree, vks)
    @assert get_param(lnr, :regression)
    return Dict(leaf.idx => leaf.label for leaf in find_leaves(lnr))
end

function deepen_to_max_depth!(mt::MIOTree)
    md = get_param(mt, :max_depth)
    queue = allnodes(mt)
    idx = maximum([nd.idx for nd in queue])
    while !isempty(queue)
        node = popfirst!(queue)
        if is_leaf(node) && depth(node) < md
            idx += 1
            leftchild(node, BinaryNode(idx))
            idx += 1
            rightchild(node, BinaryNode(idx))
            append!(queue, [node.left, node.right])
        end
    end
    return
end

function deepen_one_level!(mt::MIOTree)
    md = get_param(mt, :max_depth)
    queue = allnodes(mt)
    idx = maximum([nd.idx for nd in queue])
    while !isempty(queue)
        node = popfirst!(queue)
        if is_leaf(node) && depth(node) < md
            idx += 1
            leftchild(node, BinaryNode(idx))
            idx += 1
            rightchild(node, BinaryNode(idx))
        end
    end
    return
end


text:
module AlgThys

# Algebraic Theories
# Owen & Christian [Fall 2021]

# TODO
# DATA
# data structure that contains names for operations
# easy way of parsing a presentation

# CATEGORY
# generate the theory from the presentation
#       [obj: contexts, mor: terms, comp: substitution, 2mor: equations]

# MORPHISMS
# morphisms of presentations
#       [send op to op]
# morphisms of theories
# functorial: sends an equation to a composite equation
#       [send op to composite term]
#       [push forward a term]
# pushouts of morphisms

# REWRITING
# contexts with equations?

using Reexport

include("LabeledVectors.jl")
include("Theories.jl")
include("Parsing.jl")
include("Examples.jl")

@reexport using .LabeledVectors
@reexport using .Theories
@reexport using .Parsing

end


text:
using HPAT

#HPAT.DomainPass.set_debug_level(3)
#HPAT.DistributedPass.set_debug_level(3)

@acc hpat function cumsum_test1(n)
    A = ones(n)
    B = cumsum(A)
    return sum(B)
end

using Base.Test
@test_approx_eq cumsum_test1(5) 15.0


text:
function export_tox_consts(filename)
	consts = ASCIIString[]
	file = open(filename, "r")

	while(!eof(file))
		line = readline(file)

		if ismatch(r"const .*", line)
			# tox.h
			m = match(r"TOX_([A-Z,_])*", line)
			if m != nothing
				push!(consts, m.match)
			end

			# toxav.h
			m = match(r"Tox([A-z,_])*", line)
			if m != nothing
				push!(consts, m.match)
			end

			m = match(r"av_([A-z,_])*", line)
			if m != nothing
				push!(consts, m.match)
			end
		end
	end

	for c in consts
		println("export $(c)")
	end
end

function export_tox_typealias(filename)
	typealiases = ASCIIString[]
	file = open(filename, "r")

	while(!eof(file))
		line = readline(file)

		if ismatch(r"typealias .*", line)

			# tox.h
			m = match(r"TOX_([A-Z,_])*", line)
			if m != nothing
				push!(typealiases, m.match)
			end

			m = match(r"tox_([A-z,_])*", line)
			if m != nothing
				push!(typealiases, m.match)
			end

			# toxav.h
			m = match(r"Tox([A-z,_])*", line)
			if m != nothing
				push!(typealiases, m.match)
			end

			m = match(r"av_([A-z,_])*", line)
			if m != nothing
				push!(typealiases, m.match)
			end
		end
	end

	for t in typealiases
		println("export $(t)")
	end
end

println("# typealias")
export_tox_typealias(  Pkg.dir("Toxcore.jl", "src", "gen", "tox.h.jl") )
println("")
println("# consts")
export_tox_consts(  Pkg.dir("Toxcore.jl", "src", "gen", "tox.h.jl") )
println("")
println("# typealias")
export_tox_typealias(  Pkg.dir("Toxcore.jl", "src", "gen", "toxav.h.jl") )
println("")
println("# consts")
export_tox_consts( Pkg.dir("Toxcore.jl", "src", "gen", "toxav.h.jl") )


text:
module KDESieve

import Distributions
import KernelDensity
import LinearAlgebra

include("comparator.jl")
include("sieve.jl")

export Sieve, predict

end # module


text:
using Documenter
using LazyGroupBy

makedocs(
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.makedocs
    sitename = "LazyGroupBy",
    format = Documenter.HTML(),
    pages = ["Home" => "index.md", hide("internals.md")],
    modules = [LazyGroupBy],
    doctest = false,
)

deploydocs(
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
    repo = "github.com/JuliaFolds/LazyGroupBy.jl",
    push_preview = true,
)


text:
abstract type AbstractDLRProblem end
abstract type AbstractDLRSolution end
abstract type AbstractDLRIntegrator end
abstract type AbstractDLRAlgorithm end
abstract type AbstractDLRAlgorithm_Cache end

"""
    Problem of approximating the solution of a matrix differential equation
    dy/dt = f(t,y) with y(0) = y₀ on [t0,tf] with a low rank factorization
    y(t) ≈ u(t) = U(t)S(t)V(t)' where U(t) and V(t) are orthonormal bases approximating
    range and co-range of y(t).
"""
mutable struct MatrixDEProblem{fType, uType, tType} <: AbstractDLRProblem
    f::fType
    u0::uType
    tspan::Tuple{tType, tType}
end

"""
    Problem of tracking the low rank decomposition u(t) = U(t)S(t) V(t)' (or U(t)Z(t)') of a
    time-dependent (or streamed) matrix y(t) with t ∈ [t_0, t_f].
"""
mutable struct MatrixDataProblem{yType, uType, tType} <: AbstractDLRProblem
    y::yType
    u0::uType
    tspan::Tuple{tType, tType}
end
function MatrixDataProblem(y::AbstractArray, u0)
    return MatrixDataProblem(y, u0, (1,length(y)))
end

"""
    Problem of identifying optimal projection basis U(t) such that we get good reconstruction
    y(t) = U(t) Z(t)' where dZ/dt = f(Z,U,t).
"""
mutable struct MatrixHybridProblem{yType, fType, uType, tType} <: AbstractDLRProblem
    y::yType
    f::fType
    u0::uType
    tspan::Tuple{tType, tType}
end

"""
    Solution object that tracks the evolution of a low rank approximation
"""
mutable struct DLRSolution{solType, tType} <: AbstractDLRSolution
    Y::Vector{solType}
    t::Vector{tType}
end

"""
    Integrator computing solution to a dynamic low rank approximation problem
"""
mutable struct DLRIntegrator{uType, tType, aType, cType, pType} <: AbstractDLRIntegrator
    u::uType
    t::tType
    dt::tType
    sol::DLRSolution{uType,tType}
    alg::aType
    cache::cType
    probType::pType
    iter::Int
end

"""
    solves the given problem with the specified algorithm and step size
"""
function solve(prob::AbstractDLRProblem, alg::AbstractDLRAlgorithm, dt)
    integrator = init(prob, alg, dt)
    T = prob.tspan[2] - prob.tspan[1]
    while (prob.tspan[2]-integrator.t)/T > 1e-8
        step!(integrator, alg, dt)
        update_sol!(integrator)
    end
    return integrator.sol
end
function solve(prob::MatrixDataProblem, alg::AbstractDLRAlgorithm)
    @assert typeof(prob.y) <: AbstractArray "If the data is not provided as array, integration stepsize needs to be specified"
    return solve(prob, alg, 1)
end

function update_sol!(integrator::AbstractDLRIntegrator)
    if integrator.iter <= length(integrator.sol.Y) - 1
        integrator.sol.Y[integrator.iter + 1] = deepcopy(integrator.u)
        integrator.sol.t[integrator.iter + 1] = integrator.t
    else
        push!(integrator.sol.Y, deepcopy(integrator.u))
        push!(integrator.sol.t, integrator.t)
    end
end

function init_sol(dt, t0, tf, u0)
    n = floor(Int,(tf-t0)/dt) + 1
    # compute more sensible dt # rest will be handled via interpolation/save_at
    dt = (tf-t0)/(n-1)
    # initialize & return solution object
    return DLRSolution(Vector{typeof(u0)}(undef, n), collect(range(t0, tf, length=n)))
end

function init_sol(dt::Int, t0, tf, u0)
    # initialize solution object
    steps = t0:dt:tf
    return DLRSolution(Vector{typeof(u0)}(undef, length(steps)), collect(steps))
end


text:
module mechanism
include("modlspc.jl")
include("assign2_rk.jl")
include("euler.jl")
include("assign3a.jl")
include("chemODE.jl")
include("stiff.jl")
#using DifferentialEquations
using LinearAlgebra
import Random
export gas

function gas(t1,t2,tm,days::Int64)
#=**************************************************************************
This is the main function in an ozone NOx formaldehyde mechanism, adapted from Fortran
by Obin Sturm (UC Davis) in October 2019.  It has been edited in order to generate training
data for machine learning and restricted inverse algorithms.

Notable changes:
Oct 2019 -- Random seeding of concentrations, looping over days, incorporating
            fluctuating daylight constant (HV) values.
         -- Option to use stiff equation solver.  To use this
            integration method, set method to "stiff!" instead of "euler!"

Mar 2020 -- Implementation of a NaN in the first column of S output matrix
            when using "euler" method.  This is used as a flag for when
            concentrations are negative and a mass balance violation occurs
            (a limitation of the original Fortran mechanism).

The original comments and file description are below.
     created by: Mike Kleeman (April 1999)
             University of California, Davis
             ECI 241 Air Quality Modeling
             Assignment #3
            The purpose of this program is to act as a driver for the evaluation and
            integration of a model ozone mechanism involving NOx and formaldehyde.
            Note: the "active" array contains active species in the first "maxact"
            locations and buildup species in the remaining "maxbo" locations.
**************************************************************************=#
# names for active species.
#    name = ["O3", "NO", "NO2","HCHO","HO2.","HO2H","HO.","O","HNO3","CO","H2"]

    Random.seed!(43) # Seed, for replicability. 43 was arbitrary 😁

    constn[1] = Float32(2.09e+05) # O2 concentration in ppm

# get the reaction rate constants at 1 atm and 25C--
    tempk = Float64(298.0)
    press = Float64(1.0)
    #      tempk = 283.
    #      press = 0.7
    global rk = getrk!(tempk,press,rk)

# call the integration subroutine to advance the solution
#    t1 = Float64(0.0)
#    t2 = Float64(60.0*24.0*1.0) # 60 mins/hr * 24 hrs/day *1 day
    dt = Float64(1.0e-3)
#    tm = Float64(6.0)
    C = zeros(1,11)
    S = zeros(1,10)
    J = zeros(1,1)

    Random.seed!(42)
    for d = 1:days
        println("Running ozone/NOx/formaldehyde mechanism: day ", d , "  of ", days)
        global c = zeros(Float64,1,maxact+maxbo)
        c[1] = 0.001*10^(2*rand())         # 03 range: 0.001 - 0.1 ppm
        c[2] = 0.0015*10^(2*rand())        # NO range: 0.0015 - 0.15
        c[3] = 0.0015*10^(2*rand())        # NO2 range: 0.0015 - 0.15
        c[4] = 0.02*10^(2*rand())          # HCHO range: 0.02 - 2 ppm
        c[5] = 1.0e-05 * rand()            # HO2. range: 1 - 10 *ppt*
        c[6] = 0.01*rand()	               # HO2H range: 0.001 - 0.01 ppm
        c[7:11] = zeros(Float64,1,5)

        # c[1] = 0.0  # For testing with original Fortran mechanism
        # c[2] = 0.0421
        # c[3] = 0.0679
        # c[4] = 0.1
        # c[5] = 0.0
        # c[6] = 0.01
        # c[7:11] = zeros(Float64,1,5)

        daymax = rand(Float64)

        # set daylight HV value
        constn[2] = 0.0

        global C_day,S_day,J_day = euler!(t1,t2,dt,tm,c,daymax,maxact+maxbo,rk,diff!)

        C = vcat(C,C_day)
        S = vcat(S,S_day)
        J = vcat(J,J_day)

    end
    C = C[2:end,:]
    S = S[2:end,:]
    J = J[2:end,:]
    return C,S,J
end
end


text:
module SeisPlots

using SeisIO, RecipesBase, Dates
import DSP

include("spectrogram.jl")

@recipe function f(C::GphysChannel; max_npts=400000)
    label --> C.id
    seriescolor --> :black
    legendfontsize --> 8
    ytickfontsize --> 8
    xtickfontsize --> 6
    t = u2d.(SeisIO.t_expand(C.t,C.fs) * SeisIO.μs)
    t, C.x
end

@recipe function f(S::GphysData; max_npts=400000, max_traces=10)
    layout := (min(S.n,max_traces),1)
    seriescolor --> :black

    # get minumum and maximum times
    starttimes = u2d.(SeisIO.μs .* [SeisIO.starttime(S[i].t,S[i].fs) for i = 1:S.n])
    endtimes = u2d.(SeisIO.μs .* [SeisIO.endtime(S[i].t,S[i].fs) for i = 1:S.n])
    xlims --> (minimum(starttimes),maximum(endtimes))
    legendfontsize --> max(3,8-S.n)
    ytickfontsize --> max(3,8-S.n)
    xtickfontsize --> 6

    # plot each series
    for i in 1:min(S.n,max_traces)
        @series begin
            subplot := i
            label --> S[i].id
            t = u2d.(SeisIO.t_expand(S[i].t,S[i].fs) * SeisIO.μs)
            if i != min(S.n,max_traces)
                xticks --> Int[]
            end
            t,S[i].x
        end
    end
end

end


text:
using Test
using BiApprox

function test1()
    L(w, a, b) = a + w^2 - a*cos( π * b * w)
    f(x,y; n=length(y)-1) = abs( 4x[n+1] - x[n+2] ) * y[end]^2 + sum(L.( abs.(x[1:n]-y[1:n]), x[n+1], x[n+2]))
    F(x,y; n=length(y)-1) = abs( 4x[n+1] - x[n+2] ) / (1+y[end]^2) + sum(L.( abs.(x[1:n]-y[1:n]), x[n+1], x[n+2]))

    Ψ(X::Matrix) = begin
        n = size(X,2)-2
        Y = zeros(N, n+1)
        Y[:,1:n] = X[:,1:n]
        Y
    end

    n = 5
    N = 300
    D_ul = n+2
    D_ll = n+1

    X = -10 .+ 20rand(N, D_ul)
    X[:, n+1:n+2] = 2 .+ 8rand(N, 2)

    Y = Ψ(X)

    Fs =  [ F(X[i,:], Y[i,:]) for i =1:N ]

    X_test = -10 .+ 20rand(N, D_ul)
    X_test[:, n+1:n+2] = 2 .+ 8rand(N, 2)
    Y_test = Ψ(X_test)
    Fs_test =  [ F(X_test[i,:], Y_test[i,:]) for i =1:N ]

    F̂s_test, α = kernel_approx_ul(Fs,X_test, X)

    true

end

function test2()
    f(x) = sum(x.^2)
    N = 100
    D = 11
    X = randn(N, D)

    method = KernelInterpolation(f, X)
    train!(method)

    X_test = rand(97, D)

    ff = evaluate(X_test, method)

    method2 = KernelInterpolation(f, X)
    train!(method2)
    ff2 = [evaluate(X_test[i,:], method2) for i = 1:size(X_test, 1)]

    return sum(abs.(ff - ff2)) ≈ 0.0

end

function test3()
    f(x) = sum(x.^2)
    N = 100
    D = 11
    X = randn(N, D)
    X2 = randn(N, D)

    model = KernelInterpolation(f, X; λ = 0.1)
    train!(model)
    f_approx = approximate(model)

    model2 = KernelInterpolation(f, X2; λ = 0.1)
    train!(model2)
    f_approx2 = approximate(model)

    X_test = rand(97, D)

    ff = [f_approx( X_test[i,:] ) for i = 1:15 ]
    ff2 = [f_approx2( X_test[i,:] ) for i = 1:15 ]

    return sum(abs.(ff - ff2)) ≈ 0.0

end

function test4()
    Ψ(x) = [ sum(x.^2), x[1] - x[2].^2, sin(x[3]), abs(x[1]) ]

    X = -10 .+ 20rand(100, 3)
    Y = zeros(100, 4)

    for i = 1:size(Y,1)
        Y[i,:] = Ψ(X[i,:])
    end

    model = LL_approx(X, Y, p_train=0.8, scale_values=0.1ones(3))

    bruteforce!(model, debug=false)

    sum(model.approx_errors) >= 0

end

@test test1()
@test test2()
@test test3()
@test test4()


text:
abstract type AbstractRBArray{T, N} <: AbstractArray{T, N} end

abstract type RBEven end
abstract type RBOdd end

const rb_even = Int(2)
const rb_odd  = Int(1)

const rb_bool_even = true
const rb_bool_odd  = false

@inline function rbidx( nidx::Int )::Int
    return div(nidx, 2) + mod(nidx, 2) + 1
end

"""
    nidx( rbidx::Int, ::Val{true}, ::Val{true})::Int

first type argument:  type of the orgal point (for even points -> `Val{true}()`, else `Val{false}()`)
second type argument: is sum of other point indices even or odd -> (if sum is even -> `Val{true}()`, else `Val{false}()`)
"""
@inline function nidx( rbidx::Int, ::Val{true}, ::Val{true})::Int
   return (rbidx - 1) * 2
end
@inline function nidx( rbidx::Int, ::Val{true}, ::Val{false})::Int
   return (rbidx - 1) * 2 - 1
end
@inline function nidx( rbidx::Int, ::Val{false}, ::Val{true})::Int
   return (rbidx - 1) * 2 - 1
end
@inline function nidx( rbidx::Int, ::Val{false}, ::Val{false})::Int
   return (rbidx - 1) * 2
end

"""
    get_rbidx_right_neighbour(rbidx::Int, ::Val{true}, ::Val{true})::Int

needs docu...
"""
@inline function get_rbidx_right_neighbour(rbidx::Int, ::Val{true}, ::Val{true})::Int
    return rbidx + 1
end
@inline function get_rbidx_right_neighbour(rbidx::Int, ::Val{true}, ::Val{false})::Int
    return rbidx
end
@inline function get_rbidx_right_neighbour(rbidx::Int, ::Val{false}, ::Val{true})::Int
    return rbidx
end
@inline function get_rbidx_right_neighbour(rbidx::Int, ::Val{false}, ::Val{false})::Int
    return rbidx + 1
end

"""
    RBExtBy2Array( et::Type, g::Grid{T, N, :cylindrical} )::Array{et, N + 1} where {T, N}

Returns a RedBlack array for the grid `g`.
"""
function RBArray( et::Type, g::Grid{T, N, :cylindrical} )::Array{et, N + 1} where {T, N}
    nr, nφ, nz = size(g)
    # new ordering in memory: r, φ, z -> z, φ, r (so inner loop goes over z)
    return zeros(T, div(nz, 2) + mod(nz, 2), nφ, nr, 2)
end
function RBArray( a::Array{T, N}, grid::Grid{TG, N, :cylindrical} )::Array{T, N + 1} where {T, N, TG}
    rbarray::Array{T, N + 1} = RBArray(T, grid)
    for iz in axes(a, 3)
        irbz::Int = div(iz, 2) + mod(iz, 2)
        for iφ in axes(a, 2)
            idxsum::Int = iz + iφ
            for ir in axes(a, 1)
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                rbarray[ irbz, iφ, ir, rbi ] = a[ir, iφ, iz]
            end
        end
    end
    return rbarray
end

"""
    RBExtBy2Array( et::Type, g::Grid{T, N, :cylindrical} )::Array{et, N + 1} where {T, N}

Returns a RedBlack array for the grid `g`. The RedBlack array is extended in its size by 2 in each geometrical dimension.
"""
function RBExtBy2Array( et::Type, g::Grid{T, N, :cylindrical} )::Array{et, N + 1} where {T, N}
    nr, nφ, nz = size(g)
    # new ordering in memory: r, φ, z -> z, φ, r (so inner loop goes over z)
    return zeros(T, div(nz, 2) + mod(nz, 2) + 2, nφ + 2, nr + 2, 2)
end
function RBExtBy2Array( a::Array{T, N}, grid::Grid{TG, N, :cylindrical} )::Array{T, N + 1} where {T, N, TG}
    rbarray::Array{T, N + 1} = RBExtBy2Array(T, grid)
    for iz in axes(a, 3)
        irbz::Int = rbidx(iz)
        for iφ in axes(a, 2)
            irbφ::Int = iφ + 1
            idxsum::Int = iz + iφ
            for ir in axes(a, 1)
                irbr::Int = ir + 1
                rbi::Int = iseven(idxsum + ir) ? rb_even::Int : rb_odd::Int
                rbarray[ irbz, irbφ, irbr, rbi ] = a[ir, iφ, iz]
            end
        end
    end
    return rbarray
end

"""
    RBExtBy2Array( et::Type, g::Grid{T, 3, :cartesian} )::Array{et, 4} where {T}

Returns a RedBlack array for the grid `g`. The RedBlack array is extended in its size by 2 in each geometrical dimension.
"""
function RBExtBy2Array( et::Type, g::Grid{T, 3, :cartesian} )::Array{et, 4} where {T}
    nx, ny, nz = size(g)
    return zeros(T, div(nx, 2) + mod(nx, 2) + 2, ny + 2, nz + 2, 2)
end
function RBExtBy2Array( a::Array{T, 3}, grid::Grid{TG, 3, :cartesian} )::Array{T, 4} where {T, TG}
    rbarray::Array{T, 4} = RBExtBy2Array(T, grid)
    for iz in axes(a, 3)
        irbz::Int = iz + 1
        for iy in axes(a, 2)
            irby::Int = iy + 1
            idxsum::Int = iz + iy
            for ix in axes(a, 1)
                irbx::Int = rbidx(ix)
                rbi::Int = iseven(idxsum + ix) ? rb_even::Int : rb_odd::Int
                rbarray[ irbx, irby, irbz, rbi ] = a[ix, iy, iz]
            end
        end
    end
    return rbarray
end


text:
using CachedCalls
using Test

@testset "CachedCalls.jl" begin
    @testset "_deconstruct" begin
        @testset "f(1, 2, c)" begin
            ex = :(f(1, 2, c))
            func, args, kw_names, kw_values = CachedCalls._deconstruct(ex)
            @test func == :f
            @test args == Any[1, 2, :c]
            @test kw_names == Any[]
            @test kw_values == Any[]
        end

        @testset "f(;a=1, b=2, c)" begin
            ex = :(f(; a=1, b=2, c))
            func, args, kw_names, kw_values = CachedCalls._deconstruct(ex)
            @test func == :f
            @test args == Any[]
            @test kw_names == [:a, :b, :c]
            @test kw_values == [1, 2, :c]
        end
    end

    @testset "_extract_kwargs" begin
        @test CachedCalls._extract_kwargs(:a) == []
        @test CachedCalls._extract_kwargs(:a; keep_args=true) == [(:a, :a)]
        @test CachedCalls._extract_kwargs(Expr(:kw, :a, 1)) == [(:a, 1),]
        @test CachedCalls._extract_kwargs(Expr(:kw, :a, :variable)) == [(:a, :variable),]

        @testset "f(a=1, b=2, c=3)" begin
            fargs = Any[:($(Expr(:kw, :a, 1))), :($(Expr(:kw, :b, 2))), :($(Expr(:kw, :c, 3)))]
            @test CachedCalls._extract_kwargs(fargs) == [(:a, 1), (:b, 2), (:c, 3)]
        end
        @testset "f(;a=1, b=2, c=3)" begin
            fargs = Any[:($(Expr(:parameters, :($(Expr(:kw, :a, 1))), :($(Expr(:kw, :b, 2))), :($(Expr(:kw, :c, 3))))))]
            @test CachedCalls._extract_kwargs(fargs) == [(:a, 1), (:b, 2), (:c, 3)]
        end

        @testset "f(;a=1, b=2, c)" begin
            fargs = Any[:($(Expr(:parameters, :($(Expr(:kw, :a, 1))), :($(Expr(:kw, :b, 2))), :c)))]
            @test CachedCalls._extract_kwargs(fargs) == [(:a, 1), (:b, 2), (:c, :c)]
        end
    end

    # @cached_call calls the @hash_call macro under the hood, which means that most of the
    # functionality is tested via @cached_call tests. However, the expression passed to
    # @cached_call is escaped (macro hygiene) before being passed to @hash_call. Since
    # @hash_call is exported we also want to test it parses non-escaped expressions
    # correctly and can be used directly. This is why @hash_call only tests the type of the
    # output (it tries to catch errors rather than failed tests)
    @testset "@cached_call and @hash_call" begin
        @testset "f()" begin
            f() = 2
            call1 = @cached_call f()
            call2 = @cached_call f()
            @test f() == call1 == call2

            @test @hash_call(f()) isa UInt
        end

        @testset "f(a, b)" begin
            f(a, b) = a + b
            a = 1
            @test f(a, 2) == @cached_call f(a, 2)
            @test f(a, 2) != @cached_call f(2, 2)

            @test @hash_call(f(a, 2)) isa UInt
        end

        @testset "f(;kw=kw)" begin
            f(;kw=1) = kw
            call1 = @cached_call f(;kw=1)
            call2 = @cached_call f(kw=1)
            kw = 1

            @test @hash_call(f(;kw=1)) isa UInt
            @test @hash_call(f(kw=1)) isa UInt

            @static if VERSION >= v"1.5"
                call3 = @cached_call f(;kw)
                @test f(;kw=1) == call1 == call2 == call3
                @test @hash_call(f(;kw)) isa UInt
            else
                @test f(;kw=1) == call1 == call2
            end
        end

        @testset "f(a ;kw=kw)" begin
            f(a; kw=1) = a - kw
            a = 2.0
            kw = 1
            @test f(a; kw=kw) == @cached_call f(a; kw=kw)
            @test @hash_call(f(a; kw=kw)) isa UInt
        end

        @testset "dot access" begin
            f(a; kw=0) = a + kw
            nt = (one=1, two=2)
            @test f(1, kw=2) == @cached_call f(nt.one; kw=nt.two)
            @test f(1, kw=2) == @cached_call f(nt.one, kw=nt.two)
            @test @hash_call(f(nt.one; kw=nt.two)) isa UInt
        end

        @testset "square bracket access" begin
            f(a; kw=0) = a + kw

            array = [1, 2, 3]
            call1 = @cached_call f(array[1]; kw=array[2])
            call2 = @cached_call f(array[1], kw=array[2])
            @test f(1, kw=2) == call1 == call2
            @test @hash_call(f(array[1], kw=array[2])) isa UInt

            array = [10, 20, 30]
            call3 = @cached_call f(array[1]; kw=array[2])
            call4 = @cached_call f(array[1], kw=array[2])
            @test call1 != call3
            @test call2 != call4
        end

        @testset "splatting" begin
            f(a, b; kw1=1, kw2=2) = a + b + kw1 + kw2

            one = 1
            a = [1, 2]
            kw = (kw1=1, kw2=2)

            @testset "basic splatting" begin
                @test f(a...) == @cached_call f(a...)
                @test f(a...) == @cached_call f([1, 2]...)
                @test f(a...) == @cached_call f([one, 2]...)
                @test f(a...; kw...) == @cached_call f(a...; kw...)
                @test f(a...; kw...) == @cached_call f(a...; (kw1=1, kw2=2)...)
                @test f(a...; kw...) == @cached_call f(a...; (kw1=one, kw2=2)...)

                @test @hash_call(f(a...)) isa UInt
                @test @hash_call(f([1, 2]...)) isa UInt
                @test @hash_call(f([one, 2]...)) isa UInt
                @test @hash_call(f(a...; kw...)) isa UInt
                @test @hash_call(f(a...; (kw1=1, kw2=2)...)) isa UInt
                @test @hash_call(f(a...; (kw1=one, kw2=2)...)) isa UInt
            end

            @testset "different containers" begin
                call1 = @cached_call f(1, 2)
                call2 = @cached_call f([1, 2]...)
                call3 = @cached_call f((1, 2)...)
                @test call1 == call2 == call3

                @test @hash_call(f([1, 2]...)) isa UInt
                @test @hash_call(f((1, 2)...)) isa UInt
            end

            @testset "make sure we hash splatted values (not variable names)" begin
                a = [1, 2]
                call4 = @cached_call f(a...)
                a = [10, 20]
                call5 = @cached_call f(a...)
                @test call4 != call5

                kw = (kw1=1, kw2=2)
                call6 = @cached_call f(1, 2; kw...)
                kw = (kw1=10, kw2=20)
                call7 = @cached_call f(1, 2; kw...)
                @test call6 != call7
            end
        end
    end
end