text:
# Part of submodule Utils of BetaML - The Beta Machine Learning Toolkit
# Functions typically used for processing (manipulating) data, typically preprocessing data before running a ML model

# ------------------------------------------------------------------------------
# Various reshaping functions
import Base.reshape
""" reshape(myNumber, dims..) - Reshape a number as a n dimensional Array """
reshape(x::T, dims...) where {T <: Number} =   (x = [x]; reshape(x,dims) )
makeColVector(x::T) where {T} =  [x]
makeColVector(x::T) where {T <: AbstractArray} =  reshape(x,length(x))
makeRowVector(x::T) where {T <: Number} = return [x]'
makeRowVector(x::T) where {T <: AbstractArray} =  reshape(x,1,length(x))
"""Transform an Array{T,1} in an Array{T,2} and leave unchanged Array{T,2}."""
makeMatrix(x::AbstractArray) = ndims(x) == 1 ? reshape(x, (size(x)...,1)) : x

"""Return wheather an array is sortable, i.e. has methos issort defined"""
issortable(::AbstractArray{T,N})  where {T,N} = hasmethod(isless, Tuple{nonmissingtype(T),nonmissingtype(T)})

"""
    getPermutations(v::AbstractArray{T,1};keepStructure=false)

Return a vector of either (a) all possible permutations (uncollected) or (b) just those based on the unique values of the vector

Useful to measure accuracy where you don't care about the actual name of the labels, like in unsupervised classifications (e.g. clustering)

"""
function getPermutations(v::AbstractArray{T,1};keepStructure=false) where {T}
    if !keepStructure
        return Combinatorics.permutations(v)
    else
        classes       = unique(v)
        nCl           = length(classes)
        N             = size(v,1)
        pSet          = Combinatorics.permutations(1:nCl)
        nP            = length(pSet)
        vPermutations = fill(similar(v),nP)
        vOrigIdx      = [findfirst(x -> x == v[i] , classes) for i in 1:N]
        for (pIdx,perm) in enumerate(pSet)
            vPermutations[pIdx] = classes[perm[vOrigIdx]] # permuted specific version
        end
        return vPermutations
    end
end

""" singleUnique(x) Return the unique values of x whether x is an array of arrays, an array or a scalar"""
function singleUnique(x::Union{T,AbstractArray{T}}) where {T <: Union{Any,AbstractArray{T2}} where T2 <: Any }
    if typeof(x) <: AbstractArray{T2} where {T2 <: AbstractArray}
        return unique(vcat(unique.(x)...))
    elseif typeof(x) <: AbstractArray{T2} where {T2}
        return unique(x)
    else
        return [x]
    end
end

function oneHotEncoderRow(x::Union{AbstractArray{T},T}; d=maximum(x), factors=1:d,count = false) where {T <: Integer}
    x = makeColVector(x)
    out = zeros(Int64,d)
    for j in x
        out[j] = count ? out[j] + 1 : 1
    end
    return out
end

function oneHotEncoderRow(x::Union{AbstractArray{T},T};factors=singleUnique(x),d=length(factors),count = false) where {T}
    x = makeColVector(x)
    return oneHotEncoderRow(integerEncoder(x;factors=factors),d=length(factors);count=count)
end

"""
    oneHotEncoder(x;d,factors,count)

Encode arrays (or arrays of arrays) of categorical data as matrices of one column per factor.

The case of arrays of arrays is for when at each record you have more than one categorical output. You can then decide to encode just the presence of the factors or their counting

# Parameters:
- `x`: The data to convert (array or array of arrays)
- `d`: The number of dimensions in the output matrix [def: `maximum(x)` for integers and `length(factors)` otherwise]
- `factors`: The factors from which to encode [def: `1:d` for integer x or `unique(x)` otherwise]
- `count`: Wether to count multiple instances on the same dimension/record (`true`) or indicate just presence. [def: `false`]

# Examples
```julia
julia> oneHotEncoder(["a","c","c"],factors=["a","b","c","d"])
3×4 Matrix{Int64}:
 1  0  0  0
 0  0  1  0
 0  0  1  0
julia> oneHotEncoder([2,4,4])
3×4 Matrix{Int64}:
 0  1  0  0
 0  0  0  1
 0  0  0  1
 julia> oneHotEncoder([[2,2,1],[2,4,4]],count=true)
2×4 Matrix{Int64}:
 1  2  0  0
 0  1  0  2
```
"""
function oneHotEncoder(x::Union{T,AbstractVector{T}};factors=singleUnique(x),d=length(factors),count=false) where {T <: Union{Any,AbstractVector{T2}} where T2 <: Any  }
    if typeof(x) <: AbstractVector
        n  = length(x)
        out = zeros(Int64,n,d)
        for (i,x) in enumerate(x)
          out[i,:] = oneHotEncoderRow(x;factors=factors,count = count)
        end
        return out
    else
       out = zeros(Int64,1,d)
       out[1,:] = oneHotEncoderRow(x;factors=factors,count = count)
       return out
   end
end

function oneHotEncoder(Y::Union{Ti,AbstractVector{Ti}};d=maximum(maximum.(Y)),factors=1:d,count=false) where {Ti <: Union{Integer,AbstractVector{Ti2}} where Ti2 <: Integer  }
    n   = length(Y)
    if d < maximum(maximum.(Y))
        error("Trying to encode elements with indexes greater than the provided number of dimensions. Please increase d.")
    end
    out = zeros(Int64,n,d)
    for (i,y) in enumerate(Y)
        out[i,:] = oneHotEncoderRow(y;d=d,factors=1:d,count = count)
    end
    return out
end

findfirst(el::T,cont::Array{T};returnTuple=true) where {T<:Union{AbstractString,Number}} = ndims(cont) > 1 && returnTuple ? Tuple(findfirst(x -> isequal(x,el),cont)) : findfirst(x -> isequal(x,el),cont)
#findfirst(el::T,cont::Array{T,N};returnTuple=true) where {T,N} = returnTuple ? Tuple(findfirst(x -> isequal(x,el),cont)) : findfirst(x -> isequal(x,el),cont)
#findfirst(el::T,cont::Array{T,1};returnTuple=true) where {T} =  findfirst(x -> isequal(x,el),cont)

findall(el::T, cont::Array{T};returnTuple=true) where {T} = ndims(cont) > 1 && returnTuple ? Tuple.(findall(x -> isequal(x,el),cont)) : findall(x -> isequal(x,el),cont)

"""
    integerEncoder(x;factors=unique(x))

Encode an array of T to an array of integers using the their position in `factor` vector (default to the unique vector of the input array)

# Parameters:
- `x`: The vector to encode
- `factors`: The vector of factors whose position is the result of the encoding [def: `unique(x)`]
# Return:
- A vector of [1,length(x)] integers corresponding to the position of each element in the `factors` vector`
# Note:
- Attention that while this function creates a ordered (and sortable) set, it is up to the user to be sure that this "property" is not indeed used in his code if the unencoded data is indeed unordered.
# Example:
```
julia> integerEncoder(["a","e","b","e"],factors=["a","b","c","d","e"]) # out: [1,5,2,5]
```
"""
function integerEncoder(x::AbstractVector;factors=Base.unique(x))
    #return findfirst.(x,Ref(factors)) slower
    return  map(i -> findfirst(j -> j==i,factors) , x  )
end

"""
    integerDecoder(x,factors::AbstractVector{T};unique)

Decode an array of integers to an array of T corresponding to the elements of `factors`

# Parameters:
- `x`: The vector to decode
- `factors`: The vector of elements to use for the encoding
- `unique`: Wether `factors` is already made of unique elements [def: `true`]
# Return:
- A vector of length(x) elements corresponding to the (unique) `factors` elements at the position x
# Example:
```
julia> integerDecoder([1, 2, 2, 3, 2, 1],["aa","cc","bb"]) # out: ["aa","cc","cc","bb","cc","aa"]
```
"""
function integerDecoder(x,factors::AbstractVector{T};unique=true) where{T}
    uniqueTarget =  unique ? factors :  Base.unique(factors)
    return map(i -> uniqueTarget[i], x )
end

"""
    partition(data,parts;shuffle,dims,rng)

Partition (by rows) one or more matrices according to the shares in `parts`.

# Parameters
* `data`: A matrix/vector or a vector of matrices/vectors
* `parts`: A vector of the required shares (must sum to 1)
* `shufle`: Whether to randomly shuffle the matrices (preserving the relative order between matrices)
* `dims`: The dimension for which to partition [def: `1`]
* `copy`: Wheter to _copy_ the actual data or only create a reference [def: `true`]
* `rng`: Random Number Generator (see [`FIXEDSEED`](@ref)) [deafult: `Random.GLOBAL_RNG`]

# Notes:
* The sum of parts must be equal to 1
* The number of elements in the specified dimension must be the same for all the arrays in `data`

# Example:
```julia
julia> x = [1:10 11:20]
julia> y = collect(31:40)
julia> ((xtrain,xtest),(ytrain,ytest)) = partition([x,y],[0.7,0.3])
 ```
 """
function partition(data::AbstractArray{T,1},parts::AbstractArray{Float64,1};shuffle=true,dims=1,copy=true,rng = Random.GLOBAL_RNG) where T <: AbstractArray
        # the sets of vector/matrices
        N = size(data[1],dims)
        all(size.(data,dims) .== N) || @error "All matrices passed to `partition` must have the same number of elements for the required dimension"
        ridx = shuffle ? Random.shuffle(rng,1:N) : collect(1:N)
        return partition.(data,Ref(parts);shuffle=shuffle,dims=dims,fixedRIdx = ridx,copy=copy,rng=rng)
end

function partition(data::AbstractArray{T,Ndims}, parts::AbstractArray{Float64,1};shuffle=true,dims=1,fixedRIdx=Int64[],copy=true,rng = Random.GLOBAL_RNG) where {T,Ndims}
    # the individual vector/matrix
    N        = size(data,dims)
    nParts   = size(parts)
    toReturn = toReturn = Array{AbstractArray{T,Ndims},1}(undef,nParts)
    if !(sum(parts) ≈ 1)
        @error "The sum of `parts` in `partition` should total to 1."
    end
    ridx = fixedRIdx
    if (isempty(ridx))
       ridx = shuffle ? Random.shuffle(rng, 1:N) : collect(1:N)
    end
    allDimIdx = convert(Vector{Union{UnitRange{Int64},Vector{Int64}}},[1:i for i in size(data)])
    current = 1
    cumPart = 0.0
    for (i,p) in enumerate(parts)
        cumPart += parts[i]
        final = i == nParts ? N : Int64(round(cumPart*N))
        allDimIdx[dims] = ridx[current:final]
        toReturn[i]     = copy ? data[allDimIdx...] : @views data[allDimIdx...]
        current         = (final +=1)
    end
    return toReturn
end

"""
    getScaleFactors(x;skip)

Return the scale factors (for each dimensions) in order to scale a matrix X (n,d)
such that each dimension has mean 0 and variance 1.

# Parameters
- `x`: the (n × d) dimension matrix to scale on each dimension d
- `skip`: an array of dimension index to skip the scaling [def: `[]`]

# Return
- A touple whose first elmement is the shift and the second the multiplicative
term to make the scale.
"""
function getScaleFactors(x;skip=[])
    μ  = mean(x,dims=1)
    σ² = var(x,corrected=false,dims=1)
    sfμ = - μ
    sfσ² = 1 ./ sqrt.(σ²)
    for i in skip
        sfμ[i] = 0
        sfσ²[i] = 1
    end
    return (sfμ,sfσ²)
end

"""
    scale(x,scaleFactors;rev)

Perform a linear scaling of x using scaling factors `scaleFactors`.

# Parameters
- `x`: The (n × d) dimension matrix to scale on each dimension d
- `scalingFactors`: A tuple of the constant and multiplicative scaling factor
respectively [def: the scaling factors needed to scale x to mean 0 and variance 1]
- `rev`: Whether to invert the scaling [def: `false`]

# Return
- The scaled matrix

# Notes:
- Also available `scale!(x,scaleFactors)` for in-place scaling.
- Retrieve the scale factors with the `getScaleFactors()` function
"""
function scale(x,scaleFactors=(-mean(x,dims=1),1 ./ sqrt.(var(x,corrected=false,dims=1))); rev=false )
    if (!rev)
      y = (x .+ scaleFactors[1]) .* scaleFactors[2]
    else
      y = (x ./ scaleFactors[2]) .- scaleFactors[1]
    end
    return y
