text:
abstract type AbstractRecord{E<:AbstractEntry,N,C<:AbstractClock} <:
              ResizingTools.AbstractRNArray{E,N} end
const Record = AbstractRecord

newentry(r::AbstractRecord{E}) where {E} = Base.Fix2(E, currenttime(getclock(r)))
newentry(r::AbstractRecord{E}, v) where {E} = E(v, currenttime(getclock(r)))

function Base.getindex(r::AbstractRecord{E,N,C}, I::Vararg{Int,N}) where {E,N,C}
    Base.@_propagate_inbounds_meta
    I′ = to_entryind(r, I...)
    return getentry(r, I′...)
end
function Base.setindex!(r::AbstractRecord{E,N,C}, v, I::Vararg{Int,N}) where {E,N,C}
    Base.@_propagate_inbounds_meta
    I′ = to_entryind(r, I...)
    e = getentry!(r, I′...)
    store!(e, v, getclock(r))
    return r
end

delat!(r::AbstractRecord{<:DynamicEntry}, ::Any...) = r
function delat!(r::AbstractRecord, inds...)
    t = currenttime(getclock(r))
    _delall!(r[inds...], t)
    return r
end

_delall!(e::AbstractEntry, t::Real) = (del!(e, t); e)
function _delall!(es, t::Real)
    for entry in es
        del!(entry, t)
    end
    return es
end

const AbstractScalar{T} = AbstractArray{T,0}
struct ScalarRecord{E<:AbstractEntry,C<:AbstractClock} <: AbstractRecord{E,0,C}
    c::C
    e::E
end
Record{E}(c::AbstractClock, s::AbstractScalar) where {E<:AbstractEntry} =
    ScalarRecord(c, E(s[1], c))
Record{E}(c::AbstractClock, v::Number) where {E<:AbstractEntry} = ScalarRecord(c, E(v, c))

Base.parent(r::ScalarRecord) = r.e
ArrayInterface.parent_type(::Type{R}) where {E,R<:ScalarRecord{E}} = E
ResizingTools.getsize(::ScalarRecord) = ()
to_entryind(::ScalarRecord) = ()
getentry(r::ScalarRecord) = r.e
getentry!(r::ScalarRecord) = r.e
getclock(r::ScalarRecord) = r.c

struct VectorRecord{E<:AbstractEntry,C<:AbstractClock} <: AbstractRecord{E,1,C}
    c::C
    es::Vector{E}
    im::Vector{Int}
end
function Record{E}(c::AbstractClock, v::AbstractVector) where {E<:AbstractEntry}
    es = map(v -> E(v, c), v)
    im = collect(Base.OneTo(length(v)))
    return VectorRecord(c, es, im)
end

Base.parent(r::VectorRecord) = r.es
ArrayInterface.parent_type(::Type{R}) where {E,R<:VectorRecord{E}} = Vector{E}
ResizingTools.getsize(r::VectorRecord) = (length(r.im),)
Base.@propagate_inbounds to_entryind(r::VectorRecord, I::Int) = (r.im[I],)
Base.@propagate_inbounds getentry(r::VectorRecord, i′::Int) = r.es[i′] # i′ is parent index
Base.@propagate_inbounds getentry!(r::VectorRecord, i′::Int) = r.es[i′] # i′ is parent index
getclock(r::VectorRecord) = r.c

function ResizingTools.resize_buffer!(A::VectorRecord{E}, I) where {E<:AbstractEntry}
    len = length(A)
    nl = to_dims(I)
    pl = length(parent(A))
    diff = nl - len
    if diff > 0
        append!(A.im, pl+1:pl+diff)
        append!(A.es, map(_ -> E(), 1:diff))
    elseif diff < 0
        delat!(A, _del_ind(len, I))
        resize!(A.im, (I,))
    end
    return A
end
function ResizingTools.resize_buffer_dim!(A::VectorRecord, d::Int, I)
    @boundscheck d == 1 || throw(BoundsError())
    return ResizingTools.resize_buffer!(A, I)
end

_del_ind(len::Int, ind::Integer) = Int(ind):len
_del_ind(len::Int, ind::Base.LogicalIndex) = Base.OneTo(len)[not(ind)]

Base.sizehint!(r::VectorRecord, sz::Integer) = sizehint!(r.es, sz)
function Base.push!(r::VectorRecord, v)
    push!(r.es, newentry(r, v))
    push!(r.im, lastindex(r.es))
    return r
end
Base.append!(r::VectorRecord{E}, v::V) where {V,E<:AbstractEntry{V}} = push!(r, v)
function Base.append!(r::VectorRecord, vs)
    len = length(r.es)
    append!(r.es, map(newentry(r), vs))
    append!(r.im, len+1:len+length(vs))
    return r
end
function Base.insert!(r::VectorRecord, i::Integer, v)
    push!(r.es, newentry(r, v))
    insert!(r.im, i, lastindex(r.es))
end
function Base.deleteat!(r::VectorRecord, inds)
    delat!(r, inds)
    deleteat!(r.im, inds)
    return r
end
function Base.resize!(r::VectorRecord, nl::Integer)
    len = length(r)
    if nl > len
        E = eltype(r)
        append!(r.es, map(_ -> E(), len+1:nl))
        append!(r.im, len+1:nl)
    elseif nl != len
        nl < 0 && throw(ArgumentError("new length must be ≥ 0"))
        delat!(r, nl+1:len)
        resize!(r.im, len - nl)
    end
    return r
end

struct DOKRecord{E<:AbstractEntry,N,C<:AbstractClock} <: AbstractRecord{E,N,C}
    c::C
    dok::DOKSparseArray{E,N}
    sz::Size{N}
    im::MCIndices{N}
end
function Record{E}(c::AbstractClock, A::AbstractArray) where {E<:AbstractEntry}
    sz = Size(A)
    dok = Dict{NTuple{ndims(A),Int},return_type(E, A, c)}()
    sizehint!(dok, prod(size(A)))
    for (i, ind) in enumerate(MCIndices(A))
        dok[ind] = E(A[i], c)
    end
    dok_array = DOKSparseArray(dok, Size(A))
    im = MCIndices(A)
    return DOKRecord(c, dok_array, sz, im)
end
Base.parent(r::DOKRecord) = r.dok
ArrayInterface.parent_type(::Type{R}) where {E,N,R<:DOKRecord{E,N}} = DOKSparseArray{E,N}
ResizingTools.getsize(r::DOKRecord) = r.sz
ResizingTools.size_type(::Type{T}) where {T<:DOKRecord} = Size{ndims(T)}

to_entryind(r::DOKRecord{E,N}, I::Vararg{Int,N}) where {E,N} =
    (Base.@_propagate_inbounds_meta; r.im[I...])
getentry(r::DOKRecord{E,N}, I::Vararg{Int,N}) where {E,N} = get(parent(r), I, E())
getentry!(r::DOKRecord{E,N}, I::Vararg{Int,N}) where {E,N} = get!(parent(r), I, E())
getclock(r::DOKRecord) = r.c

function ResizingTools.pre_resize!(r::DOKRecord{E,N}, inds::NTuple{N,Any}) where {E,N}
    sz = size(r)
    nsz = to_dims(inds)
    psz = size(parent(r))
    for ind in CartesianIndices(r)
        # del entries that are no longer in the new size
        _checkindices(inds, ind.I) || delat!(r, ind)
    end
    for d in 1:N # change im
        diff = nsz[d] - sz[d]
        if diff > 0
            append!(getdim(r.im, d), psz[d]+1:psz[d]+diff)
        elseif diff < 0
            resize!(getdim(r.im, d), (inds[d],))
        end
    end
    return r
end
function ResizingTools.pre_resize!(r::DOKRecord, d::Int, I)
    n = size(r, d)
    nn = to_dims(I)
    pn = size(parent(r), d)
    for ind in CartesianIndices(r)
        # del entries that are no longer in the new size
        _checkindex(I, ind[d]) || delat!(r, ind)
    end
    diff = nn - n
    if diff > 0
        append!(getdim(r.im, d), pn+1:pn+diff)
    elseif diff < 0
        resize!(getdim(r.im, d), (I,))
    end
    return r
end

_checkindices(inds::Tuple, I::Dims) = all(map(_checkindex, inds, I))
_checkindex(ind, i::Int) = i in ind
_checkindex(n::Integer, i::Int) = 1 <= i <= n

ResizingTools.to_parentinds(r::DOKRecord, Is::Tuple) =
    map(_to_parentind, size(r), size(parent(r)), Is)
ResizingTools.to_parentinds(r::DOKRecord, d::Integer, I) =
    Int(d), _to_parentind(size(r, d), size(parent(r), d), I)

_to_parentind(::Int, pn::Int, ::Any) = pn
_to_parentind(n::Int, pn::Int, nn::Int) = (d = nn - n; d > 0 ? pn + d : pn)

Base.show(io::IO, ::MIME"text/plain", r::AbstractRecord) = summary(io, r)


