text:
module Datasets
import ..Cubes.ESDLZarr: toaxis, axname, AbstractCubeData, ZArrayCube, propfromattr, subsetcube, caxes, concatenateCubes
import Zarr: ZGroup, zopen
import ..Cubes.Axes: axsym, CubeAxis, findAxis, CategoricalAxis
import ..Cubes: AbstractCubeData, Cube
import DataStructures: OrderedDict, counter
struct Dataset
    cubes::OrderedDict{Symbol,AbstractCubeData}
    axes::Dict{Symbol,CubeAxis}
end
function Dataset(;cubesnew...)
    axesall = Set{CubeAxis}()
    foreach(values(cubesnew)) do c
        ax = caxes(c)
        foreach(a->push!(axesall,a),ax)
    end
    axesall = collect(axesall)
    axnameall = [axsym(a) in (:Lon, :Lat, :Time) ? Symbol(lowercase(string(axsym(a)))) : axsym(a) for a in axesall]
    axesnew = Dict{Symbol,CubeAxis}(axnameall[i]=>axesall[i] for i in 1:length(axesall))
    Dataset(OrderedDict(cubesnew), axesnew)
end

function Base.show(io::IO,ds::Dataset)
    println(io,"ESDL Dataset")
    println(io,"Dimensions: ")
    foreach(a->println(io,"   ", a),values(ds.axes))
    print(io,"Variables: ")
    foreach(i->print(io,i," "),keys(ds.cubes))
end
function Base.propertynames(x::Dataset, private=false)
    if private
        Symbol[:cubes; :axes; keys(x.cubes); keys(x.axes)]
    else
        Symbol[collect(keys(x.cubes)); collect(keys(x.axes))]
    end
end
function Base.getproperty(x::Dataset,k::Symbol)
    if k === :cubes
        return getfield(x,:cubes)
    elseif k === :axes
        return getfield(x,:axes)
    else
        x[k]
    end
end
Base.getindex(x::Dataset,i::Symbol) = haskey(x.cubes, i) ? x.cubes[i] : haskey(x.axes,i) ? x.axes[i] : throw(ArgumentError("$i not found in Dataset"))
function Base.getindex(x::Dataset,i::Vector{Symbol})
    cubesnew = Dict{Symbol, AbstractCubeData}(j=>x.cubes[j] for j=i)
    Dataset(;cubesnew...)
end

function fuzzyfind(s::String,comp::Vector{String})
    sl = lowercase(s)
    f = findall(i->startswith(lowercase(i),sl),comp)
    if length(f) != 1
        throw(KeyError("Name $s not found"))
    else
        f[1]
    end
end
function Base.getindex(x::Dataset,i::Vector{String})
    istr = string.(keys(x.cubes))
    ids = map(name->fuzzyfind(name,istr),i)
    syms   = map(j->Symbol(istr[j]),ids)
    cubesnew = Dict{Symbol, AbstractCubeData}(Symbol(i[j])=>x.cubes[syms[j]] for j=1:length(ids))
    Dataset(;cubesnew...)
end
Base.getindex(x::Dataset,i::String)=getproperty(x,Symbol(i))
function subsetcube(x::Dataset; var=nothing, kwargs...)
    if var ===nothing
        cc = x.cubes
        Dataset(;map(ds->ds=>subsetcube(cc[ds];kwargs...),collect(keys(cc)))...)
    elseif isa(var,String) || isa(var, Symbol)
        subsetcube(getproperty(x,Symbol(var));kwargs...)
    else
        cc = x[var].cubes
        Dataset(;map(ds->ds=>subsetcube(cc[ds];kwargs...),collect(keys(cc)))...)
    end
end
function collectdims(g::ZGroup)
  dlist = Set{Tuple{String,Int,Int}}()
  foreach(g.arrays) do ar
    k,v = ar
    for (len,dname) in zip(size(v),reverse(v.attrs["_ARRAY_DIMENSIONS"]))
      if !occursin("bnd",dname) && !occursin("bounds",dname)
        offs = get(g.arrays[dname].attrs,"_ARRAY_OFFSET",0)
        push!(dlist,(dname,offs,len))
      end
    end
  end
  outd = Dict(d[1] => (ax = toaxis(d[1],g,d[2],d[3]), offs = d[2]) for d in dlist)
  length(outd)==length(dlist) || throw(ArgumentError("All Arrays must have the same offset"))
  outd
end

function Dataset(g::ZGroup)

  isempty(g.arrays) && throw(ArgumentError("Zarr Group does not contain datasets."))
  dimlist = collectdims(g)
  dnames  = string.(keys(dimlist))
  varlist = filter(g.arrays) do ar
    upname = uppercase(ar[1])
    !occursin("BNDS",upname) && !occursin("BOUNDS",upname) && !any(i->isequal(upname,uppercase(i)),dnames)
  end
  allcubes = OrderedDict{Symbol,AbstractCubeData}()
  for iv in  varlist
    vname, zarray = iv
    vardims = reverse((zarray.attrs["_ARRAY_DIMENSIONS"]...,))
    iax = [dimlist[vd].ax for vd in vardims]
    offs = [dimlist[vd].offs for vd in vardims]
    subs = if all(iszero,offs)
      nothing
    else
      ntuple(i->(offs[i]+1):(offs[i]+length(iax[i])),length(offs))
    end
    allcubes[Symbol(vname)] = ZArrayCube{eltype(zarray),ndims(zarray),typeof(zarray),typeof(subs)}(zarray,iax,subs,true,propfromattr(zarray.attrs))
  end
  sdimlist = Dict(Symbol(k)=>v.ax for (k,v) in dimlist)
  Dataset(allcubes,sdimlist)
end
Base.getindex(x::Dataset;kwargs...) = subsetcube(x;kwargs...)
Dataset(s::String;kwargs...) = Dataset(zopen(s);kwargs...)
ESDLDataset(;kwargs...) = Dataset(ESDL.ESDLDefaults.cubedir[];kwargs...)

function Cube(ds::Dataset; joinname="Variable")

  dl = collect(keys(ds.axes))
  dls = string.(dl)
  length(ds.cubes)==1 && return first(values(ds.cubes))
  #TODO This is an ugly workaround to merge cubes with different element types,
  # There should bde a more generic solution
  eltypes = map(eltype,values(ds.cubes))
  majtype = findmax(counter(eltypes))[2]
  newkeys = Symbol[]
  for k in keys(ds.cubes)
    c = ds.cubes[k]
    if all(axn->findAxis(axn,c)!==nothing,dls) && eltype(c)==majtype
      push!(newkeys,k)
    end
  end
  if length(newkeys)==1
    return ds.cubes[first(newkeys)]
  else
    varax = CategoricalAxis(joinname, string.(newkeys))
    return concatenateCubes([ds.cubes[k] for k in newkeys], varax)
  end
end

Cube(z::ZGroup;joinname="Variable") = Cube(Dataset(z),joinname=joinname)

end


text:
module RecombinatorKMeans

using Random
using Statistics
using StatsBase
using ExtractMacro

export kmeans, reckmeans, gakmeans, rswapkmeans

mutable struct Configuration
    m::Int
    k::Int
    n::Int
    c::Vector{Int}
    cost::Float64
    costs::Vector{Float64}
    centroids::Matrix{Float64}
    active::BitVector
    nonempty::BitVector
    csizes::Vector{Int}
    function Configuration(m::Int, k::Int, n::Int, c::Vector{Int}, costs::Vector{Float64}, centroids::Matrix{Float64})
        @assert length(c) == n
        @assert length(costs) == n
        @assert size(centroids) == (m, k)
        cost = sum(costs)
        active = trues(k)
        nonempty = trues(k)
        csizes = zeros(Int, k)
        if !all(c .== 0)
            for i = 1:n
                csizes[c[i]] += 1
            end
            nonempty .= csizes .> 0
        end
        return new(m, k, n, c, cost, costs, centroids, active, nonempty, csizes)
    end
    function Base.copy(config::Configuration)
        @extract config : m k n c cost costs centroids active nonempty csizes
        return new(m, k, n, copy(c), cost, copy(costs), copy(centroids), copy(active), copy(nonempty), copy(csizes))
    end
end

function Configuration(data::Matrix{Float64}, centroids::Matrix{Float64})
    m, n = size(data)
    k = size(centroids, 2)
    @assert size(centroids, 1) == m

    c = zeros(Int, n)
    costs = fill(Inf, n)
    config = Configuration(m, k, n, c, costs, centroids)
    partition_from_centroids!(config, data)
    return config
end

function Base.copy!(config::Configuration, other::Configuration)
    @extract config : m k n c costs centroids active nonempty csizes
    @extract other : om=m ok=k on=n oc=c ocosts=costs ocentroids=centroids oactive=active ononempty=nonempty ocsizes=csizes
    @assert m == om && k == ok && n == on
    copy!(c, oc)
    config.cost = other.cost
    copy!(costs, ocosts)
    copy!(centroids, ocentroids)
    copy!(active, oactive)
    copy!(nonempty, ononempty)
    copy!(csizes, ocsizes)
    return config
end

function remove_empty!(config::Configuration)
    @extract config: m k n c costs centroids active nonempty csizes

    k_new = sum(nonempty)
    if k_new == k
        return config
    end
    centroids = centroids[:, nonempty]
    new_inds = cumsum(nonempty)
    for i = 1:n
        @assert nonempty[c[i]]
        c[i] = new_inds[c[i]]
    end
    csizes = csizes[nonempty]
    nonempty = trues(k_new)

    config.k = k_new
    config.centroids = centroids
    config.csizes = csizes
    config.nonempty = nonempty

    return config
end

Base.@propagate_inbounds function _cost(d1, d2)
    v1 = 0.0
    @simd for l = 1:length(d1)
        v1 += (d1[l] - d2[l])^2
    end
    return v1
end

# grouping method from Kaukoranta, Fränti, Nevlainen, "Reduced comparison for the exact GLA"
function partition_from_centroids!(config::Configuration, data::Matrix{Float64})
    @extract config: m k n c costs centroids active nonempty csizes
    @assert size(data) == (m, n)

    active_inds = findall(active)
    all_inds = collect(1:k)

    fill!(nonempty, false)
    fill!(csizes, 0)

    num_fullsearch = 0
    cost = 0.0
    @inbounds for i in 1:n
        ci = c[i]
        if ci > 0 && active[ci]
            old_v′ = costs[i]
            @views new_v′ = _cost(data[:,i], centroids[:,ci])
            fullsearch = (new_v′ > old_v′)
        else
            fullsearch = (ci == 0)
        end
        num_fullsearch += fullsearch

        v, x, inds = fullsearch ? (Inf, 0, all_inds) : (costs[i], ci, active_inds)
        for j in inds
            @views v′ = _cost(data[:,i], centroids[:,j])
            if v′ < v
                v, x = v′, j
            end
        end
        costs[i], c[i] = v, x
        nonempty[x] = true
        csizes[x] += 1
        cost += v
    end

    config.cost = cost
    return config
end

let centroidsdict = Dict{NTuple{2,Int},Matrix{Float64}}()

    global function centroids_from_partition!(config::Configuration, data::Matrix{Float64})
        @extract config: m k n c costs centroids active nonempty csizes
        @assert size(data) == (m, n)

        new_centroids = get!(centroidsdict, (m,k)) do
            zeros(Float64, m, k)
        end
        fill!(new_centroids, 0.0)
        fill!(active, false)
        @inbounds for i = 1:n
            j = c[i]
            for l = 1:m
                new_centroids[l,j] += data[l,i]
            end
        end
        @inbounds for j = 1:k
            z = csizes[j]
            z > 0 || continue
            @assert nonempty[j]
            for l = 1:m
                new_centroids[l,j] /= z
                new_centroids[l,j] ≠ centroids[l,j] && (active[j] = true)
                centroids[l,j] = new_centroids[l,j]
            end
        end
        return config
    end

    global function clear_cache!()
        empty!(centroidsdict)
    end
end

function check_empty!(config::Configuration, data::Matrix{Float64})
    @extract config: m k n c costs centroids active nonempty csizes
    num_nonempty = sum(nonempty)
    num_centroids = min(config.n, config.k)
    gap = num_centroids - num_nonempty
    gap == 0 && return false
    to_fill = findall(.~(nonempty))[1:gap]
    for j in to_fill
        i = rand(1:n)
        centroids[:,j] .= data[:,i]
        active[j] = true
    end
    return true
end

function init_centroid_unif(data::Matrix{Float64}, k::Int; kw...)
    m, n = size(data)
    centroids = zeros(m, k)
    for j = 1:k
        i = rand(1:n)
        centroids[:,j] .= data[:,i]
    end
    return Configuration(data, centroids)
end

function compute_costs_one!(costs::Vector{Float64}, data::AbstractMatrix{<:Float64}, x::AbstractVector{<:Float64})
    m, n = size(data)
    @assert length(costs) == n
    @assert length(x) == m

    @inbounds for i = 1:n
        @views costs[i] = _cost(data[:,i], x)
    end
    return costs
end
compute_costs_one(data::AbstractMatrix{<:Float64}, x::AbstractVector{<:Float64}) = compute_costs_one!(Array{Float64}(undef,size(data,2)), data, x)