end
function scale!(x,scaleFactors=(-mean(x,dims=1),1 ./ sqrt.(var(x,corrected=false,dims=1))); rev=false)
    if (!rev)
        x .= (x .+ scaleFactors[1]) .* scaleFactors[2]
    else
        x .= (x ./ scaleFactors[2]) .- scaleFactors[1]
    end
    return nothing
end

"""
pca(X;K,error)

Perform Principal Component Analysis returning the matrix reprojected among the dimensions of maximum variance.

# Parameters:
- `X` : The (N,D) data to reproject
- `K` : The number of dimensions to maintain (with K<=D) [def: `nothing`]
- `error`: The maximum approximation error that we are willing to accept [def: `0.05`]

# Return:
- A named tuple with:
  - `X`: The reprojected (NxK) matrix with the column dimensions organized in descending order of of the proportion of explained variance
  - `K`: The number of dimensions retieved
  - `error`: The actual proportion of variance not explained in the reprojected dimensions
  - `P`: The (D,K) matrix of the eigenvectors associated to the K-largest eigenvalues used to reproject the data matrix
  - `explVarByDim`: An array of dimensions D with the share of the cumulative variance explained by dimensions (the last element being always 1.0)

# Notes:
- If `K` is provided, the parameter `error` has no effect.
- If one doesn't know _a priori_ the error that she/he is willling to accept, nor the wished number of dimensions, he/she can run this pca function with `out = pca(X,K=size(X,2))` (i.e. with K=D), analise the proportions of explained cumulative variance by dimensions in `out.explVarByDim`, choose the number of dimensions K according to his/her needs and finally pick from the reprojected matrix only the number of dimensions needed, i.e. `out.X[:,1:K]`.

# Example:
```julia
julia> X = [1 10 100; 1.1 15 120; 0.95 23 90; 0.99 17 120; 1.05 8 90; 1.1 12 95]
6×3 Matrix{Float64}:
 1.0   10.0  100.0
 1.1   15.0  120.0
 0.95  23.0   90.0
 0.99  17.0  120.0
 1.05   8.0   90.0
 1.1   12.0   95.0
julia> X = pca(X,error=0.05).X
6×2 Matrix{Float64}:
 100.449    3.1783
 120.743    6.80764
  91.3551  16.8275
 120.878    8.80372
  90.3363   1.86179
  95.5965   5.51254
```

"""
function pca(X;K=nothing,error=0.05)
    # debug
    #X = [1 10 100; 1.1 15 120; 0.95 23 90; 0.99 17 120; 1.05 8 90; 1.1 12 95]
    #K = nothing
    #error=0.05

    (N,D) = size(X)
    if !isnothing(K) && K > D
        @error("The parameter K must be ≤ D")
    end
    Σ = (1/N) * X'*(I-(1/N)*ones(N)*ones(N)')*X
    E = eigen(Σ) # eigenvalues are ordered from the smallest to the largest
    totVar  = sum(E.values)
    explVarByDim = [sum(E.values[D-k+1:end])/totVar for k in 1:D]
    propVarExplained = 0.0
    if K == nothing
        for k in 1:D
            if explVarByDim[k] >= (1 - error)
                propVarExplained  = explVarByDim[k]
                K                 = k
                break
            end
        end
    else
        propVarExplained = explVarByDim[K]
    end

    P = E.vectors[:,end:-1:D-K+1] # bug corrected 2/9/2021

    return (X=X*P,K=K,error=1-propVarExplained,P=P,explVarByDim=explVarByDim)
end

"""
    colsWithMissing(x)

Retuyrn an array with the ids of the columns where there is at least a missing value.
"""
function colsWithMissing(x)
    colsWithMissing = Int64[]
    (N,D) = size(x)
    for d in 1:D
        for n in 1:N
            if ismissing(x[n,d])
                push!(colsWithMissing,d)
                break
            end
        end
    end
    return colsWithMissing
end

"""
    crossValidation(f,data,sampler;dims,verbosity,returnStatistics)

Perform crossValidation according to `sampler` rule by calling the function f and collecting its output

# Parameters
- `f`: The user-defined function that consume the specific train and validation data and return somehting (often the associated validation error). See later
- `data`: A single n-dimenasional array or a vector of them (e.g. X,Y), depending on the tasks required by `f`.
- sampler: An istance of a ` AbstractDataSampler`, defining the "rules" for sampling at each iteration. [def: `KFold(nSplits=5,nRepeats=1,shuffle=true,rng=Random.GLOBAL_RNG)` ]
- `dims`: The dimension over performing the crossValidation i.e. the dimension containing the observations [def: `1`]
- `verbosity`: The verbosity to print information during each iteration (this can also be printed in the `f` function) [def: `STD`]
- `returnStatistics`: Wheter crossValidation should return the statistics of the output of `f` (mean and standard deviation) or the whole outputs [def: `true`].

# Notes

crossValidation works by calling the function `f`, defined by the user, passing to it the tuple `trainData`, `valData` and `rng` and collecting the result of the function f. The specific method for which `trainData`, and `valData` are selected at each iteration depends on the specific `sampler`, whith a single 5 k-fold rule being the default.

This approach is very flexible because the specific model to employ or the metric to use is left within the user-provided function. The only thing that crossValidation does is provide the model defined in the function `f` with the opportune data (and the random number generator).

**Input of the user-provided function**
`trainData` and `valData` are both themselves tuples. In supervised models, crossValidations `data` should be a tuple of (X,Y) and `trainData` and `valData` will be equivalent to (xtrain, ytrain) and (xval, yval). In unsupervised models `data` is a single array, but the training and validation data should still need to be accessed as  `trainData[1]` and `valData[1]`.
**Output of the user-provided function**
The user-defined function can return whatever. However, if `returnStatistics` is left on its default `true` value the user-defined function must return a single scalar (e.g. some error measure) so that the mean and the standard deviation are returned.

Note that `crossValidation` can beconveniently be employed using the `do` syntax, as Julia automatically rewrite `crossValidation(data,...) trainData,valData,rng  ...user defined body... end` as `crossValidation(f(trainData,valData,rng ), data,...)`

# Example

```
julia> X = [11:19 21:29 31:39 41:49 51:59 61:69];
julia> Y = [1:9;];
julia> sampler = KFold(nSplits=3);
julia> (μ,σ) = crossValidation([X,Y],sampler) do trainData,valData,rng
                 (xtrain,ytrain) = trainData; (xval,yval) = valData
                 trainedModel    = buildForest(xtrain,ytrain,30)
                 predictions     = predict(trainedModel,xval)
                 ϵ               = meanRelError(predictions,yval,normRec=false)
                 return ϵ
               end
(0.3202242202242202, 0.04307662219315022)
```

"""
function crossValidation(f,data,sampler=KFold(nSplits=5,nRepeats=1,shuffle=true,rng=Random.GLOBAL_RNG);dims=1,verbosity=STD, returnStatistics=true)
    iterResults = []
    for (i,iterData) in enumerate(SamplerWithData(sampler,data,dims))
       iterResult = f(iterData[1],iterData[2],sampler.rng)
       push!(iterResults,iterResult)
       if verbosity > STD
           println("Done iteration $i. This iteration output: $iterResult")
       end
    end
    if returnStatistics  return (mean(iterResults),std(iterResults)) else return iterResults end
end

"""
   classCountsWithLabels(x)

Return a dictionary that counts the number of each unique item (rows) in a dataset.

"""
function classCountsWithLabels(x;classes=nothing)
    dims = ndims(x)
    if dims == 1
        T = eltype(x)
    else
        T = Array{eltype(x),1}
    end
    if classes != nothing
        counts = Dict([u=>0 for u in classes])
    else
        counts = Dict{T,Int64}()  # a dictionary of label -> count.
    end
    for i in 1:size(x,1)
        if dims == 1
            label = x[i]
        else
            label = x[i,:]
        end
        if !(label in keys(counts))
            counts[label] = 1
        else
            counts[label] += 1
        end
    end
    return counts
end

"""
   classCounts(x;classes=nothing)

Return a (unsorted) vector with the counts of each unique item (element or rows) in a dataset.

If order is important or not all classes are present in the data, a preset vectors of classes can be given in the parameter `classes`

"""
function classCounts(x; classes=nothing)
   if classes == nothing # order doesn't matter
      return values(classCountsWithLabels(x;classes=classes))
   else
       cWithLabels = classCountsWithLabels(x;classes=classes)
       return [cWithLabels[k] for k in classes]
   end
end

"""
   mode(dict::Dict{T,Float64};rng)

Return the key with highest mode (using rand in case of multimodal values)

"""
function mode(dict::Dict{T,Float64};rng = Random.GLOBAL_RNG) where {T}
    mks = [k for (k,v) in dict if v==maximum(values(dict))]
    if length(mks) == 1
        return mks[1]
    else
        return mks[rand(rng,1:length(mks))]
    end
end

"""
   mode(v::AbstractVector{T};rng)

Return the position with the highest mode (using rand in case of multimodal values)

"""
function mode(v::AbstractVector{T};rng = Random.GLOBAL_RNG) where {T <: Number}
    mpos = findall(x -> x == maximum(v),v)
    if length(mpos) == 1
        return mpos[1]
    else
        return mpos[rand(rng,1:length(mpos))]
    end
end

"""
  mode(elements,rng)

Given a vector of dictionaries whose key is numerical (e.g. probabilities), a vector of vectors or a matrix, it returns the mode of each element (dictionary, vector or row) in terms of the key or the position.

Use it to return a unique value from a multiclass classifier returning probabilities.

# Note:
- If multiple classes have the highest mode, one is returned at random (use the parameter `rng` to fix the stochasticity)

"""
function mode(dicts::AbstractArray{Dict{T,Float64}};rng = Random.GLOBAL_RNG) where {T}
    return mode.(dicts;rng=rng)
end

function mode(vals::AbstractArray{T,1};rng = Random.GLOBAL_RNG) where {T <: AbstractArray{T2,1} where T2 <: Number}
    return mode.(vals;rng=rng)
end
function mode(vals::AbstractArray{T,2};rng = Random.GLOBAL_RNG) where {T <: Number}
    return [mode(r;rng=rng) for r in eachrow(vals)]
end

"""
   meanDicts(dicts)

Compute the mean of the values of an array of dictionaries.

Given `dicts` an array of dictionaries, `meanDicts` first compute the union of the keys and then average the values.
If the original valueas are probabilities (non-negative items summing to 1), the result is also a probability distribution.

"""
function meanDicts(dicts; weights=ones(length(dicts)))
    T = eltype(keys(dicts[1]))
    allkeys = union([keys(i) for i in dicts]...)
    outDict = Dict{T,Float64}()
    ndicts = length(dicts)
    totWeights = sum(weights)
    for k in allkeys
        v = 0
        for (i,d) in enumerate(dicts)
            if k in keys(d)
                v += (d[k])*(weights[i]/totWeights)
            end
        end
        outDict[k] = v
    end

    return outDict
end

# ------------------------------------------------------------------------------
# Other mathematical/computational functions

""" LogSumExp for efficiently computing log(sum(exp.(x))) """
lse(x) = maximum(x)+log(sum(exp.(x .- maximum(x))))
""" Sterling number: number of partitions of a set of n elements in k sets """
sterling(n::BigInt,k::BigInt) = (1/factorial(k)) * sum((-1)^i * binomial(k,i)* (k-i)^n for i in 0:k)
sterling(n::Int64,k::Int64)   = sterling(BigInt(n),BigInt(k))


text:
# Functions for choosing initial centers (often called seeds)

function randseed_initialize!{T<:FloatingPoint}(x::Matrix{T}, centers::Matrix{T})
    n = size(x, 2)
    k = size(centers, 2)
    si = sample_without_replacement(1:n, k)
    centers[:,:] = x[:,si]
end

function kmeanspp_initialize!{T<:FloatingPoint}(x::Matrix{T}, centers::Matrix{T})
    n = size(x, 2)
    k = size(centers, 2)

    # randomly pick the first center
    si = rand(1:n)
    v = x[:,si]
    centers[:,1] = v

    # initialize the cost vector
    costs = colwise(SqEuclidean(), v, x)

    # pick remaining (with a chance proportional to cost)
    for i = 2 : k
        si = sample_by_weights(costs)
        v = x[:,si]
        centers[:,i] = v

        # update costs

        if i < k
            new_costs = colwise(SqEuclidean(), v, x)
            costs = min(costs, new_costs)
        end
    end