text:
"""
`deleteitems!(vector, items)` deletes a group of `items` (test it using `in`) from `vector`.
"""
function deleteitems!(vector::Vector, items)
  i = 1
  j = 0
  @inbounds while j < length(vector)
    j += 1
    value = vector[j]
    if !(value in items)
      if i != j
        vector[i] = value
      end
      i += 1
    end
  end
  resize!(vector, i-1)
end

_get_data(str::ASCIIString) = str.data
_get_data(str::UTF8String) = collect(str)

"""
`get_n_words{T <: Union{ASCIIString, UTF8String}}(line::T, n::Int)` returns a `Vector{T}` with the first `n` (possibles) words/fields (delimited by space, tab or newline).
If there is more than `n` words, the last word returned contains the finals words and the delimiters.
The length of the returned vector is `n` or less (if the number of words is less than `n`).
The last newline character is always removed.
This is used for parsing the Stockholm format.
```
julia> get_n_words("#=GR O31698/18-71 SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH\n", 3)
3-element Array{ASCIIString,1}:
 "#=GR"
 "O31698/18-71"
 "SS    CCCHHHHHHHHHHHHHHHEEEEEEEEEEEEEEEEHHH"
```
"""
function get_n_words{T <: ByteString}(line::T, n::Int)
  linedata = _get_data(line)
  CharType = eltype(linedata)
  N = length(linedata)
  coords = Array(T, n)
  spaces = CharType[' ', '\n', '\t']
  j = 0
  start = 0
  ended = false
  if N > 1
    block = !(linedata[1] in spaces)
    if block
      start = 1
    end
    for i in 2:N-1
      if linedata[i] in spaces
        if block
          j += 1
          if j == n
            last = linedata[N] == CharType('\n') ? N-1 : N
            @inbounds coords[j] = UTF8String(getindex(linedata,start:last))
            ended = true
            break
          else
            @inbounds coords[j] = UTF8String(getindex(linedata,start:i-1))
          end
        end
        block = false
      else
        if ! block
          block = true
          start = i
        end
      end
    end
    if !ended
      j += 1
      last = linedata[N] == CharType('\n') ? N-1 : N
      @inbounds coords[j] = UTF8String(getindex(linedata,start:last))
    end
  elseif N==1
    if linedata[1] in spaces
      j=0
    else
      j=1
      coords[1] = line
    end
  else
    j=0
  end
  if j != n
    resize!(coords, j)
  end
  coords
end

"""
`hascoordinates(id)` returns `true` if `id` has the format: **uniprot/start-end** i.e. O83071/192-246
"""
function hascoordinates(id)
  ismatch(r"^\w+/\d+-\d+$", id)
end

"""
Selects the first element of the vector.
This is useful for unpacking one element vectors.
Throws a warning if there are more elements.
The `element_name` for the warning is *element* by default, but the name can be changed.
"""
function select_element{T}(vector::Array{T,1}, element_name::ASCIIString="element")
  len = length(vector)
  if len == 0
    throw(ErrorException(string("There is not ", element_name)))
  elseif len != 1
    warn(string("There are more than one (", len, ") ", element_name, " using the first."))
  end
  @inbounds return(vector[1])
end

"""
Returns a vector with the `part` ("upper" or "lower") of the square matrix `mat`.
The `diagonal` is not included by default.
"""
function matrix2list{T}(mat::AbstractMatrix{T}; part="upper", diagonal::Bool=false)
  nrow, ncol = size(mat)
  if nrow != ncol
    throw(ErrorException("Should be a square matrix"))
  end
  if diagonal
    d = 0
    N = div((ncol * ncol) + ncol, 2)
  else
    d = 1
    N = div((ncol * ncol) - ncol, 2)
  end
  list = Array(T, N)
  k = 1
  if part=="upper"
    for i in 1:(ncol-d)
      for j in (i+d):ncol
        list[k] = mat[i, j]
        k += 1
      end
    end
  elseif part=="lower"
    for j in 1:(ncol-d)
      for i in (j+d):ncol
        list[k] = mat[i, j]
        k += 1
      end
    end
  else
    throw(ErrorException("part should be \"upper\" or \"lower\""))
  end
  list
end

"""
Returns a square symmetric matrix from the vector `vec`.
`side` is the number of rows/columns.
The `diagonal` is not included by default, set to `true` if the diagonal elements are in the list.
"""
function list2matrix{T}(vec::AbstractVector{T}, side::Int; diagonal::Bool=false)
  d = diagonal ? 0 : 1
  mat = zeros(T, side, side)
  k = 1
  for i in 1:(side-d)
    for j in (i+d):side
      value = vec[k]
      mat[i, j] = value
      mat[j, i] = value
      k += 1
    end
  end
  mat
end

ext: jl


text:
include("model_tests.jl")
include("states_tests.jl")
include("actions_tests.jl")
include("mdp_tests.jl")

max_stars_repo_name: ancorso/MAST

max_issues_repo_name: ancorso/MAST


text:
using LaserTag
using Plots
using POMDPToolbox
using Reel
using ProgressMeter
using PmapProgressMeter
using ParticleFilters
using ContinuousPOMDPTreeSearchExperiments
using Plots
using QMDP
using JLD
# using DESPOT
using BasicPOMCP
using ARDESPOT

file_contents = readstring(@__FILE__())

@everywhere begin
# begin
    using POMDPs
    using POMDPToolbox
    using ContinuousPOMDPTreeSearchExperiments
    using DiscreteValueIteration
    using ParticleFilters
    using POMCPOW
    using LaserTag
    using QMDP
    import DESPOT
    using BasicPOMCP
    using ARDESPOT

    N = 1000
    K = 5000
    T_max = 10.0
    P = typeof(gen_lasertag(rng=MersenneTwister(18)))

    solvers = Dict{String, Union{Policy, Solver}}(

        "qmdp" => QMDPSolver(max_iterations=1000),

        "ardespot" => begin
            DESPOTSolver(lambda=0.0,
                         max_trials=1_000_000,
                         T_max=T_max,
                         K=K,
                         bounds=LaserBounds{P}(),
                         default_action=nogap_tag,
                         rng=MersenneTwister(4))
        end,

        "despot" => begin
            DESPOT.DESPOTSolver{LTState,
                                Int,
                                DMeas,
                                LaserBounds{P},
                                DESPOT.MersenneStreamArray}(
                bounds=LaserBounds{P}(),
                time_per_move=T_max,
                max_trials=1_000_000,
                n_particles=K,
                random_streams=DESPOT.MersenneStreamArray(MersenneTwister(4)),
                next_state=LTState([-1,-1],[-1,-1], false)
            )
        end
    )
end

@show K
@show N
@show T_max

rdict = Dict{String, Any}()
for (k,sol) in solvers
    prog = Progress(N, desc="Simulating...")
    @show k
    rewards = pmap(prog, 1:N) do i
    # rewards = map(1:N) do i
        pomdp = gen_lasertag(rng=MersenneTwister(i+600_000))
        if isa(sol,Solver)
            p = solve(deepcopy(sol), pomdp)
        else
            p = sol
        end
        hr = HistoryRecorder(max_steps=100, rng=MersenneTwister(i))
        up_rng = MersenneTwister(i+100_000)
        up = ObsAdaptiveParticleFilter(deepcopy(pomdp), LowVarianceResampler(100_000), 0.05, up_rng)
        hist = simulate(hr, pomdp, p, up)
        discounted_reward(hist)
    end
    @show mean(rewards)
    @show std(rewards)/sqrt(N)
    rdict[k] = rewards
end

filename = Pkg.dir("ContinuousPOMDPTreeSearchExperiments", "data", "despot_comparison_run_$(Dates.format(now(), "E_d_u_HH_MM")).jld")
println("saving to $filename...")
@save(filename, rdict, file_contents)
println("done.")


text:
using Catlab.CategoricalAlgebra, Catlab.Graphs, Catlab.Present, Catlab.Graphics, Catlab.Theories
using Catlab.CategoricalAlgebra.FinCats: FinCatGraphEq

using Distributions: Exponential, Geometric, cdf
using Random: randsubseq
using StatsBase: sample

using Plots
using ProgressBars

include("utils.jl")

# Lotka-Volterra
@present ThLV(FreeSchema) begin
    (Rabbit,Fox)::Ob
end

@acset_type LV(ThLV)

# birth of prey
L_birth = @acset LV begin Rabbit=1 end
I_birth = LV()
R_birth = @acset LV begin Rabbit=2 end

l_birth = ACSetTransformation(I_birth, L_birth)
r_birth = ACSetTransformation(I_birth, R_birth)

# predation
L_predation = @acset LV begin Rabbit=1; Fox=1 end
I_predation = @acset LV begin Fox=1 end
R_predation = @acset LV begin Fox=2 end

l_predation = ACSetTransformation(I_predation, L_predation; Fox=[1])
r_predation = ACSetTransformation(I_predation, R_predation; Fox=[1])

