text:
"""
  Clustering.jl file

Clustering and collaborative filtering (via clustering) algorithms

- [Importable source code (most up-to-date version)](https://github.com/sylvaticus/BetaML.jl/blob/master/src/Clustering.jl) - [Julia Package](https://github.com/sylvaticus/BetaML.jl)
- [Demonstrative static notebook](https://github.com/sylvaticus/BetaML.jl/blob/master/notebooks/Clustering.ipynb)
- [Demonstrative live notebook](https://mybinder.org/v2/gh/sylvaticus/BetaML.jl/master?filepath=notebooks%2FClustering.ipynb) (temporary personal online computational environment on myBinder) - it can takes minutes to start with!
- Theory based on [MITx 6.86x - Machine Learning with Python: from Linear Models to Deep Learning](https://github.com/sylvaticus/MITx_6.86x) ([Unit 4](https://github.com/sylvaticus/MITx_6.86x/blob/master/Unit%2004%20-%20Unsupervised%20Learning/Unit%2004%20-%20Unsupervised%20Learning.md))
- New to Julia? [A concise Julia tutorial](https://github.com/sylvaticus/juliatutorial) - [Julia Quick Syntax Reference book](https://julia-book.com)
"""

"""
    Clustering module (WIP)

Provide clustering methods and missing values imputation / collaborative filtering / reccomendation systems using clustering methods as backend.

The module provides the following functions. Use `?[function]` to access their full signature and detailed documentation:

- [`initRepresentatives(X,K;initStrategy,Z₀)`](@ref initRepresentatives): Initialisation strategies for Kmean and Kmedoids
- [`kmeans(X,K;dist,initStrategy,Z₀)`](@ref kmeans): Classical KMean algorithm
- [`kmedoids(X,K;dist,initStrategy,Z₀)`](@ref kmedoids): Kmedoids algorithm
- [`gmm(X,K;p₀,mixtures,tol,verbosity,minVariance,minCovariance,initStrategy)`](@ref gmm): gmm algorithm over GMM
- [`predictMissing(X,K;p₀,mixtures,tol,verbosity,minVariance,minCovariance)`](@ref predictMissing): Impute mixing values ("matrix completion") using gmm as backbone. Note that this can be used for collaborative filtering / reccomendation systems often with better results than traditional algorithms as k-nearest neighbors (KNN)

{Spherical|Diagonal|Full}Gaussian mixtures for `gmm` / `predictMissing` are already provided. User defined mixtures can be used defining a struct as subtype of `Mixture` and implementing for that mixture the following functions:
- `initMixtures!(mixtures, X; minVariance, minCovariance, initStrategy)`
- `lpdf(m,x,mask)` (for the e-step)
- `updateParameters!(mixtures, X, pₙₖ; minVariance, minCovariance)` (the m-step)

"""
module Clustering

using LinearAlgebra, Random, Statistics, Reexport, CategoricalArrays
import Distributions

using  ForceImport
@force using ..Api
@force using ..Utils

export initRepresentatives, kmeans, kmedoids, gmm, predictMissing

abstract type AbstractMixture end
include("Mixtures.jl")

## Helper functions

"""
  initRepresentatives(X,K;initStrategy,Z₀)

Initialisate the representatives for a K-Mean or K-Medoids algorithm

# Parameters:
* `X`: a (N x D) data to clusterise
* `K`: Number of cluster wonted
* `initStrategy`: Whether to select the initial representative vectors:
  * `random`: randomly in the X space
  * `grid`: using a grid approach [default]
  * `shuffle`: selecting randomly within the available points
  * `given`: using a provided set of initial representatives provided in the `Z₀` parameter
 * `Z₀`: Provided (K x D) matrix of initial representatives (used only together with the `given` initStrategy) [default: `nothing`]
 * `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Returns:
* A (K x D) matrix of initial representatives

# Example:
```julia
julia> Z₀ = initRepresentatives([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.6 38],2,initStrategy="given",Z₀=[1.7 15; 3.6 40])
```
"""
function initRepresentatives(X,K;initStrategy="grid",Z₀=nothing,rng = Random.GLOBAL_RNG)
    X  = makeMatrix(X)
    (N,D) = size(X)
    # Random choice of initial representative vectors (any point, not just in X!)
    minX = minimum(X,dims=1)
    maxX = maximum(X,dims=1)
    Z = zeros(K,D)
    if initStrategy == "random"
        for i in 1:K
            for j in 1:D
                Z[i,j] = rand(rng,Distributions.Uniform(minX[j],maxX[j]))
            end
        end
    elseif initStrategy == "grid"
        for d in 1:D
                # same "space" for each class on each dimension
                Z[:,d] = collect(range(minX[d] + (maxX[d]-minX[d])/(K*2) , stop=maxX[d] - (maxX[d]-minX[d])/(K*2)  , length=K))
                #ex: collect(range(minX[d], stop=maxX[d], length=K))
                #collect(range(s+(e-s)/(K*2), stop=e-(e-s)/(K*2), length=K))
        end
    elseif initStrategy == "given"
        if isnothing(Z₀) error("With the `given` strategy you need to provide the initial set of representatives in the Z₀ parameter.") end
        Z₀ = makeMatrix(Z₀)
        Z = Z₀
    elseif initStrategy == "shuffle"
        zIdx = shuffle(rng,1:size(X)[1])[1:K]
        Z = X[zIdx, :]
    else
        error("initStrategy \"$initStrategy\" not implemented")
    end
    return Z
end

## Basic K-Means Algorithm (Lecture/segment 13.7 of https://www.edx.org/course/machine-learning-with-python-from-linear-models-to)

"""
  kmeans(X,K;dist,initStrategy,Z₀)

Compute K-Mean algorithm to identify K clusters of X using Euclidean distance

# Parameters:
* `X`: a (N x D) data to clusterise
* `K`: Number of cluster wonted
* `dist`: Function to employ as distance (see notes). Default to Euclidean distance.
* `initStrategy`: Whether to select the initial representative vectors:
  * `random`: randomly in the X space
  * `grid`: using a grid approach [default]
  * `shuffle`: selecting randomly within the available points
  * `given`: using a provided set of initial representatives provided in the `Z₀` parameter
* `Z₀`: Provided (K x D) matrix of initial representatives (used only together with the `given` initStrategy) [default: `nothing`]
* `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Returns:
* A tuple of two items, the first one being a vector of size N of ids of the clusters associated to each point and the second one the (K x D) matrix of representatives

# Notes:
* Some returned clusters could be empty
* The `dist` parameter can be:
  * Any user defined function accepting two vectors and returning a scalar
  * An anonymous function with the same characteristics (e.g. `dist = (x,y) -> norm(x-y)^2`)
  * One of the above predefined distances: `l1_distance`, `l2_distance`, `l2²_distance`, `cosine_distance`

# Example:
```julia
julia> (clIdx,Z) = kmeans([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.3 38; 5.1 -2.3; 5.2 -2.4],3)
```
"""
function kmeans(X,K;dist=(x,y) -> norm(x-y),initStrategy="grid",Z₀=nothing,rng = Random.GLOBAL_RNG)
    X  = makeMatrix(X)
    (N,D) = size(X)
    # Random choice of initial representative vectors (any point, not just in X!)
    minX = minimum(X,dims=1)
    maxX = maximum(X,dims=1)
    Z₀ = initRepresentatives(X,K,initStrategy=initStrategy,Z₀=Z₀,rng=rng)
    Z  = Z₀
    cIdx_prev = zeros(Int64,N)

    # Looping
    while true
        # Determining the constituency of each cluster
        cIdx      = zeros(Int64,N)
        for (i,x) in enumerate(eachrow(X))
            cost = Inf
            for (j,z) in enumerate(eachrow(Z))
               if (dist(x,z)  < cost)
                   cost    =  dist(x,z)
                   cIdx[i] = j
               end
            end
        end

        # Determining the new representative by each cluster
        # for (j,z) in enumerate(eachrow(Z))
        for j in  1:K
            Cⱼ = X[cIdx .== j,:] # Selecting the constituency by boolean selection
            #debug  =  sum(Cⱼ,dims=1)
            #debug2 = size(Cⱼ)[1]
            if size(Cⱼ)[1] > 0
                Z[j,:] = sum(Cⱼ,dims=1) ./ size(Cⱼ)[1]
            else
                # move toward the center if no costituency
                xAvg = mean(X,dims=1)'
                Z[j,:] = Z[j,:] .+ ((xAvg - Z[j,:]) .* 0.01)
                #debug  = sum(Cⱼ,dims=1)
                # debug2 = size(Cⱼ)[1]
            end
            #Z[j,:] = median(Cⱼ,dims=1) # for l1 distance
        end

        # Checking termination condition: clusters didn't move any more
        if cIdx == cIdx_prev
            return (cIdx,Z)
        else
            cIdx_prev = cIdx
        end

    end
end

## Basic K-Medoids Algorithm (Lecture/segment 14.3 of https://www.edx.org/course/machine-learning-with-python-from-linear-models-to)
"""
  kmedoids(X,K;dist,initStrategy,Z₀)

Compute K-Medoids algorithm to identify K clusters of X using distance definition `dist`

# Parameters:
* `X`: a (n x d) data to clusterise
* `K`: Number of cluster wonted
* `dist`: Function to employ as distance (see notes). Default to Euclidean distance.
* `initStrategy`: Whether to select the initial representative vectors:
  * `random`: randomly in the X space
  * `grid`: using a grid approach
  * `shuffle`: selecting randomly within the available points [default]
  * `given`: using a provided set of initial representatives provided in the `Z₀` parameter
 * `Z₀`: Provided (K x D) matrix of initial representatives (used only together with the `given` initStrategy) [default: `nothing`]
 * `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Returns:
* A tuple of two items, the first one being a vector of size N of ids of the clusters associated to each point and the second one the (K x D) matrix of representatives

# Notes:
* Some returned clusters could be empty
* The `dist` parameter can be:
  * Any user defined function accepting two vectors and returning a scalar
  * An anonymous function with the same characteristics (e.g. `dist = (x,y) -> norm(x-y)^2`)
  * One of the above predefined distances: `l1_distance`, `l2_distance`, `l2²_distance`, `cosine_distance`

# Example:
```julia
julia> (clIdx,Z) = kmedoids([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.3 38; 5.1 -2.3; 5.2 -2.4],3,initStrategy="grid")
```
"""
function kmedoids(X,K;dist=(x,y) -> norm(x-y),initStrategy="grid",Z₀=nothing,rng = Random.GLOBAL_RNG)
    X  = makeMatrix(X)
    (n,d) = size(X)
    # Random choice of initial representative vectors
    Z₀ = initRepresentatives(X,K,initStrategy=initStrategy,Z₀=Z₀,rng=rng)
    Z = Z₀
    cIdx_prev = zeros(Int64,n)

    # Looping
    while true
        # Determining the constituency of each cluster
        cIdx      = zeros(Int64,n)
        for (i,x) in enumerate(eachrow(X))
            cost = Inf
            for (j,z) in enumerate(eachrow(Z))
               if (dist(x,z) < cost)
                   cost =  dist(x,z)
                   cIdx[i] = j
               end
            end
        end

        # Determining the new representative by each cluster (within the points member)
        #for (j,z) in enumerate(eachrow(Z))
        for j in  1:K
            Cⱼ = X[cIdx .== j,:] # Selecting the constituency by boolean selection
            nⱼ = size(Cⱼ)[1]     # Size of the cluster
            if nⱼ == 0 continue end # empty continuency. Let's not do anything. Stil in the next batch other representatives could move away and points could enter this cluster
            bestCost = Inf
            bestCIdx = 0
            for cIdx in 1:nⱼ      # candidate index
                 candidateCost = 0.0
                 for tIdx in 1:nⱼ # target index
                     candidateCost += dist(Cⱼ[cIdx,:],Cⱼ[tIdx,:])
                 end
                 if candidateCost < bestCost
                     bestCost = candidateCost
                     bestCIdx = cIdx
                 end
            end
            Z[j,:] = reshape(Cⱼ[bestCIdx,:],1,d)
        end

        # Checking termination condition: clusters didn't move any more
        if cIdx == cIdx_prev
            return (cIdx,Z)
        else
            cIdx_prev = cIdx
        end
    end

end

## The gmm algorithm (Lecture/segment 16.5 of https://www.edx.org/course/machine-learning-with-python-from-linear-models-to)