end

ext: jl


text:
using brainflow

# specify logging library to use
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_current_board_data(32, board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)

data_first_channel = data[eeg_channels[1], :]
println("Original Data First Channel")
println(data_first_channel)
brainflow.perform_rolling_filter(data_first_channel, 3, brainflow.MEAN)
println("After Rolling Filter")
println(data_first_channel)

data_second_channel = data[eeg_channels[2], :]
println("Original Data Second Channel")
println(data_second_channel)
brainflow.perform_wavelet_denoising(data_second_channel, "db4", 2)
println("After Wavelet Denoising")
println(data_second_channel)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
# ---
# title: 677. Map Sum Pairs
# id: problem677
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Trie
# link: <https://leetcode.com/problems/map-sum-pairs/description/>
# hidden: true
# ---
#
# Implement the `MapSum` class:
#
#   * `MapSum()` Initializes the `MapSum` object.
#   * `void insert(String key, int val)` Inserts the `key-val` pair into the map. If the `key` already existed, the original `key-value` pair will be overridden to the new one.
#   * `int sum(string prefix)` Returns the sum of all the pairs' value whose `key` starts with the `prefix`.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["MapSum", "insert", "sum", "insert", "sum"]
#     [[], ["apple", 3], ["ap"], ["app", 2], ["ap"]]
#     **Output**
#     [null, null, 3, null, 5]
#
#     **Explanation**
#     MapSum mapSum = new MapSum();
#     mapSum.insert("apple", 3);
#     mapSum.sum("ap");           // return 3 ( _ap_ ple = 3)
#     mapSum.insert("app", 2);
#     mapSum.sum("ap");           // return 5 ( _ap_ ple + _ap_ p = 3 + 2 = 5)
#
#
#
#
# **Constraints:**
#
#   * `1 <= key.length, prefix.length <= 50`
#   * `key` and `prefix` consist of only lowercase English letters.
#   * `1 <= val <= 1000`
#   * At most `50` calls will be made to `insert` and `sum`.
#
#
## @lc code=start
using LeetCode

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


text:
include("../coin-or-common.jl")

sources = [
    GitSource(
        "https://github.com/coin-or/SHOT.git",
        SHOT_gitsha,
    ),
    ArchiveSource("https://github.com/phracker/MacOSX-SDKs/releases/download/10.15/MacOSX10.15.sdk.tar.xz",
                  "2408d07df7f324d3beea818585a6d990ba99587c218a3969f924dfcc4de93b62"),
    DirectorySource("./bundled"),
]

script = raw"""
cd $WORKSPACE/srcdir/SHOT
git submodule update --init --recursive
# Disable run_source_test in CppAD
atomic_patch -p1 ../patches/CppAD.patch
if [[ "${target}" == x86_64-apple-darwin* ]]; then
    # Work around the issue
    #     /workspace/srcdir/SHOT/src/Model/../Model/Simplifications.h:1370:26: error: 'value' is unavailable: introduced in macOS 10.14
    #                     optional.value()->coefficient *= -1.0;
    #                              ^
    #     /opt/x86_64-apple-darwin14/x86_64-apple-darwin14/sys-root/usr/include/c++/v1/optional:947:27: note: 'value' has been explicitly marked unavailable here
    #         constexpr value_type& value() &
    #                               ^
    export CXXFLAGS="-mmacosx-version-min=10.15"
    # ...and install a newer SDK which supports `std::filesystem`
    pushd $WORKSPACE/srcdir/MacOSX10.*.sdk
    rm -rf /opt/${target}/${target}/sys-root/System
    cp -ra usr/* "/opt/${target}/${target}/sys-root/usr/."
    cp -ra System "/opt/${target}/${target}/sys-root/."
    popd
elif [[ "${target}" == aarch64-apple-darwin* ]]; then
    # TODO: we need to fix this in the compiler wrappers
    export CXXFLAGS="-mmacosx-version-min=11.0"
elif [[ ${target} == *mingw* ]]; then
    export LDFLAGS="-L${libdir}"
fi

mkdir -p build
cd build

cmake \
    -DCMAKE_INSTALL_PREFIX=${prefix} \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DHAS_CBC=on \
    -DCBC_DIR=${prefix} \
    -DHAS_IPOPT=on \
    -DIPOPT_DIR=${prefix} \
    -DHAS_AMPL=on \
    -DGENERATE_EXE=on \
    ..

make -j${nproc}
make install

install_license ../LICENSE
"""

products = [
    ExecutableProduct("SHOT", :amplexe),
    LibraryProduct("libSHOTSolver", :libshotsolver),
]

dependencies = [
    Dependency("ASL_jll", ASL_version),
    Dependency("Cbc_jll", compat="$(Cbc_version)"),
    Dependency("Ipopt_jll", compat="$(Ipopt_version)"),
    Dependency("CompilerSupportLibraries_jll"),
]

build_tarballs(
    ARGS,
    "SHOT",
    SHOT_version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    julia_compat = "1.6",
    preferred_gcc_version = v"9",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module SDtoolbox

export zndsolve, cvsolve, cell_size, CJspeed, PostShock_eq, PostShock_fr

include("postshock.jl")
include("cv.jl")
include("znd.jl")

# submodules that are essentially standalone
using .Postshock
using .CV
using .ZND

# ensure all python packages are loaded at run time to make the package
# precompile safe.
using PyCall

const ct = PyNULL()

function __init__()
    # load cantera python package
    copy!(ct, pyimport_conda("cantera","cantera","cantera"))

    # use xml mech files in tests for cantera before v2.5, otherwise yaml
    if VersionNumber(ct.__version__) < v"2.5"
        global CANTERA_MECH_FILETYPE = "xml"
    else
        global CANTERA_MECH_FILETYPE = "yaml"
    end
end

include("cell_size.jl")

end


text:
export Lemma

const SPACE = ' '

struct Lemma
    word::AbstractString
    pos::Char
    tagsense_count::Int
    synset_offsets::Vector{Int}
    id::Int
    pointer_syms::Vector{AbstractString}
end

function Lemma(lexicon_line, id)
    parts = split(lexicon_line, SPACE)

    word = shift!(parts)
    pos = shift!(parts)[1]
    synset_count = parse(Int, shift!(parts))

    n_syms = parse(Int, shift!(parts))
    pointer_syms = [s for s in parts[1:n_syms]]
    parts = parts[n_syms+2:end]

    tagsense_count = parse(Int, shift!(parts))
    synset_offsets = [parse(Int, c) for c in parts[1:synset_count]]

    Lemma(word, pos, tagsense_count, synset_offsets, id, pointer_syms)
end

Base.show(io::IO, lemma::Lemma) = print(io, "$(lemma.word).$(lemma.pos)")


text:
# Julia code
# Copyright (c) 2022 the SLICOTMath.jl developers
# Portions extracted from SLICOT-Reference distribution:
# Copyright (c) 2002-2020 NICONET e.V.

function run_mb02jx(datfile, io=stdout)
    NIN = 5
    NOUT = 6
    KMAX = 20
    LMAX = 20
    MMAX = 20
    NMAX = 20
    LDR  = NMAX*LMAX
    LDQ  = MMAX*KMAX
    LDTC = MMAX*KMAX
    LDTR = KMAX
    LDWORK = ( MMAX*KMAX + NMAX*LMAX ) *( LMAX + 2*KMAX ) + 5*LMAX + MMAX*KMAX + NMAX*LMAX
    TC = Array{Float64,2}(undef, LDTC,LMAX)
    TR = Array{Float64,2}(undef, LDTR,NMAX*LMAX)
    Q = Array{Float64,2}(undef, LDQ,NMAX*LMAX)
    R = Array{Float64,2}(undef, LDR,NMAX*LMAX)
    JPVT = Array{BlasInt,1}(undef, NMAX*LMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    K = parse(BlasInt, vs[1])
    L = parse(BlasInt, vs[2])
    M = parse(BlasInt, vs[3])
    N = parse(BlasInt, vs[4])
    TOL1 = parse(Float64, replace(vs[5],'D'=>'E'))
    TOL2 = parse(Float64, replace(vs[6],'D'=>'E'))
    JOB =  vs[7][1]
if ( K<=0 || K>KMAX )
elseif ( L<=0 || L>LMAX )
elseif ( M<=0 || M>MMAX )
elseif ( N<=0 || N>NMAX )
else

    vs = String[]
    _isz,_jsz = (M*K,L)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       TC[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (K,N*L-L)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       TR[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
# interp call 1

    RNK, INFO = SLICOT.mb02jx!(JOB, K, L, M, N, TC, TR, Q, R, JPVT, TOL1, TOL2, LDWORK)
    println(io, "RNK = $RNK")
    @test INFO == 0
    INFO == 0 || return

if ( INFO!=0 )
else
if ( LSAME( JOB, 'Q' ) )
# interp output 1
    println(io,"Q:")
    _nc = RNK
    _nr = M*K
    show(io,"text/plain",Q[1:_nr,1:_nc])
    println(io,)

end # if
# interp output 2
    println(io,"R:")
    _nc = RNK
    _nr = N*L
    show(io,"text/plain",R[1:_nr,1:_nc])
    println(io,)

# interp output 3
    println(io,"JPVT:")
    _nr = min( M*K, N*L )
    show(io,"text/plain",JPVT[1:_nr])
    println(io,)

end # if
end # if
    close(f)
end # run_X()


text:
export marginalrule

@marginalrule MvNormalMeanCovariance(:out_μ_Σ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, m_Σ::PointMass) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanCovariance(mean(m_μ), mean(m_Σ)), m_out), μ = m_μ, Σ = m_Σ)
end

@marginalrule MvNormalMeanCovariance(:out_μ_Σ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, m_Σ::PointMass) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), m_μ, MvNormalMeanCovariance(mean(m_out), mean(m_Σ))), Σ = m_Σ)
end

@marginalrule MvNormalMeanCovariance(:out_μ_Σ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, m_Σ::PointMass) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_m, W_m = weightedmean_precision(m_μ)

    W_bar = cholinv(mean(m_Σ))

    xi = [ xi_out; xi_m ]
    W  = [ W_out + W_bar -W_bar; -W_bar W_m + W_bar ]

    return (out_μ = MvNormalWeightedMeanPrecision(xi, W), Σ = m_Σ)
end

@marginalrule MvNormalMeanCovariance(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::MultivariateNormalDistributionsFamily, q_Σ::Any) = begin
    xi_out, W_out = weightedmean_precision(m_out)
    xi_m, W_m = weightedmean_precision(m_μ)

    W_bar = cholinv(mean(q_Σ))

    xi = [ xi_out; xi_m ]
    W  = [ W_out + W_bar -W_bar; -W_bar W_m + W_bar ]

    return MvNormalWeightedMeanPrecision(xi, W)
end

@marginalrule MvNormalMeanCovariance(:out_μ) (m_out::PointMass, m_μ::MultivariateNormalDistributionsFamily, q_Σ::Any) = begin
    return (out = m_out, μ = prod(ProdAnalytical(), MvNormalMeanCovariance(mean(m_out), mean(q_Σ)), m_μ))
end

@marginalrule MvNormalMeanCovariance(:out_μ) (m_out::MultivariateNormalDistributionsFamily, m_μ::PointMass, q_Σ::Any) = begin
    return (out = prod(ProdAnalytical(), MvNormalMeanCovariance(mean(m_μ), mean(q_Σ)), m_out), μ = m_μ)
end


text:
using Documenter, GenomicAnnotations

makedocs(sitename = "GenomicAnnotations.jl", authors = "Karl Dyrhage",
    pages = [
        "index.md",
        "I/O" => "io.md",
        "Accessing and modifying annotations" => "accessing.md",
        "Filtering: the @genes macro" => "genes.md"
    ],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true")
)

deploydocs(
    repo = "github.com/kdyrhage/GenomicAnnotations.jl.git",
)


text:
using CUTEst, DelimitedFiles

probs= CUTEst.select()

similar(a,b;atol=10,rtol=0.5) = (abs(a-b) < atol) || (abs(a-b)/max(a,b) < rtol)