# death of predator
L_death = @acset LV begin Fox=1 end
I_death = LV()
R_death = LV()

l_death = ACSetTransformation(I_death, L_death)
r_death = ACSetTransformation(I_death, R_death)

# rules
rule_birth = Rule(l_birth, r_birth)
rule_predation = Rule(l_predation, r_predation)
rule_death = Rule(l_death, r_death)

# parameters
X = 50
Y = 10

Δt = 0.1
tmax = 365*2
steps = Int(tmax/Δt)

b = 1.0
mu = 0.01
β = 0.0025

# run model
state = @acset LV begin Rabbit=X; Fox=Y end

out = fill(-1, (steps, 2))

for t = ProgressBar(1:steps)

    # possible events
    birth = homomorphisms(L_birth, state, monic = true)
    predation = homomorphisms(L_predation, state, monic = true)
    death = homomorphisms(L_death, state, monic = true)

    # sample firings
    birth = sample_matches(birth, b, Δt)
    predation = sample_matches(predation, β, Δt)
    death = sample_matches(death, mu, Δt)

    # apply events
    events = MatchedRule[
        [MatchedRule(rule_birth, m) for m in birth]...,
        [MatchedRule(rule_predation, m) for m in predation]...,
        [MatchedRule(rule_death, m) for m in death]...
    ]
    while length(events) > 0
        ev = pop!(events)
        # apply event
        _, kg, _, kh = rewrite_match_maps(ev.rule.L, ev.rule.R, ev.match)
        state = codom(kh)
        # update the remaining matches post rewrite
        for j in 1:length(events)
            events[j].match = postcompose_partial(kg, kh, events[j].match)
        end
        events = filter((e) -> !isnothing(e.match), events)
    end
    # write output
    out[t, :] = [nparts(state, x) for x in [:Rabbit, :Fox]]
end

plot(
    (1:steps) * Δt,
    out,
    label=["Rabbit" "Fox"],
    xlabel="Time",
    ylabel="Number"
)


text:
struct Graph
    edges::Vector{Vector{Tuple{Int64, Float64}}}
    Graph(n) = new([Tuple{Int64, Float64}[] for _ in 1 : n])
end

function add_edge(G::Graph, v::Tuple{Int64, Int64, Float64})
    i, j, c = v
    push!(G.edges[i], (j, c))
    push!(G.edges[j], (i, c))
end

Base.getindex(G::Graph, i::Int) = G.edges[i]

mutable struct HDBSCANCluster
    parent::Int64
    children::Vector{Int64}
    points::Vector{Int64}
    λp::Vector{Float64}
    stability::Float64
    children_stability::Float64
    function HDBSCANCluster(noise::Bool, points::Vector{Int64})
        if noise
            return new(0, [], [], [], -1, -1)
        else
            return new(0, [], points, [], 0, 0)
        end
    end
end

Base.length(c::HDBSCANCluster) = size(c.points, 1)
join(c1::HDBSCANCluster, c2::HDBSCANCluster, id) = HDBSCANCluster(nothing, vcat(c1.points, c2.points), id, 0)
isnoise(c::HDBSCANCluster) = c.stability == -1
hasstability(c::HDBSCANCluster) = c.stability != 0
function compute_stability(c::HDBSCANCluster, λbirth)
    c.stability += sum(c.λp.-λbirth)
end

"""
    HDBSCAN(ε, minpts)
Density-Based Clustering Based on Hierarchical Density Estimates.
This algorithm performs clustering as follows.
1. generate a minimum spanning tree
2. build a HDBSCAN hierarchy
3. extract the target cluster
4. generate the list of cluster assignment for each point
The detail is so complex it is difficult to explain the detail in here. But, if you want to know more about this algorithm, you should read [this docs](https://hdbscan.readthedocs.io/en/latest/how_hdbscan_works.html).

# Parameters:
- `k`: we will define "core distance of point A" as the distance between point A and the `k` th neighbor point of point A.
- `min_cluster_size`: minimum number of points in the cluster

# Example
```jldoctest
julia> model = HDBSCAN(3, 5)
HDBSCAN(3, 5, #undef)

julia> fit!(model, data); #What is returned at this time is the minimum spanning tree of the data generated during the clustering process.

julia> model.labels |> size
(100,)
```
"""
mutable struct HDBSCAN
    k::Int64
    min_cluster_size::Int64
    labels::Union{Vector{Int64}, Nothing}
    function HDBSCAN(k::Int64, min_cluster_size::Int64)
        if min_cluster_size < 1
            throw(DomainError(min_cluster_size, "The `min_cluster_size` must be greater than or equal to 1"))
        end
        return new(k, min_cluster_size)
    end
end

function fit!(model::HDBSCAN, x; gen_mst::Bool=true, mst=nothing)
    n = size(x, 1)
    if gen_mst
        #calculate core distances for each point
        core_dists = core_dist(x, model.k)
        #calculate mutual reachability distance between any two points
        mrd = mutual_reachability_distance(core_dists, x)
        #compute a minimum spanning tree by prim method
        mst = prim(mrd, n)
    elseif mst == nothing
        throw(ArgumentError("if you set `gen_mst` to false, you must pass a minimum spanning tree as `mst`"))
    end
    #build a HDBSCAN hierarchy
    hierarchy = build_hierarchy(mst, model.min_cluster_size)
    #extract the target cluster
    extract_cluster!(hierarchy)
    #generate the list of cluster assignment for each point
    result = fill(-1, n)
    for (i, j) in enumerate(hierarchy[2n-1].children)
        c = hierarchy[j]
        for k in c.points
            result[k] = i
        end
    end
    model.labels = result
    if gen_mst
        return mst
    end
end

function core_dist(points, k)
    core_dists = Array{Float64}(undef, size(points, 1))
    for i in 1 : size(points, 1)
        p = points[i:i, :]
        dists = vec(sum((@. (points - p)^2), dims=2))
        sort!(dists)
        core_dists[i] = dists[k]
    end
    return core_dists
end

function mutual_reachability_distance(core_dists, points)
    n = size(points, 1)
    graph = Graph(div(n * (n-1), 2))
    for i in 1 : n-1
        for j in i+1 : n
            c = max(core_dists[i], core_dists[j], sum((points[i, :]-points[j, :]).^2))
            add_edge(graph, (i, j, c))
        end
    end
    return graph
end

function prim(graph, n)
    minimum_spanning_tree = Array{Tuple{Float64, Int64, Int64}}(undef, n-1)

    marked = falses(n)
    marked_cnt = 1
    marked[1] = true

    h = []

    for (i, c) in graph[1]
        heappush!(h, (c, 1, i))
    end

    while marked_cnt < n
        c, i, j = popfirst!(h)

        marked[j] == true && continue
        minimum_spanning_tree[marked_cnt] = (c, i, j)
        marked[j] = true
        marked_cnt += 1

        for (k, c) in graph[j]
            marked[k] == true && continue
            heappush!(h, (c, j, k))
        end
    end
    return minimum_spanning_tree
end

function build_hierarchy(mst, min_size)
    n = length(mst) + 1
    cost = 0
    uf = UnionFind(n)
    Hierarchy = Array{HDBSCANCluster}(undef, 2n-1)
    if min_size == 1
        for i in 1 : n
            Hierarchy[i] = HDBSCANCluster(false, [i])
        end
    else
        for i in 1 : n
            Hierarchy[i] = HDBSCANCluster(true, Int64[])
        end
    end
    sort!(mst)

    for i in 1 : n-1
        c, j, k = mst[i]
        cost += c
        λ = 1 / cost
        #child clusters
        c1 = group(uf, j)
        c2 = group(uf, k)
        #reference to the parent cluster
        Hierarchy[c1].parent = Hierarchy[c2].parent = n+i
        nc1, nc2 = isnoise(Hierarchy[c1]), isnoise(Hierarchy[c2])
        if !(nc1 || nc2)
            #compute stability
            compute_stability(Hierarchy[c1], λ)
            compute_stability(Hierarchy[c2], λ)
            #unite cluster
            unite!(uf, j, k)
            #create parent cluster
            points = members(uf, group(uf, j))
            Hierarchy[n+i] = HDBSCANCluster(false, points)
        elseif !(nc1 && nc2)
            if nc2 == true
                (c1, c2) = (c2, c1)
            end
            #record the lambda value
            append!(Hierarchy[c2].λp, fill(λ, length(Hierarchy[c1])))
            #unite cluster
            unite!(uf, j, k)
            #create parent cluster
            points = members(uf, group(uf, j))
            Hierarchy[n+i] = HDBSCANCluster(false, points)
        else
            #unite the noise cluster
            unite!(uf, j, k)
            #create parent cluster
            points = members(uf, group(uf, j))
            if length(points) < min_size
                Hierarchy[n+i] = HDBSCANCluster(true, Int64[])
            else
                Hierarchy[n+i] = HDBSCANCluster(false, points)
            end
        end
    end
    return Hierarchy