# no longer true with the numerical trick implemented
# - For mixtures with full covariance matrix (i.e. `FullGaussian(μ,σ²)`) the minCovariance should NOT be set equal to the minVariance, or if the covariance matrix goes too low, it will become singular and not invertible.
"""
  gmm(X,K;p₀,mixtures,tol,verbosity,minVariance,minCovariance,initStrategy)

Compute Expectation-Maximisation algorithm to identify K clusters of X data, i.e. employ a Generative Mixture Model as the underlying probabilistic model.

X can contain missing values in some or all of its dimensions. In such case the learning is done only with the available data.
Implemented in the log-domain for better numerical accuracy with many dimensions.

# Parameters:
* `X`  :           A (n x d) data to clusterise
* `K`  :           Number of cluster wanted
* `p₀` :           Initial probabilities of the categorical distribution (K x 1) [default: `nothing`]
* `mixtures`:      An array (of length K) of the mixture to employ (see notes) [def: `[DiagonalGaussian() for i in 1:K]`]
* `tol`:           Tolerance to stop the algorithm [default: 10^(-6)]
* `verbosity`:     A verbosity parameter regulating the information messages frequency [def: `STD`]
* `minVariance`:   Minimum variance for the mixtures [default: 0.05]
* `minCovariance`: Minimum covariance for the mixtures with full covariance matrix [default: 0]. This should be set different than minVariance (see notes).
* `initStrategy`:  Mixture initialisation algorithm [def: `kmeans`]
* `maxIter`:       Maximum number of iterations [def: `-1`, i.e. ∞]
* `rng`:           Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Returns:
* A named touple of:
  * `pₙₖ`:      Matrix of size (N x K) of the probabilities of each point i to belong to cluster j
  * `pₖ`:       Probabilities of the categorical distribution (K x 1)
  * `mixtures`: Vector (K x 1) of the estimated underlying distributions
  * `ϵ`:        Vector of the discrepancy (matrix norm) between pⱼₓ and the lagged pⱼₓ at each iteration
  * `lL`:       The log-likelihood (without considering the last mixture optimisation)
  * `BIC`:      The Bayesian Information Criterion (lower is better)
  * `AIC`:      The Akaike Information Criterion (lower is better)

 # Notes:
 - The mixtures currently implemented are `SphericalGaussian(μ,σ²)`,`DiagonalGaussian(μ,σ²)` and `FullGaussian(μ,σ²)`
 - Reasonable choices for the minVariance/Covariance depends on the mixture. For example 0.25 seems a reasonable value for the SphericalGaussian, 0.05 seems better for the DiagonalGaussian, and FullGaussian seems to prefer either very low values of variance/covariance (e.g. `(0.05,0.05)` ) or very big but similar ones (e.g. `(100,100)` ).
 - For `initStrategy`, look at the documentation of `initMixtures!` for the mixture you want. The provided gaussian mixtures support `grid`, `kmeans` or `given`. `grid` is faster (expecially if X contains missing values), but `kmeans` often provides better results.

 # Resources:
 - [Paper describing gmm with missing values](https://doi.org/10.1016/j.csda.2006.10.002)
 - [Class notes from MITx 6.86x (Sec 15.9)](https://stackedit.io/viewer#!url=https://github.com/sylvaticus/MITx_6.86x/raw/master/Unit 04 - Unsupervised Learning/Unit 04 - Unsupervised Learning.md)
 - [Limitations of gmm](https://www.r-craft.org/r-news/when-not-to-use-gaussian-mixture-model-gmm-clustering/)

# Example:
```julia
julia> clusters = gmm([1 10.5;1.5 0; 1.8 8; 1.7 15; 3.2 40; 0 0; 3.3 38; 0 -2.3; 5.2 -2.4],3,verbosity=HIGH)
```
"""
function gmm(X,K;p₀=nothing,mixtures=[DiagonalGaussian() for i in 1:K],tol=10^(-6),verbosity=STD,minVariance=0.05,minCovariance=0.0,initStrategy="kmeans",maxIter=-1,rng = Random.GLOBAL_RNG)
    if verbosity > STD
        @codeLocation
    end
    # debug:
    #X = [1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4]
    #K = 3
    #p₀=nothing; tol=0.0001; msgStep=1; minVariance=0.25; initStrategy="grid"
    #mixtures = [SphericalGaussian() for i in 1:K]
    # ---------
    X     = makeMatrix(X)
    (N,D) = size(X)
    pₖ    = isnothing(p₀) ? fill(1/K,K) : p₀

    # no longer true with the numerical trick implemented
    #if (minVariance == minCovariance)
    #    @warn("Setting the minVariance equal to the minCovariance may lead to singularity problems for mixtures with full covariance matrix.")
    #end

    msgStepMap = Dict(NONE => 0, LOW=>100, STD=>20, HIGH=>5, FULL=>1)
    msgStep    = msgStepMap[verbosity]

    # Initialisation of the parameters of the mixtures
    initMixtures!(mixtures,X,minVariance=minVariance,minCovariance=minCovariance,initStrategy=initStrategy,rng=rng)

    pₙₖ = zeros(Float64,N,K) # The posteriors, i.e. the prob that item n belong to cluster k
    ϵ = Float64[]

    # Checking dimensions only once (but adding then inbounds doesn't change anything. Still good
    # to provide a nice informative message)
    if size(pₖ,1) != K || length(mixtures) != K
        error("Error in the dimensions of the inputs. Please check them.")
    end

    # finding empty/non_empty values
    Xmask     =  .! ismissing.(X)
    #XdimCount = sum(Xmask, dims=2)

    lL = -Inf
    iter = 1
    while(true)
        oldlL = lL
        # E Step: assigning the posterior prob p(j|xi) and computing the log-Likelihood of the parameters given the set of data
        # (this last one for informative purposes and terminating the algorithm)
        pₙₖlagged = copy(pₙₖ)
        logpₙₖ = log.(pₙₖ)
        lL = 0
        for n in 1:N
            if any(Xmask[n,:]) # if at least one true
                Xu = X[n,Xmask[n,:]]
                #=if (length(ϵ) == 2)
                    println("here I am")
                    for m in mixtures[3:end]
                        println(m.μ)
                        println(m.σ²)
                        println(Xu)
                        println(Xmask[n,:])
                        lpdf(m,Xu,Xmask[n,:])
                        println("here I am partially")
                    end
                    println("here I am dead")
                end=#
                logpx = lse([log(pₖ[k] + 1e-16) + lpdf(mixtures[k],Xu,Xmask[n,:]) for k in 1:K])
                lL += logpx
                #px = sum([pⱼ[k]*normalFixedSd(Xu,μ[k,XMask[n,:]],σ²[k]) for k in 1:K])
                #println(n)
                for k in 1:K
                    logpₙₖ[n,k] = log(pₖ[k] + 1e-16)+lpdf(mixtures[k],Xu,Xmask[n,:])-logpx
                end
            else
                logpₙₖ[n,:] = log.(pₖ)
            end
        end
        pₙₖ = exp.(logpₙₖ)

        push!(ϵ,norm(pₙₖlagged - pₙₖ))

        # M step: find parameters that maximise the likelihood
        # Updating the probabilities of the different mixtures
        nₖ = sum(pₙₖ,dims=1)'
        n  = sum(nₖ)
        pₖ = nₖ ./ n
        updateParameters!(mixtures, X, pₙₖ; minVariance=minVariance,minCovariance=minCovariance)

        # Information. Note the likelihood is whitout accounting for the new mu, sigma
        if msgStep != 0 && (length(ϵ) % msgStep == 0 || length(ϵ) == 1)
            println("Iter. $(length(ϵ)):\tVar. of the post  $(ϵ[end]) \t  Log-likelihood $(lL)")
        end

        # Closing conditions. Note that the logLikelihood is those without considering the new mu,sigma
        if ((lL - oldlL) <= (tol * abs(lL))) || (maxIter > 0 && iter == maxIter)
            npars = npar(mixtures) + (K-1)
            #BIC  = lL - (1/2) * npars * log(N)
            BICv = bic(lL,npars,N)
            AICv = aic(lL,npars)
        #if (ϵ[end] < tol)
           return (pₙₖ=pₙₖ,pₖ=pₖ,mixtures=mixtures,ϵ=ϵ,lL=lL,BIC=BICv,AIC=AICv)
       else
            iter += 1
       end
    end # end while loop
end # end function

#using BenchmarkTools
#@benchmark clusters = emGMM([1 10.5;1.5 10.8; 1.8 8; 1.7 15; 3.2 40; 3.6 32; 3.3 38; 5.1 -2.3; 5.2 -2.4],3,msgStep=0)
#@benchmark clusters = emGMM([1 10.5;1.5 0; 1.8 8; 1.7 15; 3.2 40; 0 0; 3.3 38; 0 -2.3; 5.2 -2.4],3,msgStep=0,missingValue=0)
#@benchmark clusters = emGMM([1 10.5;1.5 0; 1.8 8; 1.7 15; 3.2 40; 0 0; 3.3 38; 0 -2.3; 5.2 -2.4],3,msgStep=0,missingValue=0)
#@benchmark clusters = emGMM([1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4],3,msgStep=0)
#@code_warntype gmm([1 10.5;1.5 0; 1.8 8; 1.7 15; 3.2 40; 0 0; 3.3 38; 0 -2.3; 5.2 -2.4],3,msgStep=0,missingValue=0)
#using Profile
#Juno.@profiler (for i = 1:1000 gmm([1 10.5;1.5 0; 1.8 8; 1.7 15; 3.2 40; 0 0; 3.3 38; 0 -2.3; 5.2 -2.4],3,msgStep=0,missingValue=0) end)
#Profile.clear()
#Profile.print()

#  - For mixtures with full covariance matrix (i.e. `FullGaussian(μ,σ²)`) the minCovariance should NOT be set equal to the minVariance, or if the covariance matrix goes too low, it will become singular and not invertible.
"""
  predictMissing(X,K;p₀,mixtures,tol,verbosity,minVariance,minCovariance)

Fill missing entries in a sparse matrix (i.e. perform a "matrix completion") assuming an underlying Gaussian Mixture probabilistic Model (GMM) and implementing
an Expectation-Maximisation algorithm.

While the name of the function is `predictMissing`, the function can be also used for system reccomendation / collaborative filtering and GMM-based regressions. The advantage over traditional algorithms as k-nearest neighbors (KNN) is that GMM can "detect" the hidden structure of the observed data, where some observation can be similar to a certain pool of other observvations for a certain characteristic, but similar to an other pool of observations for other characteristics.

Implemented in the log-domain for better numerical accuracy with many dimensions.

# Parameters:
* `X`  :           A (N x D) sparse matrix of data to fill according to a GMM model
* `K`  :           Number of mixtures (latent classes) to consider [def: 3]
* `p₀` :           Initial probabilities of the categorical distribution (K x 1) [default: `nothing`]
* `mixtures`:      An array (of length K) of the mixture to employ (see notes) [def: `[DiagonalGaussian() for i in 1:K]`]
* `tol`:           Tolerance to stop the algorithm [default: 10^(-6)]
* `verbosity`:     A verbosity parameter regulating the information messages frequency [def: `STD`]
* `minVariance`:   Minimum variance for the mixtures [default: 0.05]
* `minCovariance`: Minimum covariance for the mixtures with full covariance matrix [default: 0]. This should be set different than minVariance (see notes).
* `initStrategy`:  Mixture initialisation algorithm [def: `grid`]
* `maxIter`:       Maximum number of iterations [def: `-1`, i.e. ∞]
* `rng`:           Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Returns:
* A named touple of:
  * `̂X̂`    : The Filled Matrix of size (N x D)
  * `nFill`: The number of items filled
  * `lL`   : The log-likelihood (without considering the last mixture optimisation)
  * `BIC` :  The Bayesian Information Criterion (lower is better)
  * `AIC` :  The Akaike Information Criterion (lower is better)

  # Notes:
  - The mixtures currently implemented are `SphericalGaussian(μ,σ²)`,`DiagonalGaussian(μ,σ²)` and `FullGaussian(μ,σ²)`
  - For `initStrategy`, look at the documentation of `initMixtures!` for the mixture you want. The provided gaussian mixtures support `grid`, `kmeans` or `given`. `grid` is faster, but `kmeans` often provides better results.
  - The algorithm requires to specify a number of "latent classes" (mlixtures) to divide the dataset into. If there isn't any prior domain specific knowledge on this point one can test sevaral `k` and verify which one minimise the `BIC` or `AIC` criteria.

# Example:
```julia
julia>  cFOut = predictMissing([1 10.5;1.5 missing; 1.8 8; 1.7 15; 3.2 40; missing missing; 3.3 38; missing -2.3; 5.2 -2.4],3)
```
"""
function predictMissing(X,K=3;p₀=nothing,mixtures=[DiagonalGaussian() for i in 1:K],tol=10^(-6),verbosity=STD,minVariance=0.05,minCovariance=0.0,initStrategy="kmeans",maxIter=-1,rng = Random.GLOBAL_RNG)
    if verbosity > STD
        @codeLocation
    end
    emOut = gmm(X,K;p₀=p₀,mixtures=mixtures,tol=tol,verbosity=verbosity,minVariance=minVariance,minCovariance=minCovariance,initStrategy=initStrategy,maxIter=maxIter,rng=rng)
    (N,D) = size(X)
    nDim  = ndims(X)
    nmT   = nonmissingtype(eltype(X))
    #K = size(emOut.μ)[1]
    XMask = .! ismissing.(X)
    nFill = (N * D) - sum(XMask)
    X̂ = copy(X)
    for n in 1:N
        for d in 1:D
            if !XMask[n,d]
                 X̂[n,d] = sum([emOut.mixtures[k].μ[d] * emOut.pₙₖ[n,k] for k in 1:K])
            end
        end
    end
    X̂ = convert(Array{nmT,nDim},X̂)
    return (X̂=X̂,nFill=nFill,lL=emOut.lL,BIC=emOut.BIC,AIC=emOut.AIC)
end

# MLJ interface
include("Clustering_MLJ.jl")

end


text:
"""
Functions for analysing simulations

Output functions: Target the peak timing for each county, cases by each county,
timing of total peak and total cases
"""

"""
Define uncertainty of parameter estimates
"""
d_incubation = LogNormal(log(5.),0.25) #Liu et al
mean(d_incubation)
(quantile(d_incubation,0.025),median(d_incubation),quantile(d_incubation,0.975))
d_R₀ = Gamma(100,2.5/100) ##Liu et al
mean(d_R₀)
(quantile(d_R₀,0.025),median(d_R₀),quantile(d_R₀,0.975))

"""
Simulation functions
"""

"""
    function output_daily_and_final_incidence(sol,i)

This function stores the cumulative incidence and hospitalisation (summed over areas and age groups),
    and the final numbers by area and age group.
"""
function output_daily_and_final_incidence(sol,i)
    times = sol.prob.tspan[1]:1:sol.prob.tspan[end]
    z = [sum(sol(t)[:,:,9:12],dims = 2)[:,1,:]  for t in times]
    return (z,sol[end][:,:,9:12]),false
end