function is_similar_exists(bin,nlp)
    for (name,nvar,ncon,nnzh,nnzj) in bin
        if similar(nvar,nlp.meta.nvar) && similar(ncon,nlp.meta.ncon) && similar(nnzh,nlp.meta.nnzh) && similar(nnzj,nlp.meta.nnzj)
            return true
        end
    end
    return false
end

bin = Tuple{String,Int,Int,Int,Int}[]

for i=1:length(probs)
    prob = probs[i]
    nlp = CUTEstModel(prob)
    finalize(nlp)

    if is_similar_exists(bin,nlp)
        println("#$cnt Skipping $(nlp.meta.name)")
    else
        println("#$cnt Putting $(nlp.meta.name) in the bin")
        push!(bin, (nlp.meta.name,nlp.meta.nvar,nlp.meta.ncon,nlp.meta.nnzh,nlp.meta.nnzj))
    end
end

writedlm("cutest-quick-names.csv",[name for (name,nvar,ncon,nnzh,nnzj) in bin],',')


text:
using Unicode

"""

# cmdline

Recursively parse the model to construct command line.

### Method
```julia
cmdline(m)
```

### Required arguments
```julia
* `m::Stanmodel`                : Stanmodel
```

### Related help
```julia
?Stanmodel                      : Create a StanModel
```
"""
function cmdline(m)
  cmd = ``
  if isa(m, Stanmodel)
    # Handle the model name field for unix and windows
    cmd = @static Sys.isunix() ? `./$(getfield(m, :name))` : `cmd /c $(getfield(m, :name)).exe`

    # Method (sample, optimize, variational and diagnose) specific portion of the model
    cmd = `$cmd $(cmdline(getfield(m, :method)))`

    # Common to all models
    cmd = `$cmd $(cmdline(getfield(m, :random)))`

    # Init file required?
    if length(m.init_file) > 0 && isfile(m.init_file)
      cmd = `$cmd init=$(m.init_file)`
    end

    # Data file required?
    if length(m.data_file) > 0 && isfile(m.data_file)
      cmd = `$cmd id=$(m.id) data file=$(m.data_file)`
    end

    # Output options
    cmd = `$cmd output`
    if length(getfield(m, :output).file) > 0
      cmd = `$cmd file=$(string(getfield(m, :output).file))`
    end
    if length(getfield(m, :output).diagnostic_file) > 0
      cmd = `$cmd diagnostic_file=$(string(getfield(m, :output).diagnostic_file))`
    end
    cmd = `$cmd refresh=$(string(getfield(m, :output).refresh))`
  else
    # The 'recursive' part
    #println(lowercase(string(typeof(m))))
    if isa(m, SamplingAlgorithm) || isa(m, OptimizeAlgorithm)
      cmd = `$cmd algorithm=$(split(lowercase(string(typeof(m))), '.')[end])`
    elseif isa(m, Engine)
      cmd = `$cmd engine=$(split(lowercase(string(typeof(m))), '.')[end])`
    elseif isa(m, Diagnostics)
      cmd = `$cmd test=$(split(lowercase(string(typeof(m))), '.')[end])`
    else
      cmd = `$cmd $(split(lowercase(string(typeof(m))), '.')[end])`
    end
    #println(cmd)
    for name in fieldnames(typeof(m))
      #println("$(name) = $(getfield(m, name)) ($(typeof(getfield(m, name))))")
      if  isa(getfield(m, name), String) || isa(getfield(m, name), Tuple)
        cmd = `$cmd $(name)=$(getfield(m, name))`
      elseif length(fieldnames(typeof(getfield(m, name)))) == 0
        if isa(getfield(m, name), Bool)
          cmd = `$cmd $(name)=$(getfield(m, name) ? 1 : 0)`
        else
          if name == :metric || isa(getfield(m, name), DataType)
            cmd = `$cmd $(name)=$(split(lowercase(string(typeof(getfield(m, name)))), '.')[end])`
          else
            cmd = `$cmd $(name)=$(getfield(m, name))`
          end
        end
      else
        cmd = `$cmd $(cmdline(getfield(m, name)))`
      end
    end
  end
  cmd
end


text:
struct Action <: AbstractAction
    i::UInt32
end


text:
module StandardizedPredictors

export
       center,
       center!,
       Center,
       CenteredTerm,
       scale,
       scale!,
       Scale,
       ScaledTerm,
       zscore,    # from StatsBase
       zscore!,   # from StatsBase
       ZScore,
       ZScoredTerm

using StatsModels
using StatsBase
using Statistics

"""
    _standard(xs::AbstractArray, val)

Translate an abstract standardization value to a concrete one based on `xs`.

`nothing` and already concrete `Number` `val`s are passed through.
Otherwise, `val(xs)` is returned.
"""
_standard(::AbstractArray, t::Number) = t
_standard(::AbstractArray, ::Nothing) = nothing
_standard(xs::AbstractArray, t) = t(xs)

include("utils.jl")
include("centering.jl")
include("scaling.jl")
include("zscoring.jl")

end


text:
using StructJuMP, JuMP
using StructJuMPSolverInterface

include("select_solver.jl")

#############
# A sample model
#############
scen = 1
m = StructuredModel(num_scenarios=scen)
@variable(m, x[1:4])
# @variable(m, -100<=x[1:4]<=100)
@NLobjective(m, Min,  1*x[1] + 3*x[3] + 4*x[4] )

for i in 1:scen
    bl = StructuredModel(parent=m)
    # @variable(bl, -200<=y1[1:3]<=200)
    @variable(bl, y1[1:3])

    @NLconstraint(bl, -100<= 0.1*x[1] + 0.2*x[2]                          <= 100)
    @NLconstraint(bl, -101<=               1.1*x[2] + 1.2*x[3]            <= 101)
    @NLconstraint(bl, -102<=                          2.1*x[3] + 2.2*x[4] <= 102)

    # @NLconstraint(bl,  9.1*x[1]           + 9.2*x[3]            -109 == 0)
    # @NLconstraint(bl,            8.1*x[2] + 8.2*x[3]            -108 == 0)
    # @NLconstraint(bl,                       7.1*x[3] + 7.2*x[4] -107 == 0)
    @NLconstraint(bl, 1*x[3] + 2*x[4]  == 107 )

    @NLconstraint(bl, -201<= 10.1*x[1] + 10.2*x[2]                           - 10.3*y1[1] + 10.4*y1[2]               <= 201)
    @NLconstraint(bl, -202<=               11.1*x[2] + 11.2*x[3]             - 11.3*y1[1] + 11.4*y1[2] + 11.5*y1[3]  <= 202)
    @NLconstraint(bl, -203<=                           12.1*x[3] + 12.2*x[4]              + 12.4*y1[2] - 12.5*y1[3]  <= 203)

    @NLobjective(bl, Min, 12*y1[2] + 13*y1[3])
end

structJuMPSolve(m)

getVarValue(m)


text: # LLᵀ (Choleski) factorization for symmetric systems


text:
"""
    struct DielectricDistribution{T, N, S, AT} <: AbstractArray{T, N}

Dielectric distribution, or distribution of the relative permittivity, ``\\epsilon_r``,
needed to calculate the [`ElectricPotential`](@ref).
The dielectric distribution is unitless.

## Parametric types
* `T`: Element type of `data`.
* `N`: Dimension of the `grid` and `data` array.
* `S`: Coordinate system (`Cartesian` or `Cylindrical`).
* `AT`: Axes type.

## Fields
* `data::Array{T, N}`: Array containing the values of the dielectric distribution at the discrete points of the `grid`.
* `grid::Grid{T, N, S, AT}`: [`Grid`](@ref) defining the discrete points at which the electric potential is determined.

!!! note
    The `data` array contains the values of the dielectric distribution at the discrete points
    between the points defined by the axes ticks of the extended grid of `grid`.\n
    Thus, `size(data) == size(grid) .+ 1` !
"""
struct DielectricDistribution{T, N, S, AT} <: AbstractArray{T, N}
    data::Array{T, N}
    grid::Grid{T, N, S, AT}
end

@inline size(ϵ::DielectricDistribution{T, N, S}) where {T, N, S} = size(ϵ.data)
@inline length(ϵ::DielectricDistribution{T, N, S}) where {T, N, S} = length(ϵ.data)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, I::Vararg{Int, N}) where {T, N, S} = getindex(ϵ.data, I...)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, i::Int) where {T, N, S} = getindex(ϵ.data, i)
@inline getindex(ϵ::DielectricDistribution{T, N, S}, s::Symbol) where {T, N, S} = getindex(ϵ.grid, s)

function DielectricDistribution(pss::PotentialSimulationSetup{T, N, S})::DielectricDistribution{T, N, S} where {T, N, S}
    return DielectricDistribution{T, N, S}( pss.ϵ_r, pss.grid )
end

function NamedTuple(ϵ::DielectricDistribution{T}) where {T <: SSDFloat}
    return (
        grid = NamedTuple(ϵ.grid),
        values = ϵ.data,
    )
end
Base.convert(T::Type{NamedTuple}, x::DielectricDistribution) = T(x)

function DielectricDistribution(nt::NamedTuple)
    grid = Grid(nt.grid)
    T = typeof(ustrip(nt.values[1]))
    S = get_coordinate_system(grid)
    N = get_number_of_dimensions(grid)
    DielectricDistribution{T, N, S, typeof(grid.axes)}( nt.values, grid )
end
Base.convert(T::Type{DielectricDistribution}, x::NamedTuple) = T(x)


text:
ps = [37/40]
refractive_indexes = Float64[1., 1., sqrt(3)]
x, y, Ez, dpml, dsource, resolution = simulation_hole_layers_unit_cell(ps, refractive_indexes=refractive_indexes)
δ = x[2] - x[1]
Ly_pml = y[end] - y[1] + δ

# test that the unit-cell is correct
@test resolution == 40
@test Ly_pml == 16

# test that some permittivity are 2 which is the desired pixel-averaging
geometry = real.(ϵ_hole_layers(x, y, ps, refractive_indexes=refractive_indexes))
@test any(isapprox.(geometry, 2.))


text:
#!/usr/bin/env julia
#
# Run _join_ benchmarks
#
# Usage: julia runtests.jl [--single-process] <llen> <rlen> <join_type>
#     --single-process      Runs all tests in a single julia process instead of
#                           spawning multiple julia subprocesses
#     <llen>                The length of the left dataframe
#     <rlen>                The length of the right dataframe
#     <join_type>           One of 'inner', 'left', 'right', 'outer', 'semi', 'anti'
#
# Note: Either provide all of <llen>, <rlen>, <join_type> arguments or leave blank
#       to run all test cases (see `run_all_tests` below)

using DataFrames
include("join_performance.jl")

# re-enable "cat" when
# https://github.com/JuliaData/CategoricalArrays.jl/pull/327
# and
# https://github.com/JuliaData/CategoricalArrays.jl/pull/331
# are merged and tagged

const file_loc = joinpath(dirname(@__FILE__), "join_performance.jl")

function run_all_tests(single_process)
    cases = [
        (100000, 50000000, "inner"),
        (5000000, 10000000, "inner"),
        (100000, 50000000, "left"),
        (5000000, 10000000, "left"),
        (100000, 50000000, "right"),
        (5000000, 10000000, "right"),
        (100000, 50000000, "outer"),
        (5000000, 10000000, "outer"),
        (100000, 50000000, "semi"),
        (5000000, 10000000, "semi"),
        (100000, 50000000, "anti"),
        (5000000, 10000000, "anti")
    ]

    for (llen, rlen, join_type) in cases
        run_test(single_process, llen, rlen, join_type)
    end
end