end

function extract_cluster!(hierarchy::Vector{HDBSCANCluster})
    for i in 1 : length(hierarchy)-1
        if isnoise(hierarchy[i])
            c = hierarchy[i]
            push!(hierarchy[c.parent].children, i)
            hierarchy[c.parent].children_stability += c.stability
        else
            c = hierarchy[i]
            if c.stability > c.children_stability
                push!(hierarchy[c.parent].children, i)
                hierarchy[c.parent].children_stability += c.stability
            else
                append!(hierarchy[c.parent].children, c.children)
                hierarchy[c.parent].children_stability += c.children_stability
            end
        end
    end
end


text:
abstract type MassComponent{D,K} end

"""
    get_kernel(mass_component::MassComponent{MD,K}) -> K

Return the smoothing kernel of the mass component.
"""
function get_kernel end

"""
    get_particle_mass(mass_component::MassComponent) -> Float

Return the mass of the particles in the mass component.
"""
function get_particle_mass end

"""
    get_mass_densities(mass_component::MassComponent) -> Vector{Float}

Return the mass densities of the particles in the mass component.
"""
function get_mass_densities end

"""
    get_kernel_widths(mass_component::MassComponent) -> Vector{Float}

Return the kernel widths for the particles in the mass component.
"""
function get_kernel_widths end

"""
    get_mass_densities(
        mass_component::MassComponent,
        (i, j)::Tuple{Unsigned,Unsigned},
    ) -> (Float, Float)

Return the mass densities of the particles with the given indices.
"""
function get_mass_densities(
    mass_component::MassComponent,
    (i, j)::Tuple{Unsigned,Unsigned},
)
    mass_densities = get_mass_densities(mass_component)
    mass_densities[i], mass_densities[j]
end

"""
    get_mass_density(mass_component::MassComponent, i::Unsigned) -> Float

Return the mass density of the particle with the given index.
"""
function get_mass_density(mass_component::MassComponent, i::Unsigned)
    mass_densities = get_mass_densities(mass_component)
    mass_densities[i]
end

"""
    get_kernel_widths(
        mass_component::MassComponent,
        (i, j)::Tuple{Unsigned,Unsigned},
    ) -> (Float, Float)

Return the kernel widths of the particles with the given indices.
"""
function get_kernel_widths(
    mass_component::MassComponent,
    (i, j)::Tuple{Unsigned,Unsigned},
)
    kernel_widths = get_kernel_widths(mass_component)
    kernel_widths[i], kernel_widths[j]
end

"""
    get_kernel_width(mass_component::MassComponent, i::Unsigned) -> Float

Return the kernel width of the particles with the given index.
"""
function get_kernel_width(mass_component::MassComponent, i::Unsigned)
    kernel_widths = get_kernel_widths(mass_component)
    kernel_widths[i]
end

"""
    initderivatives(mass_component::MassComponent{MD}) -> MD

Creates initial derivatives for the mass component.
"""
initderivatives(::MassComponent{Nothing}) = nothing

"""
    initderivatives!(
        derivatives::MD,
        ::MassComponent{MD},
    )

Resets the given derivatives for the mass component.
"""
function initderivatives!(::Nothing, ::MassComponent{Nothing}) end

"""
    updatederivatives!(
        mass_derivatives::MD,
        (i, j)::Tuple{Unsigned,Unsigned},
        vᵣ::Number,
        m_ddr_W_r_hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
        mass_component::MassComponent{MD},
    )

Update mass component derivatives for particles `i` and `j`.
"""
function updatederivatives!(
    ::Nothing,
    ::Tuple{Unsigned,Unsigned},
    ::Number,
    ::Union{Number,Tuple{Number,Number}},
    ::MassComponent{Nothing},
) end

"""
    updatederivative!(
        mass_derivatives::MD,
        i::Unsigned,
        vᵣ::Number,
        m_ddr_W_r_hᵢ_hⱼ::Union{Number,Tuple{Number,Number}},
        mass_component::MassComponent{MD},
    )

Update mass component derivative for particle `i`.
"""
function updatederivative!(
    ::Nothing,
    ::Unsigned,
    ::Number,
    ::Union{Number,Tuple{Number,Number}},
    ::Number,
    ::MassComponent{Nothing},
) end

"""
    evolvevariables!(
        mass_component::MassComponent{MD},
        mass_derivatives::MD,
        Δt::Number,
    )

Evolve mass component with the given derivatives and time step.
"""
function evolvevariables!(::MassComponent{Nothing}, ::Nothing, ::Number) end

"""
    normalize_distance(
        ::MassComponent,
        r::Number,
        (hᵢ, hⱼ)::Tuple{Number,Number},
    ) -> (Float, Float)

Normalize the distance by the two kernel widths.
"""
normalize_distance(::MassComponent, r::Number, (hᵢ, hⱼ)::Tuple{Number,Number}) =
    r / hᵢ, r / hⱼ

"""
    normalize_distance(::MassComponent, r::Number, h::Number) -> Float

Normalize the distance by the kernel width.
"""
normalize_distance(::MassComponent, r::Number, h::Number) = r / h

"""
    estimate_initial_kernel_width(
        ρ::Number,
        m::Number,
        η::Number,
        ::Dim{N},
    ) -> Vector{Float}

Estimate initial widths of the SPH smoothing kernels.
"""
function estimate_initial_kernel_width(
    ρ::Number,
    m::Number,
    η::Number,
    ::Dim{N},
) where {N}
    η * (m / ρ)^(1 / N)
end

"""
    updatemasses!(
        mass_component::MassComponent,
        positions::AbstractVector{SVectorF{N}},
    )

Update the state of the mass component according to the given particle
positions.
"""
function updatemasses!(::MassComponent, ::AbstractVector{SVectorF{N}}) where {N} end

function compute_m_ddr_W_r_h(
    ::Union{Unsigned,Tuple{Unsigned,Unsigned}},
    ddr_W_r_h::Number,
    mass_component::MassComponent,
)
    get_particle_mass(mass_component) * ddr_W_r_h
end

function compute_m_ddr_W_r_h(
    ::Tuple{Unsigned,Unsigned},
    (ddr_W_r_hᵢ, ddr_W_r_hⱼ)::Tuple{Number,Number},
    mass_component::MassComponent,
)
    m = get_particle_mass(mass_component)
    m * ddr_W_r_hᵢ, m * ddr_W_r_hⱼ
end