function output_daily_incidence_and_hosp(sol,i)
    times = sol.prob.tspan[1]:1:sol.prob.tspan[end]
    z = [sol(t)[:,:,9:12]  for t in times]
    return z,false
end

function randomise_params(prob,i,repeat) #Remember to rescale susceptibility by the inverse leading eigenvalue
    _P = deepcopy(prob.p)
    _P.isolating_detecteds = true
    _P.τ = _P.τ_initial
    _P.σ = 1/rand(d_incubation)
    _P.β = rand(d_R₀)*_P.γ
    return remake(prob,p=_P)
end

function consensus_randomise_params(prob,i,repeat) #Remember to rescale susceptibility by the inverse leading eigenvalue
    _P = deepcopy(prob.p)
    _P.β = rand(d_R₀)
    return remake(prob,p=_P)
end

function randomise_params_and_infectiousness(prob,i,repeat)
    _P = deepcopy(prob.p)
    _P.isolating_detecteds = true
    _P.τ = _P.τ_initial
    _P.σ = 1/rand(d_incubation)
    _P.β = rand(d_R₀)*_P.γ
    _P.ϵ = rand(Uniform(0.,0.5))
    return remake(prob,p=_P)
end

"""
Scenario functions
"""
# function run_simulations(P::KenyaCoV.CoVParameters_AS,prob,n_traj,τ,ϵ_D)
#     P.τ_initial = τ
#     P.ϵ_D = ϵ_D
#     ensemble_prob = EnsembleProblem(prob,
#                                     prob_func = randomise_params,
#                                     output_func = output_daily_and_final_incidence)
#     return solve(ensemble_prob,FunctionMap(),dt = P.dt,trajectories = n_traj)
# end

function run_simulations(P::KenyaCoV.CoVParameters_AS,prob,n_traj,τ,ϵ_D,cb)
    P.τ_initial = τ
    P.ϵ_D = ϵ_D
    ensemble_prob = EnsembleProblem(prob,
                                    prob_func = randomise_params,
                                    output_func = output_daily_and_final_incidence)
    return solve(ensemble_prob,FunctionMap(),dt = P.dt,callback = cb,trajectories = n_traj)
end

function run_consensus_simulations(P::KenyaCoV.CoVParameters_AS,prob,n_traj,cb)
    ensemble_prob = EnsembleProblem(prob,
                                    prob_func = consensus_randomise_params,
                                    output_func = output_daily_incidence_and_hosp)
    return solve(ensemble_prob,FunctionMap(),dt = P.dt,callback = cb,trajectories = n_traj)
end

# function run_simulations(P::KenyaCoV.CoVParameters_AS,prob,n_traj,τ,ϵ_D,prob_func)
#     P.τ_initial = τ
#     P.ϵ_D = ϵ_D
#     ensemble_prob = EnsembleProblem(prob,
#                                     prob_func = prob_func,
#                                     output_func = output_daily_and_final_incidence)
#     return solve(ensemble_prob,FunctionMap(),dt = P.dt,trajectories = n_traj)
# end

# function run_scenario(P::KenyaCoV.CoVParameters_AS,prob,n_traj,treatment_rates)
#     results = []
#     for (τ,ϵ_D) in treatment_rates
#         sims = run_simulations(P,prob,n_traj,τ,ϵ_D)
#         analysisdata = incidence_from_sims(sims)
#         push!(results,analysisdata)
#     end
#     return results
# end

function run_scenario(P::KenyaCoV.CoVParameters_AS,prob,n_traj,treatment_rates,cb)
    results = []
    for (τ,ϵ_D) in treatment_rates
        sims = run_simulations(P,prob,n_traj,τ,ϵ_D,cb)
        analysisdata = incidence_from_sims(sims)
        push!(results,analysisdata)
    end
    return results
end

# function run_scenario(P::KenyaCoV.CoVParameters_AS,prob,n_traj,treatment_rates,prob_func)
#     results = []
#     for (τ,ϵ_D) in treatment_rates
#         sims = run_simulations(P,prob,n_traj,τ,ϵ_D,prob_func)
#         analysisdata = incidence_from_sims(sims)
#         push!(results,analysisdata)
#     end
#     return results
# end

"""
Analysis functions
"""
function incidence_from_sims(sims)
    n = length(sims)
    m = length(sims[1][1])
    inc_A_data = zeros(n_wa+1,m-1,n)
    inc_D_data = zeros(n_wa+1,m-1,n)
    final_D = zeros(n_wa,n_a,n)
    final_A = zeros(n_wa,n_a,n)
    for (k,sim) in enumerate(sims)
        for i = 1:n_wa,t = 1:(m-1)
            inc_A_data[i,t,k] = sim[1][t+1][i,1] - sim[1][t][i,1]#Daily incidence rather than cumulative
            inc_D_data[i,t,k] = sim[1][t+1][i,2] - sim[1][t][i,2]
        end
        for t = 1:(m-1)
            inc_A_data[n_wa+1,t,k] = sum(sim[1][t+1][:,1] .- sim[1][t][:,1])
            inc_D_data[n_wa+1,t,k] = sum(sim[1][t+1][:,2] .- sim[1][t][:,2])
        end
        final_A[:,:,k] .= sim[2][:,:,1]
        final_D[:,:,k] .= sim[2][:,:,2]
    end
    peak_times = zeros(n,21)
    inc_A_conf_intvs = zeros(21,m-1,3)
    inc_D_conf_intvs = zeros(21,m-1,3)
    for i = 1:(n_wa+1),t = 1:(m-1)
        inc_A_conf_intvs[i,t,1] = quantile(inc_A_data[i,t,:],0.5)
        inc_A_conf_intvs[i,t,2] = inc_A_conf_intvs[i,t,1] - quantile(inc_A_data[i,t,:],0.025) #Lower conf. int.
        inc_A_conf_intvs[i,t,3] =  quantile(inc_A_data[i,t,:],0.975) - inc_A_conf_intvs[i,t,1] #Upper conf. int.
        inc_D_conf_intvs[i,t,1] = quantile(inc_D_data[i,t,:],0.5)
        inc_D_conf_intvs[i,t,2] = inc_D_conf_intvs[i,t,1] - quantile(inc_D_data[i,t,:],0.025) #Lower conf. int.
        inc_D_conf_intvs[i,t,3] =  quantile(inc_D_data[i,t,:],0.975) - inc_D_conf_intvs[i,t,1] #Upper conf. int.
    end

    for i = 1:(n_wa+1), k =1:n
        peak_times[k,i] = argmax(inc_A_data[i,:,k])
    end

    return inc_D_conf_intvs,inc_A_conf_intvs,peak_times,final_D,final_A
end

ext: jl


text:
module TestThresholdPredictors
using Test, MLJModels, CategoricalArrays
import MLJBase
using CategoricalDistributions

const MMI = MLJModels.MLJModelInterface

X = NamedTuple{(:x1,:x2,:x3)}((rand(4), rand(4), rand(4)))
yraw = ["in", "out", "out", "out"]
y = categorical(yraw, ordered = true) # `AbstractArray{OrderedFactor{2}, 1}`
y1 = categorical(yraw)                # `AbstractArray{Multiclass{2}, 1}
y2 = categorical(yraw[2:end], ordered=true)

@testset "BinaryThresholdPredictor - Probablistic" begin

    atom = ConstantClassifier()

    @test_throws MLJModels.ERR_MODEL_UNSPECIFIED BinaryThresholdPredictor()
    model = BinaryThresholdPredictor(atom)

    # Check warning when `y` is not ordered:
    @test_logs((:warn, MLJModels.warn_classes(levels(y)...)),
                MMI.fit(model, 1, X, y1))
    # Check predictions containing two classes
    @test_throws ArgumentError BinaryThresholdPredictor(ConstantRegressor())
    @test_logs((:warn, r"`threshold` should be"),
               BinaryThresholdPredictor(atom, threshold=-1))
    @test_logs((:warn, r"`threshold` should be"),
               BinaryThresholdPredictor(atom, threshold=1))

    # Compare fitresult and fitted_params with that of model, and
    # check no warning issued:
    model_fr, model_cache, model_report = @test_logs MMI.fit(model, 1, X, y)
    atom_fr, atom_cache, atom_report =
        MMI.fit(model.model, 1, X, y)
    @test model_fr[1] == atom_fr

    # Check model update
    model_up, model_cache_up, model_report_up =
        MMI.update(model, 1, model_fr, model_cache, X, y)
    atom_up, atom_cache_up, atom_report_up =
        MMI.update(model.model, 1, atom_fr, atom_cache, X, y)
    @test model_up[1] == atom_up
    @test model_cache_up[1] == atom_cache_up
    @test model_report_up[1] == atom_report_up

    # Check fitted_params
    @test MMI.fitted_params(model, model_fr).model_fitted_params ==
         MMI.fitted_params(model.model, atom_fr)

    # Check deterministic predictions
    @test MMI.predict(model, model_fr, X) ==
        MMI.predict_mode(model.model, atom_fr, X)

    model.threshold = 0.8
    model_fr, cache, report = MMI.fit(model, 1, X, y)
    @test MMI.predict(model, model_fr, X) ==
        [y[1] for i in 1:MMI.nrows(X)]

    d = MLJModels.info_dict(model)
    @test d[:supports_weights] == MMI.supports_weights(model.model)
    @test d[:input_scitype] == MMI.input_scitype(model.model)
    @test d[:target_scitype] == AbstractVector{<:MMI.Finite{2}}
    @test d[:is_pure_julia] == MMI.is_pure_julia(model.model)
    @test d[:name] == "BinaryThresholdPredictor"
    @test d[:load_path] == "MLJModels.BinaryThresholdPredictor"
end

@testset "_predict_threshold" begin
    v1 = categorical(['a', 'b', 'a'])
    v2 = categorical(['a', 'b', 'a', 'c'])
    # Test with UnivariateFinite object
    d1 = UnivariateFinite(MMI.classes(v1), [0.4, 0.6])
    @test_throws ArgumentError MLJModels._predict_threshold(d1, 0.7)
    @test MLJModels._predict_threshold(d1, (0.7, 0.3)) == v1[2]
    @test MLJModels._predict_threshold(d1, [0.5, 0.5]) == v1[2]
    @test MLJModels._predict_threshold(d1, (0.4, 0.6)) == v1[1]
    @test MLJModels._predict_threshold(d1, [0.2, 0.8]) == v1[1]
    d2 = UnivariateFinite(MMI.classes(v2), [0.4, 0.3, 0.3])
    @test_throws ArgumentError MLJModels._predict_threshold(d2, (0.7, 0.3))
    @test MLJModels._predict_threshold(d2, (0.2, 0.5, 0.3)) == v2[1]
    @test MLJModels._predict_threshold(d2, [0.3, 0.2, 0.5]) == v2[2]
    @test MLJModels._predict_threshold(d2, (0.4, 0.4, 0.2)) == v2[4]
    @test MLJModels._predict_threshold(d2, [0.2, 0.5, 0.3]) == v2[1]

    # Test with Array{UnivariateFinite, 1} object
    d1_arr = [d1 for i in 1:3]
    @test_throws ArgumentError MLJModels._predict_threshold(d1_arr, 0.7)
    @test MLJModels._predict_threshold(d1_arr, (0.7, 0.3)) == [v1[2] for i in 1:3]
    @test MLJModels._predict_threshold(d1_arr, [0.5, 0.5]) == [v1[2] for i in 1:3]
    @test MLJModels._predict_threshold(d1_arr, (0.4, 0.6)) == [v1[1] for i in 1:3]
    @test MLJModels._predict_threshold(d1_arr, [0.2, 0.8]) == [v1[1] for i in 1:3]
    d2_arr = [d2 for i in 1:3]
    @test_throws ArgumentError MLJModels._predict_threshold(d2_arr, (0.7, 0.3))
    @test MLJModels._predict_threshold(d2_arr, (0.2, 0.5, 0.3)) == [v2[1] for i in 1:3]
    @test MLJModels._predict_threshold(d2_arr, [0.3, 0.2, 0.5]) == [v2[2] for i in 1:3]
    @test MLJModels._predict_threshold(d2_arr, (0.4, 0.4, 0.2)) == [v2[4] for i in 1:3]
    @test MLJModels._predict_threshold(d2_arr, [0.2, 0.5, 0.3]) == [v2[1] for i in 1:3]

    # Test with UnivariateFiniteArray oject
    probs1 = [0.2 0.8; 0.7 0.3; 0.1 0.9]
    unf_arr1 = UnivariateFinite(MMI.classes(v1), probs1)
    @test_throws ArgumentError MLJModels._predict_threshold(unf_arr1, 0.7)
    @test MLJModels._predict_threshold(unf_arr1, (0.7, 0.3)) == [v1[2], v1[1], v1[2]]
    @test MLJModels._predict_threshold(unf_arr1, [0.5, 0.5]) == [v1[2], v1[1], v1[2]]
    @test MLJModels._predict_threshold(unf_arr1, (0.4, 0.6)) == [v1[2], v1[1], v1[2]]
    @test MLJModels._predict_threshold(unf_arr1, [0.2, 0.8]) == [v1[1], v1[1], v1[2]]
    probs2 = [0.2 0.3 0.5;0.1 0.6 0.3; 0.4 0.0 0.6]
    unf_arr2 = UnivariateFinite(MMI.classes(v2), probs2)
    @test_throws ArgumentError MLJModels._predict_threshold(unf_arr2, (0.7, 0.3))
    @test MLJModels._predict_threshold(unf_arr2, (0.2, 0.5, 0.3)) == [v2[4], v2[2], v2[1]]
    @test MLJModels._predict_threshold(unf_arr2, [0.3, 0.2, 0.5]) == [v2[2], v2[2], v2[1]]
    @test MLJModels._predict_threshold(unf_arr2, (0.4, 0.4, 0.2)) == [v2[4], v2[2], v2[4]]
    @test MLJModels._predict_threshold(unf_arr2, [0.2, 0.5, 0.3]) == [v2[4], v2[2], v2[1]]