# This function was first written from scratch from the k-means++ paper, then improved after reading
# the corresponding scikit-learn's code at https://github.com/scikit-learn/scikit-learn, then heavily modified.
# The scikit-learn's version was first coded by Jan Schlueter as a port of some other code that is now lost.
function init_centroid_pp(pool::Matrix{Float64}, k::Int; ncandidates = nothing, w = nothing, data = nothing)
    m, np = size(pool)
    @assert np ≥ k

    ncandidates::Int = ncandidates ≡ nothing ? floor(Int, 2 + log(k)) : ncandidates
    w = (w ≡ nothing ? ones(np) : w)::Vector{Float64}
    data = (data ≡ nothing ? pool : data)::Matrix{Float64}
    dataispool = data ≡ pool
    n = size(data, 2)
    @assert size(data, 1) == m

    centr = zeros(m, k)
    i = sample(1:np, Weights(w))
    pooli = pool[:,i]
    centr[:,1] .= pooli

    costs = compute_costs_one(data, pooli)
    pcosts = dataispool ? costs : compute_costs_one(pool, pooli)

    curr_cost = sum(costs)
    c = ones(Int, n)

    new_costs, new_c = similar(costs), similar(c)
    new_costs_best, new_c_best = similar(costs), similar(c)
    for j = 2:k
        pw = Weights(pcosts .* w)
        nonz = count(pw .≠ 0)
        candidates = sample(1:np, pw, min(ncandidates,np,nonz), replace = false)
        cost_best = Inf
        y_best = 0
        for y in candidates
            pooly = pool[:,y]
            compute_costs_one!(new_costs, data, pooly)
            cost = 0.0
            @inbounds for i = 1:n
                v = new_costs[i]
                v′ = costs[i]
                if v < v′
                    new_c[i] = j
                    cost += v
                else
                    new_costs[i] = v′
                    new_c[i] = c[i]
                    cost += v′
                end
            end
            if cost < cost_best
                cost_best = cost
                y_best = y
                new_costs_best, new_costs = new_costs, new_costs_best
                new_c_best, new_c = new_c, new_c_best
            end
        end
        @assert y_best ≠ 0 && cost_best < Inf
        pooly = pool[:,y_best]
        centr[:,j] .= pooly
        costs, new_costs_best = new_costs_best, costs
        if dataispool
            pcosts = costs
        elseif j < k
            pcosts .= min.(pcosts, compute_costs_one(pool, pooly))
        end
        c, new_c_best = new_c_best, c
        dataispool && (pc = c)
    end
    return Configuration(m, k, n, c, costs, centr)
end

"""
  kmeans(data::Matrix{Float64}, k::Integer; keywords...)

Runs k-means using Lloyd's algorithm on the given data matrix (if the size is `d`×`n` then `d` is
the dimension and `n` the number of points, i.e. data is organized by column).
It returns: 1) a vector of labels (`n` integers from 1 to `k`); 2) a `d`×`k` matrix of centroids;
3) the final cost; 4) whether it converged or not

The possible keyword arguments are:

* `max_it`: maximum number of iterations (default=1000). Normally the algorithm stops at fixed
  points.
* `seed`: random seed, either an integer or `nothing` (this is the default, it means no seeding
  is performed).
* `init`: how to initialize (default=`"++"`). It can be a string or a `Matrix{Float64}`. If it's
  a matrix, it represents the initial centroids (by column). If it is a string, it can be either
  `"++"` for k-means++ initialization, or `"unif"` for uniform.
* `tol`: a `Float64`, relative tolerance for detecting convergence (default=1e-5).
* `verbose`: a `Bool`; if `true` (the default) it prints information on screen.
* `ncandidates`: if init=="++", set the number of candidates for k-means++ (the default is
  `nothing`, which means that it is set automatically)
"""
function kmeans(
        data::Matrix{Float64}, k::Integer;
        max_it::Integer = 1000,
        seed::Union{Integer,Nothing} = nothing,
        init::Union{String,Matrix{Float64}} = "++",
        tol::Float64 = 1e-5,
        verbose::Bool = true,
        ncandidates::Union{Nothing,Int} = nothing,
    )
    seed ≢ nothing && Random.seed!(seed)
    m, n = size(data)

    if init isa String
        if init == "++"
            config = init_centroid_pp(data, k; ncandidates)
        elseif init == "unif"
            config = init_centroid_unif(data, k)
        else
            throw(ArgumentError("unrecognized init string \"$init\"; should be \"++\" or \"unif\""))
        end
    else
        size(init) == (m, k) || throw(ArgumentError("invalid size of init, expected $((m,k)), found: $(size(init))"))
        centroids = copy(init)
        config = Configuration(data, centroids)
    end

    verbose && println("initial cost = $(config.cost)")
    converged = lloyd!(config, data, max_it, tol, verbose)

    return config.c, config.centroids, config.cost, converged
end

function crossover_GA(configs::Vector{Configuration}, data::Matrix{Float64}, k::Int, a::Int; kw...)
    Jc = length(configs)
    (i1, i2) = [(x,y) for x = 1:Jc for y = (x+1):Jc][a] # TODO: quite ugly...
    return merge_pnn(configs[i1], configs[i2], data, k)
end

function lloyd!(config::Configuration, data::Matrix{Float64}, max_it::Int, tol::Float64, verbose::Bool)
    cost0 = config.cost
    converged = false
    it = 0
    for outer it = 1:max_it
        centroids_from_partition!(config, data)
        old_cost = config.cost
        found_empty = check_empty!(config, data)
        partition_from_centroids!(config, data)
        new_cost = config.cost
        if new_cost ≥ old_cost * (1 - tol) && !found_empty
            verbose && println("converged cost = $new_cost")
            converged = true
            break
        end
        verbose && println("lloyd it = $it cost = $new_cost")
    end
    return converged
end

function combine_configs(config1::Configuration, config2::Configuration)
    @extract config1 : m1=m k1=k n1=n c1=c costs1=costs centroids1=centroids
    @extract config2 : m2=m k2=k n2=n c2=c costs2=costs centroids2=centroids
    @assert (m1, n1) == (m2, n2)
    m, n = m1, n1
    k = k1 + k2
    centroids_new = hcat(centroids1, centroids2)
    @assert size(centroids_new) == (m, k)
    c_new = zeros(Int, n)
    costs_new = zeros(n)
    for i = 1:n
        j1, j2 = c1[i], c2[i]
        cost1, cost2 = costs1[i], costs2[i]
        if cost1 ≤ cost2
            j = j1
            cost = cost1
        else
            j = k1 + j2
            cost = cost2
        end
        c_new[i], costs_new[i] = j, cost
    end
    return Configuration(m, k, n, c_new, costs_new, centroids_new)
end

function pairwise_nn!(config::Configuration, tgt_k::Int)
    @extract config : m k n c costs centroids csizes
    if k < tgt_k
        @assert false # TODO: inflate the config? this shouldn't normally happen
    end
    if k == tgt_k
        return config
    end

    nns = zeros(Int, k)
    nns_costs = fill(Inf, k)
    @inbounds for j = 1:k
        z = csizes[j]
        v, x = Inf, 0
        for j′ = 1:k
            j′ == j && continue
            z′ = csizes[j′]
            @views v1 = (z * z′) / (z + z′) * _cost(centroids[:,j], centroids[:,j′])
            if v1 < v
                v = v1
                x = j′
            end
        end
        nns_costs[j], nns[j] = v, x
    end

    @inbounds while k > tgt_k
        jm = findmin(@view(nns_costs[1:k]))[2]
        js = nns[jm]
        @assert nns_costs[js] == nns_costs[jm]
        @assert jm < js

        # update centroid
        zm, zs = csizes[jm], csizes[js]
        for l = 1:m
            centroids[l,jm] = (zm * centroids[l,jm] + zs * centroids[l,js]) / (zm + zs)
            centroids[l,js] = centroids[l,k]
        end

        # update csizes
        csizes[jm] += zs
        csizes[js] = csizes[k]
        csizes[k] = 0

        # update partition
        # not needed since we don't use c in this computation and
        # we call partition_from_centroids! right after this function
        # for i = 1:n
        #     ci = c[i]
        #     if ci == js
        #         c[i] = jm
        #     elseif ci == k
        #         c[i] = js
        #     end
        # end

        # update nns
        nns[js] = nns[k]
        nns[k] = 0
        nns_costs[js] = nns_costs[k]
        nns_costs[k] = Inf

        num_fullupdates = 0
        for j = 1:(k-1)
            # 1) merged cluster jm, or clusters which used to point to either jm or js
            #    perform a full update
            if j == jm || nns[j] == jm || nns[j] == js
                num_fullupdates += 1
                z = csizes[j]
                v, x = Inf, 0
                for j′ = 1:(k-1)
                    j′ == j && continue
                    z′ = csizes[j′]
                    @views v1 = (z * z′) / (z + z′) * _cost(centroids[:,j], centroids[:,j′])
                    if v1 < v
                        v, x = v1, j′
                    end
                end
                nns_costs[j], nns[j] = v, x
            # 2) clusters that did not point to jm or js
            #    only compare the old cost with the cost for the updated cluster
            else
                z = csizes[j]
                # note: what used to point at k now must point at js
                v, x = nns_costs[j], (nns[j] ≠ k ? nns[j] : js)
                j′ = jm
                z′ = csizes[j′]
                @views v′ = (z * z′) / (z + z′) * _cost(centroids[:,j], centroids[:,j′])
                if v′ < v
                    v, x = v′, j′
                end
                nns_costs[j], nns[j] = v, x
                @assert nns[j] ≠ j
            end
        end

        k -= 1
    end

    config.k = k
    config.centroids = centroids[:,1:k]
    config.csizes = csizes[1:k]
    @assert all(config.csizes .> 0)
    config.active = trues(k)
    config.nonempty = trues(k)
    fill!(config.c, 0) # reset in order for partition_from_centroids! to work

end

function merge_pnn(config1::Configuration, config2::Configuration, data::Matrix{Float64}, tgt_k::Int)
    config = combine_configs(config1, config2)
    centroids_from_partition!(config, data)
    remove_empty!(config)
    pairwise_nn!(config, tgt_k)
    partition_from_centroids!(config, data)
    return config
end

struct Results
    exit_status::Symbol
    labels::Vector{Int}
    centroids::Matrix{Float64}
    cost::Float64
end

"""
  reckmeans(data::Matrix{Float64}, k::Integer, J::Integer; keywords...)

Runs recombinator-k-means on the given data matrix (if the size is `d`×`n` then `d` is
the dimension and `n` the number of points, i.e. data is organized by column). `k` is
the number of clusters, `J` is the population size

It returns an object of type `Results`, which contains the following fields:
* exit_status: a symbol that indicates the reason why the algorithm stopped. It can take three
  values, `:collapsed`, `:maxgenerations` or `:noimprovement`.
* labels: a vector of labels (`n` integers from 1 to `k`)
* centroids: a `d`×`k` matrix of centroids
* cost: the final cost

The possible keyword arguments are:

* `Δβ`: a `Float64` (default=0.1), the reweigting parameter increment (only makes sense if
  non-negative).
* `seed`: random seed, either an integer or `nothing` (this is the default, it means no seeding
  is performed).
* `tol`: a `Float64` (default=1e-4), the relative tolerance for determining whether the solutions
  have collapsed.
* `lltol`: a `Float64` (default=1e-5), relative tolerance for Lloyd (kmeans) convergence.
* `max_it`: an `Int` (default=10), maximum number of Lloyd (kmeans) iterations.
* `max_gen`: an `Int` (default=1_000), maximum number of generations.
* `stop_if_noimprovement`: a `Bool` (default=false), stop if the new generation did not improve on
  the old one
* `verbose`: a `Bool`; if `true` (the default) it prints information on screen.

If there are multiple threads available this function will parallelize each batch.
"""
function reckmeans(
        data::Matrix{Float64}, k::Integer, J::Int;
        Δβ::Float64 = 0.1,
        seed::Union{Integer,Nothing} = nothing,
        tol::Float64 = 1e-4,
        lltol::Float64 = 1e-5,
        verbose::Bool = true,
        max_it::Int = 10,
        stop_if_noimprovement::Bool = false,
        max_gen::Int = 1_000
    )
    J ≥ 2 || throw(ArgumentError("J must be at least 2"))

    seed ≢ nothing && Random.seed!(seed)
    m, n = size(data)

    β = 0.0
    best_cost = Inf
    local best_config::Configuration

    pool = data
    w = ones(size(pool, 2))

    configs = Configuration[]

    exit_status = :running
    for generation = 0:max_gen
        verbose && @info "gen = $generation"
        @assert length(w) == size(pool, 2)
        h0 = hash(pool)
        new_configs = Vector{Configuration}(undef, J)
        Threads.@threads for a = 1:J
            h = hash((seed, a), h0)
            Random.seed!(h)  # horrible hack to ensure determinism (not really required, only useful for experiments)
            config = init_centroid_pp(pool, k; w, data)

            converged = lloyd!(config, data, max_it, lltol, false)
            verbose && println("  a = $a cost = $(config.cost)")
            new_configs[a] = config
        end
        append!(configs, new_configs)
        sort!(configs, by=config->config.cost)
        configs = configs[1:J]

        if configs[1].cost < best_cost
            best_cost = configs[1].cost
            best_config = copy(configs[1])
            improved = true
        else
            improved = false
        end
        mean_cost = mean(config.cost for config in configs)
        stddev_cost = std((config.cost for config in configs), mean = mean_cost)
        verbose && println("  best_cost = $best_cost mean cost = $mean_cost ± $stddev_cost")
        if stop_if_noimprovement && !improved
            verbose && @info "no improvement"
            exit_status = :noimprovement
            break
        end
        if mean_cost ≤ best_cost * (1 + tol)
            verbose && @info "collapsed"
            exit_status = :collapsed
            break
        end

        if generation < max_gen
            pool = hcat((config.centroids for config in configs)...)
            resize!(w, size(pool,2))
            β += Δβ
            for a = 1:(length(w) ÷ k)
                w[(1:k) .+ (a-1)*k] .= exp(-β * ((configs[a].cost - best_cost) / (mean_cost - best_cost)))
            end
            # w ./= sum(w)
        end
    end
    exit_status == :running && (exit_status = :maxgenerations)

    verbose && @info "final cost = $best_cost"

    clear_cache!()

    return Results(
        exit_status,
        best_config.c,
        best_config.centroids,
        best_config.cost,
        )