function run_test(single_process, llen, rlen, join_type)
    for
        arr_type in ["str", "int", "pool"], # "cat"],
        dup_mode in ["uniq", "dup", "manydup"],
        sort_mode in ["sort", "rand"],
        num_keys in ["1", "2"]
        dup_mode == "uniq" && arr_type in ["pool", "cat"] && continue

        if single_process
            @info "Running test: $llen $rlen $arr_type $dup_mode $sort_mode $num_keys $join_type"
            flush(stdout)

            t = run_bench(llen, rlen, arr_type, dup_mode, sort_mode, num_keys, join_type)

            @info "Timing: $t seconds"
            flush(stdout)
        else
            run(`julia $file_loc $llen $rlen $arr_type $dup_mode $sort_mode $num_keys $join_type`)
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    path = pathof(DataFrames)
    @info "DataFrames.jl at $path"

    # do some argument processing
    args = copy(ARGS)
    single_process = false
    pos = []
    while length(args) > 0
        arg = popfirst!(args)
        if arg == "--single-process"
            @info "Using single Julia process"
            global single_process = true
        elseif startswith(arg, "--")
            @error("Unknown option $arg")
            exit(1)
        else
            push!(pos, arg)
        end
    end

    @assert length(args) in [0, 3] "Wrong positional arguments. Expect [llen, rlen, join_type], got $args."

    if length(pos) == 0
        @info "Running all tests"
        run_all_tests(single_process)
    else
        llen, rlen, join_type = pos
        llen = parse(Int, llen)
        rlen = parse(Int, rlen)
        run_test(single_process, llen, rlen, join_type)
    end
end


text:
#   Unit classification.jl, part of PosDefManifold Package for julia language
#   v 0.1.3 - last update 7th of Mai 2019
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home
#
#   DESCRIPTION
#   This Unit implements classification methods and related functions
#   useful in relation to the Riemannian Geometry on the manifold
#   of Symmetric Positive Definite (SPD) or Hermitian matrices.
#
#   CONTENT
#   1. Probability
# __________________________________________________________________

#  ------------------------
## 1. Probability
#  ------------------------

"""
    softmax(χ::Vector{T}) where T<:Real

 Given a real vector of ``k`` non-negative scores ``χ=c_1,...,c_k``,
 return the vector ``π=p_1,...,p_k`` of their
 [softmax](https://en.wikipedia.org/wiki/Softmax_function) probabilities,
 as per

 `` p_i=\\frac{\\textrm{e}^{c_i}}{\\sum_{i=1}^{k}\\textrm{e}^{c_i}} ``.

 ## Examples
    χ=[1.0, 2.3, 0.4, 5.0]
    π=softmax(χ)

"""
softmax(χ::Vector{T}) where T<:Real = exp.(χ) ./ 𝚺(exp.(χ))


text:
include("../src/SOM.jl")
# using SOM

# using Test
#

using Distances
using ProgressMeter
using StatsBase
using Distributions
using NearestNeighbors
#using TensorToolbox
using LinearAlgebra

using PyPlot
using PyCall

using DataFrames
using RDatasets
using Primes

include("../src/errors.jl")

include("../src/types.jl")
include("../src/helpers.jl")
include("../src/grids.jl")
include("../src/kernels.jl")
include("../src/soms.jl")
include("../src/api.jl")

include("../src/plotTerminal.jl")

include("testFuns.jl")

include("../src/mongo.jl")
include("../src/namesPreprocessor.jl")

# include("benchmark.jl")

function generate_integer_hash_function(a,b,cardinality)
    # cardinality must be prime number
    # a,b: 1 ≤ a ≤ p−1  , 0 ≤ b ≤ p−1
    @assert a < cardinality
    @assert a > 0
    @assert b < cardinality
    @assert b >= 0
    @assert isprime(cardinality)
    m = cardinality
    x -> (x^a%m + b)%m  #TODO use modular exponentiation
end

function generate_vector_deterministic_pseudorandom(seed,hash_func,dimension)
    [hash_func(x+seed) for x in 1:dimension]
end

function generate_dataset(number,dimension,max_value,y=(x->x))
    # dimension = 16
    # max_value = 1013
    params = Dict(:a => 13, :b => 17, :hash_seed => 23, :prime => 1013)
    hash_func = generate_integer_hash_function(params[:a],params[:b],params[:prime])

    seed = params[:hash_seed]
    matrix = zeros((number,dimension))
    for i in 1:number
        vector = generate_vector_deterministic_pseudorandom(seed,hash_func,dimension)
        vector = map((x->x*max_value/params[:prime]),vector)
        vector = map((x->(x+max_value)/2),vector)
        vector = map(y,vector)
        seed = hash_func(seed+i)
        matrix[i,:] = vector
    end
    convert(DataFrame,matrix)
end

function benchmark_init(train, topol; toroidal = false,
                    normType = :zscore)

    xdim = 40
    ydim = 40

    initSOM(train, xdim, ydim, norm = normType,
                  topol = topol, toroidal = toroidal)
end

function benchmark_train(som,train, kernel = gaussianKernel)

    return trainSOM(som, train, 20000, kernelFun = kernel)
end

# query = generate_dataset(1000,4,1)
# query = generate_dataset(10000,4,3,x->x+1)
# iris = dataset("datasets", "iris")
# train = iris[:,1:4]
# #

# auto = dataset("ISLR","Auto")
# train = auto[:,1:7]

#

# train = copy(names_data)
#
# virgin_som = benchmark_init(train, :rectangular, toroidal = true)
# som = benchmark_train(virgin_som,train)

 # train = generate_dataset(1000,10,1)

 # som = benchmark_init(train, :rectangular, toroidal = true)
 # som = benchmark_train(som,train)

# test hexagonal, rectangular and spherical training:
#
# @time benchmarkTrain(train, :hexagonal, toroidal = false)
# @time benchmarkTrain(train, :rectangular, toroidal = false)
# @time benchmarkTrain(train, :spherical, toroidal = false)


text:
"""Usage:
```
For Unconstrained Optimsation:
  f = X->((x,y)=X;x^3 + 2y + 5)
  A = IntervalBox(2..4, 2..3)
  (global_min, minimisers, info) = ibc_minimise(f, A)
  (global_max, maximisers, info) = ibc_maximise(f, A)

For Constrained Optimisation:
  f = X->((x,y)=X;-(x-4)^2-(y-4)^2)
  A = IntervalBox(-4..4, -4..4)

  vars = ModelingToolkit.@variables x y
  C1 = Constraint(vars, x+y, -Inf..4)
  C2 = Constraint(vars, x+3y, -Inf..9)

  (global_min, minimisers, info) = ibc_minimise(f, A, [C1, C2])
  (global_max, maximisers, info) = ibc_maximise(f, A, [C1, C2])

ibc_minimise/ibc_maximise find the global minimum/maximum value of the function in given search space by using Interval Bound & Contract(IBC) algorithm
```
"""
function ibc_minimise(f::Function , X::IntervalBox{N,T}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, T}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), debug = false, structure = SortedVector, tol=1e-6, ibc_ind = true) where{N, T}

    vars = [Variable(Symbol("x",i))() for i in 1:length(X)]
    g(x...) = f(x)
    C = BasicContractor(vars, g)

    working = structure([(X, inf(f(X)))], x->x[2]) # list of boxes with corresponding lower bound, arranged according to selected structure :

    minimizers = IntervalBox{N,T}[]
    global_min = ∞  # upper bound

    info = Information(0, 0, 0)
    num_bisections = 0

    while !isempty(working)

        info.iterations= info.iterations + 1

        if !ibc_ind
            if isready(ibc_chnl)
                from_diff = take!(ibc_chnl)     # Receiving best individual from ibc_minimise
                if debug
                    println("Box recevied from DifferentialEvolution: ", from_diff)
                end
                global_min = min(from_diff[2], global_min)
                info.de_to_ibc = info.de_to_ibc + 1
            end
       end

        (X, X_min) = popfirst!(working)

        gra = gradient(f, X.v)
        (lbb, fcb, cb) = bauman_form(X, f, gra)      # ----- second order form

        if global_min < lbb
            continue
        else
            X_min = max(X_min, lbb)
            if fcb < global_min
                global_min = fcb
                x_best = cb
            end
        end                                        # ------

        if debug
            println("New search-space : ", X)
        end

        A = -∞..global_min
        X = invokelatest(C, A, X)                        # Contracting the box by constraint f(X) < globla_min
        X_min = inf(f(X))

        if debug
            println("Contracted search_space: ", X)
        end

        if X_min > global_min    # X_min == inf(f(X))
            continue
        end

        # find candidate for upper bound of global minimum by just evaluating a point in the interval:
        m = sup(f(Interval.(mid.(X))))   # evaluate at midpoint of current interval

        if m < global_min
            global_min = m
            x_best = SVector(mid(X))
            if !ibc_ind
                if debug
                    println("Box send to DifferentialEvolution: ", x_best )
                end
                if info.iterations % 200 == 0
                    put!(diffevol_chnl, (x_best, global_min, X))  # sending best individual to diffevol
                else
                    put!(diffevol_chnl, (x_best, global_min, nothing))
                end
                info.ibc_to_de = info.ibc_to_de + 1
            end
        end

        filter_elements!(working , (X, global_min) )   # Remove all boxes whose lower bound is greater than the current one:

        if diam(X) < tol
            push!(minimizers, X)
        else
            X1, X2 = bisect(X)
            push!( working, (X1, inf(f(X1))) )
            push!( working, (X2, inf(f(X2))) )
            num_bisections += 1
        end
    end

    if debug
        println("IBC search is ended")
    end

    if !ibc_ind
        if isready(diffevol_chnl)
            take!(diffevol_chnl)
            put!(diffevol_chnl,(SVector(mid(X)), Inf, nothing) )
        else
            put!(diffevol_chnl,(SVector(mid(X)), Inf, nothing) )
        end
        if debug
            println("DifferentialEvolution is also terminated")
        end
        if isready(ibc_chnl)
            take!(ibc_chnl)
        end
    end

    lower_bound = minimum(inf.(f.(minimizers)))

    return Interval(lower_bound,global_min), minimizers, info

end

function ibc_maximise(f::Function , X::IntervalBox{N,T}; ibc_chnl = RemoteChannel(()->Channel{Tuple{IntervalBox{N,T}, T}}(0)),
               diffevol_chnl = RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(0)), debug = false, structure = SortedVector, tol=1e-6, ibc_ind = true) where{N, T}
    bound, minimizer, info = ibc_minimise(x -> -f(x), X, ibc_chnl = ibc_chnl, diffevol_chnl = diffevol_chnl, debug = debug, structure = structure, tol = tol, ibc_ind = true)
    return -bound, minimizer, info
end

max_issues_repo_name: yashcodes/Charibde

max_forks_repo_name: yashcodes/Charibde


text:
module FoamWorld
using Gtk
using StaticArrays
using ColorTypes
using FixedPointNumbers:N0f8
# 宏变量
include("consts.jl")
include("globals.jl")

# 功能
include("lib_collide.jl")
include("lib_draw.jl")
include("lib_IB.jl")
include("lib_rand.jl")
include("lib_string.jl")
include("codes.jl")

# 内容
include("chunk.jl")
include("saves.jl")
include("tick.jl")
include("operation.jl")

# 细节
include("blk/blk.jl")
include("cmd/cmd.jl")
include("ent/ent.jl")
include("gen/gen.jl")
include("itm/itm.jl")
include("mat/mat.jl")
include("sta/sta.jl")

end # module


text:
module OpenBLAS32

using LinearAlgebra.BLAS
using OpenBLAS32_jll

__init__() = BLAS.lbt_forward(OpenBLAS32_jll.libopenblas_path)

end # module


text:
##############################
## by Qin Yu, Nov 2018
## using Julia 1.0.1
##############################

# If you have not installed these packages:
# using Pkg # Julia 1.0.1 only
# Pkg.add("GR")
# Pkg.add("Plots")
# Pkg.add("PyPlot")
# Pkg.add("Distributions")

# I will not include Pkg.add info in future .jl files,
# add it if you can't be "using" any.

############################## Starting here:
using LinearAlgebra
using Plots
using Random
using Distributions
# Load my own package:
push!(LOAD_PATH, ".")
using SuperLearn

############################## Change Basis, Run:
# Core.eval(SuperLearn, :(POLY_OR_SINE = "poly"))  # Change to Polynomial Basis (default)
# Core.eval(SuperLearn, :(POLY_OR_SINE = "sine"))  # Change to Sine Basis

############################## Another Simple Dataset:
# random function 𝑔_σ(𝑥) := sin²(2π𝑥) + ϵ, where random var ϵ ~ 𝐍(0, σ²)
# sample size n = 30, iid random variable 𝑋ᵢ ~ 𝐔(0, 1)
# srand(123)  # Julia 0.6.4
Random.seed!(777)  # Julia 1.0.1

# Define a distribution that we are sampling from, and sample:
𝑥ᵢ = Uniform(0, 1)
𝒙 = rand(𝑥ᵢ, 30)
sort!(𝒙)