end

# dummy detector always predicts outliers and inliers with equal
# probability:
struct DummyDetector <: MMI.ProbabilisticUnsupervisedDetector end
MMI.fit(::DummyDetector, verbosity, X) = nothing, nothing, nothing
MMI.predict(::DummyDetector, verbosity, X) =
    MLJBase.UnivariateFinite(["in", "out"],
                             fill(0.5, MLJBase.nrows(X)),
                             augment=true, pool=missing)
MMI.input_scitype(::Type{<:DummyDetector}) = MMI.Table

@testset "BinaryThresholdPredictor - ProbabilisticUnsupervisedDetector" begin
    detector = BinaryThresholdPredictor(DummyDetector(), threshold=0.2)
    @test_throws MLJModels.ERR_CLASSES_DETECTOR MMI.fit(detector, 1, X, y1)

    fr, _, _ = MMI.fit(detector, 0, X, y)
    @test MMI.predict(detector, fr, X) == fill("out", length(y))
    fr, _, _ = MMI.fit(detector, 0, X)
    @test MMI.predict(detector, fr, X) == fill("out", length(y))

    detector.threshold = 0.8
    fr, _, _ = MMI.fit(detector, 0, X, y)
    @test MMI.predict(detector, fr, X) == fill("in", length(y))
    fr, _, _ = MMI.fit(detector, 0, X)
    @test MMI.predict(detector, fr, X) == fill("in", length(y))

    # integration (y == ["in", "out", "out", "out"]):
    e = MLJBase.evaluate(detector, X, y,
                         resampling=MLJBase.Holdout(fraction_train=0.5),
                         measure=MLJBase.accuracy)
    @test e.measurement[1] ≈ 0
end

end
true


text:
"""
    Hypergeometric(s, f, n)

A *Hypergeometric distribution* describes the number of successes in `n` draws without replacement from a finite population containing `s` successes and `f` failures.

```math
P(X = k) = {{{s \\choose k} {f \\choose {n-k}}}\\over {s+f \\choose n}}, \\quad \\text{for } k = \\max(0, n - f), \\ldots, \\min(n, s).
```

```julia
Hypergeometric(s, f, n)  # Hypergeometric distribution for a population with
                         # s successes and f failures, and a sequence of n trials.

params(d)       # Get the parameters, i.e. (s, f, n)
```

External links

* [Hypergeometric distribution on Wikipedia](http://en.wikipedia.org/wiki/Hypergeometric_distribution)

"""
struct Hypergeometric <: DiscreteUnivariateDistribution
    ns::Int     # number of successes in population
    nf::Int     # number of failures in population
    n::Int      # sample size

    function Hypergeometric(ns::Real, nf::Real, n::Real)
        @check_args(Hypergeometric, ns >= zero(ns) && nf >= zero(nf))
        @check_args(Hypergeometric, zero(n) <= n <= ns + nf)
        new(ns, nf, n)
    end
end

@distr_support Hypergeometric max(d.n - d.nf, 0) min(d.ns, d.n)

### Parameters

params(d::Hypergeometric) = (d.ns, d.nf, d.n)

### Statistics

mean(d::Hypergeometric) = d.n * d.ns / (d.ns + d.nf)

function var(d::Hypergeometric)
    N = d.ns + d.nf
    p = d.ns / N
    d.n * p * (1.0 - p) * (N - d.n) / (N - 1.0)
end
mode(d::Hypergeometric) = floor(Int, (d.n + 1) * (d.ns + 1) / (d.ns + d.nf + 2))

function modes(d::Hypergeometric)
    if (d.ns == d.nf) && mod(d.n, 2) == 1
        [(d.n-1)/2, (d.n+1)/2]
    else
        [mode(d)]
    end
end

skewness(d::Hypergeometric) = (d.nf-d.ns)*sqrt(d.ns+d.nf-1)*(d.ns+d.nf-2*d.n)/sqrt(d.n*d.ns*d.nf*(d.ns+d.nf-d.n))/(d.ns+d.nf-2)
function kurtosis(d::Hypergeometric)
    N = d.ns + d.nf
    a = (N-1) * N^2 * (N * (N+1) - 6*d.ns * (N-d.ns) - 6*d.n*(N-d.n)) + 6*d.n*d.ns*(d.nf)*(N-d.n)*(5*N-6)
    b = (d.n*d.ns*(N-d.ns) * (N-d.n)*(N-2)*(N-3))
    a/b
end

### Evaluation & Sampling

@_delegate_statsfuns Hypergeometric hyper ns nf n

rand(d::Hypergeometric) = convert(Int, StatsFuns.RFunctions.hyperrand(d.ns, d.nf, d.n))

struct RecursiveHypergeomProbEvaluator <: RecursiveProbabilityEvaluator
    ns::Float64
    nf::Float64
    n::Float64
end

RecursiveHypergeomProbEvaluator(d::Hypergeometric) = RecursiveHypergeomProbEvaluator(d.ns, d.nf, d.n)

nextpdf(s::RecursiveHypergeomProbEvaluator, p::Float64, x::Integer) =
    ((s.ns - x + 1) / x) * ((s.n - x + 1) / (s.nf - s.n + x)) * p

Base.broadcast!(::typeof(pdf), r::AbstractArray, d::Hypergeometric, rgn::UnitRange) =
    _pdf!(r, d, rgn, RecursiveHypergeomProbEvaluator(d))

function Base.broadcast(::typeof(pdf), d::Hypergeometric, X::UnitRange)
    r = similar(Array{promote_type(partype(d), eltype(X))}, indices(X))
    r .= pdf.(d,X)
end


text:
function run(env = ENV,
             cicfg::CIService=auto_detect_ci_service(;env=env);
             merge_new_packages::Bool,
             merge_new_versions::Bool,
             new_package_waiting_period,
             new_version_waiting_period,
             registry::String,
             authorized_authors::Vector{String},
             #
             master_branch::String = "master",
             master_branch_is_default_branch::Bool = true,
             suggest_onepointzero::Bool = true)

    registry_head = directory_of_cloned_registry(cicfg; env=env)

    # Run tests on the registry (these are very quick)
    RegistryCI.test(registry_head)

    # Figure out what type of build this is
    run_pr_build = conditions_met_for_pr_build(cicfg; env=env, master_branch=master_branch)
    run_merge_build = conditions_met_for_merge_build(cicfg; env=env, master_branch=master_branch)

    if !(run_pr_build || run_merge_build)
        @info "Build not determined to be either a PR build or a merge build. Exiting."
        return nothing
    end

    # Authentication
    auth = my_retry(() -> GitHub.authenticate(env["AUTOMERGE_GITHUB_TOKEN"]))
    whoami = my_retry(() -> username(cicfg; auth=auth))
    @info("Authenticated to GitHub as \"$(whoami)\"")
    registry_repo = my_retry(() -> GitHub.repo(registry; auth=auth))

    if run_pr_build
        pr_number = pull_request_number(cicfg; env=env)
        pr_head_commit_sha = current_pr_head_commit_sha(cicfg; env=env)
        pull_request_build(pr_number,
                           pr_head_commit_sha,
                           registry_repo,
                           registry_head;
                           auth = auth,
                           authorized_authors = authorized_authors,
                           master_branch = master_branch,
                           master_branch_is_default_branch = master_branch_is_default_branch,
                           suggest_onepointzero = suggest_onepointzero,
                           whoami = whoami)
        return nothing
    elseif run_merge_build
        cron_or_api_build(registry_repo;
                          auth = auth,
                          authorized_authors = authorized_authors,
                          merge_new_packages = merge_new_packages,
                          merge_new_versions = merge_new_versions,
                          new_package_waiting_period = new_package_waiting_period,
                          new_version_waiting_period = new_version_waiting_period,
                          whoami = whoami)
        return nothing
    else
        error("This should not happen.")
    end
end


text:
"""
Evaluate value of the graph with a given root node
"""
function evaluate(policy::POMDPs.Policy, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP)
    N = 5000                    # TODO: pompd/solver config
    n = policy.updater.root
    v = 0.0
    v = @distributed (+) for i in 1:N
        simulate(sim, pomdp, policy, policy.updater, n)
    end
    v /= N
    return v
end

"""
Evaluate a batch of states
"""
function evaluate(sts::Vector{S}, policy::POMDPs.Policy, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP, n::MCVINode) where {S}
    # vs = pmap((s)->simulate(sim, pomdp, policy, policy.updater, n, s), sts)
    vs = zeros(Float64, length(sts))
    for (i,s) in enumerate(sts)
        vs[i] = simulate(sim, pomdp, policy, policy.updater, n, s)
    end
    return vs
end

"""
Evaluate belief
"""
function evaluate(b::MCVIBelief, policy::POMDPs.Policy, sim::MCVISimulator, pomdp::POMDPs.POMDP, n::MCVINode, num_eval_belief::Int64)
    val::Float64 = 0.0
    val = @distributed (+) for i in 1:num_eval_belief
        s = rand(sim.rng, b)    # This is the initial state stuff :/
        # sim.init_state = s      # TODO maybe roll this into simulate as well?
        simulate(sim, pomdp, policy, policy.updater, n, s)
    end
    # sim.init_state = nothing    # Back to being nothing
    return val/num_eval_belief
end

"""
Alpha vectors
"""
function compute(sts::Vector{S}, policy::POMDPs.Policy, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP{S,A,O}, n::MCVINode) where {S,A,O}
    a = evaluate(sts, policy, sim, pomdp, n)
    edge = AlphaEdge(a, n.id)
    return edge
end

function _fill_obs!(obs::Vector{O}, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP{S,A,O}, ba::MCVIBelief) where {S,A,O}
    for i in 1:length(obs)
        # sample observation from belief
        s = rand(sim.rng, ba)
        obs[i] = gen(DDNNode(:o), pomdp, s, sim.rng)
    end
end

function _fill_ov!(ov::Vector{Float64}, osum::Vector{Float64}, obs::Vector{O},
                   policy::POMDPs.Policy, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP,
                   n::MCVINode, ba::MCVIBelief, num_states::Int64) where {O}
    sts = [rand(sim.rng, ba) for _ in 1:num_states]
    v = evaluate(sts, policy, sim, pomdp, n)
    for i in 1:num_states
        for (j,o) in enumerate(obs)
            wt = pdf(sts[i], o)
            osum[j] += wt
            ov[j] += v[i]*wt
        end
    end
    # normalize
    for (i,p) in enumerate(osum)
        ov[i] /= p
    end
end

function _fill_X!(X::Array{Float64,2}, obs::Vector{O}, sts::Vector{S}) where {O,S}
    for (j,s) in enumerate(sts)
        for (i,o) in enumerate(obs)
            X[i,j] = pdf(s, o)
        end
    end
end

mutable struct Scratch{O}
    obs::Vector{O}
    ov::Vector{Float64}
    osum::Vector{Float64}
    X::Array{Float64,2}
end

"""
Least square computation
"""
function compute(sts::Vector{S}, policy::POMDPs.Policy, sim::POMDPs.Simulator, pomdp::POMDPs.POMDP{S,A,O}, n::MCVINode, ba::MCVIBelief, scratch::Scratch) where {S,A,O}
    _fill_obs!(scratch.obs, sim, pomdp, ba)
    _fill_ov!(scratch.ov, scratch.osum, scratch.obs, policy, sim, pomdp, n, ba, size(scratch.X,2))
    _fill_X!(scratch.X, scratch.obs, sts)

    b = scratch.X \ scratch.ov                  # Least squares

    @assert length(b) == length(sts)
    @assert !isnan(b[1])

    edge = AlphaEdge(b, n.id)
    return edge
end


text:
using Pkg.Artifacts

artifact_toml = joinpath(@__DIR__, "Artifacts.toml")
iris_hash = artifact_hash("iris", artifact_toml)

if iris_hash == nothing || !artifact_exists(iris_hash)
     iris_hash = create_artifact() do artifact_dir
        iris_url_base = "https://archive.ics.uci.edu/ml/machine-learning-databases/iris"
        download("$(iris_url_base)/iris.data",
                 joinpath(artifact_dir, "iris.csv"))
        download("$(iris_url_base)/bezdekIris.data",
                 joinpath(artifact_dir, "bezdekIris.csv"))
        download("$(iris_url_base)/iris.names",
                 joinpath(artifact_dir, "iris.names"))
    end

    bind_artifact!(artifact_toml, "iris", iris_hash)
end

iris_dataset_path = artifact_path(iris_hash)


text:
# grids.jl

using OffsetArrays
using CUDA
using KernelAbstractions
using CUDAKernels

abstract type AbstractArchitecture end
struct CPUArch <: AbstractArchitecture end
struct GPUArch <: AbstractArchitecture end

device(::CPUArch) = KernelAbstractions.CPU()
device(::GPUArch) = CUDAKernels.CUDADevice()

mutable struct Grid{AT,FT}
    # grid_type::String
    nx::Int64
    ny::Int64
    nz::Int64
    Δx::FT
    Δy::FT
    Δz::FT
    x::AT
    y::AT
    z::AT
    arch::AbstractArchitecture
end

function gridCoords(nx, ny, nz, Δx, Δy, Δz, useCUDA, useOffset, myFloat)
    # Create 3D structured grid
    # For now cell centered

    x = initField1D(0.0, nx, useCUDA, useOffset, myFloat)
    y = initField1D(0.0, ny, useCUDA, useOffset, myFloat)
    z = initField1D(0.0, nz, useCUDA, useOffset, myFloat)

    # Now fill the grid values
    for n = 0:nx + 1
        x[n] = myFloat(n - 0.5) * Δx
    end

    for n = 0:ny + 1
        y[n] = myFloat(n - 0.5) * Δy
    end

    for n = 0:nz + 1
        z[n] = myFloat(n - 0.5) * Δz
    end

    return x, y, z