end

# find c such that: (c * (c-1)) ÷ 2 ≥ J
elitist_crossset_size(J) = ceil(Int, (1 + √(1 + 8J)) / 2)

"""
  gakmeans(data::Matrix{Float64}, k::Integer, J::Integer; keywords...)

Runs ga-k-means on the given data matrix (if the size is `d`×`n` then `d` is
the dimension and `n` the number of points, i.e. data is organized by column). `k` is
the number of clusters, `J` is the population size

It returns an object of type `Results`, which contains the following fields:
* exit_status: a symbol that indicates the reason why the algorithm stopped. It can take three
  values, `:collapsed`, `:maxgenerations` or `:noimprovement`.
* labels: a vector of labels (`n` integers from 1 to `k`)
* centroids: a `d`×`k` matrix of centroids
* cost: the final cost

The possible keyword arguments are:

* `seed`: random seed, either an integer or `nothing` (this is the default, it means no seeding
  is performed).
* `init`: a `String` (default=`"++"`). It can be either `"++"` for k-means++ initialization,
  or `"unif"` for uniform+Lloyd, `"raw"` for uniform (non-optimized).
* `tol`: a `Float64` (default=1e-4), the relative tolerance for determining whether the solutions
  have collapsed.
* `lltol`: a `Float64` (default=1e-5), relative tolerance for Lloyd (kmeans) convergence.
* `max_it`: an `Int` (default=10), maximum number of Lloyd (kmeans) iterations.
* `max_gen`: an `Int` (default=1_000), maximum number of generations.
* `stop_if_noimprovement`: a `Bool` (default=true), stop if the new generation did not improve on
  the old one
* `verbose`: a `Bool`; if `true` (the default) it prints information on screen.

If there are multiple threads available this function will parallelize each batch.
"""
function gakmeans(
        data::Matrix{Float64}, k::Integer, J::Int;
        seed::Union{Integer,Nothing} = nothing,
        init::String = "++",
        tol::Float64 = 1e-4,
        lltol::Float64 = 1e-5,
        verbose::Bool = true,
        max_it::Integer = 10,
        max_gen::Integer = 1_000,
        stop_if_noimprovement::Bool = true,
    )
    init ∈ ["++", "unif", "raw"] || throw(ArgumentError("init must be one of \"++\", \"unif\" or \"raw\""))
    J ≥ 2 || throw(ArgumentError("J must be at least 2"))

    seed ≢ nothing && Random.seed!(seed)
    m, n = size(data)

    Jc = elitist_crossset_size(J)

    init_func = init == "++" ? init_centroid_pp :
                init_centroid_unif

    configs = Vector{Configuration}(undef, J)

    h0 = hash(data)

    exit_status = :running
    Threads.@threads for a = 1:J
        h = hash((seed, a), h0)
        Random.seed!(h)  # horrible hack to ensure determinism (not really required, only useful for experiments)
        config = init_func(data, k)
        init ≠ "raw" && (converged = lloyd!(config, data, max_it, lltol, false))
        configs[a] = config
    end
    sort!(configs, by=config->config.cost)
    costs = [config.cost for config in configs]
    verbose && @info "costs after init: $(costs)"

    best_cost = configs[1].cost
    best_config = copy(configs[1])

    mean_cost = mean(config.cost for config in configs)
    stddev_cost = std((config.cost for config in configs), mean = mean_cost)

    for generation = 1:max_gen
        old_costs = costs
        elite_configs = configs[1:Jc]
        new_configs = Vector{Configuration}(undef, J)
        h0 = hash(configs)
        Threads.@threads for a = 1:J
            h = hash((seed, a), h0)
            Random.seed!(h)  # horrible hack to ensure determinism (not really required, only useful for experiments)
            config::Configuration = crossover_GA(elite_configs, data, k, a)
            converged = lloyd!(config, data, max_it, lltol, false)

            verbose && println("  a = $a cost = $(config.cost)")
            new_configs[a] = config
        end
        sort!(new_configs, by=config->config.cost) # TODO ? the scaling is bad, but in practice this is negligible
        configs = new_configs
        costs = [config.cost for config in configs]

        if configs[1].cost < best_cost
            best_cost = configs[1].cost
            best_config = copy(configs[1])
            improved = true
        else
            improved = false
        end

        verbose && @info "costs after gen $generation: $costs"
        mean_cost = mean(costs)
        stddev_cost = std(costs, mean = mean_cost)
        verbose && println("  best_cost = $best_cost mean cost = $mean_cost ± $stddev_cost")

        ## same condition as in the original implementation
        if stop_if_noimprovement && !improved && generation > 2
            verbose && @info "no improvement"
            exit_status = :noimprovement
            break
        end
        if mean_cost ≤ best_cost * (1 + tol)
            verbose && @info "collapsed"
            exit_status = :collapsed
            break
        end
    end
    if exit_status == :running
        exit_status = :maxgenerations
    end

    verbose && @info "final cost = $best_cost"

    clear_cache!()

    return Results(
        exit_status,
        best_config.c,
        best_config.centroids,
        best_config.cost,
        )
end

"""
  rswapkmeans(data::Matrix{Float64}, k::Integer; keywords...)

Runs random-swap-k-means on the given data matrix (if the size is `d`×`n` then `d` is
the dimension and `n` the number of points, i.e. data is organized by column). `k` is
the number of clusters, `J` is the population size

It returns an object of type `Results`, which contains the following fields:
* exit_status: a symbol that indicates the reason why the algorithm stopped. It can take three
  values, `:solved`, `:maxswaps` or `:outoftime`.
* labels: a vector of labels (`n` integers from 1 to `k`)
* centroids: a `d`×`k` matrix of centroids
* cost: the final cost

The possible keyword arguments are:

* `seed`: random seed, either an integer or `nothing` (this is the default, it means no seeding
  is performed).
* `init`: a `String` (default=`"++"`). It can be either `"++"` for k-means++ initialization,
  or `"unif"` for uniform+Lloyd.
* `lltol`: a `Float64` (default=1e-5), relative tolerance for Lloyd (kmeans) convergence.
* `max_it`: an `Int` (default=2), maximum number of Lloyd (kmeans) iterations.
* `max_swaps`: an `Int` (default=1_000), maximum number of swaps to attempt.
* `max_time`: a `Float64` (default=Inf), maximum allowed (wall-clock) time.
* `target_cost`: a `Float64` (default=0.0), stop if the cost achieves this value (the problem is
  considered solved)
* `final_convergence`: a `Bool`; if `true` (the default) a final run of Lloyd (kmeans) with
  unbounded iterations is performed until convergence, in case it hadn't already been achieved
* `verbose`: a `Bool`; if `true` (the default) it prints information on screen.
"""
function rswapkmeans(
        data::Matrix{Float64}, k::Integer;
        seed::Union{Integer,Nothing} = nothing,
        init::String = "++",
        lltol::Float64 = 1e-5,
        max_swaps::Integer = 1_000,
        max_time::Float64 = Inf,
        target_cost::Float64 = 0.0,
        max_it::Integer = 2,
        final_converge::Bool = true,
        verbose::Bool = true,
    )
    init ∈ ["++", "unif"] || throw(ArgumentError("init must be one of \"++\" or \"unif\""))

    seed ≢ nothing && Random.seed!(seed)
    m, n = size(data)

    t0 = time()

    init_func = init == "++" ? init_centroid_pp :
                init_centroid_unif

    config = init_func(data, k)
    converged = lloyd!(config, data, max_swaps, lltol, false)

    verbose && @info "cost after init: $(config.cost)"

    exit_status = :running
    best_config = copy(config)
    best_converged = converged
    it = 0
    for outer it = 1:max_swaps
        copy!(config, best_config)

        j = rand(1:k)
        i = rand(1:n)
        config.centroids[:,j] = @view data[:,i]
        config.active[j] = true
        partition_from_centroids!(config, data)

        converged = lloyd!(config, data, max_it, lltol, false)

        if config.cost < best_config.cost
            best_config, config = config, best_config
            best_converged = converged
            t = time() - t0
            verbose && println("  it = $it cost = $(best_config.cost) time = $t")
        end
        if best_config.cost ≤ target_cost
            exit_status = :solved
            verbose && @info "target_cost achieved"
            break
        end
        if time() - t0 ≥ max_time
            exit_status = :outoftime
            verbose && @info "max_time reached"
            break
        end
    end
    if exit_status == :running
        exit_status = :maxswaps
    end
    if final_converge && !best_converged
        best_converged = lloyd!(best_config, data, typemax(Int), lltol, false)
    end
    t = time() - t0
    verbose && @info "final cost = $(best_config.cost) time = $t"

    clear_cache!()

    return Results(
        exit_status,
        best_config.c,
        best_config.centroids,
        best_config.cost
      )
end

# Centroid Index
# P. Fränti, M. Rezaei and Q. Zhao, Centroid index: cluster level similarity measure, Pattern Recognition, 2014
function CI(true_centroids::Matrix{Float64}, centroids::Matrix{Float64})
    m, tk = size(true_centroids)
    @assert size(centroids, 1) == m
    k = size(centroids, 2)

    matched = falses(tk)
    @inbounds for j = 1:k
        v = Inf
        p = 0
        for tj = 1:tk
            @views v1 = _cost(true_centroids[:,tj], centroids[:,j])
            if v1 < v
                v = v1
                p = tj
            end
        end
        matched[p] = true
    end

    return tk - count(matched)
end

# CI_sym(centroids1::Matrix{Float64}, centroids2::Matrix{Float64}) =
#     max(CI(centroids1, centroids2), CI(centroids2, centroids1))

function CI_sym(centroids1::Matrix{Float64}, centroids2::Matrix{Float64})
    m, k1 = size(centroids1)
    @assert size(centroids2, 1) == m
    k2 = size(centroids2, 2)

    a12 = zeros(Int, k1)
    a21 = zeros(Int, k2)
    v12 = fill(Inf, k1)
    v21 = fill(Inf, k2)
    @inbounds for j1 = 1:k1, j2 = 1:k2
        @views v = _cost(centroids1[:,j1], centroids2[:,j2])
        if v < v12[j1]
            v12[j1] = v
            a12[j1] = j2
        end
        if v < v21[j2]
            v21[j2] = v
            a21[j2] = j1
        end
    end

    return max(k1 - length(BitSet(a12)), k2 - length(BitSet(a21)))
end

# Variation of Information
# M. Meilă, Comparing clusterings—an information based distance, Journal of multivariate analysis, 2007

xlogx(x) = ifelse(iszero(x), zero(x), x * log(x))
entropy(p) = -sum(xlogx, p)

function VI(c1::Vector{Int}, c2::Vector{Int})
    n = length(c1)
    length(c2) == n || throw(ArgumentError("partitions bust have the same length, given: $(n) and $(length(c2))"))
    a, k1 = extrema(c1)
    a ≥ 1 || throw(ArgumentError("partitions elements must be ≥ 1, found $(a)"))
    a, k2 = extrema(c2)
    a ≥ 1 || throw(ArgumentError("partitions elements must be ≥ 1, found $(a)"))
    o = zeros(k1, k2)
    o1 = zeros(k1)
    o2 = zeros(k2)
    for i = 1:n
        j1, j2 = c1[i], c2[i]
        o[j1, j2] += 1
        o1[j1] += 1
        o2[j2] += 1
    end
    o ./= n
    o1 ./= n
    o2 ./= n
    vi = 2entropy(o) - entropy(o1) - entropy(o2)
    @assert vi > -1e-12
    return max(0.0, vi)
end

end # module

ext: jl


text:
@testset "get burgers data" begin
    xs, ys = get_burgers_data(n=1000)

    @test size(xs) == (2, 1024, 1000)
    @test size(ys) == (1024, 1000)
end

@testset "unit gaussian normalizer" begin
    dims = (3, 3, 5, 6)
    𝐱 = rand(Float32, dims)

    n = UnitGaussianNormalizer(𝐱)

    @test size(n.mean) == size(n.std)
    @test size(encode(n, 𝐱)) == dims
    @test size(decode(n, encode(n, 𝐱))) == dims
end

@testset "get darcy flow data" begin
    xs, ys, _, _ = get_darcy_flow_data()

    @test size(xs) == (1, 85, 85, 1024)
    @test size(ys) == (1, 85, 85, 1024)
end


text:
using Weave, Compat
using Base.Test

function mmtest(source, resfile, doctype)
  weave("documents/$source", out_path = "documents/multimedia/$resfile", doctype=doctype, plotlib=nothing)
  result = @compat readstring(open("documents/multimedia/$resfile"))
  ref = @compat readstring(open("documents/multimedia/$resfile.ref"))
  @test result == ref
end

mmtest("rich_output.jmd", "rich_output.md", "pandoc")
mmtest("rich_output.jmd", "rich_output.tex", "tex")