# Plot the function sin²(2πx) in the range 0 ≤ x ≤ 1
f(x) = (sin(2 * π * x))^2
plot(f, 0, 1, lab="sin²(2 pi x)")

# with the points of the above data set superimposed:
g(σ, x) = f(x) + rand(Normal(0, σ), 1)[1]  # add noise when we sample
g₀₀₇(x) = g(0.07, x)
𝒚 = vcat(g₀₀₇.(𝒙)...)
scatter!(𝒙, 𝒚, lab="S")

savefig("./graph/2.1.pdf")

############################## Train 20 regressions:
# I am using trained_regression_line() from my own package,
# which is extracted from 1.Linear_Regression_Basis_Function.jl
regression_curves = [x -> trained_regression_line(x, 𝒙, 𝒚, i) for i in 1:20]

# Plot selected learned linear regression curves (every third):
for i = 1:3:size(regression_curves, 1)
    display(plot(regression_curves[i], 0, 1, lab="k = $i"))
    display(scatter!(𝒙, 𝒚, lab="S"))
    savefig("./graph/2.2-$i.pdf")
end

############################## Training Errors for Different Sizes of Basis:
te_k(k) = training_error_k_dim_basis(𝒙, 𝒚, k)
te_k_1to20 = te_k.(1:20)
log_te_k_1to20 = log.(te_k_1to20)
plot(1:20, log_te_k_1to20, xlabel="k", xticks=0:20, ylabel="log training MSE", lab="log(te(k, S))")

############################## Testing Error with Testing Set of Size 1000:
𝒙_test = rand(𝑥ᵢ, 1000)
sort!(𝒙_test)
𝒚_test = vcat(g₀₀₇.(𝒙_test)...)
scatter(𝒙_test, 𝒚_test)

tse_k(k) = test_error_k_dim_basis(𝒙_test, 𝒚_test, 𝒙, 𝒚, k)
tse_k_1to20 = tse_k.(1:20)
log_tse_k_1to20 = log.(tse_k_1to20)
plot(1:20, log_te_k_1to20, xlabel="k", xticks=0:20, ylabel="log MSE", lab="log(te(k, S))")
plot!(1:20, log_tse_k_1to20, lab="log(tse(k, S, T))")
savefig("./graph/2.3.pdf")

############################## Obtain Previous Result 100 Times and Get Log-Average:
sum_all_100_te = fill(0, 20)
sum_all_100_tse = fill(0, 20)
for n = 1:100
    𝒙 = rand(𝑥ᵢ, 30)  # increasing the dimension of training set REDUCEs the error
    sort!(𝒙)
    𝒚 = vcat(g₀₀₇.(𝒙)...)
    te_k(k) = training_error_k_dim_basis(𝒙, 𝒚, k)
    te_k_1to20 = te_k.(1:20)
    global sum_all_100_te += te_k.(1:20)

    𝒙_test = rand(𝑥ᵢ, 1000)
    sort!(𝒙_test)
    𝒚_test = vcat(g₀₀₇.(𝒙_test)...)
    tse_k(k) = test_error_k_dim_basis(𝒙_test, 𝒚_test, 𝒙, 𝒚, k)
    tse_k_1to20 = tse_k.(1:20)
    global sum_all_100_tse += tse_k.(1:20)
end

avg_all_100_te = sum_all_100_te / 100
avg_all_100_tse = sum_all_100_tse / 100
log_all_100_te = log.(avg_all_100_te)
log_all_100_tse = log.(avg_all_100_tse)
plot(1:20, log_all_100_te, xlabel="k", xticks=0:20, ylabel="log average MSE", lab="log(avg(te(k, S)))")
plot!(1:20, log_all_100_tse, lab="log(avg(tse(k, S, T)))")  # run this section again, the larger k tail varies
savefig("./graph/2.4.pdf")

############################## To Change Polynomial Basis to Sine Basis, Run:
Core.eval(SuperLearn, :(POLY_OR_SINE = "sine"))
# Changed the basis, run everything again.

max_stars_repo_name: qin-yu/julia-regression-boston-housing

max_issues_repo_name: qin-yu/ml-julia-boston-housing

max_forks_repo_name: qin-yu/ml-julia-boston-housing


text:
function (method::GuaranteedEllipseFit)(observations::AbstractObservations)
    @unpack parametrisation = method
    @unpack optimisation_scheme = method
    @unpack seed = optimisation_scheme

    𝛈₀ = to_latent_parameters(parametrisation, seed.𝛉)
    optimisation_scheme = @set optimisation_scheme.seed = ManualEstimation(𝛈₀)
    N = length(observations.data[1])
    jacobian_matrix = JacobianMatrix(observations, zeros(N,5))

    optimisation_result = optimisation_scheme(observations, jacobian_matrix)
    𝛈₁ = optimisation_result.minimiser

    if (method.verbose)
        cost = optimisation_result.minimum
        info_str = optimisation_result.info_str
        @info "AML Cost: $cost"
        @info info_str
    end

    # TODO Remove parametrisation option and just always use "SecondLatentEllipseParametrisation".
    # The lm method below assumes "SecondLatentEllipseParametrisation".
    𝛉 = from_latent_parameters(parametrisation, 𝛈₁)
    return 𝛉 / norm(𝛉)
end