end


text:
function newpoints_for_particular_k!(new_points, new_impoints, split_inds)
    # Fill the array by looping over all simplices that we want to split
    for i = 1:length(inds_toolarge)
        # Figure out what the row indices corresponding to the ith simplex
        # must be. Marks the beginning of each of the simplex stacks in new_points
        k = string(ks[i])
        n_newverts_persplit = size(rules[k], 1)
        ind = n_newverts_persplit * (i - 1)

        # Index of the simplex we need to split
        simplex_idx = inds_toolarge[i]

        # Get the vertices of the simplex currently being splitted. Each of the
        # n_newverts_persplit new vertices will be a linear combination
        # of these vertices. Each row is a vertex.
        ith_simpl_vertices = points[simplex_inds[simplex_idx, :], :]
        ith_simpl_imvertices = impoints[simplex_inds[simplex_idx, :], :]

        # Create new vertices within parent simplex as linear combinations of its original
        # vertices such that each generated new vertex lies strictly within the simplex
        for j = 1:n_newverts_persplit
            # Rules for constructing the jth subsimplex (i. e. indices indicating
            # which of the vertices of the parent simplex to combine into the jth
            # subsimplex; could be for example [1, 1, 2, 4], meaning that the new
            # vertex will be the following linear combination of the original vertices:
            # (2*V1 + 1*V2 + 1*V4)/(splitting factor).
            jth_simpl_rules = rules[k][j, :]

            original_points = ith_simpl_vertices[jth_simpl_rules, :]
            original_impoints = ith_simpl_imvertices[jth_simpl_rules, :]

            append!(new_points, sum(original_points, dims=1) ./ ks[i])
            append!(new_impoints, sum(original_impoints, dims=1) ./ ks[i])
        end
    end

    new_points = reshape(new_points, ÷(length(new_points), dim), dim)
    new_impoints = reshape(new_impoints, ÷(length(new_points), dim), dim)
end


text:
using FileIO
import Knet

# BUGS:
# After optimization arrays in x.opt (adam statistics) get saved as well.
# RNNs have r.c, r.h etc arrays that are not Params.
# If we do Params only we miss BatchNorm parameters.

# How to save/load weights in a way that is robust to code change (except the order of weights ;)

function getweights(model; atype=Knet.atype())
    # Can't just to Params: there are non-Param weights in BatchNorm, should we call these Const?
    # On the other hand we don't want to save param.opt variables, or do it optionally
    # We could make bn mean/var Params with lr=0, but an optimizer may change that. Add a freeze flag to Param?
    w = []
    deepmap(atype, model) do a
        push!(w, convert(Array,a))
    end
    return w
end

function setweights!(model, weights; atype=Knet.atype())
    # The trouble is a newly initialized model does not have weights until first input
    n = 0
    deepmap(atype, model) do w
        n += 1
        @assert size(w) == size(weights[n]) "$(size(w)) != $(size(weights[n]))"
        copyto!(w, weights[n])
    end
    @assert n == length(weights)  "The number of weights do not match: $n != $(length(weights))"
    return model
end

function setweights!(model, file::String; atype=Knet.atype())
    setweights!(model, load(file, "weights"); atype)
end

function saveweights(file::String, model; atype=Knet.atype())
    save(file, "weights", getweights(model; atype))
end


text:
using Test
using DataStructure101
const DS = DataStructure101

@testset "Test Deque" begin

    # Initial a Deque with data '0'
    T_e = Int64
    T = DS.Deque{T_e}
    data= T_e[0]
    deque = T(0)
    @test deque.data == data

    #Type test:
    @test deque isa T

    # Push test: Push data into deque 20 times by 1:20
    passed_check_data = []
    for i = 1:20
        push!(deque, i)
        push!(data, i)
        push!(passed_check_data, deque.data == data)
    end
    @test all(passed_check_data)

    # Pushfirst test: Push data into deque 20 times by 31:35
    for i = 31:35
        pushfirst!(deque, i)
        pushfirst!(data, i)
        push!(passed_check_data, deque.data == data)
    end
    @test all(passed_check_data)

    # Pop test: Pop data from tail of deque
    passed_pop = [
        pop!(deque) == pop!(data),
        deque.data == data
    ]
    @test all(passed_pop)

    # Pop test: Pop data from head of deque
    passed_pop = [
        popfirst!(deque) == popfirst!(data),
        deque.data == data
    ]
    @test all(passed_pop)

    # Show Test
    @test repr(deque) == "Deque(["*
    "34, 33, 32, 31, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19"*
    "])"

end


text:
using SymFEL
using Test
using SymPy
using LinearAlgebra

# Testing Lagrange elements
println("# Testing Lagrange elements")
p1 = SymFEL.get_lagrange_basis()
@test SymPy.degree(p1[1], gen=x) == 1
@test SymPy.degree(p1[2], gen=x) == 1
@test SymPy.subs(p1[1], x, 0) == 1
@test SymPy.subs(p1[1], x, h) == 0
@test SymPy.subs(p1[2], x, h) == 1
@test SymPy.subs(p1[2], x, 0) == 0

p4 = SymFEL.get_lagrange_basis(4)
E = Matrix(1I, 5, 5)
for i = 1:4
    @test SymPy.degree(p4[i], gen=x) == 4
    for j = 1:5
        @test SymPy.subs(p4[i], x, (j - 1) * h / 4) == E[i, j]
    end
end

M = SymFEL.get_lagrange_em()
@test (M[1, 1] == h / 3) && ((M[1, 2] == h / 6)) &&
    (M[2, 1] == h / 6) && ((M[2, 2] == h / 3))

MVC = SymFEL.get_lagrange_em_varcoeff()
MVC = MVC.subs([(xa, 0), (xb, h)])
@test (MVC[1, 1] == h / 3) && (MVC[1, 2] == h / 6) &&
    (MVC[2, 1] == h / 6) && (MVC[2, 2] == h / 3)

# Testing Hermite elements
println("# Testing Hermite elements")
p3 = SymFEL.get_hermite_basis()
@test SymPy.subs(p3[1], x, 0) == 1
@test SymPy.subs(p3[1], x, h) == 0
@test SymPy.subs(SymPy.diff(p3[1], x), x, 0) == 0
@test SymPy.subs(SymPy.diff(p3[1], x), x, h) == 0

@test SymPy.subs(p3[2], x, 0) == 0
@test SymPy.subs(p3[2], x, h) == 0
@test SymPy.subs(SymPy.diff(p3[2], x), x, 0) == 1
@test SymPy.subs(SymPy.diff(p3[2], x), x, h) == 0

@test SymPy.subs(p3[3], x, 0) == 0
@test SymPy.subs(p3[3], x, h) == 1
@test SymPy.subs(SymPy.diff(p3[3], x), x, 0) == 0
@test SymPy.subs(SymPy.diff(p3[3], x), x, h) == 0

@test SymPy.subs(p3[4], x, 0) == 0
@test SymPy.subs(p3[4], x, h) == 0
@test SymPy.subs(SymPy.diff(p3[4], x), x, 0) == 0
@test SymPy.subs(SymPy.diff(p3[4], x), x, h) == 1

MH = SymFEL.get_hermite_em()
MH_hc = [13*h/35 11*h^2/210 9*h/70 -13*h^2/420;
         11*h^2/210 h^3/105 13*h^2/420 -h^3/140;
         9*h/70 13*h^2/420 13*h/35 -11*h^2/210;
         -13*h^2/420 -h^3/140 -11*h^2/210 h^3/105]
@test MH == MH_hc

MHVC = SymFEL.get_hermite_em_varcoeff()
MHVC = MHVC.subs([(xa, 0), (xb, h)])
@test MHVC == MH_hc

try
    SymFEL.get_hermite_basis(2)
catch y
    println(y)
end