include("mass/standard.jl")
include("mass/adaptive.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
include("../tcm_code/tcm_estimateForest.jl")

function tcm_estimateForest_L1(N, assortments, choice_probs, time_limit_overall, depth_limit, initial_left, initial_right, initial_product, initial_isLeaf)

	M = size(assortments,1)

	dim_V = sum(assortments)

	m_dual = Model(solver = GurobiSolver(OutputFlag = 0))
	@variable(m_dual, -1 <= p[1:M, 1:N] <= +1)
	@variable(m_dual, q)

	@objective(m_dual, Max, sum( p[m,i]* choice_probs[m,i] for m in 1:M, i in 1:N ) + q)

	forest_left = deepcopy(initial_left)
    forest_right = deepcopy(initial_right)
    forest_product = deepcopy(initial_product)
    forest_isLeaf = deepcopy(initial_isLeaf)

	K = length(forest_left)

	v_by_class = Array{Int64,2}[];

	for k in 1:K
		v_one_class = tcm_treeToA( forest_left[k], forest_right[k], forest_product[k], forest_isLeaf[k], assortments)
		push!(v_by_class, v_one_class)
		@constraint(m_dual, sum( v_one_class[m,i] * p[m,i] for m in 1:M, i in 1:N) + q <= 0)
	end

	tol_RC = -1e-4

	start_time = time()

	solve(m_dual)

	dual_objval = getobjectivevalue(m_dual)

	iter = 0;

    println("CG iter ", iter, " -- time elapsed = ", 0.0, " objective = ",dual_objval)

	while (time() - start_time < time_limit_overall)
        iter += 1

        p_val = getvalue(p)
        q_val = getvalue(q)

        single_left, single_right, single_product, single_isLeaf, sub_objval = tcm_solveTreeSubproblem_L1(N, p_val, assortments, depth_limit)

        # @show sub_objval
        v_one_class = tcm_treeToA( single_left, single_right, single_product, single_isLeaf, assortments)
        # @show dot(p_val, v_one_class)

        # @show p_val

        # error("Stop here")

        RC = -sub_objval - q_val

        @show RC
        @show sub_objval
        @show q_val
        @show sub_objval + q_val

        println("CG iter ", iter, " -- time elapsed = ", time() - start_time, " objective = ", dual_objval, " -- RC = ", RC)

        if (RC >= tol_RC)
        	# Stop here -- we are at optimality
        	println("\t RC = ", RC, " -- terminating ... ")
        	break;
        else
        	# Add the tree to the dual problem, and re-solve the master.
        	push!(forest_left, single_left)
        	push!(forest_right, single_right)
        	push!(forest_product, single_product)
        	push!(forest_isLeaf, single_isLeaf)

        	v_one_class = tcm_treeToA( single_left, single_right, single_product, single_isLeaf, assortments)
        	push!(v_by_class, v_one_class)
        	@constraint(m_dual, sum( v_one_class[m,i] * p[m,i] for m in 1:M, i in 1:N) + q <= 0)

        	solve(m_dual)
        	dual_objval = getobjectivevalue(m_dual)
        end

        # break;
	end

	end_time = time()
	elapsed_time = end_time - start_time

	# Have terminated; solve primal to get what the lambda's are.
	K = length(forest_left)

	m_primal = Model(solver = GurobiSolver())
	@variable(m_primal, lambda_var[1:K] >= 0)
	@variable(m_primal, eps_plus[1:M,1:N] >= 0)
	@variable(m_primal, eps_minus[1:M,1:N] >= 0)

	for m in 1:M
		for i in 1:N
			@constraint(m_primal, sum( v_by_class[k][m,i] * lambda_var[k] for k in 1:K) + eps_plus[m,i] - eps_minus[m,i] == choice_probs[m,i])
		end
	end
	@constraint(m_primal, sum(lambda_var[k] for k in 1:K) == 1)

	@objective(m_primal, Min, sum( eps_minus[m,i] for m in 1:M, i in 1:N) + sum(eps_plus[m,i] for m in 1:M, i in 1:N))

	solve(m_primal)

	lambda_val = getvalue(lambda_var)

	primal_objval = getobjectivevalue(m_primal)

	return lambda_val, forest_left, forest_right, forest_product, forest_isLeaf, primal_objval, elapsed_time

end

function tcm_treeToA(tree_left, tree_right, tree_product, tree_isLeaf, assortments)
	M = size(assortments, 1)
	v_one_class = zeros(Float64,size(assortments))
    for m = 1 : M
        cn = 1;
        while (!tree_isLeaf[cn])
        	if (assortments[m, tree_product[cn]] > 0)
        		cn = tree_left[cn]
        	else
        		cn = tree_right[cn]
        	end
        end
        v_one_class[m, tree_product[cn]] = 1.0
    end

    return v_one_class
end

function tcm_solveTreeSubproblem_L1(N, u, assortments, depth_limit)
    # @objective(m_sub, Max, sum(u[m,p] * a[m,p] for m in 1:M, p in 1:N))
    # solve(m_sub)

    # x_val = getvalue(x)
    # objval = getobjectivevalue(m_sub)

    # single_ranking = sum(x_val, 2)
 #    single_ranking = single_ranking[:]
    # single_ordering = sortperm(single_ranking, rev = true)

    # return single_ordering, objval

    M = size(assortments,1)

    single_left = Int64[0]
    single_right = Int64[0]
    single_product = Int64[N]
    single_isLeaf = Bool[true]
    single_depth = Int64[1]
    single_include_set = Array{Int64}[ [N] ];
    single_assortment_binary = Array{Int64}[ ones(M) ]
    # single_reduced_cost = Float64[ sum(u[:,N]) ] #Float64[0.0];
    single_reduced_cost = Float64[ -Inf ] #Float64[0.0];
    single_isClosed = Bool[false]

    iter = 0 ;

    while (true)
        current_reduced_cost = sum( single_reduced_cost[single_isLeaf])
        println("\t\t Subproblem Iter : ", iter, " -- reduced cost: ", current_reduced_cost, " -- num leaves: ", sum(single_isLeaf))
        candidate_leaves = find(single_isLeaf .& (single_depth .< depth_limit) .&  (.!single_isClosed) )

        if (isempty(candidate_leaves))
            println("\t\t\t No more eligible leaves to split; exiting procedure...")
            break;
        end

        num_candL = length(candidate_leaves)
        candidate_reduced_costs = -Inf* ones(num_candL)
        best_split_p_by_leaf = zeros(Int64, num_candL)
        best_left_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_right_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_left_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_right_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_left_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)
        best_right_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)

        for ell_ind in 1:num_candL

            leaf = candidate_leaves[ell_ind]

            assortment_binary = copy(single_assortment_binary[leaf])

            candidate_split_products = setdiff(collect(1:(N-1)), single_include_set[leaf])

            for split_p in candidate_split_products

                candidate_left_leaf_products = union(single_include_set[leaf], split_p)
                best_left_reduced_cost = -Inf
                best_left_p = 0;
                best_left_assortment_binary = zeros(M)
                for left_p in candidate_left_leaf_products
                    left_assortment_binary = assortment_binary .* assortments[:, split_p][:]
                    temp = dot( left_assortment_binary, u[:, left_p])
                    if (temp > best_left_reduced_cost)
                        best_left_reduced_cost = temp
                        best_left_p = left_p
                        best_left_assortment_binary = copy(left_assortment_binary)
                    end
                end

                candidate_right_leaf_products = copy(single_include_set[leaf])
                best_right_reduced_cost = -Inf
                best_right_p = 0;
                best_right_assortment_binary = zeros(M)
                for right_p in candidate_right_leaf_products
                    right_assortment_binary = assortment_binary .* (1 - assortments[:,split_p][:])
                    temp = dot(right_assortment_binary, u[:, right_p])
                    if (temp > best_right_reduced_cost)
                        best_right_reduced_cost = temp
                        best_right_p = right_p
                        best_right_assortment_binary = copy(right_assortment_binary)
                    end
                end

                new_reduced_cost = best_left_reduced_cost + best_right_reduced_cost

                if (new_reduced_cost > candidate_reduced_costs[ell_ind] )
                    candidate_reduced_costs[ell_ind] = new_reduced_cost
                    best_split_p_by_leaf[ell_ind] = split_p
                    best_left_leaf_p_by_leaf[ell_ind] = best_left_p
                    best_right_leaf_p_by_leaf[ell_ind] = best_right_p
                    best_left_reduced_cost_by_leaf[ell_ind] = best_left_reduced_cost
                    best_right_reduced_cost_by_leaf[ell_ind] = best_right_reduced_cost
                    best_left_assortment_binary_by_leaf[ell_ind] = best_left_assortment_binary # VVM: should there be a copy here?
                    best_right_assortment_binary_by_leaf[ell_ind] = best_right_assortment_binary # VVM: should there be a copy here?
                end
            end
        end

        incumbent_reduced_costs = single_reduced_cost[candidate_leaves]

        tol_RC = 1e-7

        for ell_ind in 1:num_candL
            if (incumbent_reduced_costs[ell_ind] + tol_RC < candidate_reduced_costs[ell_ind])
                # Splitting this leaf does improve things, so go ahead with

                leaf = candidate_leaves[ell_ind]

                # Grow tree.
                numNodes = length(single_left)
                push!(single_left, 0, 0)
                push!(single_right, 0, 0)
                push!(single_isLeaf, true, true)
                push!(single_product, best_left_leaf_p_by_leaf[ell_ind], best_right_leaf_p_by_leaf[ell_ind])
                push!(single_depth, single_depth[leaf]+1, single_depth[leaf]+1)
                push!(single_include_set, union(single_include_set[leaf],best_split_p_by_leaf[ell_ind]), copy(single_include_set[leaf]) )
                push!(single_assortment_binary, copy(best_left_assortment_binary_by_leaf[ell_ind]), copy(best_right_assortment_binary_by_leaf[ell_ind]))
                push!(single_reduced_cost, best_left_reduced_cost_by_leaf[ell_ind], best_right_reduced_cost_by_leaf[ell_ind])
                push!(single_isClosed, false, false)

                # Update the leaf
                single_left[leaf] = numNodes+1
                single_right[leaf] = numNodes+2
                single_product[leaf] = best_split_p_by_leaf[ell_ind]
                single_isClosed[leaf] = true
                single_isLeaf[leaf] = false
            else
                # Leaf does not improve
                leaf = candidate_leaves[ell_ind]
                single_isClosed[leaf] = true
            end
        end
    end

    optimal_reduced_cost = sum( single_reduced_cost[single_isLeaf] )

    return single_left, single_right, single_product, single_isLeaf, optimal_reduced_cost