function (lm::LevenbergMarquardt)(observations::AbstractObservations, jacobian_matrix::JacobianMatrix)
    @unpack data = observations
    @unpack γ, γ_inc, γ_dec = lm
    @unpack λ, λ_lower_bound = lm
    @unpack tol_𝛉, tol_cost, tol_Δ, tol_∇ = lm
    @unpack max_iter = lm
    @unpack seed = lm
    @unpack max_func_eval = lm

    # minimum allowable magnitude of conic determinant to prevent ellipse from
    # convering on a degenarate parabola (e.g. two parallel lines).
    tol_detD = 1e-5
    # barrier tolerance tp prevent ellipse from convering on parabola.
    tol_barrier = 15.5

    𝛈ₖ = seed.𝛉

    info_str = ""
    keep_going = true
    theta_updated = false
    func_eval = 0
    k = 1

    ∂𝐫 = jacobian_matrix
    𝛉ₖ = from_latent_parameters(SecondLatentEllipseParametrisation(), 𝛈ₖ)
    Δₖ = ones(length(𝛉ₖ))

    N = length(first(data))
    𝐫 = vector_valued_objective(observations, 𝛉ₖ)
    # Sum of squared residuals: 𝐫'*𝐫
    costₖ = sum(abs2, 𝐫)

    𝐈 = UniformScaling(1)
    𝐅₀ = SA_F64[0 0 2; 0 -1 0; 2 0 0]
    𝐅 = SA_F64[1 0; 0 0] ⊗ 𝐅₀

    was_updated = false
    while keep_going && k < max_iter

        𝛉ₖ = from_latent_parameters(SecondLatentEllipseParametrisation(), 𝛈ₖ)
        𝐫 = vector_valued_objective(observations, 𝛉ₖ)
        𝐉 = ∂𝐫(𝛈ₖ)
        𝐇 = 𝐉'*𝐉

        𝐏ₜ = UniformScaling(1) - (𝛉ₖ*𝛉ₖ')/norm(𝛉ₖ)^2
        ∂π = norm(𝛉ₖ)^-1 * 𝐏ₜ *  ∂𝛋(𝛈ₖ)
        𝐖 = (∂π'*∂π)

        # Compute two potential updates based on different weightings of 𝐖.
        λ₁ = max(λ,  λ_lower_bound)
        Δ₁ = -(𝐇 + λ₁*𝐖) \ (𝐉'*𝐫)
        𝛈₁ = 𝛈ₖ + Δ₁

        λ₂ = max(λ / γ,  λ_lower_bound)
        Δ₂ = -(𝐇 + λ₂*𝐖) \ (𝐉'*𝐫)
        𝛈₂ = 𝛈ₖ + Δ₂

        𝛉₁ = from_latent_parameters(SecondLatentEllipseParametrisation(), 𝛈₁)
        𝛉₂ = from_latent_parameters(SecondLatentEllipseParametrisation(), 𝛈₂)

        # Compute new residuals and costs based on these updates.
        𝐫₁ = vector_valued_objective(observations, 𝛉₁)
        cost₁ = sum(abs2, 𝐫₁)
        𝐫₂ = vector_valued_objective(observations, 𝛉₂)
        cost₂ = sum(abs2, 𝐫₂)

        if cost₁ >= costₖ && cost₂ >= costₖ
            # Neither potential update reduced the cost.
            was_updated = false
            𝛉ₖ₊₁ = 𝛉ₖ
            𝛈ₖ₊₁ = 𝛈ₖ
            Δₖ₊₁ = Δₖ
            costₖ₊₁ = costₖ
            λ = λ * γ # In the next iteration add more of the identity matrix.
            func_eval = func_eval + 1
        elseif cost₂ < costₖ
            # Update (2) reduced the cost function.
            was_updated = true
            𝛈ₖ₊₁ = 𝛈₂
            𝛉ₖ₊₁ = 𝛉₂
            Δₖ₊₁ = Δ₂
            costₖ₊₁ = cost₂
            λ = λ / γ  # In the next iteration add less of the identity matrix.
        else
            # Update (1) reduced the cost function.
            was_updated = true
            𝛈ₖ₊₁ = 𝛈₁
            𝛉ₖ₊₁ = 𝛉₁
            Δₖ₊₁ = Δ₁
            costₖ₊₁ = cost₁
            λ = λ  # Keep the same damping for the next iteration.
        end

        barrier = (𝛉ₖ₊₁'*𝐈*𝛉ₖ₊₁)/(𝛉ₖ₊₁'*𝐅*𝛉ₖ₊₁)
        D = SA_F64[𝛉ₖ₊₁[1] 𝛉ₖ₊₁[2]/2 𝛉ₖ₊₁[4]/2 ;
                   𝛉ₖ₊₁[2]/2 𝛉ₖ₊₁[3] 𝛉ₖ₊₁[5]/2 ;
                   𝛉ₖ₊₁[4]/2 𝛉ₖ₊₁[5]/2 𝛉ₖ₊₁[6]]

        detD = det(D)

        # Since 𝛈 is a projective entity this converge criterion will have
        # to change to take into account the scale/sign ambiguity.
        if min(norm(𝛈ₖ₊₁ - 𝛈ₖ), norm(𝛈ₖ₊₁ + 𝛈ₖ)) < tol_𝛉 && was_updated
            info_str = "Breaking because of tolerance."
            keep_going = false
        elseif abs(costₖ₊₁ - costₖ) < tol_cost && was_updated
            info_str = "Breaking because of cost."
            keep_going = false
        elseif norm(Δₖ₊₁) < tol_Δ
            info_str = "Breaking because of update norm."
            keep_going = false
        elseif  norm(𝐉'*𝐫, Inf) < tol_∇
            info_str = "Breaking because of gradient norm."
            keep_going = false
        elseif func_eval > max_func_eval
            info_str = "Breaking because maximum func evaluations reached."
            keep_going = false
        elseif log(barrier) > tol_barrier || abs(detD) < tol_detD
            info_str = "Breaking because approaching degenerate ellipse."
            keep_going = false
        end

        𝛈ₖ = 𝛈ₖ₊₁
        Δₖ = Δₖ₊₁
        costₖ = costₖ₊₁
        k = was_updated ? k + 1 : k
    end

    return OptimisationResult(𝛈ₖ, costₖ, info_str)
end

function vector_valued_objective(observations::AbstractObservations, 𝛉::AbstractVector)
    @unpack data = observations
    ℳ = data[1]
    N = length(ℳ)
    𝐫 = zeros(N)
    for n = 1:N
        𝐦 = ℳ[n]
        𝐮ₙ = SA_F64[𝐦[1]^2, 𝐦[1]*𝐦[2], 𝐦[2]^2, 𝐦[1], 𝐦[2], 1]
        ∂𝐮ₙ = SA_F64[2*𝐦[1]  𝐦[2]  0  1  0  0; 0 𝐦[1] 2*𝐦[2] 0 1 0]'
        𝐀ₙ = 𝐮ₙ * 𝐮ₙ'
        𝚲ₙ =  SA_F64[1 0 ; 0 1]
        𝐁ₙ = ∂𝐮ₙ * 𝚲ₙ * ∂𝐮ₙ'
        𝐫[n] = sqrt(abs((𝛉' * 𝐀ₙ * 𝛉)/(𝛉' * 𝐁ₙ * 𝛉)))
    end
    return 𝐫
end

function to_latent_parameters(::FirstLatentEllipseParametrisation, 𝛉::AbstractVector)
    p = 𝛉[2] / (2*𝛉[1])
    q = ((𝛉[3] / 𝛉[1]) - (𝛉[2]/(2*𝛉[1]))^2)^(-0.5)
    r = (𝛉[4] / 𝛉[1])
    s = (𝛉[5] / 𝛉[1])
    t = (𝛉[6] / 𝛉[1])
    return SVector(p, q, r, s, t)
end

function from_latent_parameters(::FirstLatentEllipseParametrisation, 𝛋::AbstractVector)
    p, q, r, s, t = 𝛋
    a = 1
    b = 2*p
    c = p^2 + q^(-2)
    d = r
    e = s
    f = t
    𝛉 = SVector(a, b, c, d, e ,f)
    return 𝛉 / norm(𝛉)
end

function to_latent_parameters(::SecondLatentEllipseParametrisation, 𝛉::AbstractVector)
    p = 𝛉[2] / (2*𝛉[1])
    q = ((𝛉[3] / 𝛉[1]) - (𝛉[2]/(2*𝛉[1]))^2)^(0.5)
    r = (𝛉[4] / 𝛉[1])
    s = (𝛉[5] / 𝛉[1])
    t = (𝛉[6] / 𝛉[1])
    return SVector(p, q, r, s, t)
end

function from_latent_parameters(::SecondLatentEllipseParametrisation, 𝛋::AbstractVector)
    p, q, r, s, t = 𝛋
    a = 1
    b = 2*p
    c = p^2 + q^(2)
    d = r
    e = s
    f = t
    𝛉 = SVector(a, b, c, d, e ,f)
    return 𝛉 / norm(𝛉)
end

function ∂𝛋(𝛈)
    return SA_F64[0      0       0 0 0 ;
                  2      0       0 0 0 ;
                  2*𝛈[1] 2*𝛈[2]  0 0 0 ;
                  0      0       1 0 0 ;
                  0      0       0 1 0 ;
                  0      0       0 0 1]
end

function (jacobian_matrix::JacobianMatrix)(𝛈::AbstractVector)
    @unpack observations = jacobian_matrix
    @unpack data = observations
    ℳ = data[1]
    N = length(ℳ)

    𝛉 = from_latent_parameters(SecondLatentEllipseParametrisation(), 𝛈)

    𝐏ₜ = UniformScaling(1) - (𝛉*𝛉')/norm(𝛉)^2
    ∂π = norm(𝛉)^-1 * 𝐏ₜ *  ∂𝛋(𝛈)
    # TODO overwrite the pre-allocated array instead (for performance).
    ∂𝐫′ = zeros(N, 5)
    for n = 1:N
        𝐦 = ℳ[n]
        𝐮ₙ = SA_F64[𝐦[1]^2, 𝐦[1]*𝐦[2], 𝐦[2]^2, 𝐦[1], 𝐦[2], 1]
        ∂𝐮ₙ = SA_F64[2*𝐦[1]  𝐦[2]  0  1  0  0; 0 𝐦[1] 2*𝐦[2] 0 1 0]'
        𝐀ₙ = 𝐮ₙ * 𝐮ₙ'
        𝚲ₙ =  SA_F64[1 0 ; 0 1]
        𝐁ₙ = ∂𝐮ₙ * 𝚲ₙ * ∂𝐮ₙ'
        𝐌ₙ = 𝐀ₙ / (𝛉' * 𝐁ₙ * 𝛉)
        𝐗ₙ =  𝐌ₙ - 𝐁ₙ * ((𝛉' * 𝐀ₙ * 𝛉)/ (𝛉' * 𝐁ₙ * 𝛉)^2)
        ∂𝐫 = (𝐗ₙ*𝛉 / sqrt(abs((𝛉' * 𝐀ₙ * 𝛉) / (𝛉' * 𝐁ₙ * 𝛉)) + eps()))'
        ∂𝐫′[n,:] = ∂𝐫 * ∂π
     end
    return ∂𝐫′
end


text:
# ---------------------------------------------------------------------------------------------
#   Solve covariance estimation problem
#           min_{X \in \R^{d x r}}   F(X)   :=    E_{a,b} | <XX^T, aa^T> - b |
#   using either 1) subgradient method or
#                       2) mirror descent (with a polynomial mirror map)
#
#   Inputs:     X0 = initialization (d x r matrix)
#                   Xtrue = true matrix (used to plot dist. to solution)
#                   stochErr = variance of stochastic errors in residuals <XX^T,aa^T> - b
#                   maxIter = maximum no. iterations
#                   stepSizes = one-dim array of stepsizes (length maxIter)
#                   method = "subgradient" or "mirror"
#                   clipped = whether to use clipping (true/false)
#                   verbose = whether to print output (true/false)
#
#   Outputs:   (err_hist, fun_hist) = (array of normalized distances to solution,
#                                       array of absolute errors in empirical function values)
# ---------------------------------------------------------------------------------------------
include("func.jl");
using Printf

function solve_cov_est(X0::Array{Float64,2},
                                     Xtrue::Array{Float64,2},
                                     stochErr::Float64,
                                     maxIter::Int64,
                                     stepSizes::Array{Float64,1};
                                     method::String="subgradient",
                                     clipped::Bool=false,
                                     verbose::Bool=true
                                    )
    # Extract basic data
    (d,r) = size(X0);
    XTtrue = Xtrue';
    sqnrmXtrue = sum(abs2, Xtrue);

    #  Initialize
    X = copy(X0);
    a = zeros(d);
    XTa = zeros(r);
    b = 0;
    res = 0;  # residual <XX^T, aa^T> - b
    G = zeros(d,r);  # subgradient at X
    η = 0;  # stepsize

    if method=="mirror"
        # Set coefficients for the Bregman divergence polynomials
        #       p(u) = a0 + a1 * u + a2 * u²
        #       Φ(x) = c0 ||x||₂² + c1 ||x||₂³ + c2 ||x||₂⁴   (mirror map)
        (a0, a1, a2) = (1.0, 0.0, 1.0);  # chosen by us
        (c0, c1, c2) = (a0*7/2,  a1*10/3,  a2*13/4);  # set according to theorem

        # Initialize more quantities
        Y = copy(X0);  # need an extra copy of current iterate
        V = zeros(d,r);  # V = ∇Φ(X) - η * G
        λ = 0;  # λ = ||X||
    end

    # Draw stochastic a, b for all iterations at once
    (A,B) = get_ab(XTtrue, stochErr, maxIter);
    true_empirical_value = compute_empirical_function(Xtrue, A, B);

    # Initialize vectors for tracking progress
    err = NaN;  # distance to solution
    err_hist =  fill(NaN, maxIter);  # distances to solution
    fun_val = compute_empirical_function(X, A, B);  # empirical function value
    fun_hist =  fill(NaN, maxIter);  # empirical function values

    #   Run the method (stochatic subgradient or stochastic mirror)
    for k=1:maxIter
        # Get stochastic a, b
        a = A[:,k];
        b = B[k];
        XTa = X'*a;

        # Update subgradient in place
        subgrad!(G, XTa, a, b);
        thresh = dot(G, X) - fun_val;  # clipping threshold

        # Update X in place
        if method=="subgradient"
            # Get step size
            if clipped
                η = min(stepSizes[k], fun_val/sum(abs2,G));
            else
                η = stepSizes[k];
            end
            BLAS.axpy!(-η,G,X);
        elseif method=="mirror"
            # Update V = ∇Φ(X) - η * G
            η = stepSizes[k];
            gradPhi = ( 2*c0 + 3*c1*norm(X,2) + 4*c2*sum(abs2, X) ) * X;
            V = gradPhi  - η*G;
            nrmV = norm(V,2);

            # Note that direction of X is identical to direction of V.
            # Root finding problem to find λ = ||X|| is
            #    (2*c0) * λ + (3*c1) * λ^2 + (4*c2) * λ^3  = nrmV
            λ = get_root([-nrmV, 2*c0, 3*c1, 4*c2]);

            # Record intermediate matrix Y
            lmul!(0.0, Y)  # reset to zero
            BLAS.axpy!(λ/nrmV, V, Y)

            if (! clipped) | (dot(G,Y) >= thresh)
                X = Y;
            else  # Need to find a better Y.
                # First try: ∇ϕ(Y) = ∇ϕ(X_k), i.e. Lagrange multiplier is zero
                nrmGrad= norm(gradPhi,2);
                λ = get_root([-nrmGrad, 2*c0, 3*c1, 4*c2]);
                lmul!(0.0, Y)  # reset to zero
                BLAS.axpy!(λ/nrmGrad, gradPhi, Y)

                if (dot(G,Y) <= thresh)
                    X = Y;
                else  # Y still doesn't satisfy clipped constraint
                    # Second try: use gradPhi = W + α*G,  W ⟂ G,  and ∇Φ(Y) = W + (α-μ)*G
                    sqnormG = sum(abs2,G);
                    α = dot(G, gradPhi) / sqnormG;
                    W = gradPhi - α*G;

                    #  Inner product with G + norm squared gives coefficients
                    d0, d1, d2 = 2*c0, 3*c1, 4*c2;
                    d3 = thresh^2/sqnormG;

                    # Now solve for λ=||Y|| using (a + bλ + cλ^2)^2 (λ^2 - d) - ||W||^2 = 0
                    λ = get_root([ -d0^2*d3 - sum(abs2,W),
                                        -(2*d0*d1)*d3,
                                        d0^2 - (2*d0*d2-d1^2)*d3,
                                        (2*d0*d1) - (2*d1*d2)*d3,
                                        (2*d0*d2 + d1^2) - d2^2*d3,
                                        2*d1*d2,
                                        d2^2
                                        ]);

                    # Get μ from λ
                    μ = α  - (d0+d1*λ+d2*λ^2) * thresh / sqnormG ;
                    if μ > 0
                        V = W + (α-μ)*G;
                        Y = (λ/nrmV) * V;
                        X = Y;
                    else
                        # Can't find a good Y, so don't take a step at all. Keep X as is.
                    end
                end
            end
        end

        # record error status and print to console
        err = sqnrmXtrue + sum(abs2, X) - 2 * sum(svdvals(XTtrue * X));
        normalized_err = err / sqnrmXtrue;
        err_hist[k] = normalized_err;

        fun_val = compute_empirical_function(X, A, B);
        abs_value_error = fun_val - true_empirical_value ;
        fun_hist[k] = abs_value_error;

        if verbose
            @printf("iter %3d: function err = %1.2e, dist err = %1.2e, stepsize = %1.2e\n",
                            k, value_error, normalized_err, η);
        end
    end

    return (err_hist, fun_hist)
end

max_stars_repo_name: kelliemac/stochastic-model-based

max_issues_repo_name: kelliemac/stochastic-model-based

max_forks_repo_name: kelliemac/stochastic-model-based


text:
function initialize!(integrator, cache::DImplicitEulerConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.du
end

function initialize!(integrator, cache::DImplicitEulerCache)
  integrator.kshortsize = 2
  @unpack k₁,k₂ = cache
  resize!(integrator.k, integrator.kshortsize)
  integrator.k .= [k₁,k₂]
  integrator.k[1] .= integrator.du
  nothing
end

@muladd function perform_step!(integrator, cache::DImplicitEulerConstantCache, repeat_step=false)
  @unpack t,dt,uprev,u,f,p = integrator
  alg = unwrap_alg(integrator, true)
  @unpack nlsolver = cache

  nlsolver.z = zero(u)
  nlsolver.tmp = zero(u)
  nlsolver.γ = 1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  u = uprev + z

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    tmp = r*integrator.opts.internalnorm.((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    atmp = calculate_residuals(tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  integrator.u = u
  integrator.du = (u-uprev)/dt

  if integrator.opts.calck
    integrator.k[2] = integrator.k[1]
    integrator.k[1] = integrator.du
  end
end

@muladd function perform_step!(integrator, cache::DImplicitEulerCache, repeat_step=false)
  @unpack t,dt,uprev,du,u,f,p = integrator
  @unpack atmp,nlsolver = cache
  @unpack tmp = nlsolver
  alg = unwrap_alg(integrator, true)

  @. nlsolver.z = false
  @. nlsolver.tmp = false
  nlsolver.γ = 1
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return
  @.. u = uprev + z
  @.. du = z * inv(dt)

  if integrator.opts.adaptive && integrator.success_iter > 0
    # local truncation error (LTE) bound by dt^2/2*max|y''(t)|
    # use 2nd divided differences (DD) a la SPICE and Shampine

    # TODO: check numerical stability
    uprev2 = integrator.uprev2
    tprev = integrator.tprev

    dt1 = dt*(t+dt-tprev)
    dt2 = (t-tprev)*(t+dt-tprev)
    c = 7/12 # default correction factor in SPICE (LTE overestimated by DD)
    r = c*dt^2 # by mean value theorem 2nd DD equals y''(s)/2 for some s

    @.. tmp = r*integrator.opts.internalnorm((u - uprev)/dt1 - (uprev - uprev2)/dt2,t)
    calculate_residuals!(atmp, tmp, uprev, u, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  else
    integrator.EEst = 1
  end

  if integrator.opts.calck
    integrator.k[2] .= integrator.k[1]
    integrator.k[1] .= du
  end
end

function initialize!(integrator, cache::DABDF2ConstantCache)
  integrator.kshortsize = 2
  integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
  integrator.k[1] = integrator.du
end

@muladd function perform_step!(integrator, cache::DABDF2ConstantCache, repeat_step=false)
  @unpack t,f,p = integrator
  @unpack dtₙ₋₁,nlsolver = cache
  alg = unwrap_alg(integrator, true)
  dtₙ, uₙ, uₙ₋₁, uₙ₋₂ = integrator.dt, integrator.u, integrator.uprev, integrator.uprev2

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    integrator.fsalfirst = @.. (integrator.u - integrator.uprev) / dtₙ
    cache.fsalfirstprev = integrator.fsalfirst
    return
  end

  # precalculations
  ρ = dtₙ/dtₙ₋₁
  c1 = ρ^2/(1+2ρ)

  nlsolver.γ = (1+ρ)/(1+2ρ)
  nlsolver.α = 1//1

  nlsolver.z = zero(uₙ)

  nlsolver.tmp = -c1 * uₙ₋₁ + c1 * uₙ₋₂
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  uₙ = uₙ₋₁ + z
  integrator.fsallast = @.. z/dtₙ

  if integrator.opts.adaptive
    tmp = integrator.fsallast - (1+dtₙ/dtₙ₋₁)*integrator.fsalfirst + (dtₙ/dtₙ₋₁)*cache.fsalfirstprev
    est = (dtₙ₋₁+dtₙ)/6 * tmp
    atmp = calculate_residuals(est, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end

  integrator.u = uₙ
  integrator.du = du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dtₙ)

  if integrator.opts.calck
    integrator.k[2] = integrator.k[1]
    integrator.k[1] = integrator.du
  end
  return
end

function initialize!(integrator, cache::DABDF2Cache)
  integrator.fsalfirst = cache.fsalfirst
  integrator.fsallast = du_alias_or_new(cache.nlsolver, integrator.fsalfirst)

  integrator.kshortsize = 2
  @unpack k₁,k₂ = cache.eulercache
  resize!(integrator.k, integrator.kshortsize)
  integrator.k .= [k₁,k₂]
  integrator.k[1] .= integrator.du
  nothing
end

@muladd function perform_step!(integrator, cache::DABDF2Cache, repeat_step=false)
  @unpack t,dt,du,f,p = integrator
  @unpack atmp,dtₙ₋₁,nlsolver = cache
  @unpack z,tmp = nlsolver
  alg = unwrap_alg(integrator, true)
  uₙ,uₙ₋₁,uₙ₋₂,dtₙ = integrator.u,integrator.uprev,integrator.uprev2,integrator.dt

  if integrator.iter == 1 && !integrator.u_modified
    cache.dtₙ₋₁ = dtₙ
    perform_step!(integrator, cache.eulercache, repeat_step)
    @.. integrator.fsalfirst = (uₙ - uₙ₋₁) / dt
    cache.fsalfirstprev .= integrator.fsalfirst
    return
  end

  # precalculations
  ρ = dtₙ/dtₙ₋₁
  c1 = ρ^2/(1+2ρ)

  nlsolver.γ = (1+ρ)/(1+2ρ)
  nlsolver.α = 1//1
  @.. nlsolver.tmp = -c1 * uₙ₋₁ + c1 * uₙ₋₂
  nlsolver.z .= zero(eltype(z))
  z = nlsolve!(nlsolver, integrator, cache, repeat_step)
  nlsolvefail(nlsolver) && return

  @.. uₙ = uₙ₋₁ + z
  @.. du = (nlsolver.α * z + nlsolver.tmp) * inv(nlsolver.γ * dt)

  @.. integrator.fsallast = du
  integrator.destats.nf += 1
  if integrator.opts.adaptive
    btilde0 = (dtₙ₋₁+dtₙ)*1//6
    btilde1 = 1+ρ
    btilde2 = ρ
    @.. tmp = btilde0*(integrator.fsallast - btilde1*integrator.fsalfirst + btilde2*cache.fsalfirstprev)
    calculate_residuals!(atmp, tmp, uₙ₋₁, uₙ, integrator.opts.abstol, integrator.opts.reltol,integrator.opts.internalnorm,t)
    integrator.EEst = integrator.opts.internalnorm(atmp,t)
  end

  ################################### Finalize

  if integrator.EEst < one(integrator.EEst)
    @.. cache.fsalfirstprev = integrator.fsalfirst
    cache.dtₙ₋₁ = dtₙ
  end
  return
end


text:
using Base.Test
import WaterData

"Resources for testing thermal expansivity"
module test_thermalexpansivity_resources

import WaterData

""" A simple EOS to check thermal expansivity determination.

    ρ(P, T) = P^n * T^m so that αᵥ(P, T) = -m/T """
struct SimpleEOS <: WaterData.EOS
    n::Int
    m::Int
end
(e::SimpleEOS)(P, T) = P^(e.n) * T^(e.m)

end # module test_thermalexpansivity_resources

res = test_thermalexpansivity_resources

@testset "Thermal expansivity" begin
    e1 = res.SimpleEOS(1, 1)
    e2 = res.SimpleEOS(2, 3)

    @testset "Get gradient of an EOS" begin
        @test WaterData.gradientlog(e1, 1., 1.) ≈ [1/1, 1/1]
        @test WaterData.gradientlog(e1, 4., 1.) ≈ [1/4, 1/1]
        @test WaterData.gradientlog(e2, 1., 1.) ≈ [2/1, 3/1]
        @test WaterData.gradientlog(e2, 4., 2.) ≈ [2/4, 3/2]
    end

    @testset "Get directional derivative of an EOS" begin
        @test WaterData.partiallnT(e1, 1., 1.) ≈ 1/1
        @test WaterData.partiallnT(e1, 4., 1.) ≈ 1/1
        @test WaterData.partiallnT(e2, 1., 1.) ≈ 3/1
        @test WaterData.partiallnT(e2, 4., 2.) ≈ 3/2
    end
end


text:
using LightNLP
using LightNLP.NER
using JSON

config = JSON.parsefile(ARGS[1])

if config["training"]
    ner = NER.Model(config)
    #LightNLP.save("ner.jld2", ner)
else
    ner = LightNLP.load("ner.jld2")
    testfile = ".data/eng.testb.BIOES"
    ner(testfile)
end
println("Finish.")


text:
function fit end
function loglikelihood end


text:
import Base.getindex

## Implement some edge buffering/connector

# TODO how to declare/document interface?

struct Latest
    state
end

push!(aggr::Latest, v) = aggr.state[] = v

generate(aggr::Latest) = aggr.state[]

# TODO `generate` could imply `reset` (one interface only)
reset!(aggr::Latest) = nothing

"""
    latest create a connector that always return last seen value.
"""
latest(st) = Latest(Ref(st))

## Implement combination of node values

# TODO can probably be replaced by NamedTuple
"""
    BVal

Value created by a Builder (combining multiple node into Tuple)
"""
struct BVal{T <: Tuple}
    fields
    val::T
end

eltype(::Type{BVal{T}}) where {T} = T

function getindex(bval::BVal, i::String)
    bval.val[bval.fields[i]]
end
function getindex(bval::BVal, i)
    bval.val[i]
end

# NTicker: ref to a node inside the Graph (reducer, triggering, name)
struct NTicker
    name
    node
    buf
    triggering
end

"""
    Builder combine multiple node inputs into one output.

Internal component used by `combine`.
"""
struct Builder
    fields
    bufs
    triggering
end

function Builder(nticks)
    # names must be unique
    @assert length(Set([ntkr.name for ntkr in nticks])) == length(nticks)
    # one edge must be triggering
    @assert any([ntkr.triggering for ntkr in nticks])

    Builder(
        Dict(ntkr.name => i for (i, ntkr) in enumerate(nticks)),
        [ntkr.buf for ntkr in nticks],
        [ntkr.triggering for ntkr in nticks],
    )
end

function push!(b::Builder, i, v)
    push!(b.bufs[i], v)
    return b.triggering[i]
end

function generate(b::Builder)  # TODO Nullable
    # TODO even if triggering, maybe can ask each aggregator if
    #      they are ready (maybe some don't want the universe to be
    #      created if they have no value)
    BVal(
        b.fields,  # dict are mutable, pass by ref
        Tuple(generate.(b.bufs)),
    )
end

reset!(b::Builder) = reset!.(b.bufs)

function feed(b::Builder, i, v)
    trig = push!(b, i, v)
    if trig
        # TODO generate could still fail (e.g. missing required input)
        bv = generate(b)
        reset!(b)
        return Nullable(bv)
    else
        return Nullable()
    end
end

# TODO multiple combine version could handle other type than BVal, e.g. Array
#      or even splatting small number of parent values
function combine!(d::Dag, parents)
    # TODO NTicker is not needed ... but parents could benefit type for doc
    uticks = map(t->NTicker(t...), parents)
    builder = Builder(uticks)

    TX = map(ut->eltype(ut.node), uticks)
    OT = BVal{Tuple{TX...}}

    bparents = [
        (ntkr.node, v -> feed(builder, i, v))
        for (i, ntkr) in enumerate(uticks)
    ]

    return node!(d::Dag, OT, bparents)
end

# combine assuming Latest/init and apply a function to all parts
function apply!(d::Dag, fun, init, nodes...)
    CT = mapreduce(eltype, promote_type, nodes)
    bufs = [latest(init) for n in nodes]
    return node!(d, CT, [
        (n, v -> begin
            push!(bufs[i], v)
            Nullable(fun(map(generate, bufs)))
        end)
        for (i, n) in enumerate(nodes)
    ])
end

"""
    tr!(d::Dag, parent::Node, fun::Function)

Create a new node transforming any parent input through `fn`.
"""
function tr!(d::Dag, parent, fn)
    NT = return_type_fn1(fn, eltype(parent))
    node!(d, NT, [(parent, x->Nullable(fn(x)))])
end


text:
#Parameter file for instiating equation 31 in Kloden & Platen
include("problems.jl")

const avec = [0.8, 0.4, 0.9, 0.4, 0.2, 0.21, 0.6, 0.43, 0.21, -0.1]
const prob = Problem31(avec)

#Must be between -1 and 1!
const x0 = [-0.5, 0.5, 0.3, -0.2, 0.7, 0.85, 0.0, -0.3, -0.34, 0.62]
t = linspace(0., 1., 5001)

max_stars_repo_name: onoderat/SDEs

max_issues_repo_name: onoderat/SDEs

max_forks_repo_name: onoderat/SDEs


text:
mutable struct LinearEliminationMap
    a::Vector{Float64}                        = Float64[]
    b::Vector{Float64}                        = Float64[]
    c::Vector{Float64}                        = Float64[]
    xi::Vector{MOI.VariableIndex}             = MOI.VariableIndex[]
    xj::Vector{MOI.VariableIndex}             = MOI.VariableIndex[]
    eliminated::OrderedDict{CI{SAF,ET}, Bool} = Tuple{SAF,ET}()
end

UnitUpperTriangular

Ax = b -> (UU)x = b*


text:
@inline default_distance(r₁, r₂) = r₁ - r₂

const R_MIN = 2 ^ (1.0 / 6)