text:
# ---
# title: 659. Split Array into Consecutive Subsequences
# id: problem659
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Heap, Greedy
# link: <https://leetcode.com/problems/split-array-into-consecutive-subsequences/description/>
# hidden: true
# ---
#
# Given an array `nums` sorted in ascending order, return `true` if and only if
# you can split it into 1 or more subsequences such that each subsequence
# consists of consecutive integers and has length at least 3.
#
#
#
# **Example 1:**
#
#
#
#     Input: [1,2,3,3,4,5]
#     Output: True
#     Explanation:
#     You can split them into two consecutive subsequences :
#     1, 2, 3
#     3, 4, 5
#
#
# **Example 2:**
#
#
#
#     Input: [1,2,3,3,4,4,5,5]
#     Output: True
#     Explanation:
#     You can split them into two consecutive subsequences :
#     1, 2, 3, 4, 5
#     3, 4, 5
#
#
# **Example 3:**
#
#
#
#     Input: [1,2,3,4,4,5]
#     Output: False
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 10000`
#
#
#
#
## @lc code=start
using LeetCode

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


text:
using Test, CatBoost, DataFrames, PyCall
using Aqua

EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")

@testset "`to_pandas` and `pandas_to_df`" begin
    df = DataFrame(; floats=0.5:0.5:3.0, ints=1:6)
    pd = CatBoost.to_pandas(df)
    @test pd isa PyObject
    df2 = pandas_to_df(pd)
    @test df2 == df
end

@testset "Examples" begin
    for ex in readdir(EXAMPLES_DIR)
        @testset "$ex" begin
            # Just check the examples run, for now.
            include(joinpath(EXAMPLES_DIR, ex))
        end
    end
end

Aqua.test_all(CatBoost; ambiguities=false)


text:
function foo()
    lines = readlines(open("data/12.txt", "r"))
    #= lines = readlines(open("h12", "r")) =#
    idx = findfirst(isequal(':'), lines[1]) + 2
    x = lines[1][idx:end]
    offset = 100
    pots = zeros(Int, length(x) + 2 * offset)
    for (i, c) in enumerate(x)
        if c == '#'
            pots[i + offset] = 1
        end
    end
    maps = Dict{Array{Int, 1}, Int}()
    for line in lines[3:end]
        println(line)
        z = zeros(5)
        label = 0
        for (i, c) in enumerate(line[1:5])
            if c == '#'
                z[i] = 1
            end
        end
        if line[end] == '#'
            label = 1
        end
        maps[z] = label
    end

    for q in 1:20
        new_pots = copy(pots)
        for i in 3:length(pots)-2
            key = pots[i-2:i+2]
            if key in keys(maps)
                new_pots[i] = maps[key]
            else
                new_pots[i] = 0
            end
        end
        pots = new_pots
    end

    summed = 0
    for (i, v) in enumerate(pots)
        if v == 1
            summed += i - offset - 1
        end
    end
    return summed
end

println(foo())


text:
niters() = niters(30)
function niters(niters)
    return (τ,n) -> begin
        return n >= niters
    end
end

tolerance_reached() = at_tolerance(1e-3)
function tolerance_reached(tol)
    return (τ,n) -> begin
        return τ <= tol
    end
end


text:
__precompile__(false)
module LCIO_Cxx
using Cxx
using StaticArrays
using Libdl

import Base: getindex, iterate, length, convert, +

const lcio_path="/home/jstrube/.julia/dev/LCIO/deps/usr/"
addHeaderDir(lcio_path*"include", kind=C_System)
Libdl.dlopen(lcio_path*"lib/liblcio.so", Libdl.RTLD_GLOBAL)
cxx"""
#include "lcio.h"
#include "IO/LCReader.h"
#include "IOIMPL/LCFactory.h"
#include "UTIL/CellIDDecoder.h"
#include "EVENT/LCEvent.h"
#include "EVENT/LCCollection.h"
#include "EVENT/MCParticle.h"
#include "EVENT/SimCalorimeterHit.h"
#include "EVENT/TrackerHit.h"
#include "EVENT/SimTrackerHit.h"
#include "EVENT/Track.h"
#include "EVENT/TrackerRawData.h"
#include "EVENT/ReconstructedParticle.h"
#include "EVENT/Vertex.h"
#include "EVENT/LCRelation.h"
#include "EVENT/LCGenericObject.h"
#include "EVENT/Cluster.h"
#include "EVENT/CalorimeterHit.h"
#include <vector>
#include <iostream>
#include <string>
"""

function __init__()
end

const ClusterVec = cxxt"const std::vector<EVENT::Cluster*>*"
const CalorimeterHitVec = cxxt"const std::vector<EVENT::CalorimeterHit*>*"
const TrackVec = cxxt"const std::vector<EVENT::Track*>*"
const StringVec = cxxt"const std::vector<std::string>*"
const MCParticleVec = cxxt"const std::vector<EVENT::MCParticle*>*"

# iteration over std vectors
const StdVecs = Union{ClusterVec, CalorimeterHitVec, TrackVec, StringVec, MCParticleVec}

# uses Julia counting, 1..n
iterate(it::StdVecs, i=1) = i <= length(it) ? (it[i], i+1) : nothing
length(it::StdVecs) = icxx"$(it)->size();"
# 'at' uses C counting, 0..n-1
getindex(it::StdVecs, i) = icxx"$(it)->at($(i)-1);"
eltype(::Type{ClusterVec}) = cxxt"EVENT::Cluster*"
eltype(::Type{CalorimeterHitVec}) = cxxt"EVENT::CalorimeterHit*"
eltype(::Type{TrackVec}) = cxxt"EVENT::Track*"
eltype(::Type{StringVec}) = cxxt"std::string"
eltype(::Type{MCParticleVec}) = cxxt"EVENT::MCParticle*"

const LCReader = cxxt"IO::LCReader*"
function iterate(it::LCReader)
    event = icxx"$(it)->readNextEvent();"
    event == C_NULL && return nothing
    return (event, nothing)
end
iterate(it::LCReader, state) = iterate(it)
length(it::LCReader) = icxx"$(it)->getNumberOfEvents();"
eltype(::Type{LCReader}) = LCEvent

function open(f::Function, fn::AbstractString)
    reader = icxx"IOIMPL::LCFactory::getInstance()->createLCReader();"
    reader == C_NULL && return nothing
    try
        icxx"""$(reader)->open($(fn));"""
        f(reader)
    finally
        icxx"""$(reader)->close();"""
        icxx"""delete $(reader);"""
    end
end

# The iterators use the events as state
# the lcio reader returns NULL at the end of the file
# Because the LCIO getNext function could re-use the memory for the next event,
# the next function should not hold the current and next event at the same time.
# Returning current and reading nextEvent as the new state causes memory corruption
const LCEvent = cxxt"EVENT::LCEvent*"

function iterate(it::LCEvent)
    event = icxx"$(reader)->readNextEvent();"
    event != C_NULL ? (event, nothing) : nothing
end
iterate(it::LCEvent, state) = iterate(it)
length(it::LCEvent) = icxx"$(reader)->getNumberOfEvents();"

#type LCStdHepRdr
#    r::_LCStdHepRdrCpp
#    e::LCEventImpl
#end
#LCStdHepRdr(filename) = LCStdHepRdr(_LCStdHepRdrCpp(filename), LCEventImpl())

# the LCStdHepRdr implementation in C++ is not consistent with the LCIO reader
# this is me trying to make things a bit better
#start(it::LCStdHepRdr) = length(it.r)
#next(it::LCStdHepRdr, state) = readNextEvent(it.r), state-1
#done(it::LCStdHepRdr, state) = state < 1
#length(it::LCStdHepRdr) = getNumberOfEvents(it.r)

#function openStdhep(f::Function, fn::AbstractString)
#    reader = LCStdHepRdr(string(fn))
#    try
#       f(reader)
#    end
#end

# stdhep files only contain one specific type of collection: MCParticle
#function readNextEvent(r::LCStdHepRdr)
#    updateNextEvent(r.r, r.e)
#    getCollection(r.e, "MCParticle")
#end
# open file with reader, returns iterator

# We would like to have a typed collection, but what getCollection returns is unfortunately untyped
# The type is established by reading its name from the collection and mapping it in the LCIOTypemap
struct LCCollection{T}
        coll::cxxt"EVENT::LCCollection*"
end

const CalorimeterHit = cxxt"EVENT::CalorimeterHit*"
const Cluster = cxxt"EVENT::Cluster*"
const LCGenericObject = cxxt"EVENT::LCGenericObject*"
const LCRelation = cxxt"EVENT::LCRelation*"
const MCParticle = cxxt"EVENT::MCParticle*"
const ReconstructedParticle = cxxt"EVENT::ReconstructedParticle*"
const SimCalorimeterHit = cxxt"EVENT::SimCalorimeterHit*"
const SimTrackerHit = cxxt"EVENT::SimTrackerHit*"
const Track = cxxt"EVENT::Track*"
const TrackerHit = cxxt"EVENT::TrackerHit*"
const TrackerRawData = cxxt"EVENT::TrackerRawData*"
const Vertex = cxxt"EVENT::Vertex*"

# map from names stored in collection to actual types
const LCIOTypemap = Dict(
    "CalorimeterHit" => CalorimeterHit,
    "Cluster" => Cluster,
    "LCGenericObject" => LCGenericObject,
    "LCRelation" => LCRelation,
    "MCParticle" => MCParticle,
    "ReconstructedParticle" => ReconstructedParticle,
    "SimCalorimeterHit" => SimCalorimeterHit,
    "SimTrackerHit" => SimTrackerHit,
    "Track" => Track,
    "TrackerHit" => TrackerHit,
    "TrackerRawData" => TrackerRawData,
    "Vertex" => Vertex,
)

const CellIDDecoder{T} = cxxt"UTIL::CellIDDecoder<$(T)>"
CellIDDecoder(t::LCCollection{T}) where {T} = icxx"UTIL::CellIDDecoder<$(T)>($(t.coll));"

# decode(iddecoder) = icxx"$(iddecoder)->decode();"

iterate(it::LCCollection{T}) where {T} = length(it) > 0 ? (it[1], 2) : nothing
iterate(it::LCCollection{T}, i) where {T} = i <= length(it) ? (it[i], i+1) : nothing
length(it::LCCollection) = icxx"$(it.coll)->getNumberOfElements();"
getindex(it::LCCollection{T}, i) where {T} = icxx"static_cast<$(T)>($(it.coll)->getElementAt($(i)-1));"

function getCollection(event, collectionName)
    collection = icxx"""$(event)->getCollection($(pointer(collectionName)));"""
    collectionType = icxx"$(collection)->getTypeName();"
    return LCCollection{LCIOTypemap[String(collectionType)]}(collection)
end

getTypeName(collection::LCCollection) = String( icxx"$(collection.coll)->getTypeName();" )
getCollectionNames(event) = icxx"$(event)->getCollectionNames();"

getEnergy(particle) = icxx"$(particle)->getEnergy();"
function getMomentum(particle)
        p3 = icxx"$(particle)->getMomentum();"
        SVector{3, Float64}(unsafe_load(p3, 1), unsafe_load(p3, 2), unsafe_load(p3, 3))
end
function getPosition(hit)
        pos = icxx"$(hit)->getPosition();"
        SVector{3, Float64}(unsafe_load(pos, 1), unsafe_load(pos, 2), unsafe_load(pos, 3))
end

function getP4(particle)
        p3 = icxx"$(particle)->getMomentum();"
        e = icxx"$(particle)->getEnergy();"
        SVector{4, Float64}(unsafe_load(p3, 1), unsafe_load(p3, 2), unsafe_load(p3, 3), e)
end

getType(particle) = icxx"$(particle)->getType();"

getDetectorName(event) = String(icxx"$(event)->getDetectorName();")

include("MCParticle.jl")
include("CaloHit.jl")
export CalHit, getP4, getPosition, CellIDDecoder,
    getEventNumber, getRunNumber, getDetectorName, getCollection, getCollectionNames, # LCEvent
    getType, getEnergy, getMomentum, getPDG, getParents, getTypeName,
    decode
#    getTypeName, # LCCollection
#    getEnergy, getParents, getDaughters, getPDG, getGeneratorStatus, getSimulatorStatus, isCreatedInSimulation, isBackScatter, vertexIsNotEndpointOfParent, #isDecayedInCalorimeter, hasLeftDetector, isStopped, isOverlay, getVertex, getTime, getEndpoint, getMomentum, getMomentumAtEndpoint, getMass, getCharge, # MCParticle
#    getCalorimeterHits, # Cluster
#    getClusters, getType, isCompound, getMass, getCharge, getReferencePoint, getParticleIDs, getParticleIDUsed, getGoodnessOfPID, getParticles, getClusters, getTracks, #getStartVertex, getEndVertex # ReconstructedParticle

#const MCPARTICLE = "MCParticle"
#const WRITE_NEW = 0
#const WRITE_APPEND = 1

# iteration over std vectors
#const StdVecs = Union{ClusterVec, CalorimeterHitVec, TrackVec, StringVec, MCParticleVec}

# uses Julia counting, 1..n
#start(it::StdVecs) = convert(UInt64, 1)
#next(it::StdVecs, i) = (it[i], i+1)
#done(it::StdVecs, i) = i > length(it)
#length(it::StdVecs) = size(it)
# 'at' uses C counting, 0..n-1
#getindex(it::StdVecs, i) = at(it, i-1)

#start(it::LCReader) = getNumberOfEvents(it)
#next(it::LCReader, state) = readNextEvent(it), state-1
#done(it::LCReader, state) = state < 1
#length(it::LCReader) = getNumberOfEvents(it)

# This version of the iteration runs length() multiple times during the iteration
# if this becomes a speed problem, the length could be memoized, or iteration order could be inverted
#start(it::TypedCollection) = convert(UInt64, 1)
#done(it::TypedCollection, i) = i > length(it)
#next{T}(it::TypedCollection{T}, i) = it[i], i+1
#length(it::TypedCollection) = getNumberOfElements(it)
# getindex uses Julia counting, getElementAt uses C counting
#getindex(it::TypedCollection, i) = getElementAt(it, convert(UInt64, i-1))

#getTypeName{T}(coll::TypedCollection{T}) = "$T"

# to get the typed collection, one needs to read the typename
# then we can return the right type from the LCIOTypemap
#function getCollection(event, collectionName)
#       collection = getEventCollection(event, collectionName)
#       collectionType = getTypeName(collection)
#       return TypedCollection{LCIOTypemap[collectionType]}(collection)
#end

#function getPosition(hit)
#    p3 = Array{Float64,1}(3)
#    valid = getPosition3(hit, p3)
#    return p3
#end

#function getMomentum(particle)
#    p3 = Array{Float64,1}(3)
#    valid = getMomentum3(particle, p3)
#    return p3
#end

#function getVertex(particle)
#    p3 = Array{Float64,1}(3)
#    valid = getVertex3(particle, p3)
#    return valid, p3
#end

#function getEndpoint(particle)
#    p3 = Array{Float64,1}(3)
#    valid = getEndpoint3(particle, p3)
#    return valid, p3
#end

#function getMomentumAtEndpoint(particle)
#    p3 = Array{Float64,1}(3)
#    valid = getMomentumAtEndpoint3(particle, p3)
#    return valid, p3
#end

# the navigator gets initialized with a collection
# it defers the actual work to the C++ implementation
#immutable LCRelationNavigator
#    relnav
#    fromType
#    toType
#    LCRelationNavigator(coll::TypedCollection) = _completNavigator(new(LCRelNav(coll.coll)))
#end
#function _completNavigator(nav)
#    nav.fromType = LCIOTypemap[nav.relnav.getFromType()]
#    nav.toType = LCIOTypemap[nav.relnav.getToType()]
#    nav
#end

# this ensures that the types are appropriately cast
#function getRelatedToObjects(nav::LCRelationNavigator, obj)
#    [CastOperator{nav.toType}.cast(x) for x in getRelatedToObjects(nav.relnav)]
#end

# this ensures that the types are appropriately cast
#function getRelatedFromObjects(nav::LCRelationNavigator, obj)
#    [CastOperator{nav.fromType}.cast(x) for x in getRelatedFromObjects(nav.relnav)]
#end

#getP3(x) = getPosition(x)

# converters to keep older code working
#const CalHits = Union{SimCalorimeterHit, CalorimeterHit, RawCalorimeterHit}
#
#function CalHit(h::CalHits)
#    p = getPosition(h)
#    E = getEnergy(h)
#    return CalHit(p[1], p[2], p[3], E)
#end
#function convert(::Type(CalHit), h::CalHits)
#    p = getPosition(h)
#    E = getEnergy(h)
#    return CalHit(p[1], p[2], p[3], E)
#end

#function printParameters(p::LCParameters)
#    println("strings:")
#    for k in getStringKeys(p, StringVec())
#        println(k, "\t", getStringVal(p, k))
#    end
#    println("floats:")
#    for k in getFloatKeys(p, StringVec())
#        println(k, "\t", getFloatVal(p, k))
#    end
#    println("ints:")
#    for k in getIntKeys(p, StringVec())
#        println(k, "\t", getIntVal(p, k))
#    end
#end

end # module


text:
module ClusterManagers

using Distributed
using Sockets

export launch, manage, kill, init_worker, connect
import Distributed: launch, manage, kill, init_worker, connect

worker_arg() = `--worker=$(Distributed.init_multi(); cluster_cookie())`

# PBS doesn't have the same semantics as SGE wrt to file accumulate,
# a different solution will have to be found
include("qsub.jl")
include("scyld.jl")
include("condor.jl")
include("slurm.jl")
include("affinity.jl")
include("elastic.jl")
include("lsf.jl")

end


text:
"""
# Stipple.Typography

Typography utilities for Stipple apps.
"""
module Typography

import Genie

include(joinpath("typography", "header.jl"))

end


text:
"""Synopsis
-index space functions: middle(a,b)
oldSchoolSwap(a, b)

-value space functions:

"""

#index space & value space function
function divideConquer(arr; start = 1) # 1 D vector of Float32 (64) #cross-reference with a binary-tree
    # check length
    # l is total actual length (of a list vector)
    l = length(arr) - 1 >= 0 ? copy(length(arr)) - 1 : -1  # actual size = length() -1 # l takes 0, 1, or more (cover them )
    if l == 0  # 3,3 length =0 , singlton , scalar
        #store leaf node #TODO:
        return 0 #scalar
    elseif l == 1 # 2,3 length = 1
        #find: first last
        aContent = arr[l] #unsure if first is list[1]
        bContent = arr[l+l]
        if aContent > bContent
            aContent, bContent = oldSchoolSwap(aContent, bContent)
            #store node with their values
        end
    else # l  > 1 # there are still ranges to be explored
        #isEven()
        #your best bet
        a = start #list()
        b = list
        #isEven(a,b)
        m1, m2, whole = middle(a, b) #doesn't matter if mid , mid+1 or below, above  # how much important is whole (in tree building wanna ensure 1 node for middle , or 2 actual middles)
        #got following ranges
        (a, m1) # 1 # right #index-space
        (m1 + 1, m2) #2 #middle range # index-space
        (m2 + 1, b)  # 3 # left #index-space

        #not good enough : need a method translating indices above to actual array values (sothat they can be called by this function recursively)
    end

end
#-----------
function isPositive(number)
    try
        if number >= 0
            return true
        elseif number < 0
            return false
        else
            throw(error("Unexpected error occured"))
        end
    catch UnexpectedError
        @error "ERROR: Unexpected value error" exception = (UnexpectedError, catch_backtrace())
    end
end
#----------
#euclidean distance, for values between a & b : a+b + 1 = a+ b + start -> this changes range end accordingly
euclidDist(a::Int64, b::Int64; s = 1::Int64) =
    0 <= a && 0 <= b && 0 <= s ? max(a, b) - min(a, b) + s : -1 #both a,b > 1 positive
euclidDist(2, 4) # = 2,3,4 = 3 numbers  (2)
euclidDist(2, 5) # 2,3,4,5 = 4 numbers (2 & 5 )
euclidDist(1, 1) #correct  = 1 number (between 1 & 1 )
euclidDist(1, 2) #this should be good to go

function calcLength(a = 2, b = 4; s = 1)
    start = euclidDist(a, b; s)
    return start
end
s = 1
calcLength(2, 4; s)

#----------------------------------

#sumInterval(1, 2) # 2 + 1 -1 = 2 #should be 2 #got 3 #ERROR  # use only when intention to calculate the mid
euclidDist(1,2)
#v = fillvector(1, 2, [3, 2, 5, 6])
#typeof(v) # got 3 (expected 3,2) not passing through upper bound #increase but wich ?

# an index space function
""" checks whether input arguments a & b are positive , calculates Interval length, anc check if it's positive
```inputs:
a: lower bound
b: upper bound
```
```outputs:
true:  if Interval length is even, meaning there is only 1 unique middle point
false: if Interval length is odd, meaning there is an exact fractional middle point, that needs further processing
```
"""
function isEven(a, b)
    try
        number = -1
        if a > 0 && b > 0
            number = b + a - 1 # calculates Interval length
            number > 0 && number % 2 == 0 ? true : false   #always exists (if conditions apply)

        else
            throw(error("Unexpected value error"))
        end
    catch UnexpectedError
       @error "ERROR: Unexpected value error" exception = (UnexpectedValueError, catch_backtrace())

        #passing function pointer to catch_backtrace
    end
end

isEven(number) = number > 0 && number % 2 == 0 ? true : false   #always exists (if conditions apply)

function oldSchoolSwap(a, b)
    tmp = b #store b
    b = a #a in b #a's ready
    a = tmp
    return a, b
end

sumInterval(a, b) = a > 0 && b > 0 ? abs(a + b) - 1 : -1 #formula: b + a - 1

#index space function
function middle(a, b) # b  + a -1
    try
        _sum = copy(sumInterval(a, b)) #    b + a - 1
        isItEven = copy(iseven(_sum)) # TODO: surround by a copy
        mid = copy(_sum / 2)  #  -1  #Idea: precalculate mid here (_sum /2 )
        isWhole = nothing
        below = -1
        above = -1
        if isItEven
            # 1 middle calculate
            mid = Int(mid)
            isWhole = true
            return mid, mid + 1, isWhole
        elseif !isItEven
            # calculate fractionalMid
            below = floor(mid)
            above = ceil(mid)
            isWhole = false
            return below, above, isWhole # the differenece is still 1, only way to discriminate is by using isWhole
        else
            throw(error("Unexpected error occured"))
        end
    catch UnexpectedError
        @error "ERROR: Unexpected error occured" exception = (UnexpectedError, catch_backtrace())

    end
end

middle(1, 2)

#= this won't make it divideConquer anymore
function fillvector(a, b, arr; start = 1)
    loopstart = start # if starting index is 1
    len_ab = euclidDist(a, b) # copy(sumInterval(a, b)) #b - a # i thought b + a - 1 i.e. 1 , 3: 3+1 - 1 =3 #problem here when changing range , problem discovered need distance between 2 values:
    condition = isPositive(len_ab)
    vector = []
    try
        if condition
            n = len_ab # b - a # wrong range #should be 2  #got 3 # got 2
            upperbound = copy(n + loopstart) # add 1 if starting index is 1
            for i = 1:n #stays n (whether starts at 0 or 1) #in every for loop, ensure you're # tried upperbound = n (original + 1) it prints 3 (instead of 2)#infer n is independent of range start
                pushfirst!(vector, arr[i]) # fills corresponding range with vector values
            end
        elseif !condition
            return -1
        else
            throw(error("Unexpected error occured"))
        end
    catch UnexpectedError
        @error "ERROR: UnexpectedError occured " exception = (UnexpectedError, catch_backtrace())
    end
    return vector
end

fillvector(1, 2, [3, 2, 1, 5]) # working

fillvector(2, 3, [3, 2, 1, 5])
=#

#----

#index & value space
"""creates an indices between points a & b """
function indexify(a = 1, b = 3,arr[1,3,4,5]; s = 1) # arr is not reading it
lenA = length(arr)
    rangeIndicies = []
    #checks #s is an important input argument, its failure may result in failure of every formula including it
    #s != 0 || s != 1 ? s = 1 : s = 1 #it's forces s value to take 1 (if it was not set by user to either 0 or 1 )
    try
        if s != 0 || s != 1 # s can either be 0 or 1 , if not, force s =1
            s = 1 # set it to 1
            for i ∈ 1:lenA+s-1  # idea= keep fixating at start i=1 ,last length(arr)+s-1 i=1 #strict condition: requires s=i=1 to work well
                rangeIndicies = push!(rangeIndicies, i) #println(i) #prints range correctly 1,2,3 (accordingly)
            end
        else
            throw(error("Unexpected error"))
        end
    catch UnexpectedError
        @error "ERROR: UnexpectedError: " exception = (UnexpectedError, catch_backtrace())
    end
    rangeIndicies
end

indexify(1,3)

#index space: a=1, b=3
#value space arr
"""creates an indices between points a & b """
function indexify(a = 1, b = 3,arr=[1,3,4,5]; s = 1)
    lenA =  length(arr)
    lenV = euclidDist(a,b)

    rangeIndicies = []
    #checks #s is an important input argument, its failure may result in failure of every formula including it
    #s != 0 || s != 1 ? s = 1 : s = 1 #it's forces s value to take 1 (if it was not set by user to either 0 or 1 )
    try
          if s != 0 || s != 1 # s can either be 0 or 1 , if not, force s =1
            s = 1 # set it to 1  #thik 1 move ahead
            for i ∈ 1:lenA+s-1  # idea= keep fixating at start i=1 ,last length(arr)+s-1 i=1 #strict condition: needs s=i=1 to work well
                rangeIndicies = push!(rangeIndicies, i) #println(i) #prints range correctly 1,2,3 (accordingly)
            end
        else
            throw(error("Unexpected value for s"))
        end
    catch UnexpectedValueError
        @error "ERROR: Unexpected value for s" exception = (UnexpectedValueError, catch_backtrace())
    end
    rangeIndicies
end

euclidDist(1,3)

# index space value
function indexify(a = 1, b = 3,arr=[1,3,4,5]; s = 1)
    lenA =  length(arr)
    lenV = euclidDist(a,b)
    v = [lenV]
    #levV = []
    rangeIndicies = []
    #checks #s is an important input argument, its failure may result in failure of every formula including it
    #s != 0 || s != 1 ? s = 1 : s = 1 #it's forces s value to take 1 (if it was not set by user to either 0 or 1 )
    try
        if s != 0 || s != 1 # s can either be 0 or 1 , if not, force s =1
            s = 1 # set it to 1  #thik 1 move ahead
            for i ∈ 1:lenA+s-1  # idea= keep fixating at start i=1 ,last length(arr)+s-1 i=1 #strict condition: needs s=i=1 to work well
                rangeIndicies = push!(rangeIndicies, i) #println(i) #prints range correctly 1,2,3 (accordingly)
            end
        else
            throw(error("Unexpected value for s"))
        end
    catch UnexpectedValueError
        @error "ERROR: Unexpected value for s" exception = (UnexpectedValueError, catch_backtrace())

    end
    rangeIndicies
end

#-------------
res = indexify(1, 2)
typeof(res)
#what's distance between: 2,1 = 2 -1 + 1

function range2vector(range)
    lenA = copy(length(range)) # 3
     i = 1
    next = i + 1
    v = []

    for i ∈ 1:lenA
        if next <= lenA && i <= lenA - 1
            v = push!(v, collect(i:next))
            continue # once reached m gi ti function definition
        else
            @error UnexpectedError exception = (UnexpectedError, back_trace) #throw(error("UnexpectedError"))
        end

    end

    return v
end

res = range2vector(1:3)
typeof(res)
length(res)

#----
# last function #uses
function indicies2values(v = [], a = [1, 4, 5, 6]) # v can be a pair v[1](1)
    lenV = copy(length(v))
    next = -1
    heap = []
    subheap = []
    for i ∈ 1:lenV
        vi = copy(length(v[i]))
        if vi == 1
            push!(heap, a[v[i][1]])# got index, transofrmed into an array #TODO check isnumeric, isWhole
            _set = []
        elseif vi > 1 # 2, 3 or more range  #my bread & butter
            for j ∈ 1:vi
                next = j + 1 #warning : potential out of bounds
                if next <= vi
                    genVector = collect(a[j]:a[next])
                    Union!(_set, _set)
                    Union!(subheap, _set)

                    _set = Union!(_set, [[], genVector])

                    _set = Union!(_set, [subheap, genVector]) # j=1,2

                end
            end

        else
            @error UnexpectedError exception = (UnexpectedError, back_trace)
        end
    end
end

v = [[1, 2], [2]]
lenVi = copy(length(v))
lenVi = copy(length(enumerate(v)))
vi= []
vi = [lenVi]
typeof(vi)

function innerloop(a = [3, 2, 1, 4, 5], v = [[1, 2], [2]])
    lenVi = copy(length(v))
    vi=[lenVi]; #pushfirst!(vi,[1,2]); pushfirst!(vi,[2]) the goal
    length(vi)
    #vi =

   # for j ∈ 1:lenVi # vi = [[1, 2], [2]]
        #if   next <=
     j=1
        next = j + 1 #warning : potential out of bounds
       if next <=lenVi #vi
            genVector = collect(a[j]:a[next])
          #  Union!(_set, _set)
          #  Union!(subheap, _set)
            _set= UnionAll{_set,genVector}
          #  _set = Union!(_set, [[], genVector])
         #   _set = Union!(_set, [subheap, genVector]) # j=1,2
        return _set
      #  end
 #   end

end

innerloop()

#-----------------------------
genVector = collect(a[j]:a[next])
Union!(_set, _set)
Union!(subheap, _set)
_set = Union!(_set, [[], genVector])

#credits: @tpapp
#Yes, unfortunately array syntax ([...]) is one of the few places where whitespace has special syntactic significance.

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner


text: # print("ok")


text:
# Included in Common.jl

# for parsing
function createfeatures(sentences, featdict)

	function createedgefeat(head::Word, mod::Word)
		# Feature Template
		p_word = 	"p-word=$(head.form)"
		p_pos = 	"p-pos=$(head.postag)"
		c_word = 	"c-word=$(mod.form)"
		c_pos = 	"c-pos=$(mod.postag)"

		feature_template = [
			# Basic Features
			join([p_word, p_pos, c_word, c_pos], " "),
			join([p_pos, c_word, c_pos], " "),
			join([p_word, c_word, c_pos], " "),
			join([p_word, p_pos, c_pos], " "),
			join([p_word, p_pos, c_word], " "),
			join([p_word, c_word], " "),
			join([p_pos, c_pos], " "),
			join([p_word, p_pos], " "),
			join([c_word, p_pos], " "),
			p_word,
			p_pos,
			c_word,
			c_pos
		]

		# one-hot features
		# simply keep a list of feat indices, instead of k/v pairs
		fv = Int[]
		for feat in feature_template
			if haskey(featdict, feat)
				push!(fv, featdict[feat])
			end
		end

		return fv
	end

	root = Word(
			0, 			# id
			"<ROOT>", 	# form
			"<ROOT>", 	# lemma
			"ROOT", 	# cpostag
			"ROOT", 	# postag
			[], 		# feats
			-1, 		# head
			"ROOT" 		# deprel
		)

	edge_feats_list = Dict{(Int,Int), Vector{Int}}[] # a list of `edge_feats` (set of feature vectors per sentence)

	info("Create features and vectorize sentences ...")
	for sent in sentences
		n = length(sent)
		# edge_feats[(index_head, index_mod)] = "a list of feat indices"
		edge_feats = Dict{(Int,Int), Vector{Int}}()
		sizehint(edge_feats, n*n)

		# edges from ROOT
		for index_mod = 1:n
			edge_feats[(0, index_mod)] = createedgefeat(root, sent[index_mod])
		end

		# all other edges
		for index_head = 1:n
			for index_mod = 1:n
				if index_head != index_mod # no self-loop edge
					edge_feats[(index_head, index_mod)] = createedgefeat(sent[index_head], sent[index_mod])
				end
			end
		end

		push!(edge_feats_list, edge_feats)
	end

	return edge_feats_list
end

# for training
function createfeatures(sentences)

	function getfeatindex!(feat::String)
		if haskey(featdict, feat)
			return featdict[feat]
		else
			n_feats += 1
			featdict[feat] = n_feats
			return n_feats
		end
	end

	function createedgefeat!(head::Word, mod::Word)
		# Feature Template
		p_word = 	"p-word=$(head.form)"
		p_pos = 	"p-pos=$(head.postag)"
		c_word = 	"c-word=$(mod.form)"
		c_pos = 	"c-pos=$(mod.postag)"

		feature_template = [
			# Basic Features
			join([p_word, p_pos, c_word, c_pos], " "),
			join([p_pos, c_word, c_pos], " "),
			join([p_word, c_word, c_pos], " "),
			join([p_word, p_pos, c_pos], " "),
			join([p_word, p_pos, c_word], " "),
			join([p_word, c_word], " "),
			join([p_pos, c_pos], " "),
			join([p_word, p_pos], " "),
			join([c_word, p_pos], " "),
			p_word,
			p_pos,
			c_word,
			c_pos
		]

		# one-hot features
		# simply keep a list of feat indices, instead of k/v pairs
		fv = Int[]
		for feat in feature_template
			index = getfeatindex!(feat)
			push!(fv, index)
		end

		return fv
	end

	root = Word(
			0, 			# id
			"<ROOT>", 	# form
			"<ROOT>", 	# lemma
			"ROOT", 	# cpostag
			"ROOT", 	# postag
			[], 		# feats
			-1, 		# head
			"ROOT" 		# deprel
		)

	featdict = Dict{String, Int}() # feature index dictionary
	n_feats = 0 # number of features
	edge_feats_list = Dict{(Int,Int), Vector{Int}}[] # a list of `edge_feats` (set of feature vectors per sentence)

	info("Create features and vectorize sentences ...")
	for sent in sentences
		n = length(sent)
		# edge_feats[(index_head, index_mod)] = "a list of feat indices"
		edge_feats = Dict{(Int,Int), Vector{Int}}()
		sizehint(edge_feats, n*n)

		# edges from ROOT
		for index_mod = 1:n
			edge_feats[(0, index_mod)] = createedgefeat!(root, sent[index_mod])
		end

		# all other edges
		for index_head = 1:n
			for index_mod = 1:n
				if index_head != index_mod # no self-loop edge
					edge_feats[(index_head, index_mod)] = createedgefeat!(sent[index_head], sent[index_mod])
				end
			end
		end

		push!(edge_feats_list, edge_feats)
	end
	info("Total number of features: $n_feats")

	return edge_feats_list, featdict, n_feats
end


text:
#!/bin/bash
# -*- mode: julia -*-
#=
exec "${JULIA:-julia}" "$@" ${BASH_SOURCE[0]}
=#

using Pkg
Pkg.activate(@__DIR__)
Pkg.add("MacroTools")
#Pkg.develop(PackageSpec(name="PyCall", path=dirname(@__DIR__)))
#ENV["PYTHON"] = ""
#Pkg.build("PyCall")
Pkg.add("PyCall")
Pkg.build("PyCall")
Pkg.add("PyPlot")
#Pkg.add("Plots")
pkg"add Plots#master"
pkg"add PotentialFlow#master"
Pkg.activate()

using PackageCompiler
#
sysout, _curr_syso = compile_incremental(
    joinpath(@__DIR__, "Project.toml"),
    joinpath(@__DIR__, "precompile.jl")
)
#
#sysout, _curr_syso = compile_incremental(
#    joinpath(@__DIR__, "precompile.jl"),
#)

pysysout = joinpath(@__DIR__, basename(sysout))
cp(sysout, pysysout, force=true)

write(joinpath(@__DIR__, "_julia_path"), Base.julia_cmd().exec[1])

@info "System image: $pysysout"

max_stars_repo_name: jdeldre/BookExamples

max_issues_repo_name: jdeldre/BookExamples

max_forks_repo_name: jdeldre/BookExamples


text:
module BraidChains

using DemeNet: ID, Signer, Deme, Signer
#using PeaceVote.Plugins: AbstractChain

using PeaceVote: AbstractChain

include("Types/Types.jl")
include("Braiders/Braiders.jl")
include("Ledgers/Ledgers.jl")
include("Analysis/Analysis.jl") ### There should not be any problems with this
include("Recorders/Recorders.jl")

import .Types: Proposal, Vote
import .Braiders: braid!, BraiderConfig, Braider, Mixer
import .Recorders: RecorderConfig, Recorder, record
import .Ledgers: Ledger, record!, getrecord, readbootrecord, writebootrecord, load
import .Analysis: normalcount, voters, attest, members, proposals ### A counting strategy depends on the proposal, thus normalcount could be replaceed just with count

struct BraidChainConfig{T} # Perhaps BraidChainRemote
    server::ID
    mixerport::T
    syncport::T
    braider::BraiderConfig{T}
    recorder::RecorderConfig{T}
end

struct BraidChain <: AbstractChain
    config::BraidChainConfig # I don't need a parameter for the BraidChain itself
    deme::Deme
    ledger
end

BraidChain(config::BraidChainConfig,deme::Deme) = BraidChain(config,deme,Ledger(deme.spec.uuid))

import PeaceVote: record
record(config::BraidChainConfig,data) = record(config.recorder,data)
record(chain::BraidChain,data) = record(chain.config,data)

Recorder(chain::BraidChain,braider::Braider,signer::Signer) = Recorder(chain.config.recorder,chain.deme,chain.ledger,braider,signer)

struct BraidChainServer
    mixer
    synchronizer
    braider
    recorder
end

function BraidChainServer(chain::BraidChain,server::Signer)

    config = chain.config

    mixer = Mixer(config.mixerport,chain.deme,server)
    synchronizer = @async Ledgers.serve(config.syncport,chain.ledger)
    braider = Braider(config.braider,chain.deme,server)
    recorder = Recorder(chain,braider,server)

    return BraidChainServer(mixer,synchronizer,braider,recorder)
end

load(chain::BraidChain) = load(chain.ledger)

import PeaceVote: sync!
#sync!(chain::BraidChain,syncport) = Ledgers.sync!(chain.ledger,syncport)
#sync!(chain::BraidChain,config::BraidChainConfig) = sync!(chain,config.syncport)

sync!(chain::BraidChain) = Ledgers.sync!(chain.ledger,chain.config.syncport)

import PeaceVote: count
count(index::Int,chain::BraidChain) = normalcount(index,chain.deme,chain.ledger)

braid!(chain::BraidChain,newvoter::Signer,oldvoter::Signer) = braid!(chain.config.braider,chain.deme,newvoter,oldvoter)

attest(chain::BraidChain) = attest(load(chain.ledger),chain.deme.notary)

end


text:
module CellArrays

using Gridap
using Gridap.Helpers
using Gridap.CellValues: _test_iter_cell_value
using Gridap.CellValues: _test_index_cell_value

export CellArray
export CellMatrix
export CellVector

export IterCellArray
export IterCellMatrix
export IterCellVector

export IndexCellArray
export IndexCellMatrix
export IndexCellVector

export test_iter_cell_array
export test_index_cell_array

import Base: collect

# Iterable cell Arrays

const IterCellArray{T,N,A<:AbstractArray{T,N}} = IterCellValue{A}

const IterCellVector{T} = IterCellArray{T,1}

const IterCellMatrix{T} = IterCellArray{T,2}

# Indexable cell arrays

const IndexCellArray{T,N,A<:AbstractArray{T,N},D} = IndexCellValue{A,D}

const IndexCellVector{T} = IndexCellArray{T,1}

const IndexCellMatrix{T} = IndexCellArray{T,2}

# Cell Arrays

const CellArray{T,N} = Union{IterCellArray{T,N},IndexCellArray{T,N}}

const CellVector{T} = CellArray{T,1}

const CellMatrix{T} = CellArray{T,2}

collect(a::CellArray) = [ copy(ai) for ai in a ]

# Testers

function test_iter_cell_array(
  icv::CellArray{T,N},
  a::AbstractArray{<:AbstractArray{T,N}}) where {T,N}
  _test_iter_cell_value(icv,a)
end

function test_index_cell_array(
  icv::IndexCellArray{T,N},
  a::AbstractArray{<:AbstractArray{T,N}}) where {T,N}
  _test_index_cell_value(icv,a)
end

end # module


text:
# Adapted from:
# SOSDEMO2 --- Lyapunov Function Search
# Section 3.2 of SOSTOOLS User's Manual

for solver in sdp_solvers

  @polyvar x[1:3]

  # Constructing the vector field dx/dt = f
  f = [-x[1]^3-x[1]*x[3]^2,
      -x[2]-x[1]^2*x[2],
      -x[3]+3*x[1]^2*x[3]-3*x[3]/(x[3]^2+1)]

  m = Model(solver = solver)

  # The Lyapunov function V(x):
  Z = x.^2
  @polyvariable m V Z

  @polyconstraint m V >= sum(x.^2)

  # dV/dx*(x[3]^2+1)*f <= 0
  P = dot(differentiate(V, x), f)*(x[3]^2+1)
  @polyconstraint m P <= 0

  status = solve(m);

  flyapunov = getvalue(V)
  println("Lyapunov function using $(typeof(solver))")

  println(flyapunov)
  println(status)

  #status --> :Optimal

  #removemonomials(getvalue(V), Z) --> zero(Polynomial{true, Float64})
end

max_stars_repo_name: polislizarralde/jsos

max_issues_repo_name: polislizarralde/jsos

max_forks_repo_name: polislizarralde/jsos


text:
using JLD2
using FileIO

# tensor for classical 2-d model
const isingβc = log(1+sqrt(2))/2

"""
    tensorfromclassical(h::Matrix)

given a classical 2-body hamiltonian `h`, return the corresponding tensor
for use in e.g. `trg` for a two-dimensional square-lattice.

# Example
```julia
julia> model_tensor(Ising(),β) ≈ tensorfromclassical([β -β; -β β])