end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
################################################################################
#
#    This file contains methods to calculate the bispectrum of a given Configuration
#       Requires Configuration (see utils.jl). All of the command functions use LAMMPS.jl
#
#       The following methods are useful but require fitted snap.β
#       get_bispectrum calculates the bispectrum components using LAMMPS
#       get_dbispectrum calculuates the derivative of the bispectrum components (used for forces)
#       get_vbispectrum calculuates the virial contribution of the bispectrum components (used for virial stress tensor)
#
#       get_snap produces the A matrix in the SNAP potential fitting paradigm
#           A β = b
#           where β is the coefficients of the bispectrum components to produce energies, forces, stresses.
#           For detailed information, see http://dx.doi.org/10.1016/j.jcp.2014.12.018 (Thompson et al. 2014)
#       get_snap is instrumental in the methods that calculate potential_energy, force, virial for the snap potential.
################################################################################
using LAMMPS
function get_bispectrum(num_atoms, num_elements, file::String)
    num_coeff = num_coeffs
    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, num_coeff, num_atoms)
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units metal")
            command(lmp, "boundary f f f")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block -4 4 -4 4 -4 4")
            command(lmp, "read_data $(file)")

            # command(lmp, "mass 1 1")

            # Setup Forcefield
            command(lmp, "pair_style zero 6.0")
            command(lmp, "pair_coeff * *")

            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            radii_command = ""
            for r in radii
                radii_command *= "$r "
            end
            weight_command = ""
            for w in weight
                weight_command *= "$w "
            end
            string_command = "compute b all sna/atom $rcutfac $rcut0 $twojmax "
            string_command *= radii_command
            string_command *= weight_command
            string_command *= "rmin0 0.00 bnormflag $(Int(bnorm_flag)) bzeroflag $(Int(bzero_flag)) switchflag 0"

            if chem_flag
                chem_command = "chem $(num_elements) "
                for i = 0:(num_elements-1)
                    chem_command *= "$(i) "
                end
                string_command *= chem_command
            end
            command(lmp, string_command)
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "b", LAMMPS.API.LMP_STYLE_ATOM, LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A
end

function get_dbispectrum(num_atoms, num_elements, file::String)
    num_coeff = 3 * num_coeffs * num_elements
    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, num_coeff, num_atoms)
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units metal")
            command(lmp, "boundary f f f")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block -4 4 -4 4 -4 4")
            command(lmp, "read_data $(file)")

            # command(lmp, "mass 1 1")

            # Setup Forcefield
            command(lmp, "pair_style zero 6.0")
            command(lmp, "pair_coeff * *")

            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            radii_command = ""
            for r in radii
                radii_command *= "$r "
            end
            weight_command = ""
            for w in weight
                weight_command *= "$w "
            end
            string_command = "compute db all snad/atom $rcutfac $rcut0 $twojmax "
            string_command *= radii_command
            string_command *= weight_command
            string_command *= "rmin0 0.00 bnormflag $(Int(bnorm_flag)) bzeroflag $(Int(bzero_flag)) switchflag 0"

            if chem_flag
                chem_command = "chem $(num_elements) "
                for i = 0:(num_elements-1)
                    chem_command *= "$(i) "
                end
                string_command *= chem_command
            end
            command(lmp, string_command)
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "db", LAMMPS.API.LMP_STYLE_ATOM, LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A
end

function get_vbispectrum(num_atoms, num_elements, file::String)
    num_coeff = 6 * num_coeffs * num_elements
    A = LMP(["-screen", "none"]) do lmp
        bispectrum = Array{Float64}(undef, num_coeff, num_atoms)
        for i = 1
            # Initialize
            command(lmp, "log none")
            command(lmp, "units metal")
            command(lmp, "boundary f f f")
            command(lmp, "atom_style atomic")
            command(lmp, "atom_modify map array")

            # Setup box
            command(lmp, "region mybox block -4 4 -4 4 -4 4")
            command(lmp, "read_data $(file)")

            # command(lmp, "mass 1 1")

            # Setup Forcefield
            command(lmp, "pair_style zero 6.0")
            command(lmp, "pair_coeff * *")

            command(lmp, "compute PE all pe")
            command(lmp, "compute S all pressure thermo_temp")
            radii_command = ""
            for r in radii
                radii_command *= "$r "
            end
            weight_command = ""
            for w in weight
                weight_command *= "$w "
            end
            string_command = "compute vb all snav/atom $rcutfac $rcut0 $twojmax "
            string_command *= radii_command
            string_command *= weight_command
            string_command *= "rmin0 0.00 bnormflag $(Int(bnorm_flag)) bzeroflag $(Int(bzero_flag)) switchflag 0"

            if chem_flag
                chem_command = "chem $(num_elements) "
                for i = 0:(num_elements-1)
                    chem_command *= "$(i) "
                end
                string_command *= chem_command
            end
            command(lmp, string_command)
            command(lmp, "thermo_style custom pe")
            command(lmp, "run 0")

            ## Extract bispectrum
            bs = extract_compute(lmp, "vb", LAMMPS.API.LMP_STYLE_ATOM, LAMMPS.API.LMP_TYPE_ARRAY)
            bispectrum[:, :] = bs
            command(lmp, "clear")
        end
        return bispectrum
    end
    return A
end

A = get_bispectrum(num_atoms, num_elements, file)
dA = get_dbispectrum(num_atoms, num_elements, file)
vA = get_vbispectrum(num_atoms, num_elements, file)


text:
using ..common
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

# service methods for ApplicationClientProtocolService
const _ApplicationClientProtocolService_methods = MethodDescriptor[
        MethodDescriptor("getNewApplication", 1, GetNewApplicationRequestProto, GetNewApplicationResponseProto),
        MethodDescriptor("getApplicationReport", 2, GetApplicationReportRequestProto, GetApplicationReportResponseProto),
        MethodDescriptor("submitApplication", 3, SubmitApplicationRequestProto, SubmitApplicationResponseProto),
        MethodDescriptor("forceKillApplication", 4, KillApplicationRequestProto, KillApplicationResponseProto),
        MethodDescriptor("getClusterMetrics", 5, GetClusterMetricsRequestProto, GetClusterMetricsResponseProto),
        MethodDescriptor("getApplications", 6, GetApplicationsRequestProto, GetApplicationsResponseProto),
        MethodDescriptor("getClusterNodes", 7, GetClusterNodesRequestProto, GetClusterNodesResponseProto),
        MethodDescriptor("getQueueInfo", 8, GetQueueInfoRequestProto, GetQueueInfoResponseProto),
        MethodDescriptor("getQueueUserAcls", 9, GetQueueUserAclsInfoRequestProto, GetQueueUserAclsInfoResponseProto),
        MethodDescriptor("getDelegationToken", 10, GetDelegationTokenRequestProto, GetDelegationTokenResponseProto),
        MethodDescriptor("renewDelegationToken", 11, RenewDelegationTokenRequestProto, RenewDelegationTokenResponseProto),
        MethodDescriptor("cancelDelegationToken", 12, CancelDelegationTokenRequestProto, CancelDelegationTokenResponseProto),
        MethodDescriptor("moveApplicationAcrossQueues", 13, MoveApplicationAcrossQueuesRequestProto, MoveApplicationAcrossQueuesResponseProto),
        MethodDescriptor("getApplicationAttemptReport", 14, GetApplicationAttemptReportRequestProto, GetApplicationAttemptReportResponseProto),
        MethodDescriptor("getApplicationAttempts", 15, GetApplicationAttemptsRequestProto, GetApplicationAttemptsResponseProto),
        MethodDescriptor("getContainerReport", 16, GetContainerReportRequestProto, GetContainerReportResponseProto),
        MethodDescriptor("getContainers", 17, GetContainersRequestProto, GetContainersResponseProto),
        MethodDescriptor("submitReservation", 18, ReservationSubmissionRequestProto, ReservationSubmissionResponseProto),
        MethodDescriptor("updateReservation", 19, ReservationUpdateRequestProto, ReservationUpdateResponseProto),
        MethodDescriptor("deleteReservation", 20, ReservationDeleteRequestProto, ReservationDeleteResponseProto),
        MethodDescriptor("getNodeToLabels", 21, GetNodesToLabelsRequestProto, GetNodesToLabelsResponseProto),
        MethodDescriptor("getLabelsToNodes", 22, GetLabelsToNodesRequestProto, GetLabelsToNodesResponseProto),
        MethodDescriptor("getClusterNodeLabels", 23, GetClusterNodeLabelsRequestProto, GetClusterNodeLabelsResponseProto)
    ] # const _ApplicationClientProtocolService_methods
const _ApplicationClientProtocolService_desc = ServiceDescriptor("ApplicationClientProtocolService", 1, _ApplicationClientProtocolService_methods)

ApplicationClientProtocolService(impl::Module) = ProtoService(_ApplicationClientProtocolService_desc, impl)

type ApplicationClientProtocolServiceStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    ApplicationClientProtocolServiceStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_ApplicationClientProtocolService_desc, channel))
end # type ApplicationClientProtocolServiceStub

type ApplicationClientProtocolServiceBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    ApplicationClientProtocolServiceBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_ApplicationClientProtocolService_desc, channel))
end # type ApplicationClientProtocolServiceBlockingStub

getNewApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetNewApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[1], controller, inp, done)
getNewApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetNewApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[1], controller, inp)

getApplicationReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetApplicationReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[2], controller, inp, done)
getApplicationReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetApplicationReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[2], controller, inp)

submitApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::SubmitApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[3], controller, inp, done)
submitApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::SubmitApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[3], controller, inp)

forceKillApplication(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::KillApplicationRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[4], controller, inp, done)
forceKillApplication(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::KillApplicationRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[4], controller, inp)