##
t = range(0, stop=1,length=11)
t = convert(Array{Float64}, t)
ti = range(0, stop=1, length=101)
xv = sin.(pi * t')
yv = pi * cos.(pi * t')
fdf = [xv; yv]
f = SymFEL.interpolate(fdf, t, ti)
@test norm(f[1:10:101]' - xv) <= 1e-15

## Testing get_em() function.
println("Testing get_em() function.")
@test SymFEL.get_em(1, 1, 0, 0) == SymFEL.get_lagrange_em(1, 0, 0)
@test SymFEL.get_em(3, 3, 1, 0; fe1 = "Hermite", fe2 = "Hermite") == SymFEL.get_hermite_em(3, 1, 0)
try
    SymFEL.get_em(1, 1, 0, 0; fe1 = "a", fe2 = "Lagrange")
catch er
    println(er)
end

try
    SymFEL.get_em(1, 1, 2, 0)
catch er
    println(er)
end

try
    SymFEL.get_em(1, 1, 1, 2)
catch er
    println(er)
end

try
    SymFEL.get_em(3, 3, 3, 0; fe1 = "Hermite", fe2="Hermite")
catch er
    println(er)
end

try
    SymFEL.get_em(3, 3, 2, -1; fe1 = "Hermite", fe2="Hermite")
catch er
    println(er)
end

##
n = 1001
xv = range(0, stop=1, length=2*n-1).^2
dx = 1 / (n - 1)

elem_M = SymFEL.get_lagrange_em(2, 0, 0)
elem_M = elem_M.subs(h, dx)
elem_M = convert(Matrix{Float64}, elem_M)
M = SymFEL.assemble_1d_FE_matrix(elem_M, n, intNodes1 = 1, intNodes2 = 1, dof1 = 1, dof2 = 1)
@test abs(convert(Float64, (xv' * M * xv)[1]) - 1 / 5) < 1e-15

elem_M_p1 = SymFEL.get_lagrange_em(1, 0, 0)
elem_M_p1 = elem_M_p1.subs(h, dx)
elem_M_p1 = convert(Matrix{Float64}, elem_M_p1)
M_p1 = SymFEL.assemble_1d_FE_matrix(elem_M_p1, n)
xv = range(0, stop=1, length=n).^2
@test abs(convert(Float64, (xv' * M_p1 * xv)[1]) - 1 / 5) < 1e-6

elem_M_h3 = SymFEL.get_hermite_em(3, 0, 0)
elem_M_h3 = elem_M_h3.subs(h, dx)
elem_M_h3 = convert(Matrix{Float64}, elem_M_h3)
M_h3 = SymFEL.assemble_1d_FE_matrix(elem_M_h3, n, intNodes1 = 0, intNodes2 = 0, dof1 = 2, dof2 = 2)
xv = zeros(2 * n, 1)
xv[1:2:(2*(n-1)+1)] = convert(Array{Float64}, range(0, stop=1, length=n).^3)
xv[2:2:(2*(n-1)+2)] = 3*convert(Array{Float64}, range(0, stop=1, length=n).^2)
@test abs(convert(Float64, (xv' * M_h3 * xv)[1]) - 1 / 7) < 1e-15

nodes = convert(Array{Float64, 1}, range(0, stop=1, length=n))
elem_M_p1 = SymFEL.get_lagrange_em(1, 0, 0)
M_p1_nu = SymFEL.assemble_1d_nu_FE_matrix(elem_M_p1, nodes; intNodes1 = 0, intNodes2 = 0, dof1 = 1, dof2 = 1)
A = M_p1 - M_p1_nu
@test norm(A[:]) < 1e-15
##

nodes = [convert(Array{Float64}, range(0, stop=0.99, length=100));
         convert(Array{Float64}, range(0.99, stop=1, length=101))]
M_p1_nu2 = SymFEL.assemble_1d_nu_FE_matrix(elem_M_p1, nodes; intNodes1 = 0, intNodes2=0, dof1 = 1, dof2 = 1)
xv = nodes.^2
@test abs((xv' * M_p1_nu2 * xv)[1] - 1 / 5) < 1e-4

# lagrange - hermite
n = 1001
xv = range(0, stop=1, length=2*n-1).^4
dx = 1 / (n - 1)

elem_M = SymFEL.get_em(2, 3, 1, 0, fe1="Lagrange", fe2="Hermite")
elem_M = elem_M.subs(h, dx)
elem_M = convert(Matrix{Float64}, elem_M)
M = SymFEL.assemble_1d_FE_matrix(elem_M, n, intNodes1 = 1, intNodes2 = 0, dof1 = 1, dof2 = 2)

xr = zeros(2*n)
xvr2 = range(0, stop=1, length=n)
i = 1:n
xr[2*(i .- 1) .+ 1] = xvr2
xr[2*(i .- 1) .+ 2] = ones(n)

@test abs((xv' * M * xr)[1] - 4 / 5) < 1e-14
println("All tests passed.")


text:
#=
This file is auto-generated. Do not edit.
=#

#! format: off

"""
    mutable struct AVRSimple <: AVR
        Kv::Float64
        V_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        states_types::Vector{StateTypes}
        internal::InfrastructureSystemsInternal
    end

Parameters of a simple proportional AVR in the derivative of EMF
i.e. an integrator controller on EMF

# Arguments
- `Kv::Float64`: Proportional Gain, validation range: `(0, nothing)`
- `V_ref::Float64`: Reference Voltage Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states are:
	Vf: field voltage
- `n_states::Int`: Fixed AVR has 1 states
- `states_types::Vector{StateTypes}`: Simple AVR has 1 differential states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRSimple <: AVR
    "Proportional Gain"
    Kv::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    ext::Dict{String, Any}
    "The states are:
	Vf: field voltage"
    states::Vector{Symbol}
    "Fixed AVR has 1 states"
    n_states::Int
    "Simple AVR has 1 differential states"
    states_types::Vector{StateTypes}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRSimple(Kv, V_ref=1.0, ext=Dict{String, Any}(), )
    AVRSimple(Kv, V_ref, ext, [:Vf], 1, [StateTypes.Differential], InfrastructureSystemsInternal(), )
end

function AVRSimple(; Kv, V_ref=1.0, ext=Dict{String, Any}(), states=[:Vf], n_states=1, states_types=[StateTypes.Differential], internal=InfrastructureSystemsInternal(), )
    AVRSimple(Kv, V_ref, ext, states, n_states, states_types, internal, )
end

# Constructor for demo purposes; non-functional.
function AVRSimple(::Nothing)
    AVRSimple(;
        Kv=0,
        V_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`AVRSimple`](@ref) `Kv`."""
get_Kv(value::AVRSimple) = value.Kv
"""Get [`AVRSimple`](@ref) `V_ref`."""
get_V_ref(value::AVRSimple) = value.V_ref
"""Get [`AVRSimple`](@ref) `ext`."""
get_ext(value::AVRSimple) = value.ext
"""Get [`AVRSimple`](@ref) `states`."""
get_states(value::AVRSimple) = value.states
"""Get [`AVRSimple`](@ref) `n_states`."""
get_n_states(value::AVRSimple) = value.n_states
"""Get [`AVRSimple`](@ref) `states_types`."""
get_states_types(value::AVRSimple) = value.states_types
"""Get [`AVRSimple`](@ref) `internal`."""
get_internal(value::AVRSimple) = value.internal

"""Set [`AVRSimple`](@ref) `Kv`."""
set_Kv!(value::AVRSimple, val) = value.Kv = val
"""Set [`AVRSimple`](@ref) `V_ref`."""
set_V_ref!(value::AVRSimple, val) = value.V_ref = val
"""Set [`AVRSimple`](@ref) `ext`."""
set_ext!(value::AVRSimple, val) = value.ext = val
"""Set [`AVRSimple`](@ref) `states_types`."""
set_states_types!(value::AVRSimple, val) = value.states_types = val


text:
using MccDaqHats
using Test

@testset "MccDaqHats.jl" begin
    # Write your tests here.
end

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
using Zygote, LinearAlgebra
using IterativeSolvers: lobpcg

function test_eigmin(A::AbstractMatrix)
    results = lobpcg(A, false, 1)
    ev = minimum(results.λ)
    return ev
end

A = rand(4, 4); A = A + A';

@show eigmin(A)
@show test_eigmin(A)
@show Zygote.gradient(test_eigmin, A)

max_stars_repo_name: ericphanson/Attempt-diff-eigmin

max_issues_repo_name: ericphanson/Attempt-diff-eigmin

max_forks_repo_name: ericphanson/Attempt-diff-eigmin


text:
function plot_results(;paperplots=false)
    if paperplots
        close("all")
        figure(3, figsize=(4,3))
    end

    figure(3)
    clf()

    #N = 1./hgrid
    N = alpha_list

    loglog(N, grad_max_rel_err, "*-", label=L"||\nabla \tilde u-\nabla u||_{\infty}")
    loglog(N, grad_l2_rel_err, "+-", label=L"||\nabla \tilde u-\nabla u||_{2}")
    loglog(N, max_rel_err, ".-", label=L"||\tilde u-u||_{\infty}")
    loglog(N, l2_rel_err, "x-", label=L"||\tilde u-u||_{2}")

    legend()
    xlabel(L"\alpha")
    grid("on")
    ylim(1e-15, 1)
    ylabel("Relative error")

    tight_layout(0)
end


text:
__precompile__()
module FinancialToolbox

using SpecialFunctions: erf
using Dates, Requires
include("dates.jl")
function __init__()
    @require DualNumbers = "fa6b7ba4-c1ee-5f82-b5fc-ecf0adba8f74" include("financial_dual.jl")
    @require HyperDualNumbers = "50ceba7f-c3ee-5a84-a6e8-3ad40456ec97" include("financial_hyper.jl")
end
include("financial.jl")
export normcdf,
    normpdf,
    blsprice,
    blkprice,
    blsdelta,
    blsgamma,
    blsvega,
    blsrho,
    blstheta,
    blslambda,
    blsimpv,
    blkimpv,
    ##	ADDITIONAL Function
    blspsi,
    blsvanna,
    #Dates
    fromExcelNumberToDate,
    daysact,
    yearfrac

end#End Module


text:
module Altro

import TrajectoryOptimization
import RobotDynamics
using StaticArrays
using BenchmarkTools
using Interpolations
using SolverLogging
using Crayons

using SparseArrays
using LinearAlgebra
using Logging
using Statistics
using TimerOutputs
using ForwardDiff
using FiniteDiff
import Octavian

const TO = TrajectoryOptimization
const RD = RobotDynamics

using TrajectoryOptimization:
    num_constraints, get_trajectory

import TrajectoryOptimization: rollout!, get_constraints, get_model, get_objective,
    evaluate_constraints!, constraint_jacobians!
import RobotDynamics: discrete_dynamics, dynamics, dynamics!, evaluate, evaluate!

using TrajectoryOptimization:
    Problem,
    ConstraintList,
    AbstractObjective, Objective, #QuadraticObjective,
    # SampledTrajectory,
    # DynamicsExpansion, # TODO: Move to ALTRO
    # ALConstraintSet,
    # DynamicsConstraint,
    states, controls,
    Equality, Inequality, SecondOrderCone,
    cost

using RobotDynamics:
    AbstractModel, DiscreteDynamics, DiscreteLieDynamics,
    QuadratureRule, Implicit, Explicit,
    FunctionSignature, InPlace, StaticReturn,
    DiffMethod, ForwardAD, FiniteDifference, UserDefined,
    AbstractKnotPoint, KnotPoint, StaticKnotPoint,
    state_dim, control_dim, output_dim, dims,
    state, control, SampledTrajectory

# types
export
    ALTROSolverOld,
    ALTROSolver,
    # iLQRSolverOld,
    # AugmentedLagrangianSolver,
    SolverStats,
    SolverOptions

export
    solve!,
    benchmark_solve!,
    iterations,
    set_options!,
    max_violation,
    status

# modules
export
    Problems

const ColonSlice = Base.Slice{Base.OneTo{Int}}
const SparseView{T,I} = SubArray{T, 2, SparseMatrixCSC{T, I}, Tuple{UnitRange{I}, UnitRange{I}}, false}
const VectorView{T,I} = SubArray{T, 1, Vector{T}, Tuple{UnitRange{I}}, true}

# Select the matix multiplication kernel
const USE_OCTAVIAN = parse(Bool, get(ENV, "ALTRO_USE_OCTAVIAN", "true"))
@static if USE_OCTAVIAN
    const matmul! = Octavian.matmul!
else
    const matmul! = mul!
end

# Include the QDLDL wrapper
include("qdldl.jl")
using .Cqdldl

# High level structs
include("utils.jl")
include("infeasible_model.jl")
include("solvers.jl")
include("solver_opts.jl")

# iLQR Solver
include("ilqr/cost_expansion.jl")
include("ilqr/dynamics_expansion.jl")
include("ilqr/ilqr_solver.jl")
include("ilqr/backwardpass.jl")
include("ilqr/forwardpass.jl")
include("ilqr/ilqr_solve.jl")

# Augmented Lagrangian Solver
include("augmented_lagrangian/alcon.jl")
include("augmented_lagrangian/alconset.jl")
include("augmented_lagrangian/al_objective.jl")
include("augmented_lagrangian/al_solver.jl")
include("augmented_lagrangian/al_solve.jl")
include("direct/sparseblocks.jl")

# Projected Newton Solver
include("direct/pncon.jl")
include("direct/pnconset.jl")
include("direct/pn_solver.jl")
include("direct/pn_solve.jl")

# ALTRO Solver
include("altro/altro_solver.jl")
include("altro/altro_solve.jl")

# Example problems submodule
include("problems.jl")

end # module


text:
abstract type InitialVelocityDistribution{N} end

"""
    (::InitialVelocityDistribution{N})(
        positions::AbstractVector{SVectorF{N}},
        mass_densities::AbstractVector{Float},
    ) -> Vector{SVectorF{N}}

Compute the fluid particle velocities corresponding to this initial velocity
distribution.
"""
function (::InitialVelocityDistribution) end

include("velocity/static.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
using MultivariateStats
using LinearAlgebra
using Test
import Statistics: mean, cov
import Random
import SparseArrays

Random.seed!(34568)

## PCA with zero mean

X = randn(5, 10)
Y = randn(3, 10)

P = qr(randn(5, 5)).Q[:, 1:3]
pvars = [5., 4., 3.]
M = PCA(Float64[], P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == zeros(5)
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test transform(M, X[:,1]) ≈ P'X[:,1]
@test transform(M, X) ≈ P'X

@test reconstruct(M, Y[:,1]) ≈ P * Y[:,1]
@test reconstruct(M, Y) ≈ P * Y

## PCA with non-zero mean

mval = rand(5)
M = PCA(mval, P, pvars, 15.0)

@test indim(M) == 5
@test outdim(M) == 3
@test mean(M) == mval
@test projection(M) == P
@test principalvars(M) == pvars
@test principalvar(M, 2) == pvars[2]
@test tvar(M) == 15.0
@test tprincipalvar(M) == 12.0
@test tresidualvar(M) == 3.0
@test principalratio(M) == 0.8

@test transform(M, X[:,1]) ≈ P' * (X[:,1] .- mval)
@test transform(M, X) ≈ P' * (X .- mval)

@test reconstruct(M, Y[:,1]) ≈ P * Y[:,1] .+ mval
@test reconstruct(M, Y) ≈ P * Y .+ mval

## prepare training data

d = 5
n = 1000

R = collect(qr(randn(d, d)).Q)
@test R'R ≈ Matrix(I, 5, 5)
rmul!(R, Diagonal(sqrt.([0.5, 0.3, 0.1, 0.05, 0.05])))

X = R'randn(5, n) .+ randn(5)
mval = vec(mean(X, dims=2))
Z = X .- mval

C = cov(X, dims=2)
pvs0 = sort(eigvals(C); rev=true)
tv = sum(diag(C))

## pcacov (default using cov when d < n)

M = fit(PCA, X)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mval
@test P'P ≈ Matrix(I, 5, 5)
@test C*P ≈ P*Diagonal(pvs)
@test issorted(pvs; rev=true)
@test pvs ≈ pvs0
@test tvar(M) ≈ tv
@test sum(pvs) ≈ tvar(M)
@test reconstruct(M, transform(M, X)) ≈ X

M = fit(PCA, X; mean=mval)
@test projection(M) ≈ P

M = fit(PCA, Z; mean=0)
@test projection(M) ≈ P

M = fit(PCA, X; maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ Matrix(I, 3, 3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ Matrix(I, 3, 3)
@test issorted(pvs; rev=true)

## pcastd

M = fit(PCA, X; method=:svd)
P = projection(M)
pvs = principalvars(M)

@test indim(M) == 5
@test outdim(M) == 5
@test mean(M) == mval
@test P'P ≈ Matrix(I, 5, 5)
@test isapprox(C*P, P*Diagonal(pvs), atol=1.0e-3)
@test issorted(pvs; rev=true)
@test isapprox(pvs, pvs0, atol=1.0e-3)
@test isapprox(tvar(M), tv, atol=1.0e-3)
@test sum(pvs) ≈ tvar(M)
@test reconstruct(M, transform(M, X)) ≈ X

M = fit(PCA, X; method=:svd, mean=mval)
@test projection(M) ≈ P

M = fit(PCA, Z; method=:svd, mean=0)
@test projection(M) ≈ P

M = fit(PCA, X; method=:svd, maxoutdim=3)
P = projection(M)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ Matrix(I, 3, 3)
@test issorted(pvs; rev=true)

M = fit(PCA, X; method=:svd, pratio=0.85)

@test indim(M) == 5
@test outdim(M) == 3
@test P'P ≈ Matrix(I, 3, 3)
@test issorted(pvs; rev=true)

# Different data types
# --------------------
# test that fit works with Float32 values
X2 = convert(Array{Float32,2}, X)
# Float32 input, default pratio
M = fit(PCA, X2; maxoutdim=3)
# Float32 input, specified Float64 pratio
M = fit(PCA, X2, pratio=0.85)
# Float32 input, specified Float32 pratio
M = fit(PCA, X2, pratio=0.85f0)
# Float64 input, specified Float32 pratio
M = fit(PCA, X, pratio=0.85f0)

# views
M = fit(PCA, view(X, :, 1:500), pratio=0.85)
# sparse
@test_throws AssertionError fit(PCA, SparseArrays.sprandn(100d, n, 0.6))


text:
"""

Structure containing the options and parameters of the simulation.

### Examples
```julia-repl
julia> opt = Options()

julia> opt = Options(nsteps=2000)

```

"""
@with_kw struct Options
  # Time
  dt::Float64 = 0.05
  nsteps::Int = 2000

  # System parameters
  eps::Float64 = 1.
  sig::Float64 = 2.

  # Temperature bath parameters
  initial_velocities::Symbol = :normal
  kT::Float64 = 0.6
  ibath::Int = 10
  iequil::Int = 10
  tau::Int = 10
  lambda::Float64 = 0.1
  alpha::Float64 = 0.1

  # Output and printing options
  printxyz::Bool = true
  iprintxyz::Int = 10
  iprint::Int = div(nsteps,100)
  minprint::Bool = false
  trajectory_file::String = "traj.xyz"
  velocity_file::String = "vel.xyz"

end


text:
function myfunction(n)
    for i = 1:n
        A = randn(100,100,20)
        m = maximum(A)
        Am = mapslices(sum, A; dims=2)
        B = A[:,:,5]
        Bsort = mapslices(sort, B; dims=1)
        b = rand(100)
        C = B.*b
    end
end

# # Do the following in the Juno IDE
# using Profile
# Profile.clear()
# myfunction(1) # trigger compilation
# @profile myfunction(10)
# # up to here, everything is the same as for using ProfileView.jl

# # Now comes the cool part. Visualize your profile results:
# Juno.profiler()
# # or simply
# @profiler myfunction(10)

max_stars_repo_path: 2 Tips and tricks/2_profiling/myprogram.jl

max_stars_repo_name: crstnbr/julia-workshop


text:
# diversity / dissimilarity index calculation and plots:

using JLD2, DataFrames, CSV, Dates, Query
using LinearAlgebra, Clustering
using MAT
#using StatsPlots
using Gadfly

include("/home/kristen/Documents/coda_utilities/coda_ops_transforms.jl")

#load in the data:
@load "/home/kristen/Documents/V6V8_analysis/analysis_products/oligo_df_zero_plus.jld2"

## compute distance of each sample to all others, and plot by season:

y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)])) #rows are samples, columns are components
m=size(y,1)
oligo_dist = zeros(Float64,m,m)#Array{Float64,2}(undef,m,m)

for i=1:m
    for j=i:m
       oligo_dist[i,j] = ait_dist(y[i,:],y[j,:])
       oligo_dist[j,i] = ait_dist(y[i,:],y[j,:])
    end
end

## hmmm...perhaps, easier to organize by yearday?

ind=sortperm(dayofyear.(oligo_df_plus[!,:date])) #rev=true) #returns indicies...
# dayofyear.(oligo_df_plus[ind,:date]) #to check....
y_yrdy = y[ind,:]
oligo_byday = zeros(Float64,m,m)#Array{Float64,2}(undef,m,m)

for i=1:m
    for j=i:m
       oligo_byday[i,j] = ait_dist(y_yrdy[i,:],y_yrdy[j,:])
       oligo_byday[j,i] = ait_dist(y_yrdy[i,:],y_yrdy[j,:])
    end
end

## NEED TO WORK ON COLOR CODING, ETC...
file = matopen("/home/kristen/Documents/V6V8_analysis/paper_V6V8/matlab_figure_scripts/exports_from_julia/dissimilarity_matrix.mat", "w")
write(file, "oligo_dissim", oligo_dist)
write(file, "oligo_dissim_yrdy", oligo_byday)
write(file,"ind", ind)
write(file,"sample_yearday", dayofyear.(oligo_df_plus[!,:date]))
close(file)

## easier in a dataframe for plotting?
y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)])) #rows are samples, columns are components
m=size(y,1)
df_dist = DataFrame(x=Int64[], y=Int64[], dist=Float64[], date=DateTime[])

for i=1:m
    for j=i:m
        push!(df_dist,[i,j,ait_dist(y[i,:],y[j,:]),oligo_df_plus[i,:date]])
    end
end

## #yes - definitely easier for plotting! Very cool!

df_dist.yearday=dayofyear.(df_dist[!,:date])
plot(sort(df_dist,:yearday),x="x",y="y",color="dist", Geom.rectbin) #oops, no - this gives the same thing!

#CSV.write("/home/kristen/Documents/V6V8_analysis/paper_V6V8/matlab_figure_scripts/exports_from_julia/etc...", clade_array, delim = ',')

#what about a dendrogram?
#visualize! but must be with Plots.jl!
hc = hclust(Symmetric(oligo_dist), linkage=:single, uplo=:U)
#-> a few other linkage options, check hclust help

w=plot(hc,xticks=false,linewidth=2)
#node_labels=["O1-I","O2-CB5","O3-II","O4-III/IV","O5-I","O6-I*"] #oligotypes 1-6
#node_labels[hc.order]
#xticks!(w,([1,2,3,4,5,6],node_labels[hc.order]),fontsize=16) #maps from hc and matching up oligos to clades
#ylabel!(w,"Aitchision distance")
#savefig(w, "/home/kristen/Documents/V6V8_analysis/paper_V6V8/V6V8_Syn_CoDa_paper/figures/dendrogram_single.pdf")

## very cool! now, let's check out the aitchsion square norm for inqueality:

y=closure(convert(Array,oligo_df_plus[!,Between(:csO1,:csO6)])) #rows are samples, columns are components
m=size(y,1)
oligo_ineq = zeros(Float64,m,1)#Array{Float64,2}(undef,m,m)

for i=1:m
    oligo_ineq[i] = square_ait_norm(y[i,:])
end

## export it!
file = matopen("/home/kristen/Documents/V6V8_analysis/paper_V6V8/matlab_figure_scripts/exports_from_julia/ineq_metric.mat", "w")
write(file, "oligo_ineq", oligo_ineq)
rr=1 .- exp.(-oligo_ineq)
write(file,"scaled_ineq", rr) #should already ahve the smaples dates somewhere else...
close(file)

## plot it up!

Gadfly.set_default_plot_size(30cm,10cm)

xticks=DateTime("2010-01-01"):Year(1):DateTime("2019-01-01")
xyearticks = Dates.value.(xticks)
xyearticks_label=map(x-> Dates.format(x,"yyyy"),xticks)
labels = Dict(zip(xyearticks, xyearticks_label))

fontsize=Theme(minor_label_font_size=16pt,major_label_font_size=16pt,plot_padding=[5mm, 5mm, 2mm, 2mm])

l1=layer(x=Dates.value.(oligo_df_plus[!,:date]),y=oligo_ineq,Geom.line,Geom.point)
p2=plot(l1,Guide.xticks(ticks=xyearticks),Guide.xlabel(""),Scale.x_continuous(labels = x -> labels[x]),fontsize)

## can also scale it from 0-1:

rr=1 .- exp.(-oligo_ineq)
l2=layer(x=Dates.value.(oligo_df_plus[!,:date]),y=rr,Geom.line,Geom.point)
p2=plot(l2,Guide.xticks(ticks=xyearticks),Guide.xlabel(""),Scale.x_continuous(labels = x -> labels[x]),fontsize)

        #Guide.yticks(ticks=[0 0.25 0.5 0.75 1])
        # Coord.Cartesian(ymin=0, ymax=1),
        # Guide.ylabel("Relative abundance"),fontsize,
        # Guide.annotation(compose(context(), text(Dates.value.(DateTime(2010,1,10)), 0.9, "B"))))


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã ``2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License. You
# may obtain a copy of the License at
# http://www.apache.org/licenses/LICENSE-2.0. Unless required by applicable
# law or agreed to in writing, software distributed under the License is
# distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
# KIND, either express or implied. See the License for the specific language
# governing permissions and limitations under the License.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

"""
Derivation tree (aka parse tree or syntax tree) based on GrammaticalEvolution.jl grammar rules.
Warning: not all rules are supported
"""
module DerivationTrees

export DerivTreeParams, DerivationTree, DerivTreeNode, DecisionRule, get_expr,
    maxlength, get_children, get_sym
export initialize!, actionspace, iscomplete, isleaf, expand_node!, is_decision,
    is_terminal, swap_children!, max_depth, rm_tree!, rm_node, count_leafs,
    count_nonleafs
export IncompleteException

using RLESUtils, MemPools, TreeIterators, TreeUtils
using Reexport
@reexport using GrammaticalEvolution
using AbstractTrees
using JLD

import AbstractTrees: children, printnode
import Base: length, copy!, string, convert

typealias DecisionRule Union{OrRule, RangeRule, RepeatedRule} #rules that require a decision

type DerivTreeParams
    grammar::Grammar
end

type DerivTreeNode
    cmd::String
    rule::Union{Rule,Void}
    action::Int64
    depth::Int64
    children::Vector{DerivTreeNode}
end
function DerivTreeNode(rule::Union{Rule,Void}=nothing, depth::Int64=0,
                       cmd::AbstractString="", action::Int64=-1)
    return DerivTreeNode(cmd, rule, action, depth, DerivTreeNode[])
end

const NODEPOOL = MemPool(DerivTreeNode,10000,1000000)

type DerivationTree
    params::DerivTreeParams
    root::DerivTreeNode
    nopen::Int64
    maxactions::Int64 #maximum size of actionspace
    nodepool::MemPool
end

type DerivationTreeSerial
    params::DerivTreeParams
    root::DerivTreeNode
    nopen::Int64
    maxactions::Int64 #maximum size of actionspace
end

function DerivationTree(p::DerivTreeParams, nodepool::MemPool=NODEPOOL)
    root = DerivTreeNode(p.grammar.rules[:start])
    maxactions = maxlength(p.grammar)
    tree = DerivationTree(p, root, 0, maxactions, nodepool)
    tree
end

convert(::Type{DerivationTree}, ds::DerivationTreeSerial) =
    DerivationTree(ds.params, ds.root, ds.nopen, ds.maxactions, NODEPOOL)

function JLD.writeas(t::DerivationTree)
    DerivationTreeSerial(t.params, t.root, t.nopen, t.maxactions)
end

include("formatter.jl")

immutable IncompleteException <: Exception end

const EMPTYARRAY = DerivTreeNode[]

is_decision(node::DerivTreeNode) = isa(node.rule, DecisionRule)
is_terminal(node::DerivTreeNode) = isempty(node.children)

#reset tree, children are deallocated, root is kept
function initialize!(tree::DerivationTree)
    rm_node(tree, tree.root.children)
    initialize!(tree.root)
    tree.root.rule = tree.params.grammar.rules[:start]
    tree.nopen = 1
    tree.root
end

#reset node, children are emptied
function initialize!(node::DerivTreeNode)
    node.rule = nothing
    node.depth = 0
    node.cmd = ""
    node.action = -1
    empty!(node.children)
    node
end

function mk_node(tree::DerivationTree)
    node = checkout(tree.nodepool)
    initialize!(node)
    node
end

function rm_tree!(tree::DerivationTree)
    rm_node(tree, [tree.root])
end

function rm_node(tree::DerivationTree, nodes::Vector{DerivTreeNode})
    return_to_pool(tree.nodepool, nodes)
    nothing
end

#returns children to pool before emptying
function return_to_pool(nodepool::MemPool, nodes::Vector{DerivTreeNode})
    for node in nodes
        return_to_pool(nodepool, node)
    end
    empty!(nodes)
end

function return_to_pool(nodepool::MemPool, node::DerivTreeNode)
    return_to_pool(nodepool, node.children) #return children first
    checkin(nodepool, node) #return self
end

function iscomplete(tree::DerivationTree)
    @assert tree.nopen >= 0 #sanity check
    tree.nopen == 0
end

expand_node!(tree::DerivationTree, node::DerivTreeNode) = expand_node!(tree, node, node.rule)
expand_node!(tree::DerivationTree, node::DerivTreeNode, a::Int64) =
    expand_node!(tree, node, node.rule, a)

# Open nodes are allocated (and returned), but not populated until they are processed
###########################
### expand_node! nonterminals
function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::OrRule, a::Int64)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    node.action = a
    node.cmd = rule.name
    idx = ((a - 1) % length(rule.values)) + 1 #1-indexed
    child = mk_node(tree)
    child.rule = rule.values[idx]
    child.depth = node.depth + 1
    push!(node.children, child)
    tree.nopen += length(node.children)
    node.children
end

function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::ReferencedRule)
  #don't create a child node for reference rules, shortcut through
    node.rule = tree.params.grammar.rules[rule.symbol]
    DerivTreeNode[node]
end

function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::RepeatedRule, a::Int64)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    node.cmd = rule.name
    reps = ((a - 1) % length(rule.range)) + rule.range.start
    for i = 1:reps
        child = mk_node(tree)
        child.rule = rule.value
        child.depth = node.depth + 1
        push!(node.children, child)
    end
    tree.nopen += length(node.children)
    node.children
end

function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::AndRule)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    node.cmd = rule.name
    for subrule in rule.values
        child = mk_node(tree)
        child.rule = subrule
        child.depth = node.depth + 1
        push!(node.children, child)
    end
    tree.nopen += length(node.children)
    node.children
end

function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::ExprRule)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    node.cmd = rule.name
    for arg in rule.args
        if isa(arg, Rule)
            child = mk_node(tree)
            child.rule = arg
            child.depth = node.depth + 1
            push!(node.children, child)
        end
    end
    tree.nopen += length(node.children)
    node.children
end

###########################
### Terminals
function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::RangeRule, a::Int64)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    node.action = a
    EMPTYARRAY #use a pre-allocated empty array
end

function expand_node!(tree::DerivationTree, node::DerivTreeNode, rule::Terminal)
    tree.nopen -= 1
    @assert tree.nopen >= 0 #sanity check
    EMPTYARRAY #use a pre-allocated empty array
end

###########################
### get_expr
function get_expr(tree::DerivationTree) #entry
    iscomplete(tree) ? get_expr(tree.root) : throw(IncompleteException())
end

get_expr(node::DerivTreeNode) = get_expr(node, node.rule)
get_expr(node::DerivTreeNode, rule::Terminal) = rule.value

function get_expr(node::DerivTreeNode, rule::RangeRule)
    value = ((node.action - 1) % length(rule.range)) + rule.range.start
    if rule.action != nothing
        value = rule.action(value)
    end
    value
end

function get_expr(node::DerivTreeNode, rule::OrRule)
    child = node.children[1]
    value = get_expr(child, child.rule)
    if rule.action != nothing
        value = rule.action(value)
    end
    value
end

function get_expr(node::DerivTreeNode, rule::RepeatedRule)
    values = Any[]
    for child in node.children
        push!(values, get_expr(child, child.rule))
    end
    if rule.action !== nothing
        values = rule.action(values)
    end
    values
end

function get_expr(node::DerivTreeNode, rule::AndRule)
    values = Any[]
    for child in node.children
        push!(values, get_expr(child, child.rule))
    end
    if rule.action !== nothing
        values = rule.action(values)
    end
    values
end

function get_expr(node::DerivTreeNode, rule::ExprRule)
    xs = Any[]
    child_i = 1
    for arg in rule.args
        if isa(arg, Rule)
            child = node.children[child_i]
            push!(xs, get_expr(child, child.rule))
            child_i += 1
        else
            push!(xs, arg)
        end
    end
    Expr(xs...)
end

###########################
###
"""
Range of valid action numbers at next open node
"""
actionspace(node::DerivTreeNode) = actionspace(node, node.rule)
actionspace(node::DerivTreeNode, rule::OrRule) = 1:length(rule.values)
actionspace(node::DerivTreeNode, rule::RangeRule) = 1:length(rule.range)
actionspace(node::DerivTreeNode, rule::RepeatedRule) = 1:length(rule.range)

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

#entry
"""
Returns the highest length of actionspace of a grammar
"""
maxlength(grammar::Grammar) = reduce(max, 0, map(length, values(grammar.rules)))

#decision rules
Base.length(rule::OrRule) = length(rule.values)
Base.length(rule::RangeRule) = length(rule.range)
#other rules
Base.length(rule::Rule) = 1

function copy!(dst::DerivationTree, src::DerivationTree)
    initialize!(dst)
    dst.params = src.params
    copy!(dst.root, src.root, dst.nodepool)
    dst.nopen = src.nopen
    dst.maxactions = src.maxactions
    dst.nodepool = src.nodepool
end

function copy!(dst::DerivTreeNode, src::DerivTreeNode, nodepool::MemPool{DerivTreeNode})
    dst.rule = src.rule
    dst.depth = src.depth
    dst.cmd = src.cmd
    dst.action = src.action
    for srcchild in src.children
        dstchild = checkout(nodepool)
        copy!(dstchild, srcchild, nodepool)
        push!(dst.children, dstchild)
    end
end

isleaf(node::DerivTreeNode) = isempty(node.children)
get_children(node::DerivTreeNode) = node.children

function swap_children!(x::DerivTreeNode, y::DerivTreeNode)
    ch = x.children
    x.children = y.children
    y.children = ch
end

max_depth(tree::DerivationTree) = max_depth(tree.root)
max_depth(node::DerivTreeNode) = traverse(x->x.depth, (x,y)->max(x,y), node)

count_leafs(tree::DerivationTree) = count_nodes(isleaf, tree.root)
count_nonleafs(tree::DerivationTree) = count_nodes(x->!isleaf(x), tree.root)
length(tree::DerivationTree) = count_nodes(tree.root)

get_sym(node::DerivTreeNode) = Symbol(node.rule.name)

########################################
# Interface for AbstractTrees.jl.  Enables print_tree(STDOUT, tree)
children(tree::DerivationTree) = tree.root.children
children(node::DerivTreeNode) = node.children
function printnode(io::IO, tree::DerivationTree)
    Base.print_with_color(:blue, io, string(tree.root))
end
function printnode(io::IO, node::DerivTreeNode)
    print(io, string(node))
end
function string(node::DerivTreeNode)
    !isempty(node.cmd) ? node.cmd : string(node, node.rule)
end
string(node::DerivTreeNode, rule::RangeRule) = get_expr(node, rule)
string(node::DerivTreeNode, rule::Terminal) = get_expr(node, rule)
string(node::DerivTreeNode, rule) = string(typeof(rule))

end #module

max_stars_repo_name: sisl/ExprSearch

max_issues_repo_name: sisl/ExprSearch


text:
# Generate 3 Gaussian peak at position : 5, 10, 20 but transform true
# X into 1.1 * X + 0.2 before saving
#
using DelimitedFiles,Random
using NLS_Fit

Random.seed!(1234)

model = Gaussian_Peak() + Gaussian_Peak() + Gaussian_Peak()
θ1 = Float64[1,5,1]
θ2 = Float64[1.5,10,2]
θ3 = Float64[0.75,20,3]
θ = vcat(θ1,θ2,θ3)

X=Float64[1:0.25:30;]
n=length(X)
Y=eval_y(model,X,θ) + 0.1*(rand(n) .- 0.5)

@. X = 1.1*X + 0.2 # inverse map is: 0.91 * X - 0.18

rootDir = joinpath(dirname(pathof(NLS_Fit)), "..")
dataDir = joinpath(rootDir,"data")
dataFile = joinpath(dataDir,"recalibration.txt")
writedlm(dataFile,hcat(X,Y))


text:
@assert begin
    msh=FinVolMesh{2}()

    newptind1=add_point!(msh, [0.0, 0.0])
    newptind2=add_point!(msh, [1.0, 0.0])
    newptind3=add_point!(msh, [0.0, 1.0])

    @assert newptind1==1
    @assert newptind2==2
    @assert newptind3==3

    newfaceind1=add_face!(msh, [newptind1, newptind2])
    newfaceind2=add_face!(msh, [newptind2, newptind3])
    newfaceind3=add_face!(msh, [newptind3, newptind1])

    @assert newfaceind1==1
    @assert newfaceind2==2
    @assert newfaceind3==3

    newvol=add_volume!(msh, [newfaceind1, newfaceind2, newfaceind3])

    @assert newvol==1

    process_geometry!(msh)
    @assert isapprox(msh.volume_volumes[1], 0.5)

    mesh_smoothing!(msh)
    @assert isapprox(msh.volume_volumes[1], 0.0)

    true
end


text:
println("Db) Test of the PhotoExcitation module with ASF from an internally generated initial- and final-state multiplet.")
#
setDefaults("print summary: open", "zzz-PhotoExcitation.sum")
wa = Atomic.Computation("xx",  Nuclear.Model(36.);
                        initialConfigs=[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                        finalConfigs  =[Configuration("1s 2s^2"), Configuration("1s 2s 2p"), Configuration("1s 2p^2")],
                        process = JAC.PhotoExc,
                        processSettings=PhotoExcitation.Settings([E1, M1], [JAC.UseCoulomb, JAC.UseBabushkin], true, true, true, true,
                                                                 false, Tuple{Int64,Int64}[], 0., 0., 1.0e6, JAC.ExpStokes(0., 0., 0.) ) )

wb = perform(wa)
setDefaults("print summary: close", "")

#     [Configuration("[Ne] 3s^2 3p^5"), Configuration("[Ne] 3s 3p^6"), Configuration("[Ne] 3s^2 3p^4 3d")],


text:
"""
Drawings and calculations useful for research in the theory of Kleinian groups
and other groups of complex functions on the Riemann sphere.
"""
module SDDKleinianGroups

using Reexport

@reexport using SDDCore, SDDGeometry, SDDGraphics

include("traversegroup.jl")

include("orbits.jl")

export
    drawpointorbit,
    drawpointssetorbit,
    drawcircleorbit,
    drawcirclessetorbit

include("limitsets.jl")

export
    drawlimitset, drawΛ,
    drawfixedpointslimitset, drawfixedpointsΛ,
    drawchaosgamelimitset, drawchaosgameΛ,
#=    drawdensitylimitset, drawdensityΛ, =#
    drawtrappedpointslimitset, drawtrappedpointsΛ

#include("connectedlimitset.jl")

#=export
    drawconnectedlimitset, drawconnectedΛ,
=#

#include("fundamentalregions.jl")

#=export
=#

end


text:
module Day07
using AdventOfCode.Computer
export part1, part2

function foreach_perm(fun, n)
    v = collect(0:n-1)
    c = ones(Int, n)
    foreach_perm!(fun, v, c)
end

function foreach_perm!(fun, A::Vector, c::Vector)
    fun(A)

    n = length(c)
    i = 1
    @inbounds while i <= n
        if  c[i] < i
            if i%2 == 1
                (A[i], A[1]) = (A[1], A[i])
            else
                (A[i], A[c[i]]) = (A[c[i]], A[i])
            end
            c[i] += 1
            i = 1

            fun(A)
        else
            c[i] = 1
            i += 1
        end
    end
end

test_phases(prog, phases) = test_phases!(deepcopy(prog), prog, phases)

function test_phases!(prog_tmp, prog, phases)
    input = 0
    for phase in phases
        reset!(prog_tmp, prog)
        run_program!(prog_tmp, (phase, input))
        input = prog_tmp.output
    end
    input
end

function max_signal(prog)
    max = 0
    prog_tmp = deepcopy(prog)

    foreach_perm(5) do phases
        m = test_phases!(prog_tmp, prog, phases)
        if m > max
            max = m
        end
    end
    max
end

part1() = max_signal(read_program("input07.dat"))

function test_feedback(prog, phases)
    progs = Tuple(deepcopy(prog) for i in 1:5)
    test_feedback!(progs, phases)
end

function test_feedback!(progs, phases)
    for i in 1:5
        progs[i].input = phases[i]
        step_program!(progs[i])
        @assert progs[i].code == 3
    end

    input = 0
    while progs[5].code != 99
        for i in 1:5
            progs[i].input = input
            step_program!(progs[i])
            input = progs[i].output
        end
    end
    input
end

function max_feedback(prog)
    max = 0
    progs = Tuple(deepcopy(prog) for i in 1:5)

    foreach_perm(5) do phases
        for i in 1:5
            reset!(progs[i], prog)
        end
        m = test_feedback!(progs, 5 .+ phases)
        if m > max
            max = m
        end
    end
    max
end

part2() = max_feedback(read_program("input07.dat"))

end


text:
export @snoopi, @snoopi_deep

const __inf_timing__ = Tuple{Float64,MethodInstance}[]

if isdefined(Core.Compiler, :Params)
    function typeinf_ext_timed(linfo::Core.MethodInstance, params::Core.Compiler.Params)
        tstart = time()
        ret = Core.Compiler.typeinf_ext(linfo, params)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end
    function typeinf_ext_timed(linfo::Core.MethodInstance, world::UInt)
        tstart = time()
        ret = Core.Compiler.typeinf_ext(linfo, world)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end
    @noinline stop_timing() = ccall(:jl_set_typeinf_func, Cvoid, (Any,), Core.Compiler.typeinf_ext)
else
    function typeinf_ext_timed(interp::Core.Compiler.AbstractInterpreter, linfo::Core.MethodInstance)
        tstart = time()
        ret = Core.Compiler.typeinf_ext_toplevel(interp, linfo)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end
    function typeinf_ext_timed(linfo::Core.MethodInstance, world::UInt)
        tstart = time()
        ret = Core.Compiler.typeinf_ext_toplevel(linfo, world)
        tstop = time()
        push!(__inf_timing__, (tstop-tstart, linfo))
        return ret
    end
    @noinline stop_timing() = ccall(:jl_set_typeinf_func, Cvoid, (Any,), Core.Compiler.typeinf_ext_toplevel)
end

@noinline start_timing() = ccall(:jl_set_typeinf_func, Cvoid, (Any,), typeinf_ext_timed)

function sort_timed_inf(tmin)
    data = __inf_timing__
    if tmin > 0.0
        data = filter(tl->tl[1] >= tmin, data)
    end
    return sort(data; by=tl->tl[1])
end

"""
    inf_timing = @snoopi commands
    inf_timing = @snoopi tmin=0.0 commands

Execute `commands` while snooping on inference. Returns an array of `(t, linfo)`
tuples, where `t` is the amount of time spent inferring `linfo` (a `MethodInstance`).

Methods that take less time than `tmin` will not be reported.
"""
macro snoopi(args...)
    tmin = 0.0
    if length(args) == 1
        cmd = args[1]
    elseif length(args) == 2
        a = args[1]
        if isa(a, Expr) && a.head == :(=) && a.args[1] == :tmin
            tmin = a.args[2]
            cmd = args[2]
        else
            error("unrecognized input ", a)
        end
    else
        error("at most two arguments are supported")
    end
    return _snoopi(cmd, tmin)
end

function _snoopi(cmd::Expr, tmin = 0.0)
    return quote
        empty!(__inf_timing__)
        start_timing()
        try
            $(esc(cmd))
        finally
            stop_timing()
        end
        $sort_timed_inf($tmin)
    end
end

function __init__()
    # typeinf_ext_timed must be compiled before it gets run
    # We do this in __init__ to make sure it gets compiled to native code
    # (the *.ji file stores only the inferred code)
    if isdefined(Core.Compiler, :Params)
        @assert precompile(typeinf_ext_timed, (Core.MethodInstance, Core.Compiler.Params))
        @assert precompile(typeinf_ext_timed, (Core.MethodInstance, UInt))
    else
        @assert precompile(typeinf_ext_timed, (Core.Compiler.NativeInterpreter, Core.MethodInstance))
        @assert precompile(typeinf_ext_timed, (Core.MethodInstance, UInt))
    end
    precompile(start_timing, ())
    precompile(stop_timing, ())

    nothing
end


text:
# module Ecto.Schema

type Assoc
    vector::Vector{Tuple{Symbol,Any}}
    Assoc(; kw...) = new(Vector(kw))
    Assoc(vector) = new(vector)
end

function Base.getindex(assoc::Assoc, key::Symbol)
    for (k,v) in assoc.vector
        k==key && return v
    end
    throw(KeyError("key $key not found"))
end

Base.push!(assoc::Assoc, tup::Tuple{Symbol,Any}) = push!(assoc.vector, tup)

Base.keys(assoc::Assoc) = map(first, assoc.vector)
Base.values(assoc::Assoc) = map(last, assoc.vector)
Base.haskey(assoc::Assoc, key::Symbol) = key in keys(assoc)

Base.start(assoc::Assoc) = start(assoc.vector)
Base.next(assoc::Assoc, i::Int) = next(assoc.vector, i)
Base.done(assoc::Assoc, i::Int) = done(assoc.vector, i)

function ==(lhs::Assoc, rhs::Assoc)
    ==(lhs.vector, rhs.vector)
end

function Base.setdiff(lhs::Assoc, rhs::Assoc)::Assoc
    Assoc(setdiff(lhs.vector, rhs.vector))
end