true
```
"""
function tensorfromclassical(ham::Matrix)
    wboltzmann = exp.(ham)
    q = sqrt(wboltzmann)
    ein"ij,ik,il,im -> jklm"(q,q,q,q)
end

"""
    model_tensor(::Ising,β)

return the isingtensor at inverse temperature `β` for a two-dimensional
square lattice tensor-network.
"""
function model_tensor(::Ising, β::Real)
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 1], 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    mag_tensor(::Ising,β)

return the operator for the magnetisation at inverse temperature `β`
at a site in the two-dimensional ising model on a square lattice in tensor-network form.
"""
function mag_tensor(::Ising, β)
    a = reshape(Float64[1 0 0 0; 0 0 0 0; 0 0 0 0; 0 0 0 -1] , 2,2,2,2)
    cβ, sβ = sqrt(cosh(β)), sqrt(sinh(β))
    q = 1/sqrt(2) * [cβ+sβ cβ-sβ; cβ-sβ cβ+sβ]
    ein"abcd,ai,bj,ck,dl -> ijkl"(a,q,q,q,q)
end

"""
    energy_tensor(::Ising,β)

return the operator for the energy at inverse temperature `β`
at a site in the two-dimensional ising model on a square lattice in tensor-network form.
"""
function energy_tensor(::Ising, β)
    ham = [-1 1;1 -1]
    wboltzmann = exp.(-β .* ham)
    wenergy = ham .* wboltzmann
    qb = sqrt(wboltzmann)
    ein"as,si,bi,ci,di -> abcd"(qb^(-1),wenergy,qb,qb,qb)
end

"""
    vumps_env(model<:HamiltonianModel, β, D)

return the vumps environment of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `D` as calculated with
vumps. Save `env` in file `./data/model_β_D.jld2`. Requires that `model_tensor` are defined for `model`.
"""
function vumps_env(model::MT, β, D; tol=1e-10, maxiter=20, verbose = false) where {MT <: HamiltonianModel}
    M = model_tensor(model, β)
    mkpath("./data/")
    chkp_file = "./data/$(model)_β$(β)_D$(D).jld2"
    if isfile(chkp_file)
        rt = SquareVUMPSRuntime(M, chkp_file, D; verbose = verbose)
    else
        rt = SquareVUMPSRuntime(M, Val(:random), D; verbose = verbose)
    end
    env = vumps(rt; tol=tol, maxiter=maxiter, verbose = verbose)
    # save(chkp_file, "env", env) # if forward steps is too small, backward go to this way will make mistake!!!
    return env
end

"""
    Z(env::SquareVUMPSRuntime)

return the partition function of the `env`.
"""
function Z(env::SquareVUMPSRuntime)
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    AC = ein"asc,cb -> asb"(AL,C)
    z = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL,AC,M,FR,conj(AC))[]
    λ = ein"αcβ,βη,ηcγ,αγ ->"(FL,C,FR,conj(C))[]
    return z/λ
end

"""
    magnetisation(env::SquareVUMPSRuntime, model::MT, β)

return the magnetisation of the `model`. Requires that `mag_tensor` are defined for `model`.
"""
function magnetisation(env::SquareVUMPSRuntime, model::MT, β) where {MT <: HamiltonianModel}
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    Mag = mag_tensor(model, β)
    AC = ein"asc,cb -> asb"(AL,C)
    mag = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL,AC,Mag,FR,conj(AC))[]
    λ = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL,AC,M,FR,conj(AC))[]

    return abs(mag/λ)
end

"""
    energy(env::SquareVUMPSRuntime, model::MT, β)

return the energy of the `model` as a function of the inverse
temperature `β` and the environment bonddimension `D` as calculated with
vumps. Requires that `model_tensor` are defined for `model`.
"""
function energy(env::SquareVUMPSRuntime, model::MT, β) where {MT <: HamiltonianModel}
    M,AL,C,FL,FR = env.M,env.AL,env.C,env.FL,env.FR
    Ene = energy_tensor(model, β)
    AC = ein"asc,cb -> asb"(AL,C)
    energy = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL,AC,Ene,FR,conj(AC))[]
    λ = ein"αcβ,βsη,cpds,ηdγ,αpγ ->"(FL,AC,M,FR,conj(AC))[]

    return energy/λ*2 # factor 2 for counting horizontal and vertical links
end

"""
    magofβ(::Ising,β)

return the analytical result for the magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofβ(::Ising, β) = β > isingβc ? (1-sinh(2*β)^-4)^(1/8) : 0.

"""
    magofdβ(::Ising,β)

return the analytical result for the derivative of magnetisation at inverse temperature
`β` for the 2d classical ising model.
"""
magofdβ(::Ising, β) = β > isingβc ? (coth(2*β)*csch(2*β)^4)/(1-csch(2*β)^4)^(7/8) : 0.


text:
mutable struct MovingCells
    track_id::Int64
    coordinates::Array{Float64,2}
    time::Array{Float64,1}
    smooth_n_points::Int64
    smooth_coordinates::Array{Float64,2}
    smooth_time::Array{Float64,1}
    T::Array{Float64,2}
    N::Array{Float64,2}
    B::Array{Float64,2}
    curvature::Array{Float64,1}
    torsion::Array{Float64,1}
    speed::Array{Float64,1}
    original_shapes::Array{Mesh,1}
    reoriented_shapes::Array{Mesh,1}
    Spharm_Cx::Array{Float64,2}
    Spharm_Cy::Array{Float64,2}
    Spharm_Cz::Array{Float64,2}

    function MovingCells()
        return new()
    end

    function MovingCells(
        track_id::Int64,
        coordinates::Array{Float64,2},
        time::Array{Float64,1},
    )
        cell = new()
        cell.track_id = track_id
        cell.coordinates = coordinates .- coordinates[1, :]'
        cell.time = time

        return cell
    end
end

function get_cell(
    df::DataFrame,
    track_id::Int;
    t0 = -1.0,
    t1 = -1.0,
)::MovingCells

    if t0 == t1 == -1.0
        selected_df = filter(row -> row[:trackid] == track_id, df)
    else
        @assert t1 > t0 >= 1.0
        selected_df =
            filter(row -> (row.trackid == track_id) & (t0 <= row.t <= t1), df)
    end
    coordinates = convert(Matrix, selected_df[!, [:x, :y, :z]])
    time = collect(1.0:size(coordinates, 1))
    cell = MovingCells(track_id, coordinates, time)
    return cell
end

function add_moving_frame_data!(cell::MovingCells; step = 0.01)
    coordinates = cell.coordinates
    smooth_coordinates = smoothing(coordinates, step = step)
    cell.smooth_coordinates = smooth_coordinates
    cell.smooth_n_points = size(smooth_coordinates, 1)
    cell.smooth_time =
        range(minimum(cell.time), stop = maximum(cell.time), step = step)
    T, N, B, curvature, torsion, speed =
        parallel_transport(smooth_coordinates)
    cell.T, cell.N, cell.B = T, N, B
    cell.curvature, cell.torsion, cell.speed = curvature, torsion, speed
end

function extract_shape_data(df; data_dir = pwd())
    @assert length(unique(df.trackid)) == 1
    cell_id = df.cellid
    gc_cell_id = df.gc_cellid
    t = convert.(Int, df.t)
    shapes = Vector{Mesh}(undef, length(gc_cell_id))
    for i in eachindex(gc_cell_id)
        fname = joinpath(
            data_dir,
            "t" * lpad(t[i], 4, "0"),
            "cell" * lpad(gc_cell_id[i], 5, "0") * ".obj",
        )
        # @show fname
        shape = load(fname) |> mesh_centering |> volume_normalizing

        shapes[i] = load(fname)
    end
    return shapes
end

function add_shape_data!(cell::MovingCells,df;data_dir = pwd())
    shapes = extract_shape_data(df,data_dir=data_dir)
    cell.original_shapes = shapes
    smooth_time = cell.smooth_time
    step = smooth_time[2] - smooth_time[1]
    selected_time = [
                    1 + round(Int, 1 / step) * x
                    for
                    x in range(
                        0,
                        stop = div(length(smooth_time), round(Int, 1 / step)),
                    )
                ]
    T = cell.T[selected_time, :]
    N = cell.N[selected_time, :]
    B = cell.B[selected_time, :]
    reoriented_shapes = [mesh_reorientation(shapes[i], T[i, :], N[i, :], B[i, :]) for i in eachindex(shapes)]
    cell.reoriented_shapes = reoriented_shapes
end

function add_spharm_coef_data!(cell::MovingCells)
    reoriented_shape = cell.reoriented_shapes
    spheres = sh.spherization.(reoriented_shape, 100)
    spharm_coefs = [
        sh.spharm_descriptor(reoriented_shape[i], spheres[i].vertices)
        for i in eachindex(reoriented_shape)
    ]
    cell.Spharm_Cx =
        sh.array_of_array_to_matrix([coef[:Cx] for coef in spharm_coefs])
    cell.Spharm_Cy =
        sh.array_of_array_to_matrix([coef[:Cy] for coef in spharm_coefs])
    cell.Spharm_Cz =
        sh.array_of_array_to_matrix([coef[:Cz] for coef in spharm_coefs])
end


text:
#using Revise

using Adaptive
using Distributed

@info "Testing ask! and tell! methods"
function simpledriver!(f, learner,step)
    while !(learner.loss() < step)
        xi = ask!(learner,true)
        yi = f(xi)
        tell!(learner,(xi,yi))
    end
end

learner1d = AdaptiveLearner1D((0,1))
simpledriver!(x->exp(-x^2),learner1d,0.1)

learner2d = AdaptiveLearner2D([(0,1),(0,1)])
simpledriver!(x->exp(-x[1]^2 - x[2]^2),learner2d,0.1)

@info "Testing atribute access"

@show learner1d.loss
@show learner1d.x
@show learner1d.y
@show learner1d.data

@show learner2d.loss
@show learner2d.points
@show learner2d.vertices
@show learner2d.values
@show learner2d.ip


text:
# This file was generated, do not modify it. # hide
using LinearAlgebra, Random
Random.seed!(555)
a = randn(5)
round(norm(a), sigdigits=4)

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
#
# function cholesky_U_deravitive(A::Symmetric{T,MMatrix{N,N,T,N2}}, Ȧ::Symmetric{T,MMatrix{N,N,T,N2}}) where {N,T,N2}
#     # Consider the factorization of a PD matrix A: A = U' * U
#     # The derivative of this factorization is Ȧ = U̇' * U + U' * U̇
#     # This function finds U̇ given Ȧ using the algorithm from: https://homepages.inf.ed.ac.uk/imurray2/pub/16choldiff/choldiff.pdf
#
#     cf = cholesky(A)
#     U = cf.U
#     U⁻¹ = inv(U)
#     term_in = U⁻¹' * Ȧ * U⁻¹
#     @inbounds begin
#         for k = 1:(N+1):N2
#             term_in[k] *= 0.5
#         end
#     end
#     U̇ = UpperTriangular(term_in) * U
#     return U, U⁻¹, U̇
# end


text:
"""
    SplineWrapper{S <: Spline}

Abstract type representing a type that wraps a [`Spline`](@ref).

Such a type implements all common operations on splines, including evaluation,
differentiation, etc…
"""
abstract type SplineWrapper{S <: Spline} end

"""
    spline(w::SplineWrapper) -> Spline

Returns the [`Spline`](@ref) wrapped by the object.
"""
spline(w::SplineWrapper) = w.spline

Base.eltype(::Type{<:SplineWrapper{S}}) where {S} = eltype(S)

(I::SplineWrapper)(x) = spline(I)(x)
Base.diff(I::SplineWrapper, etc...) = diff(spline(I), etc...)
Base.:*(op::Derivative, I::SplineWrapper) = op * spline(I)

for f in (:basis, :order, :knots, :coefficients, :integral)
    @eval $f(I::SplineWrapper) = $f(spline(I))
end

max_stars_repo_name: jipolanco/BasisSplines

max_issues_repo_name: jipolanco/BasisSplines

max_forks_repo_name: jipolanco/BasisSplines


text:
#
# Least square solver
#
# author: Atsushi Sakai
#

export solve_least_square
export solve_multi_objective_least_square
export solve_constrained_least_square
export solve_nonlinear_least_square_with_newton_raphson
export solve_nonlinear_least_square_with_gauss_newton
export solve_nonlinear_least_square_with_levenberg_marquardt
export solve_constrained_nonlinear_least_square_with_augmented_lagragian

eye(T::Type, n) = Diagonal{T}(I, n)
eye(n) = eye(Float64, n)

"""
	solve least square
	xhat = argmin(|Ax = b|^2)

	All are same solution
	- xhat = inv(A’*A)*(A’*b)
	- xhat = pinv(A)*b
	- Q,R = qr(A); xhat = inv(R)*(Q’*b)
	- xhat = A\b
"""
function solve_least_square(A, b)
	return A\b
end

"""
	solve constrained least_square
	xhat = argmin(|Ax = b|^2) s.t. Cx = d
"""
function solve_least_square(A,b,C,d)
	m, n = size(A)
    p, n = size(C)
    G = A'*A  # Gram matrix
    KKT = [2*G C'; C zeros(p,p)]  # KKT matrix
    xzhat = KKT \ [2*A'*b; d]
    return xzhat[1:n,:]
end

"""
	solve multi objective least square
	xhat = argmin(λ_1|Ax = b|^2+λ_2|Ax = b|^2...)
"""
function solve_multi_objective_least_square(As, bs, lambdas)
   k = length(lambdas);
   Atil = vcat([sqrt(lambdas[i])*As[i] for i=1:k]...)
   btil = vcat([sqrt(lambdas[i])*bs[i] for i=1:k]...)
   return solve_least_square(Atil, btil)
end

"""
	solve nonlinear least square with newton-raphson method

	The inputs have to be length(x) == length(f(x)).
	If it is not, you cant use solve_nonlinear_least_square_with_gauss_newton

	xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_newton_raphson(
		f, Df, x1; kmax = 20, tol = 1e-6)

	x=x1
	@assert length(x) == length(f(x))
	for k = 1:kmax
		fk = f(x)
		dx = fk / Df(x)
		if norm(dx) < tol break end;
		x = x - dx
	end

	return x
end

"""
	solve nonlinear least square with gauss-newton method

	xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_gauss_newton(
		f, Df, x0;
	   	kmax = 20, tol = 1e-6)
	x=x0
   	obj = zeros(0,1)
   	residuals = zeros(0,1)
	for k=1:kmax
		fx = f(x)
		dfx = Df(x)
		obj = [obj; norm(fx)^2]
		res = norm(2*dfx'*fx)
		residuals = [residuals; res]
		if res < tol break end;

		tmp = (dfx'*fx)
		if isa(tmp, Array)
			dx = (dfx'*dfx) \ tmp
		else
			dx = (dfx'*dfx) \ [tmp]
		end

		x = x - dx
	end

	return x, Dict([ ("objectives", obj), ("residuals", residuals)])
end

"""
	solve nonlinear least square with levenberg marquardt

	xhat = argmin(|f(x)|^2)
"""
function solve_nonlinear_least_square_with_levenberg_marquardt(
		f, Df, x0, lambda0;
	   	kmax = 20, tol = 1e-6, lamr_n = 0.8, lamr_p = 2.0)

	n = length(x0)
	x=x0
	lambda = lambda0
	obj = zeros(0,1)
	residuals = zeros(0,1)
	xhist = x0'
	for k = 1:kmax
		fx = f(x)
		Dfk = Df(x)
		obj = [obj; norm(fx)^2]
		res = norm(2*Dfk'*fx)
		residuals = [residuals; res]
		if res < tol break end;
		tmp = Dfk'*fx
		if isa(tmp, Array)
			dx = (Dfk'*Dfk+lambda*eye(n)) \ tmp
		else
			# dx = (Dfk'*Dfk+lambda*eye(n)) \ [tmp]
			dx = (Dfk'*Dfk+lambda*1.0) \ [tmp]
		end
		xt = x - dx
		if norm(f(xt)) < norm(fx)
			lambda = lamr_n*lambda
			x = xt
		else
			lambda = lamr_p*lambda
		end
		xhist = [xhist; x']
	end

	return x, Dict([ ("objectives", obj), ("residuals", residuals), ("x_history", xhist)])
end

"""
Solve constrained nonlinear least square with augmentedl lagrangian method

xhat = argmin(|f(x)|^2) s.t g(x) = 0

"""
function solve_constrained_nonlinear_least_square_with_augmented_lagragian(
		f, Df, g, Dg, x1, lambda1;
	   	kmax = 100, feas_tol = 1e-4, oc_tol = 1e-4)
	x=x1
	z = zeros(length(g(x)))
	mu=1.0
	feas_res = [norm(g(x))]
	oc_res = [norm(2*Df(x)'*f(x) + 2*mu*Dg(x)'*z)]
	lm_iters = zeros(Int64,0,1)

	for k=1:kmax
		F(x) = [f(x); sqrt(mu)*(g(x) + z/(2*mu))]
		DF(x) = [Df(x); sqrt(mu)*Dg(x)]
		x, hist = solve_nonlinear_least_square_with_levenberg_marquardt(F, DF, x, lambda1, tol=oc_tol)
        z = z + 2*mu*g(x)
        feas_res = [feas_res; norm(g(x))]
        oc_res = [oc_res; hist["residuals"][end]]
        lm_iters = [lm_iters; length(hist["residuals"])]
        if norm(g(x)) < feas_tol
			break
		end
        mu = (norm(g(x)) < 0.25*feas_res[end-1]) ? mu : 2*mu
    end

    return x, z, Dict([ ("lm_iterations", lm_iters),
         ("feas_res", feas_res), ("oc_res", oc_res)])
end


text:
"""
    Exo1Real{T <: Real} = Parabolic{T}

A 1-exo-real number is a parabolic Cayley-Dickson construct with real numbers.
"""
const Exo1Real{T <: Real} = Parabolic{T}

function Exo1Real(a::T, b::T) where T <: Real
    Exo1Real{T}(a, b)
end

function Exo1Real(a::T) where T <: Real
    Exo1Real{T}(a, zero(T))
end

function Exo1Real(a::Real, b::Real)
    Exo1Real(promote(a, b)...)
end

function selfstar(z::Exo1Real)
    (z + star(z)) / 2
end

function anti_selfstar(z::Exo1Real)
    (z - star(z)) / 2
end

"""
    unreal(z::Exo1Real)

The unreal part of a 1-exo-real number is a 1-vector.
"""
function unreal(z::Exo1Real{T}) where T <: Real
    T[z.r]
end

function show(io::IO, z::Exo1Real)
    print(io, "[1: ")
    print(io, z.l)
    print(io, ", a: ")
    print(io, z.r)
    print(io, "]")
end

function random(::Type{Exo1Real{T}}) where T <: Real
    Exo1Real{T}(rand(T), rand(T))
end

function random(::Type{Exo1Real{BigInt}})
    Exo1Real{BigInt}(randomBigInt(), randomBigInt())
end

function random(::Type{Exo1Real{BigFloat}})
    Exo1Real{BigFloat}(randomBigFloat(), randomBigFloat())
end

function random(::Type{Exo1Real{Rational{BigInt}}})
    Exo1Real{Rational{BigInt}}(randomBigRational(), randomBigRational())
end

function convert(::Type{Exo1Real{T}}, a::Real) where T <: Real
    Exo1Real{T}(a, zero(T))
end

function convert(::Type{Exo1Real{T}}, z::Exo1Real) where T <: Real
    Exo1Real{T}(z.l, z.r)
end

function promote_rule(::Type{Exo1Real{T}}, ::Type{S}) where {T <: Real, S <: Real}
    Exo1Real{promote_type(T, S)}
end

function promote_rule(::Type{Exo1Real{T}}, ::Type{Exo1Real{S}}) where {T <: Real, S <: Real}
    Exo1Real{promote_type(T, S)}
end


text:
"""
    linearized model-predictive control policy
"""

@with_kw mutable struct LinearizedMPCOptions{T}
	altitude_update::Bool = false
	altitude_impact_threshold::T = 1.0
	altitude_verbose::Bool = false
    ip_max_time::T = 1e5     # maximum time allowed for an InteriorPoint solve
    live_plotting::Bool=false # Use the live plotting tool to debug
end

mutable struct LinearizedMPC <: Policy
	traj
	ref_traj
	im_traj
	H
	stride
	altitude
	κ
	newton
	newton_mode
	s
	q0
	N_sample
	cnt
	opts
end

function linearized_mpc_policy(traj, s, obj;
	H_mpc = traj.H,
	N_sample = 1,
	κ_mpc = traj.κ[1],
	ip_type = :interior_point,
	mode = :configurationforce,
	newton_mode = :direct,
	n_opts = NewtonOptions(
		r_tol = 3e-4,
		max_iter = 5,
		verbose = false,
		live_plotting = false),
	mpc_opts = LinearizedMPCOptions(),
	ip_opts = eval(interior_point_options(ip_type))(
				γ_reg = 0.1,
				undercut = 5.0,
				κ_tol = κ_mpc,
				r_tol = 1.0e-8,
				diff_sol = true,
				solver = :empty_solver,
				max_time = mpc_opts.ip_max_time,)
	)

	traj = deepcopy(traj)
	ref_traj = deepcopy(traj)

	im_traj = ImplicitTraj(traj, s,
		ip_type = ip_type,
		κ = κ_mpc,
		max_time = mpc_opts.ip_max_time,
		opts=ip_opts,
		mode = mode)

	stride = get_stride(s.model, traj)
	altitude = zeros(s.model.dim.c)
	if newton_mode == :direct
		newton = Newton(s, H_mpc, traj.h, traj, im_traj, obj = obj, opts = n_opts)
	elseif newton_mode == :structure
		newton = NewtonStructure(s, H_mpc, traj, obj, κ_mpc, opts = n_opts)
	else
		@error "invalid Newton solver specified"
	end

	LinearizedMPC(traj, ref_traj, im_traj, H_mpc, stride, altitude, κ_mpc, newton, newton_mode, s, copy(ref_traj.q[1]),
		N_sample, N_sample, mpc_opts)
end

function policy(p::LinearizedMPC, x, traj, t)
	# reset
	if t == 1
		p.cnt = p.N_sample
		p.q0 = copy(p.ref_traj.q[1])
	end

    if p.cnt == p.N_sample
		(p.opts.altitude_update && t > 1) && (update_altitude!(p.altitude, p.s,
									traj, t, p.N_sample,
									threshold = p.opts.altitude_impact_threshold,
									verbose = p.opts.altitude_verbose))
		# update!(p.im_traj, p.traj, p.s, p.altitude, κ = p.κ) #@@@ keep the altitude update here
		set_altitude!(p.im_traj, p.altitude) #@@@ keep the altitude update here
		newton_solve!(p.newton, p.s, p.im_traj, p.traj,
			warm_start = t > 1, q0 = copy(p.q0), q1 = copy(x))
		update!(p.im_traj, p.traj, p.s, p.altitude, κ = p.κ) #@@@ only keep the rotation stuff not the altitude update.
		p.opts.live_plotting && live_plotting(p.s.model, p.traj, traj, p.newton, p.q0, copy(x), t)

		rot_n_stride!(p.traj, p.stride)
		p.q0 .= copy(x)
		p.cnt = 0
    end

    p.cnt += 1

	if p.newton_mode == :direct
    	return p.newton.traj.u[1] / p.N_sample # rescale output
	elseif p.newton_mode == :structure
		return p.newton.u[1] / p.N_sample
	else
		@error "newton mode specified not available"
	end
end


text:
using Abaco
using Test

Abaco.loginit()

# My time zone: Saturday, October 30, 2021 3:43:35 PM GMT+02:00 DST
window = 900
metric_ts = 1635600600
metric_sn = "Civetta"

expected_triggers = 10
actual_triggers = 0

#1
abaco = abaco_init(interval=window, ages=10) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[identity] age [$ts]: scope: [$ne] $name = $value"
    @test ts == metric_ts
    @test ne == metric_sn
    @test value == 100.0
    actual_triggers += 1
end

formula(abaco, "r = x")
ingest(abaco, metric_ts, metric_sn, "x", 100.0)

#2
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[negation] age [$ts]: scope: [$ne] $name = $value"
    @test value == -100
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = -x")
ingest(abaco, metric_ts, metric_sn, "x", 100)

#3
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[add] age [$ts]: scope: [$ne] $name = $value"
    @test value == 6.2
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = x + y")
ingest(abaco, metric_ts, metric_sn, "x", 1.0)
ingest(abaco, metric_ts, metric_sn, "y", 5.2)

#4
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[mult] age [$ts]: scope: [$ne] $name = $value"
    @test value == 3.1*2.1
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = x * y")
ingest(abaco, metric_ts, metric_sn, "x", 3.1)
ingest(abaco, metric_ts, metric_sn, "y", 2.1)

#5
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[log10] age [$ts]: scope: [$ne] $name = $value"
    @test value == 2
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = log10(x)")
ingest(abaco, metric_ts, metric_sn, "x", 100.0)

#6
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[exp] age [$ts]: scope: [$ne] $name = $value"
    @test value == exp(1)
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = exp(x)")
ingest(abaco, metric_ts, metric_sn, "x", 1)

#7
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[div] age [$ts]: scope: [$ne] $name = $value"
    @test value == 1
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = div(x, y)")
ingest(abaco, metric_ts, metric_sn, "x", 5)
ingest(abaco, metric_ts, metric_sn, "y", 3.1)

#8
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[x^(y-1)] age [$ts]: scope: [$ne] $name = $value"
    @test value == 4
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = x^(y-1)")
ingest(abaco, metric_ts, metric_sn, "x", 2)
ingest(abaco, metric_ts, metric_sn, "y", 3)

#9
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[x^2 + y^3 +3*x^4] age [$ts]: scope: [$ne] $name = $value"
    @test value == 60
    actual_triggers += 1
end
metric_ts += window
formula(abaco, "r = x^2 + x^3 + 3*x^4")
ingest(abaco, metric_ts, metric_sn, "x", 2)
ingest(abaco, metric_ts, metric_sn, "y", 3)

#10
oncomplete(abaco) do ts, ne, name, value, inputs
    global actual_triggers
    @debug "[factorial] age [$ts]: scope: [$ne] $name = $value"
    @test value == factorial(6)
    actual_triggers += 1
end
metric_ts += window
# if a function requires an integer argument try to convert the arg value to Int
formula(abaco, "r = factorial(Int(x+x))")
ingest(abaco, metric_ts, metric_sn, "x", 3)

@test actual_triggers === expected_triggers