getClusterMetrics(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetClusterMetricsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[5], controller, inp, done)
getClusterMetrics(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetClusterMetricsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[5], controller, inp)

getApplications(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetApplicationsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[6], controller, inp, done)
getApplications(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetApplicationsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[6], controller, inp)

getClusterNodes(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetClusterNodesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[7], controller, inp, done)
getClusterNodes(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetClusterNodesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[7], controller, inp)

getQueueInfo(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetQueueInfoRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[8], controller, inp, done)
getQueueInfo(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetQueueInfoRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[8], controller, inp)

getQueueUserAcls(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetQueueUserAclsInfoRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[9], controller, inp, done)
getQueueUserAcls(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetQueueUserAclsInfoRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[9], controller, inp)

getDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[10], controller, inp, done)
getDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[10], controller, inp)

renewDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::RenewDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[11], controller, inp, done)
renewDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::RenewDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[11], controller, inp)

cancelDelegationToken(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::CancelDelegationTokenRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[12], controller, inp, done)
cancelDelegationToken(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::CancelDelegationTokenRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[12], controller, inp)

moveApplicationAcrossQueues(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::MoveApplicationAcrossQueuesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[13], controller, inp, done)
moveApplicationAcrossQueues(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::MoveApplicationAcrossQueuesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[13], controller, inp)

getApplicationAttemptReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetApplicationAttemptReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[14], controller, inp, done)
getApplicationAttemptReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetApplicationAttemptReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[14], controller, inp)

getApplicationAttempts(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetApplicationAttemptsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[15], controller, inp, done)
getApplicationAttempts(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetApplicationAttemptsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[15], controller, inp)

getContainerReport(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetContainerReportRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[16], controller, inp, done)
getContainerReport(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetContainerReportRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[16], controller, inp)

getContainers(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetContainersRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[17], controller, inp, done)
getContainers(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetContainersRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[17], controller, inp)

submitReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::ReservationSubmissionRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[18], controller, inp, done)
submitReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::ReservationSubmissionRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[18], controller, inp)

updateReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::ReservationUpdateRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[19], controller, inp, done)
updateReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::ReservationUpdateRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[19], controller, inp)

deleteReservation(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::ReservationDeleteRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[20], controller, inp, done)
deleteReservation(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::ReservationDeleteRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[20], controller, inp)

getNodeToLabels(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetNodesToLabelsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[21], controller, inp, done)
getNodeToLabels(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetNodesToLabelsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[21], controller, inp)

getLabelsToNodes(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetLabelsToNodesRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[22], controller, inp, done)
getLabelsToNodes(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetLabelsToNodesRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[22], controller, inp)

getClusterNodeLabels(stub::ApplicationClientProtocolServiceStub, controller::ProtoRpcController, inp::GetClusterNodeLabelsRequestProto, done::Function) = call_method(stub.impl, _ApplicationClientProtocolService_methods[23], controller, inp, done)
getClusterNodeLabels(stub::ApplicationClientProtocolServiceBlockingStub, controller::ProtoRpcController, inp::GetClusterNodeLabelsRequestProto) = call_method(stub.impl, _ApplicationClientProtocolService_methods[23], controller, inp)

export ApplicationClientProtocolService, ApplicationClientProtocolServiceStub, ApplicationClientProtocolServiceBlockingStub, getNewApplication, getApplicationReport, submitApplication, forceKillApplication, getClusterMetrics, getApplications, getClusterNodes, getQueueInfo, getQueueUserAcls, getDelegationToken, renewDelegationToken, cancelDelegationToken, moveApplicationAcrossQueues, getApplicationAttemptReport, getApplicationAttempts, getContainerReport, getContainers, submitReservation, updateReservation, deleteReservation, getNodeToLabels, getLabelsToNodes, getClusterNodeLabels


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdmodal

"""
    'feffery'_antdmodal(;kwargs...)
    'feffery'_antdmodal(children::Any;kwargs...)
    'feffery'_antdmodal(children_maker::Function;kwargs...)

An AntdModal component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional)
- `id` (String; optional)
- `cancelButtonProps` (Dict; optional)
- `cancelCounts` (Real; optional)
- `cancelText` (String; optional)
- `centered` (Bool; optional)
- `className` (String; optional)
- `closable` (Bool; optional)
- `closeCounts` (Real; optional)
- `keyboard` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `locale` (a value equal to: 'zh-cn', 'en-us'; optional)
- `mask` (Bool; optional)
- `maskClosable` (Bool; optional)
- `okButtonProps` (Dict; optional)
- `okCounts` (Real; optional)
- `okText` (String; optional)
- `renderFooter` (Bool; optional)
- `style` (Dict; optional)
- `title` (optional): . title has the following type: String | lists containing elements 'content', 'prefixIcon'.
Those elements have the following types:
  - `content` (String; optional)
  - `prefixIcon` (String; optional)
- `visible` (Bool; optional)
- `width` (Real; optional)
"""
function 'feffery'_antdmodal(; kwargs...)
        available_props = Symbol[:children, :id, :cancelButtonProps, :cancelCounts, :cancelText, :centered, :className, :closable, :closeCounts, :keyboard, :loading_state, :locale, :mask, :maskClosable, :okButtonProps, :okCounts, :okText, :renderFooter, :style, :title, :visible, :width]
        wild_props = Symbol[]
        return Component("'feffery'_antdmodal", "AntdModal", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdmodal(children::Any; kwargs...) = 'feffery'_antdmodal(;kwargs..., children = children)
'feffery'_antdmodal(children_maker::Function; kwargs...) = 'feffery'_antdmodal(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
using Base:Integer, UInt
cd(@__DIR__)  # this sets the directory pointer to this file location
cd("..")
using Pkg  # Pkg.<command> is equivalent of `]` in command line julia
Pkg.activate(".")  # this specifies the environment being used and activates it

using Plots
# using TrackingHeaps
using Random
import Base.length
using BenchmarkTools
include("..\\utils\\ModularPlane.jl")
using Main.ModularPlane

"""
fundamentally in evolution every pattern (organism) needs a source of information (matter and energy) to succeed

self, discrete and deterministic, reasonable way of gaining information to allow reproduction
    - as soon as they have enough bits to replicate, they replicate (into a new grid)

    - or an incomplete copy of them idly travels around the next grid and is built up as the pattern finds more bits (in the form of other incomplete patterns)
        - motion: one cell is moved to a edge location which rotates counterclockwise each frame (so motion is quasi-circular)
        - or, the most recently moved cell is the starting point for counting (+x or cckmost margin), count n cckwise, where n is size of self

fractal evolution (replication is at smaller scales)? Then info is exponential
"""

seed = 633
Random.seed!(seed)

struct Pattern
    coords::Set{GridIndex}
end

length(p::Pattern) = length(p.coords)

### frame to frame updates
function updategrid!(patterns::Set{Pattern}, cells::Matrix{Int8}, tmp::Matrix{Int8})
    xdim, ydim = size(cells)
    tmp .= (cells .> 0) .+ cells
    newcells = Set{GridIndex}()
    println("sum ", sum(length.(patterns)), " ", sum(cells .> 0))
    for pattern in patterns
        pat_size = length(pattern.coords)
        if pat_size > 0
            pat_sum = sum(x -> cells[x], pattern.coords)
            newcoords = Set{GridIndex}()
            if pat_sum % 4 < 2
                extremum = argmax(x -> x.I[pat_sum % 2 + 1], pattern.coords)
            else
                extremum = argmin(x -> x.I[pat_sum % 2 + 1], pattern.coords)
            end
            start = GridIndex(extremum.I[1], extremum.I[2] + 1)
            for coord in pattern.coords
                if cells[coord] == AGE_CAP
                    prev = extremum
                    current = start
                    step = 1
                    cutoff = 4 * pat_size  # max possible surface area for squares
                    while (step <= pat_size - 1 || cells[current] >= AGE_MIN || current in newcells) && step <= cutoff
                        step += 1
                        prev, current = makestep(prev, current, pattern)
                    end
                    if step == cutoff + 1
                        current = coord  # if you can't find anywhere to grow, fill the space of the cell the just died
                    end
                    new = current
                    tmp[coord] = 0
                    push!(newcells, new)
                    tmp[new] = AGE_MIN
                    push!(newcoords, new)
                    start = new
                else
                    push!(newcoords, coord)
                end
            end
            updatecoords!(pattern, newcoords)
        else
            # println("dead pattern")
        end
    end
    patterns, tmp, cells
end

function indexof(vec::Vector{T}, val::T) where T
    for i in 1:length(vec)
        if vec[i] == val
            return i
        end
    end
    -1
end

function makestep(prev::GridIndex, current::GridIndex, pattern::Pattern)
    neighs = GridIndex.([   (current.I[1], current.I[2] + 1)
                            (current.I[1] + 1, current.I[2] + 1)
                            (current.I[1] + 1, current.I[2])
                            (current.I[1] + 1, current.I[2] - 1)
                            (current.I[1], current.I[2] - 1)
                            (current.I[1] - 1, current.I[2] - 1)
                            (current.I[1] - 1, current.I[2])
                            (current.I[1] - 1, current.I[2] + 1)])
    done = false
    started = false
    # println(current)
    # println("neighs ", neighs, "\n")
    i = indexof(neighs, prev)
    prev = current
    while i <= 16  # 2 len neighs
        i += 1
        idx = i % 8 + 1  # increment idx cyclically
        if !started
            if !(neighs[idx] in pattern.coords)  # is the location part of my body
                started = true
            end
        else
            if neighs[idx] in pattern.coords
                break  # we let current be the previous neighbor once we encounter a live cell
            end
        end
        current = neighs[idx]
    end
    prev, current
end

function updatecoords!(pattern::Pattern, coords::Set{GridIndex})
    for _ in 1:length(pattern.coords)
        pop!(pattern.coords)
    end
    for _ in 1:length(coords)
        push!(pattern.coords, pop!(coords))
    end
end

const AGE_MIN = 4
const AGE_CAP = 8
xdim, ydim = 100, 140
x = 1:xdim
y = 1:ydim

# initialize patterns
pattern_locs = [GridIndex(MInt(i, xdim), MInt(j, ydim)) for i in 0:20:xdim - 1 for j in 0:20:ydim - 1]
# pattern_locs = CartesianIndex.([(5, 10), (35, 10), (40, 25), (20, 30), (10, 55), (5, 50), (35, 40), (40, 60)])
# pattern_locs = CartesianIndex.([(5, 10)])
patterns = Set{Pattern}()
for ploc in pattern_locs
    xd = rand(1:12)
    yd = rand(1:12)
    n = rand((2xd * yd + 2) ÷ 3:xd * yd)
    possible_coords = Set{GridIndex}([GridIndex(i + ploc.I[1], j + ploc.I[2]) for i in 1:xd for j in 1:yd])
    coords = Set{GridIndex}(rand(possible_coords, n))
    push!(patterns, Pattern(coords))
end

cells = zeros(Int8, xdim, ydim)
for pattern in patterns
    for c in pattern.coords
        cells[c] = rand(AGE_MIN:AGE_CAP)
    end
end

iterations = 50
name = string("main ", xdim, "x", ydim, " seed=", seed, " iters=", iterations)
tmp = similar(cells)
bsize = 0

anim = @animate for i in 1:iterations
    # NOTE: x is vertical, y is horizontal, so swap the arguments for reasonable results
    heatmap(bsize + 1:ydim - bsize, bsize + 1:xdim - bsize, view(cells, bsize + 1:xdim - bsize, bsize + 1:ydim - bsize),
    colorbar=false, grid=false, ticks=false,
    showaxis=false, title=name, aspect_ratio=1, dpi=150)
    patterns, cells, tmp = updategrid!(patterns, cells, tmp)
end

gif(anim, string("gifs\\", name, ".gif"), fps=15)

max_stars_repo_name: AlexTMallen/the-multiverse

max_issues_repo_name: AlexTMallen/the-multiverse

max_forks_repo_name: AlexTMallen/the-multiverse


text:
# Type conversion from Vlasiator to Makie

"Conversion for 1D plots"
function Makie.convert_arguments(P::PointBased, meta::MetaVLSV, var)
   data = readvariable(meta, var)
   x = LinRange(meta.coordmin[1], meta.coordmax[1], meta.ncells[1])

   ([Point2f(i, j) for (i, j) in zip(x, data)],)
end

"Conversion for 2D plots."
function Makie.convert_arguments(P::SurfaceLike, meta::MetaVLSV, var,
   axisunit=EARTH, comp=0, normal=:y, origin=0.0)

   if meta.maxamr > 0
      pArgs = Vlasiator.set_args(meta, var, axisunit; normal, origin)
      data = Vlasiator.prep2dslice(meta, var, normal, comp, pArgs)
   else
      pArgs = Vlasiator.set_args(meta, var, axisunit)
      data = Vlasiator.prep2d(meta, var, comp)
   end

   x, y = Vlasiator.get_axis(axisunit, pArgs.plotrange, pArgs.sizes)

   (x, y, data)
end

"Conversion for 3D plots."
function Makie.convert_arguments(P::VolumeLike, meta::MetaVLSV, var,
   axisunit=EARTH, comp=1)
   (;ncells, coordmin, coordmax) = meta

   # Scale the sizes to the highest refinement level
   sizes = ncells .<< meta.maxamr # data needs to be refined later

   if axisunit == EARTH
      x = LinRange(coordmin[1], coordmax[1], sizes[1]) ./ Vlasiator.RE
      y = LinRange(coordmin[2], coordmax[2], sizes[2]) ./ Vlasiator.RE
      z = LinRange(coordmin[3], coordmax[3], sizes[3]) ./ Vlasiator.RE
   else
      x = LinRange(coordmin[1], coordmax[1], sizes[1])
      y = LinRange(coordmin[2], coordmax[2], sizes[2])
      z = LinRange(coordmin[3], coordmax[3], sizes[3])
   end

   if startswith(var, "fg")
      data = meta[var]
      if comp == 0
         data = [√(data[1,i,j,k]^2 + data[2,i,j,k]^2 + data[3,i,j,k]^2)
            for k in axes(data, 4), j in axes(data, 3), i in axes(data, 2)]
      else
         data = data[comp,:,:,:]
      end
   else
      data, _ = Vlasiator.fillmesh(meta, var)
      # Select on the finest refinement level
      if comp == 0
         data = data[1][end][1,:,:,:]
      else
         data = data[1][end][comp,:,:,:]
      end
   end

   (x, y, z, data)
end


text:
# Encoder structure : a downsampling is applied then two convolutions are done
struct Encoder
    dspl
    conv
end

@Flux.functor Encoder

(e::Encoder)(x) = x |> e.dspl |> e.conv

# Decoder structure : two convolutions are done then an upsampling is applied
struct Decoder
    conv
    uspl
end

@Flux.functor Decoder

(d::Decoder)(x) = x |> d.conv |> d.uspl

# Bridge structure :
#   - downsampling
#   - two convolutions
#   - upsampling
struct Bridge
    dspl
    conv
    uspl
end

@Flux.functor Bridge

(b::Bridge)(x) = x |> b.dspl |> b.conv |> b.uspl

# Encoder builders
function encoder(lvl)
    enc = []
    if lvl == 1
        push!(enc, identity)
    else
        push!(enc, MaxPool((2, 2)))
    end
    push!(enc, Conv((3, 3), 1=>1, pad = 1))
    push!(enc, Conv((3, 3), 1=>1, pad = 1))
    enc
end

encoder(t::Symbol, lvl) = encoder(Val(t), lvl)

encoder(t::Val{:array}, lvl) = encoder(lvl)

encoder(t::Val{:chain}, lvl) = Chain(encoder(lvl)...)

encoder(t::Val{:struct}, lvl) = begin
    enc = encoder(lvl)
    Encoder(enc[1], Chain(enc[2], enc[3]))
end

# Decoder builders
function decoder(lvl)
    dec = Any[Conv((3, 3), 2=>1, pad = 1), Conv((3, 3), 1=>1, pad = 1)]
    if lvl == 1
        push!(dec, identity)
    else
        push!(dec, Upsample(:bilinear, scale = 2))
    end
    dec
end

decoder(t::Symbol, lvl) = decoder(Val(t), lvl)

decoder(t::Val{:array}, lvl) = decoder(lvl)

decoder(t::Val{:chain}, lvl) = Chain(decoder(lvl)...)

decoder(t::Val{:struct}, lvl) = begin
    dec = decoder(lvl)
    Decoder(Chain(dec[1], dec[2]), dec[3])
end

# Bridge builders
bridge() = [MaxPool((2, 2)),
            Conv((3, 3), 1=>1, pad = 1), Conv((3, 3), 1=>1, pad = 1),
            Upsample(:bilinear, scale = 2)]

bridge(t::Symbol) = bridge(Val(t))

bridge(t::Val{:array}) = bridge()

bridge(t::Val{:chain}) = Chain(bridge()...)

bridge(t::Val{:struct}) = begin
    bdg = bridge()
    Bridge(bdg[1], Chain(bdg[2], bdg[3]), bdg[4])
end

@testset "uchain" begin
    blkt = [:array :chain :struct]
    for enct ∈ blkt, dect ∈ blkt, bdgt ∈ blkt
        enc, dec = [], []
        for l ∈ 1:2
            push!(enc, encoder(enct, l))
            push!(dec, decoder(dect, l))
        end
        bdg = bridge(bdgt)
        model = uchain(encoders = enc, decoders = dec, bridge = bdg, connection = chcat)
        x = rand(Float32, 32, 32, 1, 1)
        @test (model(x) |> size) == (32, 32, 1, 1)
    end
end