text:
function updateXP(M::Model)
    fvar = Future(M.variables)
    A = hcat([M[v.target.name, 0] for v in fvar]...)
    M.G(A, M.SP, M.XP)
    for i in 1:length(fvar)
        lb = fvar[i].target.bounds[1]::Float64
        ub = fvar[i].target.bounds[2]::Float64
        @inbounds @simd for j = 1:size(M.XP, 1)
            if M.XP[j,i] < lb
                M.XP[j,i] = lb
            elseif M.XP[j,i] > ub
                M.XP[j,i] = ub
            end
        end
    end
end

function forceboundsU(M::Model)
    ns = length(State(M.variables))
    for v ∈ Policy(M.variables)
        lb = v.bounds[1]::Float64
        ub = v.bounds[2]::Float64
        vid = findnext(M.variables, v, ns + 1)
        @inbounds @simd for i = 1:length(M)
            if isnan(M.X[i,vid])
                error("NaN in $(v.name) at $i")
            elseif M.X[i,vid] < lb
                M.X[i,vid] = lb
            elseif M.X[i,vid] > ub
                M.X[i,vid] = ub
            end
        end
    end
end

function solve(M::Model, n::Int = 1000, ϕ::Float64 = 0.8;disp::Int = 10000, crit = 1e-6)
    MaxError = maximum(abs.(M.Fval), 1)
    SError = sum(abs.(M.Fval) .> crit * 5, 1) / length(M)
    updateallD(M)
    SPids  = [(find(State(M.variables), s.name), find(State(M.variables, true), s.name))::Tuple{Int, Int} for s in Endogenous(M.variables)]
    ns = length(M.G.L)
    nG = length(M)

    ne = length(MaxError)
    Fval = zeros(ne)
    Jval = zeros(ne, ne)

    for iter = 1:n
        (iter <= 5 || mod(iter, 20) == 0 || iter < 5) && (M.F(M);MaxError = maximum(abs.(M.Fval), 1);SError = sum(abs.(M.Fval) .> crit * 5, 1) / length(M))
        for s in SPids
            for j = 1:size(M.ProbWeights, 2)
                for i = 1:nG
                    @inbounds M.SP[i + (j - 1) * nG, s[1]] = M.X[i, ns + s[2]]
                end
            end
        end
        updateXP(M)
        # M.Fall(M,ϕ)
        @threadsfixed [Fval, Jval] for i = 1:nG
            M.Finner(M, i, nG, Fval, Jval, ϕ)
        end
        forceboundsU(M)
        mod(iter, disp) == 0 && println(round(log10.(MaxError), 2), "  ", round(SError, 2))
        iter > min(n, 5) && (maximum(MaxError) < crit || all(SError .< 0.025)) && (print(iter);break)
    end
    # displaysolve(M,"done",n,MaxError,SError,time()-tstart)
    updateallD(M)
end

function updateallD(M::Model)
    for v in Dependant(M.variables)
        v.update(M)
    end
end

max_stars_repo_name: ZacLN/HADSGE

max_issues_repo_name: ZacLN/HADSGE

max_forks_repo_name: ZacLN/HADSGE


text:
# get face centroids for a single coordinate array
function coordinate_face_centroids(xf, md)
    Nfaces = size(md.FToF, 1)
    Nfp = size(md.xf, 1) ÷ Nfaces
    xc = reshape(xf, Nfp, Nfaces * md.K)
    return vec(typeof(xf)(sum(xc, dims=1) / size(xc, 1)))
end

"""
    function boundary_face_centroids(md)

Returns face centroids and `boundary_face_ids` on the boundaries of the domain given by md::MeshData.
"""
function boundary_face_centroids(md)
    compute_face_centroids(md) = map(x->coordinate_face_centroids(x, md), md.xyzf)
    xyzc = compute_face_centroids(md)
    boundary_face_ids = findall(vec(md.FToF) .== 1:length(md.FToF))

    # compute coordinates of face centroids on the boundary
    xyzb = map(x->x[boundary_face_ids], xyzc)
    return xyzb, boundary_face_ids
end

"""
    function tag_boundary_faces(md, boundary_name::Symbol = :entire_boundary)
    function tag_boundary_faces(md, boundary_list::Dict{Symbol, <:Function})

When called without arguments, just returns `Dict(:entire_boundary => boundary_faces)``.

Example usage:
```julia
julia> rd = RefElemData(Tri(), N=1)
julia> md = MeshData(uniform_mesh(Tri(), 2)..., rd)
julia> on_bottom_boundary(x, y, tol = 1e-13) = abs(y+1) < tol
julia> on_top_boundary(x, y, tol = 1e-13) = abs(y-1) < tol
julia> tag_boundary_faces(Dict(:bottom => on_bottom_boundary,
                               :top    => on_top_boundary), md)
```
"""
tag_boundary_faces(md, ::Nothing) = tag_boundary_faces(md)

function tag_boundary_faces(md, boundary_name::Symbol = :entire_boundary)
    return Dict(boundary_name => findall(vec(md.FToF) .== 1:length(md.FToF)))
end

function tag_boundary_faces(md, boundary_list::Dict{Symbol, <:Function})
    xyzb, boundary_face_ids = boundary_face_centroids(md)
    boundary_face_ids_list = _tag_boundary_faces(boundary_face_ids, boundary_list, xyzb)
    return Dict(Pair.(keys(boundary_list),boundary_face_ids_list))
end

# todo: should this take in boundary_face_flag((x,y,z)) instead of boundary_face_flag(x,y,z)?
function _tag_boundary_faces(boundary_face_ids, boundary_list, xyzb)
    boundary_face_ids_list = Vector{Int}[]
    for boundary_face_flag in values(boundary_list)
        push!(boundary_face_ids_list, boundary_face_ids[boundary_face_flag.(zip(xyzb...))])
    end
    return boundary_face_ids_list
end

# todo: should I make this version with NamedTuples the default?
function tag_boundary_faces(md, boundary_list::NamedTuple)
    xyzb, boundary_face_ids = boundary_face_centroids(md)
    boundary_face_ids_list = _tag_boundary_faces(boundary_face_ids, boundary_list, xyzb)
    return NamedTuple(Pair.(keys(boundary_list), boundary_face_ids_list))
end

ext: jl


text:
# import Optim  # broken
import AbstractFFTs
import FFTW

abstract type AbstractWeeks <: AbstractILt end

function Base.show(io::IO, w::AbstractWeeks)
    print(io, string(typeof(w)), "(Nterms=", w.Nterms, ",sigma=", w.sigma, ",b=", w.b,')')
end

### Weeks

mutable struct Weeks{T} <: AbstractWeeks
    func::Function
    Nterms::Int
    sigma::Float64
    b::Float64
    coefficients::Array{T,1}
end

function _get_coefficients(func, Nterms, sigma, b, ::Type{T}) where T <:Number
    a0 = _wcoeff(func, Nterms, sigma, b, T)
    return a0[Nterms+1:2*Nterms]
end

_get_coefficients(w::Weeks{T}) where T <: Number =  _get_coefficients(w.func, w.Nterms, w.sigma, w.b, T)
_set_coefficients(w::Weeks) = (w.coefficients = _get_coefficients(w))

const weeks_default_num_terms = 64

"""
   w::Weeks{datatype} = Weeks(func::Function, Nterms::Integer=64, sigma=1.0, b=1.0; datatype=Float64)

Return `w`, which estimates the inverse Laplace transform of `func` with
the Weeks algorithm. `w(t)` evaluates the transform at `t`. The accuracy depends on the choice
of `sigma` and `b`, with the optimal choices depending on `t`. `datatype` should agree with
the `DataType` returned by `func`. For convenience, `datatype=Complex` is equivalent to
`datatype=Complex{Float64}`

The call to `Weeks` that creates `w` is expensive relative to evaluation via `w(t)`.

# Example

Compute the inverse transform of the transform of `cos` at argument `pi/2`.
```
julia> ft = Weeks(s -> s/(s^2+1), 80);

julia> ft(pi/2)
0.0
```
"""
function Weeks(func::Function, Nterms::Integer=weeks_default_num_terms,
               sigma=1.0, b=1.0; datatype=Float64)
    outdatatype = datatype == Complex ? Complex{Float64} : datatype  # allow `Complex` as abbrev for Complex{Float64}
    return Weeks{outdatatype}(func, Nterms, sigma, b, _get_coefficients(func, Nterms, sigma, b, outdatatype))
end

function eval_ilt(w::Weeks, t)
    L = _laguerre(w.coefficients, 2 * w.b * t)
    return  L * exp((w.sigma - w.b) * t)
end

function eval_ilt(w::Weeks, t::AbstractArray)
    L = _laguerre(w.coefficients, 2 * w.b * t)
    return [L1 * exp((w.sigma - w.b) * t1) for (t1, L1) in zip(t, L)]
end

## FIXME: remove magic numbers below
function optimize(w::Weeks, t)
    (w.sigma, w.b) = _optimize_sigma_and_b(w.func, t, w.Nterms, 0.0, 30, 30)
    w.coefficients = _get_coefficients(w)
    return w
end

"""
    optimize(w::AbstractWeeks, t, Nterms)

Optimize the parameters of the inverse Laplace transform `w` at the
argument `t`. If `Nterms` is ommitted, the current value of `w.Nterms`
is retained.

The accuracy of the Weeks algorithm depends strongly on `t`. For some ranges
of `t`, the accuracy is relatively insensitive to the parameters. For other
values of `t`, even using optimized parameters results in estimates of
the inverse transform that are extremely inaccurate.

`optimize` is expensive in CPU time and allocation, it performs nested single-parameter
optimization over two parameterss.
"""
function optimize(w::AbstractWeeks, t, Nterms)
    w.Nterms = Nterms
    return optimize(w, t)
end

"""
    opteval(w::AbstractWeeks, t, Nterms)

estimate an inverse Laplace transform at argument `t` using `w` after
optimizing the parameters for `t`. If `Nterms` is omitted, then the
current value of `w.Nterms` is used.

Use `Weeks` or `WeeksErr` to create `w`.
"""
function opteval(w::AbstractWeeks, t, Nterms)
    optimize(w, t, Nterms)
    return w(t)
end

opteval(w::AbstractWeeks, t) = opteval(w, t, w.Nterms)

"""
    setparameters(w::AbstractWeeks, sigma, b, Nterms)

Set the parameters for the inverse Laplace transform object `w` and recompute
the internal data. Subsequent calls `w(t)` will use these parameters. If `Nterms`
or both `Nterms` and `b` are omitted, then their current values are retained.
"""
function setparameters(w::AbstractWeeks, sigma, b, Nterms)
    (w.sigma, w.b, w.Nterms) = (sigma, b, Nterms)
    _set_coefficients(w)
   return w
end

setparameters(w,sigma,b) = setparameters(w,sigma,b,w.Nterms)
setparameters(w,sigma) = setparameters(w,sigma,w.b)

function setNterms(w::AbstractWeeks, Nterms::Integer)
    w.Nterms = Nterms
    _set_coefficients(w)
    return w
end

#### WeeksErr

mutable struct WeeksErr{T} <: AbstractWeeks
    func::Function
    Nterms::Int
    sigma::Float64
    b::Float64
    coefficients::Array{T,1}
    sa1::Float64
    sa2::Float64
end

function _get_coefficients_and_params(func, Nterms, sigma, b, ::Type{T}) where T
    M = 2 * Nterms  # why 2 * Nterms ?
    a0 = _wcoeff(func, M, sigma, b, T)
    a1 = a0[2*Nterms+1:3*Nterms]
    sa1 = sum(abs.(a1))
    sa2 = sum(abs.(@view a0[3*Nterms+1:4*Nterms]))
    return (a1,sa1,sa2)
end

_get_coefficients_and_params(w::WeeksErr{T}) where T = _get_coefficients_and_params(w.func, w.Nterms, w.sigma, w.b, T)
_get_coefficients(w::WeeksErr{T}) where T = _get_coefficients_and_params(w.func, w.Nterms, w.sigma, w.b, T)
_set_coefficients(w::WeeksErr) = (w.coefficients, w.sa1, w.sa2) = _get_coefficients_and_params(w)

# FIXME: magic numbers here
function optimize(w::WeeksErr, t)
    (w.sigma, w.b) = _optimize_sigma_and_b(w.func, t, w.Nterms, 0.0, 30, 30)
    _set_coefficients(w)
    return w
end

"""
   w::WeeksErr{datatype} = WeeksErr(func::Function, Nterms::Integer=64, sigma=1.0, b=1.0; datatype=Float64)

Return `w`, which estimates the inverse Laplace transform of `func` via the Weeks algorithm.
`w(t)` returns a tuple containing the inverse transform at `t` and an error estimate. The accuracy of the
inversion depends on the choice of `sigma` and `b`. See the documentation for `Weeks` for a
description of the parameter `datatype`.

# Example

Compute the inverse transform of the transform of `cos`, and an error estimate, at argument `pi/2` using `80` terms.
```
julia> ft = WeeksErr(s -> s/(s^2+1), 80);

julia> ft(pi/2)
(0.0,3.0872097665938698e-15)
```
This estimate is more accurate than `cos(pi/2)`.
```
julia> ft(pi/2)[1] - cos(pi/2)
-6.123233995736766e-17

julia> ft(pi/2)[1] - 0.0         # exact value
0.0

julia> ft(pi/2)[1] - cospi(1/2)  # cospi is more accurate
0.0
```
"""
function WeeksErr(func::Function, Nterms::Integer=weeks_default_num_terms, sigma=1.0, b=1.0; datatype=Float64)
    outdatatype = datatype == Complex ? Complex{Float64} : datatype  # allow `Complex` as abbrev for Complex{Float64}
    params = _get_coefficients_and_params(func, Nterms, sigma, b, outdatatype)
    return WeeksErr{outdatatype}(func, Nterms, sigma, b, params...)
end

function eval_ilt(w::WeeksErr, t)
    L = _laguerre(w.coefficients,2*w.b*t)
    f = L * exp((w.sigma-w.b)*t)
    est = exp(w.sigma*t)*(w.sa2+eps()*w.sa1)
    return (f,est)
end

function eval_ilt(w::WeeksErr, t::AbstractArray)
    L = _laguerre(w.coefficients,2*w.b*t)
    f = L .* exp.((w.sigma-w.b)*t)
    est = exp.(w.sigma*t)*(w.sa2+eps()*w.sa1)
    return (f,est)
end

for w in (:Weeks, :WeeksErr)
    @eval (w::$(w))(t) = eval_ilt(w,t)
end

##### internal functions

_wcoeff(F, N, sig, b, ::Type{T}) where T <: Real = real(_wcoeff(F, N, sig, b))
_wcoeff(F, N, sig, b, ::Type{T}) where T <: Complex = _wcoeff(F, N, sig, b)

function _wcoeff(F, N, sig, b)
    n = -N:N-1  # FIXME: remove 1 and test
    h = pi / N # FIXME: what data type ?
    th = h .* (n .+ 1//2)
    y = b .* cot.(th / 2)
    imaginary_unit = Complex(zero(eltype(y)), one(eltype(y)))
    s = sig .+ imaginary_unit * y
    FF0 = map(F, s)
    FF = [FF1 * (b + imaginary_unit * y1) for (FF1,y1) in zip(FF0,y)]
    a = (FF |> AbstractFFTs.fftshift |> FFTW.fft |> AbstractFFTs.fftshift) / (2 * N)
    return exp.(Complex(zero(h),-one(h)) * n*h/2) .* a
end

function _laguerre(a::AbstractVector,x::AbstractArray)
    N = length(a) - 1
    unp1 = zero(x)
    un = a[N+1] * one(x)
    local unm1
    for n in N:-1:1
#        unm1 =  (1//n)*(2*n-1-x) .* un - n/(n+1)*unp1 + a[n]
        unm1 =  [(1//n)*(2*n-1-x0) * un0 - n/(n+1)*unp10 + a[n] for (x0,un0,unp10) in zip(x,un,unp1)]
        unp1 = un
        un = unm1
    end
    return unm1
end

function _laguerre(a::AbstractVector,x)
    n = length(a) - 1
    unp1 = zero(x)
    un = a[n+1]*one(x)
    local unm1
    for n in n:-1:1
        unm1 = (1//n)*(2*n-1-x) * un - n/(n+1)*unp1 + a[n]
        unp1 = un
        un = unm1
    end
    return unm1
end

function _optimize_sigma_and_b(F, t, N, sig0, sigmax, bmax)
    sigma_opt = Optim.minimizer(Optim.optimize(sig -> _werr2e(sig, F,t,N, sig0, sigmax, bmax), sig0, sigmax))
    b_opt = Optim.minimizer(Optim.optimize(b -> _werr2t(b, F, N, sigma_opt), 0, bmax))
    return (sigma_opt, b_opt)
end

function _werr2e(sig, F, t, N, sig0, sigmax, bmax)
    b = Optim.minimizer(Optim.optimize((b) -> _werr2t(b, F,N, sig) , 0.0, bmax))
    M = 2*N
    a = _wcoeff(F, M, sig, b)
    a1 = @view a[2*N+1:3*N]
    sa1 = sum(abs.(a1))
    a2 = @view a[3*N+1:4*N]
    sa2 = sum(abs.(a2))
    return sig*t + log(sa2+eps()*sa1)
end

function _werr2t(b, F, N, sig)
    M = 2*N
    a = _wcoeff(F,M,sig,b)
    sa2 = sum(abs.(@view a[3*N+1:4*N]))
    return log(sa2)
end


text:
using DifferentialEquations
using Plots
using Markdown
using LinearAlgebra
using StaticArrays

description=md"""

solve the matrix system
``u' = Au``
for matrix A and i.c. u₀

"""
display(description)

# use immutable arrays that are stack-allocated
A = @SMatrix [1.0 0 0 -5
              4 -2 4 -3
              -4 0 0 1
              5 -2 2 3]

u0 = rand(4, 2)

tspan = (0.0, 1.0)

f(du, u, p, t) = mul!(du, A, u)  # in-place version of matrix multiply

prob = ODEProblem(f, u0, tspan)
sol = solve(prob)
plot(sol)

max_stars_repo_name: john-waczak/SciML-examples

max_issues_repo_name: john-waczak/SciML-examples

max_forks_repo_name: john-waczak/SciML-examples


text:
module NicheMap

using RCall
using DataFrames
using Unitful
using Microclimate

export nichemap_global

"""
    nichemap_global(location; [years, runmoist, timeinterval])

Runs the micro_global function formm the R package NicheMapR with passed in args.
Returns a nested dataframe of the results with correct column names.

# TODO allow vararg keyword splats
"""
function nichemap_global(location; years = 1, runmoist=1, timeinterval=365)
    @rput years location runmoist timeinterval
    R"""
    library(NicheMapR)
    micro <- micro_global(timeinterval=timeinterval, nyears=years, runmoist=runmoist, loc=location)

    # Convert R matrices to dataframes to preserve column names.
    for (name in names(micro)) {
        if (is.matrix(micro[[name]])) {
          micro[[name]] <- as.data.frame(micro[[name]])
        }
    }
    """
    @rget micro

    return MicroclimateData(values(micro)...)
end

end # module


text:
export viagroup

"""
    viagroup(radius, height, n, [parameters::SegmentParameters [,topequiv [,botequiv]]])
    viagroup(<keyword arguments>)

Returns a `Group` of `Segment`s for the barrel of a via.

**Keyword Arguments**

- `radius`          -- radius of via
- `height`          -- height of via
- `h`               -- thickness of copper plating
- `n`               -- number of segment around via
- `sigma`,`rho`     -- specify conductivity or resistivity
- `nhinc`           -- number of filaments per segment
- `rh`              -- ratio of filament width
- `topequiv`        -- include .equiv for all top nodes.  default = true
- `botequiv`        -- include .equiv for all bot nodes.  default = true

Terminals are in the center labeled :top and :bot.  Only `SegmentParameters`
which match keyword arguments are used.  Specify n=0 to bypass via (.equiv top bot).
Terminals :alltop and :allbot are arrays of all top and bottom nodes.  :alltop, :allbot,
`topequiv`, and `botequiv` are intended for connections to planes.
"""
viagroup(;radius=NaN, height=NaN, n=8, h=NaN, sigma=NaN, rho=NaN,
          nhinc=5, rh=NaN, topequiv = true, botequiv = true) =
  viagroup(radius, height, n,
          SegmentParameters(h=h,sigma=sigma,rho=rho,nhinc=nhinc,rh=rh),
          topequiv,botequiv)
function viagroup(radius, height, n=8,
                  sp::SegmentParameters=SegmentParameters(),
                  topequiv=true, botequiv=true)
  if isnan(radius)
    throw(ArgumentError("radius must be specified"))
  end
  if isnan(height)
    throw(ArgumentError("height must be specified"))
  end
  if isnan(h(sp))
    throw(ArgumentError("wall thickness (height of segment) must be specified"))
  end
  if n!=0 && n<2
    throw(ArgumentError("must have zero or at least 2 segments"))
  end
  group = Group()
  terminals = terms(group)
  viaelements = elements(group)
  centertop = Node(0,0,0)
  centerbot = Node(0,0,-height)
  terminals[:top] = centertop
  terminals[:bot] = centerbot
  if n!=0
    angle = range(0.0, step=2*π/n, length=n)
    x = radius .* map(cos,angle)
    y = radius .* map(sin,angle)
    wx = map(cos,angle.+(pi/2))
    wy = map(sin,angle.+(pi/2))
    w = sqrt((x[1]-x[2])^2+(y[1]-y[2])^2) * 1/cos(pi/n)
    topnodes = Array{Node}(undef,n)
    botnodes = Array{Node}(undef,n)
    for i in 1:n
      topnodes[i] = Node(x[i],y[i],0.0)
      botnodes[i] = Node(x[i],y[i],-height)
    end
    eqtop = Equiv([centertop;topnodes])
    eqbot = Equiv([centerbot;botnodes])
    segments = Array{Segment}(undef,n)
    for i in 1:n
      segments[i] = Segment(topnodes[i],botnodes[i],
        SegmentParameters(sp, w=w, wx=wx[i], wy=wy[i], wz=0.0, nwinc=1))
    end
    terminals[:alltop] = topnodes
    terminals[:allbot] = botnodes
    append!(viaelements, [topnodes;botnodes;centertop;centerbot])
    if topequiv
      push!(viaelements, eqtop)
    end
    if botequiv
      push!(viaelements, eqbot)
    end
    append!(viaelements, segments)
  else # n=0 for bypass
    terminals[:alltop] = Array{Node}([centertop])
    terminals[:allbot] = Array{Node}([centerbot])
    eq = Equiv([centertop;centerbot])
    append!(viaelements,[centertop;centerbot;eq])
  end
  return group
end


text:
# transpose, majortranspose, minortranspose
"""
    transpose(::Vec)
    transpose(::SecondOrderTensor)
    transpose(::FourthOrderTensor)

Compute the transpose of a tensor.
For a fourth order tensor, the transpose is the minor transpose.

# Examples
```jldoctest
julia> A = rand(Tensor{2,2})
2×2 Tensor{2,2,Float64,4}:
 0.590845  0.566237
 0.766797  0.460085

julia> A'
2×2 Tensor{2,2,Float64,4}:
 0.590845  0.766797
 0.566237  0.460085
```
"""
@inline function Base.transpose(S::Tensor{2, dim}) where {dim}
    Tensor{2, dim}(@inline function(i, j) @inbounds S[j,i]; end)
end

@inline Base.transpose(S::SymmetricTensor{2}) = S

"""
    minortranspose(::FourthOrderTensor)

Compute the minor transpose of a fourth order tensor.
"""
@inline function minortranspose(S::Tensor{4, dim}) where {dim}
    Tensor{4, dim}(@inline function(i, j, k, l) @inbounds S[j,i,l,k]; end)
end

@inline minortranspose(S::SymmetricTensor{4}) = S
@inline Base.transpose(S::FourthOrderTensor) = minortranspose(S)

"""
    majortranspose(::FourthOrderTensor)

Compute the major transpose of a fourth order tensor.
"""
@inline function majortranspose(S::FourthOrderTensor{dim}) where {dim}
    Tensor{4, dim}(@inline function(i, j, k, l) @inbounds S[k,l,i,j]; end)
end

@inline Base.adjoint(S::AllTensors) = transpose(S)


text:
checkDir()
cdw = pwd()

#fix the seed
Random.seed!(1)

if nprocs() <= 2
    p = addprocs(2)
end
@everywhere using DistributedArrays
@everywhere using GigaSOM
@everywhere using Distances

# only use lineageMarkers for clustering
(lineageMarkers,)= getMarkers(panel)
cc = map(Symbol, lineageMarkers)
dfSom = daf.fcstable[:,cc]

# concatenate the dataset for performance testing
n = 0
for i in 1:n
    global dfSom
    dfSom = vcat(dfSom, dfSom)
end

som2 = initGigaSOM(dfSom, 10, 10)

@testset "Dimensions - parallel" begin
    @test size(som2.codes) == (100,10)
    @test som2.xdim == 10
    @test som2.ydim == 10
    @test som2.numCodes == 100
end

som2 = trainGigaSOM(som2, dfSom, epochs = 2, rStart = 6.0)

winners = mapToGigaSOM(som2, dfSom)

embed = embedGigaSOM(som2, dfSom, k=10, smooth=0.0, adjust=0.5)

#test parallel
@testset "Parallel" begin
    codes = som2.codes
    @test size(codes) == (100,10)

    dfCodes = DataFrame(codes)
    names!(dfCodes, Symbol.(som2.colNames))
    dfEmbed = DataFrame(embed)
    CSV.write(genDataPath*"/parallelDfCodes.csv", dfCodes)
    CSV.write(genDataPath*"/parallelWinners.csv", winners)
    CSV.write(genDataPath*"/parallelEmbedded.csv", dfEmbed)

    # load the ref data
    refParallelDfCodes = CSV.File(refDataPath*"/refParallelDfCodes.csv") |> DataFrame
    refParallelWinners = CSV.File(refDataPath*"/refParallelWinners.csv") |> DataFrame
    refParallelEmbedded = CSV.File(refDataPath*"/refParallelEmbedded.csv") |> DataFrame

    # load the generated data
    parallelDfCodes = CSV.File(genDataPath*"/parallelDfCodes.csv") |> DataFrame
    parallelDfCodesTest = first(parallelDfCodes, 10)
    parallelWinners = CSV.File(genDataPath*"/parallelWinners.csv") |> DataFrame
    parallelWinnersTest = first(parallelWinners, 10)
    parallelEmbedded = CSV.File(genDataPath*"/parallelEmbedded.csv") |> DataFrame
    parallelEmbeddedTest = first(parallelEmbedded, 10)

    # test the generated data against the reference data
    @test refParallelWinners == parallelWinnersTest
    @test refParallelDfCodes == parallelDfCodesTest
    @test Array{Float64,2}(refParallelEmbedded) ≈ Array{Float64,2}(parallelEmbeddedTest) atol=1e-4

    #test parallel
    for (i, j) in zip(parallelDfCodes[:,1], refParallelDfCodes[:,1])
        @test isapprox(i, j; atol = 0.001)
    end

end

rmprocs(workers())

cd(cdw)


text:
"""
**Adamax Constructor**

```julia
Adamax(;α=0.002, β₁=0.9, β₂=0.999, ϵ=10e-8)
```
Algorithm

```math
\\begin{align*}
m_t =& \\beta_1 m_{t-1} + (1-\\beta_1)g_t\\\\
u_t =& \\max(\\beta_2 u_{t}, |g_t|)\\\\
\\Delta x_t =& \\frac{\\alpha}{(1-\\beta_1^t)}\\frac{m_t}{u_t}
\\end{align*}
```

[Algorithm Reference](https://arxiv.org/abs/1412.6980)
"""
mutable struct Adamax <: Optimizer
    opt_type::String
    t::Int64
    α::Float64
    β₁::Float64
    β₂::Float64
    m_t::AbstractArray
    u_t::AbstractArray
end

function Adamax(;α::Real=0.002, β₁::Real=0.9, β₂::Real=0.998)
    @assert α > 0.0 "α must be greater than 0"
    @assert β₁ > 0.0 "β₁ must be greater than 0"
    @assert β₂ > 0.0 "β₂ must be greater than 0"

    Adamax("Adamax", 0, α, β₁, β₂, [], [])
end

params(opt::Adamax) = "α=$(opt.α), β₁=$(opt.β₁), β₂=$(opt.β₂)"

function update(opt::Adamax, g_t::AbstractArray{T}) where {T<:Real}
    # resize biased moment estimates if first iteration
    if opt.t == 0
        opt.m_t = zero(g_t)
        opt.u_t = zero(g_t)
    end

    # update timestep
    opt.t += 1

    # update biased first moment estimate
    opt.m_t = opt.β₁ * opt.m_t + (one(T) - opt.β₁) * g_t

    # update the exponentially weighted infinity norm
    opt.u_t = max.(opt.β₂ * opt.u_t, abs.(g_t))

    # update parameters
    ρ = (opt.α / (one(T) - opt.β₁^opt.t)) * opt.m_t ./ opt.u_t

    return ρ
end


text:
struct ConstCovariate{T<:AbstractFloat} <: AbstractCovariate{T}
    name::String
    length::Integer
    value::T
end

function ConstCovariate(value::T, length::Integer) where {T<:AbstractFloat}
    ConstCovariate{T}("", length, value)
end

Base.length(covariate::ConstCovariate{T}) where {T<:AbstractFloat} = covariate.length

function slice(covariate::ConstCovariate{T}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:AbstractFloat}
    if fromobs > toobs
        EmptySeq{SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}}()
    else
        fromobs = max(1, fromobs)
        toobs = min(toobs, length(covariate))
        slicelength = verifyslicelength(fromobs, toobs, slicelength)
        buffer = Vector{T}(slicelength)
        fill!(buffer, covariate.value)
        map(Seq(Tuple{Int64, Int64}, (fromobs, toobs, slicelength), nextslice), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do rng
            view(buffer, 1:(rng[2] - rng[1] + 1))
        end
    end
end

function Base.map(covariate::ConstCovariate, dataframe::AbstractDataFrame)
    ConstCovariate(covariate.name, length(dataframe), covariate.value)
end

max_stars_repo_name: vikashdwd/JuML

max_issues_repo_name: vikashdwd/JuML

max_forks_repo_name: vikashdwd/JuML


text:
using Documenter, AnonymousTypes

makedocs(
    modules = AnonymousTypes,
    clean   = false,
)

deploydocs(
    deps = Deps.pip("pygments", "mkdocs", "mkdocs-material"),
    repo = "github.com/MichaelHatherly/AnonymousTypes.jl.git",
)

max_stars_repo_name: MichaelHatherly/AnonymousTypes

max_issues_repo_name: MichaelHatherly/AnonymousTypes

max_forks_repo_name: MichaelHatherly/AnonymousTypes


text:
module LibVulkan

include("CEnum.jl")
using .CEnum

import Libdl

@static if !isempty(get(ENV, "JULIA_VULKAN_LIBNAME", ""))
    const libvulkan = ENV["JULIA_VULKAN_LIBNAME"]
elseif Sys.iswindows()
    const libvulkan = "vulkan-1.dll"
elseif Sys.isapple()
    const libvulkan = "libvulkan.dylib"
elseif Sys.islinux()
    const libvulkan = "libvulkan.so.1"
else
    const libvulkan = "libvulkan"
end

const libvulkan_handle = Ref{Ptr{Cvoid}}(0)

function __init__()
    libname = Libdl.find_library(libvulkan)
    if isempty(libname)
            error("""
            Failed to retrieve a valid Vulkan library called '$libvulkan'.
            If you configure the `JULIA_VULKAN_LIBNAME` environment variable before precompiling VulkanCore, it will be used instead of '$libvulkan'. You may also manually add search paths by appending them to Lidbl.DL_LOAD_PATH, but note that this may have repercussions beyond this package.
            """)
    end
    libvulkan_handle[] = Libdl.dlopen(libname)
end

VK_MAKE_VERSION(major, minor, patch) = (Cuint(major) << 22) | (Cuint(minor) << 12) | patch
VK_MAKE_API_VERSION(variant, major, minor, patch) = (Cuint(variant) << 29) | (Cuint(major) << 22) | (Cuint(minor) << 12) | Cuint(patch)

VK_VERSION_MAJOR(version) = Cuint(version) >> 22
VK_VERSION_MINOR(version) = (Cuint(version) >> 12) & 0x3ff
VK_VERSION_PATCH(version) = Cuint(version) & 0xfff

const IS_LIBC_MUSL = occursin("musl", Base.BUILD_TRIPLET)
if Sys.isapple() && Sys.ARCH === :aarch64
    include("../lib/aarch64-apple-darwin20.jl")
elseif Sys.islinux() && Sys.ARCH === :aarch64 && !IS_LIBC_MUSL
    include("../lib/aarch64-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :aarch64 && IS_LIBC_MUSL
    include("../lib/aarch64-linux-musl.jl")
elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && !IS_LIBC_MUSL
    include("../lib/armv7l-linux-gnueabihf.jl")
elseif Sys.islinux() && startswith(string(Sys.ARCH), "arm") && IS_LIBC_MUSL
    include("../lib/armv7l-linux-musleabihf.jl")
elseif Sys.islinux() && Sys.ARCH === :i686 && !IS_LIBC_MUSL
    include("../lib/i686-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :i686 && IS_LIBC_MUSL
    include("../lib/i686-linux-musl.jl")
elseif Sys.iswindows() && Sys.ARCH === :i686
    include("../lib/i686-w64-mingw32.jl")
elseif Sys.islinux() && Sys.ARCH === :powerpc64le
    include("../lib/powerpc64le-linux-gnu.jl")
elseif Sys.isapple() && Sys.ARCH === :x86_64
    include("../lib/x86_64-apple-darwin14.jl")
elseif Sys.islinux() && Sys.ARCH === :x86_64 && !IS_LIBC_MUSL
    include("../lib/x86_64-linux-gnu.jl")
elseif Sys.islinux() && Sys.ARCH === :x86_64 && IS_LIBC_MUSL
    include("../lib/x86_64-linux-musl.jl")
elseif Sys.isbsd() && !Sys.isapple()
    include("../lib/x86_64-unknown-freebsd11.1.jl")
elseif Sys.iswindows() && Sys.ARCH === :x86_64
    include("../lib/x86_64-w64-mingw32.jl")
else
    error("Unknown platform: $(Base.BUILD_TRIPLET)")
end

# exports
const PREFIXES = ["VK_", "Vk", "vk"]
for name in names(@__MODULE__; all=true), prefix in PREFIXES
    if startswith(string(name), prefix)
        @eval export $name
    end
end

end # module


text:
function generic_longstaff_schwartz_regression!(simulationData::Vector{Vector{NodeData}}, basisCoefficients::Vector{Vector{Float64}})
  steps = length(simulationData)
  resize!(basisCoefficients, steps-1)

  for i = steps:-1:2
    exerciseData = simulationData[i]

    # 1) Find the covariance matrix of basis function values and deflated cash-flows
    N = length(exerciseData[1].values)
    temp = Vector{Float64}(N+1)
    stats = GenericSequenceStats(N+1, length(exerciseData))

    for j in eachindex(exerciseData)
      if exerciseData[j].isValid
        temp[1:end-1] = exerciseData[j].values
        temp[end] = exerciseData[j].cumulatedCashFlows - exerciseData[j].controlValue

        add_sample!(stats, temp, j)
      end
    end

    means = stats_mean(stats)
    covMat = stats_covariance(stats)

    C = Matrix{Float64}(N, N)
    target = Vector{Float64}(N)

    for k in eachindex(target)
      target[k] = covMat[k, N+1] + means[k] * means[N+1]
      for l = 1:k
        C[k, l] = C[l, k] = covMat[k,l] + means[k] * means[l]
      end
    end

    # 2) Solve for least squares regression
    alphas = svdfact(C) \ target
    # resize!(basisCoefficients[i-1], N)
    basisCoefficients[i-1] = copy(alphas) # copy

    # 3) use exercise strategy to divide paths into exercise and non-exercise domains
    for j in eachindex(exerciseData)
      if exerciseData[j].isValid
        exerciseValue = exerciseData[j].exerciseValue
        continuationValue = exerciseData[j].cumulatedCashFlows
        estimatedContinuationValue = dot(exerciseData[j].values, alphas)

        # for exercise paths, add deflated rebate to
        # deflated cash-flows at previous time frame;
        # for non-exercise paths, add deflated cash-flows to
        # deflated cash-flows at previous time frame
        value = estimatedContinuationValue <= exerciseValue ? exerciseValue : continuationValue

        simulationData[i-1][j].cumulatedCashFlows += value
      end
    end
  end

  # the value of the prodcut can now be estimated by averaging over all paths
  estimatedData = simulationData[1]
  estimate = gen_RiskStatistics(length(estimatedData))

  for j in eachindex(estimatedData)
    add_sample!(estimate, estimatedData[j].cumulatedCashFlows, j)
  end

  return stats_mean(estimate)
end


text:
using ThreeBodyProblem
using DifferentialEquations
using Plots

T = 1

sys = sun_earth()
Lpts = computeLpts(sys)
L1 = Lpts[1]
tf = 1*JY/sys.TUNIT
rv₀ = vcat(L1, zeros(3))
@time Wsp, Wsn, Wup, Wun, Φₜ = invariant_manifolds(sys, rv₀, T, tf)
# Juno.@enter invariant_manifolds(sys, rv₀, T, tf)
plot(sys)
plot!(Wup,vars=(1,2),label="Wu+",linecolor=:red)
plot!(Wun,vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(Wsp,vars=(1,2),label="Ws+",linecolor=:blue)
plot!(Wsn,vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,ylims=[-0.04,0.04],xlims=[0.93,1.06],legend=:outerright,flip=false)

W78 = []
C = []
a = []
e = []
t = []
systems = (sun_mercury(), sun_venus(), sun_earth(), sun_mars(), sun_jupiter(), sun_saturn(), sun_uranus(), sun_neptune())
for sys in systems
    L1, L2 = computeLpts(sys)
    C1 = []
    a1 = []
    e1 = []
    t1 = []
    for Lpt = (L1, L2)
        tf = 100*JY/sys.TUNIT # integrate for 100 years
        rv₀ = vcat(Lpt, zeros(3))
        Wsp, Wsn, Wup, Wun, Φₜ = invariant_manifolds(sys, rv₀, T, tf)
        push!(W78,Wsp,Wsn,Wup,Wun)

        C2 = []
        a2 = []
        e2 = []
        t2 = []
        for W0 = (Wsp, Wsn, Wup, Wun)
            push!(C2, [computeC(W0.u[j], sys) for j = 1:length(W0.u)])
            rv_inert = [rot2inert(W0.u[j], W0.t[j], sys.μ, origin=:prim) for j = 1:length(W0.u)]
            rv_dim = [dimensionalize(rv_inert[i], sys) for i = 1:length(rv_inert)]
            push!(a2, [cart2oe(rv_dim[i], sys.μ₁)[1] for i = 1:length(rv_dim)])
            push!(e2, [cart2oe(rv_dim[i], sys.μ₁)[2] for i = 1:length(rv_dim)])
            push!(t2, [W0.t[j]*sys.TUNIT/JY for j = 1:length(W0.u)])
        end
        push!(C1, C2)
        push!(a1, a2)
        push!(e1, e2)
        push!(t1, t2)
    end
    push!(C,C1)
    push!(a,a1)
    push!(e,e1)
    push!(t,t1)

        # plot(sys)
        # plot!(Wsp,vars=(1,2),label="Ws+",linecolor=:blue)
        # plot!(Wsn,vars=(1,2),label="Ws-",linecolor=:cyan)
        # plot!(Wup,vars=(1,2),label="Wu+",linecolor=:red)
        # plot!(Wun,vars=(1,2),label="Wu-",linecolor=:magenta)
        # plot!(aspect_ratio=1,legend=:outerright,flip=false)
end

# Mercury L1 SP
plot(t[1][1][1], C[1][1][1],label="Wₛ⁺",xlabel="time [years]", ylabel="Jacobi Constant, C", title="Mercury L1 RK78",color=:blue)
# Mercury L2 SP
plot(t[1][2][1], C[1][2][1],label="Wₛ⁺",xlabel="time [years]", ylabel="Jacobi Constant, C", title="Mercury L2 RK78",color=:blue)
# Mercury L1 SN
plot(t[1][1][2], C[1][1][2],label="Wₛ⁻",xlabel="time [years]", ylabel="Jacobi Constant, C", title="Mercury L1 RK78",color=:cyan)
# Mercury L2 SN
plot(t[1][2][2], C[1][2][2],label="Wₛ⁻",xlabel="time [years]", ylabel="Jacobi Constant, C", title="Mercury L2 RK78",color=:cyan)
plot(sys)
plot!(W78[1],vars=(1,2),label="Ws+",linecolor=:blue)
plot!(W78[2],vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(W78[5],vars=(1,2),label="Ws+",linecolor=:blue)
plot!(W78[6],vars=(1,2),label="Ws-",linecolor=:cyan)

scatter(a[1][1][1]/AU,e[1][1][1], color=:blue, label="Mercury", markerstrokewidth=0)
scatter!(a[1][1][2]/AU,e[1][1][2], color=:cyan, markerstrokewidth=0)
scatter!(a[1][1][3]/AU,e[1][1][3], color=:red, markerstrokewidth=0)
scatter!(a[1][1][4]/AU,e[1][1][4], color=:magenta, markerstrokewidth=0)

scatter(a[1][2][1]/AU,e[1][2][1], color=:blue, label="Mercury", markerstrokewidth=0)
scatter!(a[1][2][2]/AU,e[1][2][2], color=:cyan, markerstrokewidth=0)
scatter!(a[1][2][3]/AU,e[1][2][3], color=:red, markerstrokewidth=0)
scatter!(a[1][2][4]/AU,e[1][2][4], color=:magenta, markerstrokewidth=0)
#L1 SN UP
#L2 SP UN

scatter(a[1][1][1]/AU,e[1][1][1], color=:grey, label="Mercury", markerstrokewidth=0)
scatter!(a[2][1][1]/AU,e[2][1][1], color=:orange, label="Venus", markerstrokewidth=0)
scatter!(a[3][1][1]/AU,e[3][1][1], color=:blue, label="Earth", markerstrokewidth=0)
scatter!(a[4][1][1]/AU,e[4][1][1], color=:red, label="Mars", markerstrokewidth=0)
xlabel!("a [AU]")
ylabel!("e")
scatter!(a[5][1][1]/AU,e[5][1][1], color=:blue, label="Mercury", markerstrokewidth=0)
scatter!(a[6][1][1]/AU,e[6][1][1], color=:blue, label="Mercury", markerstrokewidth=0)
scatter!(a[7][1][1]/AU,e[7][1][1], color=:blue, label="Mercury", markerstrokewidth=0)
scatter!(a[8][1][1]/AU,e[8][1][1], color=:blue, label="Mercury", markerstrokewidth=0)

for (i, sys) in enumerate(systems)
    display(i)
    display(sys)
    plot(sys)
    plot!(W[(i-1)*8+1],vars=(1,2),label="Wu+",linecolor=:red)
    plot!(W[(i-1)*8+2],vars=(1,2),label="Wu-",linecolor=:magenta)
    plot!(W[(i-1)*8+3],vars=(1,2),label="Ws+",linecolor=:blue)
    plot!(W[(i-1)*8+4],vars=(1,2),label="Ws-",linecolor=:cyan)
    plot!(W[(i-1)*8+5],vars=(1,2),label="Wu+",linecolor=:red)
    plot!(W[(i-1)*8+6],vars=(1,2),label="Wu-",linecolor=:magenta)
    plot!(W[(i-1)*8+7],vars=(1,2),label="Ws+",linecolor=:blue)
    plot!(W[(i-1)*8+8],vars=(1,2),label="Ws-",linecolor=:cyan)
    plot!(aspect_ratio=1,legend=:outerright,flip=false)
end

plot(sys)
plot!(W[17],vars=(1,2),label="Wu+",linecolor=:red)
plot!(W[18],vars=(1,2),label="Wu-",linecolor=:magenta)
plot!(W[19],vars=(1,2),label="Ws+",linecolor=:blue)
plot!(W[20],vars=(1,2),label="Ws-",linecolor=:cyan)
plot!(aspect_ratio=1,legend=:outerright,flip=false)


text:
# TODO: Split this up into multiple files.

module MatroidChowRings

export direct_sum_decomp, matroid_chow_ring, apply_homorphism
export augmented_direct_sum_decomp, augmented_matroid_chow_ring

using Oscar;
const pm = Polymake;

"""
    struct MChowRing

Represents a matroid Chow ring. The result of a matroid_chow_ring call.

Fields
- chow_ring:      The Chow ring of the matroid, represented as an affine algebra.
                  It may be == nothing, in that case chow ring is QQ.
- indeterminates: The indeterminates of the chow ring. These are of the form
                  x_F, where F is a nonempty proper flat. The name of F
                  corresponds to elements inside that flat. E.g. x__1_5_7 would
                  be the flat consisting of the elements 1, 5, and 7. Another
                  way to see all the flats would be to call
                  matroid.LATTICE_OF_FLATS.FACES.
                  This vector may be empty, in that case the chow ring is trivial.
- matroid:        The matroid whose chow ring has been computed. It is
                  represented as a polymake Matroid.
"""
struct MChowRing
    indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}
    matroid::pm.BigObject
    chow_ring::MPolyQuo{MPolyElem_dec{fmpq, fmpq_mpoly}}

    function MChowRing(indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, matroid::pm.BigObject)
        if pm.type_name(matroid) != "Matroid"
            error("BigObject is not a matroid")
        end

        new(indeterminates, matroid)
    end

    function MChowRing(indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, matroid::pm.BigObject, chow_ring)
        if pm.type_name(matroid) != "Matroid"
            error("BigObject is not a matroid")
        end
        new(indeterminates, matroid, chow_ring)
    end
end

"""
    struct MAugChowRing

Represents am augmented matroid Chow ring. The result of a augmented_matroid_chow_ring call.

Fields
- chow_ring:              The augmented Chow ring of the matroid, represented as an affine algebra.
                          It may be == nothing, in that case chow ring is QQ.
- flat_indeterminates:    The flat indeterminates of the augmented Chow ring. These are of the form
                          x_F, where F is a proper flat. The name of F corresponds to elements inside
                           that flat. E.g. x__1_5_7 would be the flat consisting of the
                          elements 1, 5, and 7. Another way to see all the flats would be to call
                          matroid.LATTICE_OF_FLATS.FACES.
                          This vector may be empty, in that case the chow ring is QQ.
- element_indeterminates: The indeterminates of the augmented Chow ring. These are of the form y_i,
                          where i is an element of the groundset of the matroid. The number of these
                          corresponds to matroid.N_ELEMENTS, ranging from 0 to N_ELEMENTS-1.
- matroid:                The matroid whose augmented Chow ring has been computed. It is
                          represented as a polymake Matroid.
"""
struct MAugChowRing
    flat_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}
    element_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}
    matroid::pm.BigObject
    chow_ring::MPolyQuo{MPolyElem_dec{fmpq, fmpq_mpoly}}

    function MAugChowRing(flat_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, element_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, matroid::pm.BigObject)
        if pm.type_name(matroid) != "Matroid"
            error("BigObject is not a matroid")
        end

        new(flat_indeterminates, element_indeterminates, matroid)
    end

    function MAugChowRing(flat_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, element_indeterminates::Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}, matroid::pm.BigObject, chow_ring)
        if pm.type_name(matroid) != "Matroid"
            error("BigObject is not a matroid")
        end
        new(flat_indeterminates, element_indeterminates, matroid, chow_ring)
    end
end

"""
    struct MChowRingHomorphism

Represents a homorphism which sends the RHS of a Chow ring direct sum
decomposition to the LHS of the decompositon.

Note that the right hand side of the decomposition is not represented as
subalgebras of CH(M), but as Chow rings which are isomorphic to the subalgebras
of CH(M). If you wish to turn an element on the RHS to one which lives in the
domain of the LHS, you need to use this structure.

This structure represents the morphism:

first_term_morphism + coloop_term_morphism + sum_((projection_1, projection_2, term)) term * projection_1 * projection 2

Where g_f,h_f are projections from CH(M_F+1) and CH(M^F), respectively, to CH(M-i).
The input to g_f and h_f would then be values in the corresponding tuple in MChowRingDecomp.

Fields:
- first_term_morphism:  Sends elements of the first_term into the LHS.
- second_term_morphism: A vector of tuples of morphisms and a term x_(F+i), written
                        as (projection_1, projection_2, x_(F+i)). Each tuple represents
                        the psi described in Propositon 2.21 of Tom Braden et. al.
                        It is used in as isomorphism in the direct sum decompsotion, see
                        Propositon 3.5. It can be seen as the product of two projections, a theta_i
                        along with a constant term. I.e. as term * projection_1 * projection_2.
                        Here projection_2 already contains the theta_i.
- coloop_term_morphism: An optional morphism, that sends members of the
                        coloopterm into the LHS. Equal to X_(E-i) * theta_i
"""
struct MChowRingHomorphism
    first_term_morphism
    second_term_morphism
    coloop_term_morphism
end

"""
    struct MAugChowRingHomorphism

Represents a homorphism which sends the RHS of an augmentedChow ring direct sum
decomposition to the LHS of the decompositon.

Note that the right hand side of the decomposition is not represented as
subalgebras of CH_aug(M), but as (augmented) Chow rings which are isomorphic to the
subalgebras of CH_aug(M). If you wish to turn an element on the RHS to one
which lives in the domain of the LHS, you need to use this structure.

This structure represents the morphism:

first_term_morphism + coloop_term_morphism + sum_((projection_1, projection_2, term)) term * projection_1 * projection 2

Where g_f,h_f are projections from CH(M_F+1) and CH_aug(M^F), respectively, to CH_aug(M-i).
The input to g_f and h_f would then be values in the corresponding tuple in MAugChowRingDecomp.

Fields:
- first_term_morphism:  Sends elements of the first_term into the LHS.
- second_term_morphism: A vector of tuples of morphisms and a term x_(F+i), written
                        as (projection_1, projection_2, x_(F+i)). Each tuple represents
                        the psi described in Propositon 2.24 of Tom Braden et. al.
                        It is used in as isomorphism in the direct sum decompsotion, see
                        Propositon 3.5. It can be seen as the product of two projections, a theta_i
                        along with a constant term. I.e. as term * projection_1 * projection_2.
                        Here projection_2 already contains the theta_i.
- coloop_term_morphism: An optional morphism, that sends members of the
                        coloopterm into the LHS. Equal to X_(E-i) * theta_i
"""
struct MAugChowRingHomorphism
    first_term_morphism
    second_term_morphism
    coloop_term_morphism
end

"""
    struct MChowRingDecomp
Represents a direct sum decomposition of the Chow ring of a matroid. See
"A semi-small decomposition of the Chow ring of a matroid by Braden et. al" for
details. The result of a direct_sum_decomp call.

Note that the right hand side of the decomposition is not represented as
subalgebras of CH(M), but as Chow rings which are isomorphic to the subalgebras
of CH(M). If you wish to turn an element on the RHS to one which lives in the
domain of the LHS, you need to use the homomorphism field.

Fields:
- deleted_element: The element which was deleted on the RHS of the decomposition.
- chow_ring_LHS:   The chow ring which was decomposed.
- first_term:      Always present part of the decomposition, equal to the chow
                   ring of M-i.
- second_term:     A vector of tuples of chow rings. This has varying length,
                   and each tuple is of the form (CH(M_(F+i)), CH(M^F), where
                   i is the deleted element, and F is a flat. This is the part
                   of the decomposition indexed by the S_i in the paper.
- coloopterm:      A term that appears in the decomposition if deleted_element
                   is a coloop. Equal to the chow ring of M-i if present.
- homomorphism:    Represents the homomorphism which takes a combination of
                   elements from the first/second_term/coloop_term and maps it
                   to the LHS of the decomposition.
"""
struct MChowRingDecomp #TODO Add types here.
    deleted_element::Int64
    chow_ring_LHS::MChowRing
    first_term::MChowRing
    second_term
    coloopterm
    homomorphism::MChowRingHomorphism
end

"""
    struct MAugChowRingDecomp
Represents a direct sum decomposition of the augmented Chow ring of a matroid. See
"A semi-small decomposition of the Chow ring of a matroid by Braden et. al" for
details. The result of a direct_sum_decomp call.

Note that the right hand side of the decomposition is not represented as
subalgebras of CH_aug(M), but as Chow rings which are isomorphic to the subalgebras
of CH_aug(M). If you wish to turn an element on the RHS to one which lives in the
domain of the LHS, you need to use the homomorphism field.

Fields:
- deleted_element: The element which was deleted on the RHS of the decomposition.
- chow_ring_LHS:   The augmented Chow ring which was decomposed.
- first_term:      Always present part of the decomposition, equal to the augmented Chow
                   ring of M-i.
- second_term:     A vector of tuples of chow rings. This has varying length,
                   and each tuple is of the form (CH(M_(F+i)), CH_aug(M^F), where
                   i is the deleted element, and F is a flat. This is the part
                   of the decomposition indexed by the S_i in the paper.
- coloopterm:      A term that appears in the decomposition if deleted_element
                   is a coloop. Equal to the chow ring of M-i if present.
- homomorphism:    Represents the homomorphism which takes a combination of
                   elements from the first/second_term/coloop_term and maps it
                   to the LHS of the decomposition.
"""
struct MAugChowRingDecomp #TODO Add types here.
    deleted_element::Int64
    chow_ring_LHS::MAugChowRing
    first_term::MAugChowRing
    second_term
    coloopterm
    homomorphism::MAugChowRingHomorphism
end

"""
    function direct_sum_decomp(matroid::pm.BigObject, matroid_element::Int64)

Convinience function to compute the direct sum decomposition of the Chow ring of a matroid without
first explicitely computing the chow ring itself. See documentation of the same function applied
to a chow ring of type MChowRing for more details.

The matroid may be any loopless matroid.
Matroid element is a number from 0 to matroid.NR_ELEMENTS-1, inclusive.
The matroid element is the one which is removed in the decomposition.

# Example
```julia-repl
julia> using Oscar
julia> using MatroidChowRings
julia> f8 = Polymake.matroid.f8_matroid();
julia> decomp = direct_sum_decomp(f8, 3);
julia> length(gens(decomp.chow_ring_LHS.chow_ring))
56
julia> length(gens(decomp.first_term.chow_ring))
45
julia> length(decomp.second_term)
10
julia> decomp.coloopterm == nothing
true
```

Note that the length of the gens corresponds to the number of nonempty proper flats.
The 56 and 45 above say that the matroid corresponding to f8 delete 3 has 13 less flats.

Each element of decomp.second_term corresponds to a flat F, such that F+3 is also a flat.
Thus the length 10 says that there are 10 such choices for F.

f8 has no coloops, so the coloopterm must be empty.
"""
function direct_sum_decomp(matroid::pm.BigObject, matroid_element::Int64)
    direct_sum_decomp(matroid_chow_ring(matroid), matroid_element)
end

# TODO: Implement an augmented version of this.
"""
    function direct_sum_decomp(chow_ring::MChowRing, matroid_element::Int64)

Computes the (semi-small) direct sum decomposition of the Chow ring of a matroid, as first described
in 'A SEMI-SMALL DECOMPOSITION OF THE CHOW RING OF A MATROID' by Tom Braden, June Huh et. al.

Due to technical reasons the decomposition is not represented as
subalgebras of CH(M), but as Chow rings which are isomorphic to the subalgebras
of CH(M). Meaning that instead of a direct sum decomposition, this function instead
returns all the (isomorphic) components of the decomposition.

If you wish to instead recover the associated subalgebras, the output of this function also
contains the relevant pushforward map, which takes an element from the product of all the
isomoprhic components of the decomposition, and maps it to an element of the matroid chow ring
in the image. By then mapping a set of generators in the isomorphic component, you can get the
corresponding generators of the subalgebra. See the function apply_homorphism for details.

Note that because polymake reindexes matroid elements to always run from 0 to n-1, an indeterminate
of the form e.g x_1_3_5 need not correspond to x_1_3_5 in the original matroid. In particular, if
3 was removed in the construction of the smaller matroid. The function apply_homorphism accounts for
the reindexing.

chow_ring         the result of a direct_sum_decomp call
matroid_element   the matroid element to be removed in the decomposition.

# Example

```julia-repl
julia> using Oscar
julia> using MatroidChowRings
julia> fano_matroid = Polymake.matroid.fano_matroid();
julia> chow_ring =  matroid_chow_ring(fano_matroid);
julia> x = chow_ring.indeterminates;
julia> decomp = direct_sum_decomp(chow_ring, 0);
julia> length(decomp.first_term.indeterminates)
13
julia> length(decomp.second_term)
0
julia> decomp.coloopterm == nothing
true
julia> morphism = decomp.homomorphism.first_term_morphism;
julia> isbijective(morphism)
true
julia> morphism.image
13-element Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}:
 x__6 - x__1_2_3 - x__1_4_5 + x__2_4_6 + x__3_5_6
 x__6 - x__0_2_5 + x__0_1_6 - x__1_2_3 + x__3_5_6
 x__6 - x__0_3_4 + x__0_1_6 - x__1_2_3 + x__2_4_6
 x__6 - x__0_3_4 + x__0_1_6 - x__1_4_5 + x__3_5_6
 x__6 - x__0_2_5 + x__0_1_6 - x__1_4_5 + x__2_4_6
 x__6
 x__1_2_3
 x__1_4_5
 x__0_1_6
 x__0_2_5
 x__2_4_6
 x__0_3_4
 x__3_5_6
 julia> for img_gen in morphism.image
           if img_gen in x
               print(true,"\n")
           else
               print(false,"\n")
           end
       end
true
true
true
true
true
true
true
true
true
true
true
true
true
```

The above code says that the direct sum decomposition of the fano matroid consists of only
one component which is isomorphic to matroid you get by deleting any element from the fano
matroid. This means that the pushforward map is very simple - it is composed of a single
algebraic homomorphism. One can directly verify that this is indeed a decomposition by
checking that this algerbraic homoprhism is isomorphic.

The for loop proves that generators in the RHS of the decomposition get mapped to generators
in the Chow ring of the Fano matroid. This paired with the fact that the morphism is isomorphic,
shows that the natural set of generators for the Fano matroid contain redudant elements, in particular
there is a generating set of length 13.
"""
function direct_sum_decomp(chow_ring::MChowRing, matroid_element::Int64)
    if matroid_element < 0 || matroid_element >= chow_ring.matroid.N_ELEMENTS
        print("Matroid element out of bounds, can not decompose.\n")
        return nothing
    end
    matroid = chow_ring.matroid

    first_term = matroid_chow_ring(pm.matroid.deletion(matroid, matroid_element))
    first_term_morphism = create_theta_i(first_term, chow_ring, matroid_element)

    ground_set = Polymake.Set(range(0, matroid.N_ELEMENTS - 1, step=1))
    flats = matroid.LATTICE_OF_FLATS.FACES
    flats = pm.@pm common.convert_to{IncidenceMatrix}(flats)
    proper_flats = flats[2:pm.size(flats, 1)-1, 1:pm.size(flats,2)]
    n_proper_flats = pm.size(proper_flats, 1)

    flat_canidates = []
    for flat_index in true_indices_in_col(proper_flats, matroid_element + 1)
        proper_flat = extract_flat(proper_flats, flat_index)
        push!(flat_canidates, proper_flat)
    end
    flat_canidates = Set(flat_canidates)
    second_term = []
    second_term_morphism = []

    for flat_index in false_indices_in_col(proper_flats, matroid_element + 1)
        proper_flat = extract_flat(proper_flats, flat_index)
        proper_flat_copy = copy(proper_flat)
        push!(proper_flat_copy, matroid_element)
        if proper_flat_copy in flat_canidates
            to_remove_1 = proper_flat_copy
            to_remove_2 = set_complement(ground_set, proper_flat)

            matroid_1 = pm.matroid.contraction(matroid, to_remove_1)
            matroid_2 = pm.matroid.deletion(matroid, to_remove_2)

            chow_1 = matroid_chow_ring(matroid_1)
            chow_2 = matroid_chow_ring(matroid_2)

            push!(second_term, (chow_1, chow_2))

            projection_1 = create_projection(chow_1, chow_ring, to_remove_1, matroid_element, true)
            projection_2 = create_projection(chow_2, chow_ring, to_remove_2, matroid_element, false)
            term = find_flat_variable(chow_ring, proper_flat_copy)
            push!(second_term_morphism, (projection_1, projection_2, term))
        end
    end

    coloops = pm.matroid.coloops(matroid)
    is_coloop = pm.in(matroid_element, coloops)
    coloop_term = nothing
    coloop_term_morphism = nothing
    if is_coloop
        coloop_term = first_term
        flat = copy(ground_set)
        delete!(flat, matroid_element)
        term = find_flat_variable(chow_ring, flat)
        coloop_term_morphism = create_theta_i(coloop_term, chow_ring, matroid_element, term)
    end

    chow_hom = MChowRingHomorphism(first_term_morphism, second_term_morphism, coloop_term_morphism)

    MChowRingDecomp(matroid_element, chow_ring, first_term, second_term, coloop_term, chow_hom)

end

function augmented_direct_sum_decomp(matroid::pm.BigObject, matroid_element::Int64)
    augmented_direct_sum_decomp(augmented_matroid_chow_ring(matroid), matroid_element)
end

function augmented_direct_sum_decomp(chow_ring::MAugChowRing, matroid_element::Int64)
    if matroid_element < 0 || matroid_element >= chow_ring.matroid.N_ELEMENTS
        print("Matroid element out of bounds, can not decompose.\n")
        return nothing
    end

    matroid = chow_ring.matroid
    first_term = augmented_matroid_chow_ring(pm.matroid.deletion(matroid, matroid_element))
    first_term_morphism = create_aug_theta_i(first_term, chow_ring, matroid_element)

    ground_set = Polymake.Set(range(0, matroid.N_ELEMENTS - 1, step=1))
    flats = matroid.LATTICE_OF_FLATS.FACES
    flats = pm.@pm common.convert_to{IncidenceMatrix}(flats)

    top_node = matroid.LATTICE_OF_FLATS.TOP_NODE
    if top_node == 0
        proper_flats = flats[2:pm.size(flats, 1), 1:pm.size(flats,2)]
    else
        proper_flats = flats[1:pm.size(flats, 1)-1, 1:pm.size(flats,2)]
    end

    flat_canidates = []
    for flat_index in true_indices_in_col(proper_flats, matroid_element + 1)
        proper_flat = extract_flat(proper_flats, flat_index)
        push!(flat_canidates, proper_flat)
    end

    second_term = []
    second_term_morphism = []

    for flat_index in false_indices_in_col(proper_flats, matroid_element + 1)
        proper_flat = extract_flat(proper_flats, flat_index)
        proper_flat_copy = copy(proper_flat)
        push!(proper_flat_copy, matroid_element)
        if proper_flat_copy in flat_canidates
            to_remove_1 = proper_flat_copy
            to_remove_2 = set_complement(ground_set, proper_flat)

            matroid_1 = pm.matroid.contraction(matroid, to_remove_1)
            matroid_2 = pm.matroid.deletion(matroid, to_remove_2)

            chow_1 = matroid_chow_ring(matroid_1)
            chow_2 = augmented_matroid_chow_ring(matroid_2)

            push!(second_term, (chow_1, chow_2))

            projection_1 = create_projection(chow_1, chow_ring, to_remove_1, matroid_element, true)
            projection_2 = create_projection(chow_2, chow_ring, to_remove_2, matroid_element, false)
            term = find_flat_variable(chow_ring, proper_flat_copy)
            push!(second_term_morphism, (projection_1, projection_2, term))
        end
    end

    coloops = pm.matroid.coloops(matroid)
    is_coloop = pm.in(matroid_element, coloops)
    coloop_term = nothing
    coloop_term_morphism = nothing
    if is_coloop
        coloop_term = first_term
        flat = copy(ground_set)
        delete!(flat, matroid_element)
        term = find_flat_variable(chow_ring, flat)
        coloop_term_morphism = create_aug_theta_i(coloop_term, chow_ring, matroid_element, term)
    end

    chow_hom = MAugChowRingHomorphism(first_term_morphism, second_term_morphism, coloop_term_morphism)

    MAugChowRingDecomp(matroid_element, chow_ring, first_term, second_term, coloop_term, chow_hom)
end

"""
Creates theta_i as described on page 2 of
"A semi-small decomposition of the Chow ring of a matroid" by Braden et. al,
times a given factor in the image. E.g. x__1_3_5 * theta_i

These modified theta_i are the building blocks of the isomoprhisms of terms
inside the direct sum decomposition of a Chow ring with Chow Rings of smaller
matroids.
"""
function create_theta_i(domain::MChowRing, image::MChowRing, deleted_element::Int64, factor=1)
    image_gens = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}(undef, length(gens(domain.chow_ring)))
    domain_gens = gens(domain.chow_ring)

    for i = 1:length(gens(domain.chow_ring))
        domain_gen = domain_gens[i]
        image_gens[i] = factor * find_gen_image(domain_gen, image, deleted_element)
    end

    hom(domain.chow_ring, image.chow_ring, image_gens)
end

function create_aug_theta_i(domain::MAugChowRing, image::MAugChowRing, deleted_element::Int64, factor=1)
    image_gens = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}(undef, length(gens(domain.chow_ring)))
    domain_ys = domain.element_indeterminates
    domain_xs = domain.flat_indeterminates
    image_ys = image.element_indeterminates

    for i = 1:length(domain_ys)
        image_index = i
        if i > deleted_element
            image_index += 1
        end
        image_gens[i] = factor * image_ys[image_index]
    end

    offset = length(domain_ys)
    for i = 1:length(domain_xs)
        domain_gen = domain_xs[i]
        image_gens[i + offset] = factor * find_gen_image(domain_gen, image, deleted_element)
    end

    hom(domain.chow_ring, image.chow_ring, image_gens)
end

# TODO: This works... but it's ugly string hacking.
# One should instead pass around a dict of set ints to variables instead.
"""
Helper function for create_theta_i. Given e.g. domain_gen = x_1_3_5_7 and
deleted element = 5, it creates the polynomial

x_1_3_6_8 + x_1_3_5_6_8 in the image, where a term is set to 0 if the
corresponding variable does not exist in the image.
"""
function find_gen_image(domain_gen, image, deleted_element::Int64)
     elements_in_flat = split(string(domain_gen), "__")[2]

     first_canidate = "x_"
     second_canidate = "x_"

     previous_element = nothing
     second_canidate_adjusted = false

     for element_string in split(elements_in_flat, "_")
         if element_string == ""
             break
         end

         element = parse(Int, element_string)
         if element >= deleted_element
             element += 1
         end
         if element > deleted_element &&
            (previous_element == nothing || previous_element < deleted_element)
              second_canidate = second_canidate * "_" * string(deleted_element)
              second_canidate_adjusted = true
         end

         first_canidate = first_canidate * "_" * string(element)
         second_canidate = second_canidate * "_" * string(element)
         previous_element = element
     end

     if !second_canidate_adjusted
         second_canidate = second_canidate * "_" * string(deleted_element)
     end
     first_term = image.chow_ring()
     second_term = image.chow_ring()
     for image_gen in gens(image.chow_ring)
         image_gen_string = string(image_gen)
         if image_gen_string == first_canidate
             first_term = image_gen
         end
         if image_gen_string == second_canidate
             second_term = image_gen
         end
     end

     first_term + second_term
end

"""
Creates a projection function from a chow ring to another chow ring, where the
corresponding domain matroid is a contraction or deletion of the codomain matroid.

This is a building block of the pushforward map psi described in Proposition 2.18
and used in Proposition 3.5 to describe the isomorphisms behind the Chow ring
decomposition.

For contraction: The removed elements correspond to a flat F+i. In this case,
X_(F') gets mapped to X_(F+i union F')

For deletion: The removed elements correspond to the complement of a flat F.
There still exists a corresponding i. In this case, X_(F') gets mapped to
X_(F') + X_(F'+i).
"""
function create_projection(domain::MChowRing, image::MChowRing, removed_elements, i::Int64, is_contraction::Bool)
    if !isdefined(domain, :chow_ring) # Edge case where ring is trivial.
        return nothing
    end

    image_gens = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}(undef, length(gens(domain.chow_ring)))
    domain_gens = gens(domain.chow_ring)
    int_mapping = create_int_mapping(domain.matroid.N_ELEMENTS, removed_elements)

    for j = 1:length(gens(domain.chow_ring))
        domain_gen = domain_gens[j]
        image_gens[j] = project_gen(domain_gen, image, removed_elements, i, is_contraction, int_mapping)
    end

    hom(domain.chow_ring, image.chow_ring, image_gens)
end

function create_projection(domain::MAugChowRing, image::MAugChowRing, removed_elements, i::Int64, is_contraction::Bool)
    if !isdefined(domain, :chow_ring) # Edge case where ring is trivial.
        return nothing
    end

    image_gens = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}(undef, length(gens(domain.chow_ring)))
    domain_gens = gens(domain.chow_ring)
    int_mapping = create_int_mapping(domain.matroid.N_ELEMENTS, removed_elements)

    for j = 1:length(domain.element_indeterminates)
        domain_gen = domain.element_indeterminates[j]
        image_gens[j] = image.element_indeterminates[int_mapping[j] + 1]
    end

    offset = length(domain.element_indeterminates)
    for j = 1:length(domain.flat_indeterminates)
        domain_gen = domain.flat_indeterminates[j]
        image_gens[j + offset] = project_gen(domain_gen, image, removed_elements, i, is_contraction, int_mapping)
    end

    hom(domain.chow_ring, image.chow_ring, image_gens)
end

function create_projection(domain::MChowRing, image::MAugChowRing, removed_elements, i::Int64, is_contraction::Bool)
    if !isdefined(domain, :chow_ring) # Edge case where ring is trivial.
        return nothing
    end

    image_gens = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}(undef, length(gens(domain.chow_ring)))
    domain_gens = gens(domain.chow_ring)
    int_mapping = create_int_mapping(domain.matroid.N_ELEMENTS, removed_elements)

    for j = 1:length(gens(domain.chow_ring))
        domain_gen = domain_gens[j]
        image_gens[j] = project_gen(domain_gen, image, removed_elements, i, is_contraction, int_mapping)
    end

    hom(domain.chow_ring, image.chow_ring, image_gens)
end

function create_int_mapping(n_elements::Int64, removed_elements)
    # Relates matroid elements of the domain with the image
    # This is the reindexing of contraction/deletion in reverse.
    int_mapping = Array{Int64}(undef, n_elements)
    next_free_element = 0
    while Polymake.in(next_free_element, removed_elements)
        next_free_element += 1
    end
    for j = 1:length(int_mapping)
        int_mapping[j] = next_free_element
        next_free_element += 1
        while Polymake.in(next_free_element, removed_elements)
            next_free_element += 1
        end
    end

    int_mapping
end

function project_gen(domain_gen::MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}, image, removed_elements, i::Int64, is_contraction::Bool, int_mapping)
    elements_in_gen = split(split(string(domain_gen), "__")[2], "_")
    if elements_in_gen[1] == ""
        elements_in_gen = []
    end

    elements_in_gen = [parse(Int, e) for e in elements_in_gen]

    # +1 here because matroid elements start from 0, and arrays start from 1.
    contents_image = [int_mapping[e+1] for e in elements_in_gen]

    if is_contraction
        contents_image = vcat(contents_image, Vector(removed_elements))
        sort!(contents_image)
        return find_flat_variable(image, contents_image)
    else
        x1 = find_flat_variable(image, contents_image)
        append!(contents_image, i)
        sort!(contents_image)
        x2 = find_flat_variable(image, contents_image)
        # Due to the structure of S_i, x1 and x2 should always be found.
        return x1 + x2
    end

end

# TODO: Find some way to remove this ugly string hacking
"""
Finds the generator in chow_ring with exactly the elements as in contents.
"""
function find_flat_variable(chow_ring, contents)
    for gen in gens(chow_ring.chow_ring)
        if split(string(gen), "__")[1] != "x"
            continue
        end

        elements_in_flat = split(split(string(gen), "__")[2], "_")
        if elements_in_flat[1] == ""
            elements_in_flat = []
        end

        if length(elements_in_flat) != length(contents)
            continue
        end

        sets_equal = true
        for element_string in elements_in_flat
            element = parse(Int, element_string)
            if !Polymake.in(element, contents)
                sets_equal = false
                break
            end
        end

        if sets_equal
            return gen
        end
    end
end

# TODO: Add types here. Test it. Add example.
"""
This applies the MChowRingHomorphism, sending an element of the RHS of a
direct sum decomposition into the LHS.

Note that the right hand side of the decomposition is not represented as
subalgebras of CH(M), but as Chow rings which are isomorphic to the subalgebras
of CH(M). If you wish to turn an element on the RHS to one which lives in the
domain of the LHS, you need to use this structure.

So using this function you can instead recover an element which truly does
belong to the original CH(M).

chow_ring_hom must be the homorphism which comes as one of the return values
of direct_sum_decomp.

The first_term field must be an element inside the first_term of MChowRingDecomp.
The second_term must be a vector of 2-tuples, containing elements of the
second_term of the corresponding second_term of a MChowRingDecomp, in the same
order as they are in MChowRingDecomp.

# Example

```julia-repl
julia> using Oscar
julia> using MatroidChowRings
julia> u45 = Polymake.matroid.uniform_matroid(4,5);
julia> decomp = direct_sum_decomp(u45, 2);
julia> length(decomp.second_term)
10
julia> decomp.coloopterm == nothing
true
julia> p1 = decomp.first_term.indeterminates[2] + decomp.first_term.indeterminates[1]
x__0 + x__1
julia> p2s = [];
julia> for (t1, t2) in decomp.second_term
           x1 = 1
           x2 = 1
           if length(t1.indeterminates) > 0
               x1 = t1.indeterminates[1]
           end
           if length(t2.indeterminates) > 0
               x2 = t2.indeterminates[2]
           end

           push!(p2s, (x1, x2))
       end
julia> p2s
10-element Vector{Any}:
(x__0, 1)
(x__0, 1)
(x__0, 1)
(x__0, 1)
(1, x__1)
(1, x__1)
(1, x__1)
(1, x__1)
(1, x__1)
(1, x__1)
julia> p3 = nothing;
julia> apply_homorphism(decomp.homomorphism, p1, p2s, p3)
2*x__4 - x__0_1*x__0_1_4 - 2*x__0_1 + x__0_2*x__0_2_4 - x__0_3*x__0_3_4 - x__0_3 - x__0_4*x__0_3_4 + x__0_4 + x__1_2*x__1_2_4 - x__1_3*x__1_3_4 - x__1_3 - x__1_4*x__1_3_4 + x__1_4 + x__2_3*x__2_3_4 + x__2_4*x__2_3_4 + 2*x__2_4 - x__3_4*x__2_3_4 + 2*x__3_4 - x__0_1_2^2 - 2*x__0_1_2 - 2*x__0_1_3 - x__0_2_3^2 - x__0_2_3 - x__0_2_4^2 + x__0_2_4 + x__0_3_4 - x__1_2_3^2 - x__1_2_3 - x__1_2_4^2 + x__1_2_4 + x__1_3_4 - x__2_3_4^2 + 2*x__2_3_4
```

"""
function apply_homorphism(chow_ring_hom::MChowRingHomorphism, first_term, second_term, coloop_term)
    result = chow_ring_hom.first_term_morphism(first_term)

    i = 1
    for (contracted_term, deleted_term) in second_term
        projection_1, projection_2, term = chow_ring_hom.second_term_morphism[i]

        mapped_contracted_term = contracted_term
        mapped_deleted_term = deleted_term
        if projection_1 != nothing
            mapped_contracted_term = projection_1(contracted_term)
        end
        if projection_2 != nothing
            mapped_deleted_term = projection_2(deleted_term)
        end

        result +=  term * mapped_contracted_term * mapped_deleted_term
        i += 1
    end

    if coloop_term != nothing
        result += chow_ring_hom.coloop_term_morphism(coloop_term)
    end

    result
end

"""
    matroid_chow_ring(matroid::pm.BigObject)::MChowRing

Computes the Chow ring of a matroid. It follows the convention of the chow ring as defined
in 'A SEMI-SMALL DECOMPOSITION OF THE CHOW RING OF A MATROID' by Tom Braden, June Huh et. al.
The chow ring is represented over the rationals, and only nonempty proper flats enter as variables.

This function accepts any loopless matroid as input, where the matroid is a polymake matroid.
It returns the Chow ring described as a quitotent ring, and it does not precompute more than
it needs to. In particular, if you do any calculations on the ring, then you will likely need
to wait a bit the first time as it will likely compute a Gröbner basis

# Example
```julia-repl
julia> using Oscar
julia> using MatroidChowRings
julia> fano_matroid = Polymake.matroid.fano_matroid();
julia> chow_ring =  matroid_chow_ring(fano_matroid);
julia> x = chow_ring.indeterminates
14-element Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}:
 x__0
 x__1
 x__2
 x__3
 x__4
 x__5
 x__6
 x__0_3_4
 x__0_2_5
 x__0_1_6
 x__1_2_3
 x__1_4_5
 x__2_4_6
 x__3_5_6
 julia> x[1]
 x__0
 julia> x[8]
 x__0_3_4
 julia> x[11]
 x__1_2_3
 julia> x[1] * x[8]
 -x__3_5_6^2
 julia> x[1] * x[11]
 0
```

"""
function matroid_chow_ring(matroid::pm.BigObject)::MChowRing
    if pm.type_name(matroid) != "Matroid"
        throw(ArgumentError("BigObject is not a matroid."))
    end

    flats = matroid.LATTICE_OF_FLATS.FACES;
    #TODO: Is there a way to avoid this conversion?
    flats = pm.@pm common.convert_to{IncidenceMatrix}(flats)

    # Technically non-empty proper flats. First element is either empty or the whole set, last element is the other of the two.
    proper_flats = flats[2:pm.size(flats, 1)-1, 1:pm.size(flats,2)]

    if length(proper_flats) == 0
        empty_vector = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}()
        return MChowRing(empty_vector, matroid)
    end

    base_ring, indeterminates = generate_base_ring(proper_flats)

    type_i_polynomials = generate_type_i_ideal(base_ring, proper_flats, indeterminates, matroid)
    type_j_polynomials = generate_type_j_ideal(proper_flats, indeterminates)
    generators = vcat(type_i_polynomials, type_j_polynomials)

    chow_modulus = ideal(base_ring, generators)
    # Add std(ideal) to compute standard basis?
    chow_ring, projection = quo(base_ring, chow_modulus)
    projected_indeterminates = [projection(indeterminate) for indeterminate in indeterminates]

    MChowRing(projected_indeterminates, matroid, chow_ring)
end

"""
    augmented_matroid_chow_ring(matroid::pm.BigObject)::MAugChowRing

Computes the augmented Chow ring of a matroid. It follows the convention of the chow ring as defined
in 'A SEMI-SMALL DECOMPOSITION OF THE CHOW RING OF A MATROID' by Tom Braden, June Huh et. al.
The chowVector ring is represented over the rationals, and proper flats and matroid groundset elements
enter as variables.

This function accepts any loopless matroid as input, where the matroid is a polymake matroid.
It returns the Chow ring described as a quitotent ring, and it does not precompute more than
it needs to. In particular, if you do any calculations on the ring, then you will likely need
to wait a bit the first time as it will likely compute a Gröbner basis.

# Example
```julia-repl
julia> using Oscar
julia> using MatroidChowRings
julia> fano_matroid = Polymake.matroid.fano_matroid();
julia> chow_ring = augmented_matroid_chow_ring(fano_matroid);
julia> x = chow_ring.flat_indeterminates
15-element Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}:
 x__
 x__0
 x__1
 x__2
 x__3
 x__4
 x__5
 x__6
 x__0_3_4
 x__0_2_5
 x__0_1_6
 x__1_2_3
 x__1_4_5
 x__2_4_6
 x__3_5_6

julia> y = chow_ring.element_indeterminates
7-element Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}:
 y_0
 y_1
 y_2
 y_3
 y_4
 y_5
 y_6
julia> y[1] * y[1]
0
julia> y[1] * y[2]
-x__6^2 - 4*x__6*x__3_5_6 - x__2_4_6^2 - x__3_5_6^2
julia> x[9]
x__0_3_4
julia> y[1] * x[9]
-x__0*x__0_1_6 - x__0_3_4^2
julia> y[2] * x[9]
0
```

"""
function augmented_matroid_chow_ring(matroid::pm.BigObject)::MAugChowRing
    if pm.type_name(matroid) != "Matroid"
        throw(ArgumentError("BigObject is not a matroid."))
    end

    if matroid.N_ELEMENTS == 0
        empty_vector = Vector{MPolyQuoElem{MPolyElem_dec{fmpq, fmpq_mpoly}}}()
        return MAugChowRing(empty_vector, empty_vector, matroid)
    end

    n_elements = matroid.N_ELEMENTS
    flats = matroid.LATTICE_OF_FLATS.FACES
    flats = pm.@pm common.convert_to{IncidenceMatrix}(flats)

    top_node = matroid.LATTICE_OF_FLATS.TOP_NODE
    if top_node == 0
        proper_flats = flats[2:pm.size(flats, 1), 1:pm.size(flats,2)]
    else
        proper_flats = flats[1:pm.size(flats, 1)-1, 1:pm.size(flats,2)]
    end

    base_ring, matroid_element_vars, flat_vars = generate_augmented_base_ring(proper_flats, n_elements)
    type_i_polynomials = generate_augmented_type_i_ideal(proper_flats, matroid_element_vars, flat_vars)
    type_j_polynomials = generate_augmented_type_j_ideal(proper_flats, matroid_element_vars, flat_vars)
    generators = vcat(type_i_polynomials, type_j_polynomials)

    chow_modulus = ideal(base_ring, generators)
    chow_ring, projection = quo(base_ring, chow_modulus)

    projected_element_vars = [projection(matroid_element_var) for matroid_element_var in matroid_element_vars]
    projected_flat_vars = [projection(flat_var) for flat_var in flat_vars]

    MAugChowRing(projected_flat_vars, projected_element_vars, matroid, chow_ring)
end

function generate_base_ring(proper_flats)
    variable_names = create_flat_variables_names(proper_flats)
    GradedPolynomialRing(QQ, variable_names);
end

function generate_augmented_base_ring(proper_flats, n_elements)
    matroid_element_variable_names = Array{String}(undef, n_elements)
    for i in 0:n_elements-1
        matroid_element_variable_names[i+1] = "y_" * string(i)
    end

    flat_variable_names = create_flat_variables_names(proper_flats)
    variable_names = vcat(matroid_element_variable_names, flat_variable_names)

    ring, variables = GradedPolynomialRing( QQ, variable_names);

    ring, variables[1:n_elements], variables[n_elements + 1: length(variables)]
end

function create_flat_variables_names(flats)
    n_flats = pm.size(flats, 1)
    variable_names = Array{String}(undef, n_flats)

    for i in 1:n_flats
        row = extract_flat(flats, i)
        # TODO: Is there some better way to extract the polymake string without the type?
        set_name = split(string(extract_flat(flats, i)), "\n")[2]
        variable_descriptor = replace(set_name[2:length(set_name)-1], " " => "_")
        variable_names[i] = "x__" * variable_descriptor
    end

    variable_names
end

function generate_type_i_ideal(base_ring, proper_flats, indeterminates, matroid)
    n = matroid.N_ELEMENTS
    n_proper_flats = pm.size(proper_flats, 1)
    ideal_polynomials = Vector{MPolyElem_dec{fmpq, fmpq_mpoly}}()

    for i in 0:n-1
        for j in i+1:n-1
            ij_polynomial = base_ring() # Create zero.
            ij_set = Set{Int64}([i,j])

            for i_index in true_indices_in_col(proper_flats, i + 1)
                ij_polynomial += indeterminates[i_index]
            end

            for j_index in true_indices_in_col(proper_flats, j + 1)
                ij_polynomial -= indeterminates[j_index]
            end
            push!(ideal_polynomials, ij_polynomial)
        end
    end

    ideal_polynomials
end

function generate_augmented_type_i_ideal(proper_flats, matroid_element_vars, flat_vars)
    n_proper_flats = pm.size(proper_flats, 1)
    ideal_polynomials = Array{MPolyElem_dec{fmpq, fmpq_mpoly}}(undef, length(matroid_element_vars))

    ground_set = Polymake.Set(range(1, length(matroid_element_vars), step=1))

    for i in 1:length(matroid_element_vars)
        ideal_polynomials[i] = matroid_element_vars[i]

        for flat_index_no_i in false_indices_in_col(proper_flats, i)
            ideal_polynomials[i] -= flat_vars[flat_index_no_i]
        end
    end

    ideal_polynomials
end

function generate_type_j_ideal(proper_flats, indeterminates)
    n_proper_flats = pm.size(proper_flats, 1)
    ideal_polynomials = Vector{}()

    for i in 1:n_proper_flats
        for j in i:n_proper_flats
            i_flat = extract_flat(proper_flats, i)
            j_flat = extract_flat(proper_flats, j)
            if are_sets_incomparable(i_flat, j_flat)
                polynomial = indeterminates[i] * indeterminates[j]
                push!(ideal_polynomials, polynomial)
            end
        end
    end

    ideal_polynomials

end

function generate_augmented_type_j_ideal(proper_flats, matroid_element_vars, flat_vars)
    n_proper_flats = pm.size(proper_flats, 1)
    incomparable_polynomials = generate_type_j_ideal(proper_flats, flat_vars)
    xy_polynomials = Vector{}()

    for i in 1:length(matroid_element_vars)
        for j in false_indices_in_col(proper_flats, i)
            push!(xy_polynomials, matroid_element_vars[i] * flat_vars[j])
        end
    end

    vcat(incomparable_polynomials, xy_polynomials)

end

function are_sets_incomparable(a, b)
    return !is_subset(a,b) && !is_subset(b,a)
end

function are_sets_equal(a,b)
    is_subset(a,b) && is_subset(b,a)
end

function is_subset(a, b)
    for element in a
        if !Polymake.in(element, b)
            return false
        end
    end
    true
end

function set_complement(ground_set, to_complement)
    result = copy(ground_set)
    for i in to_complement
        delete!(result, i)
    end
    result
end

function extract_flat(flats, index)
    flat = Polymake.row(flats, index)

    # This is a hacky fix for the conversion to IncidenceMatrix adding 1 to every index.
    result = Array{Int64}(undef, length(flat))

    i = 1
    for element in flat
        result[i] = element - 1
        i += 1
    end

    Polymake.Set(result)
end

function true_indices_in_col(matrix, index)
     Polymake.SparseArrays.nonzeroinds(matrix[:,index])
end

function false_indices_in_col(matrix, index)
     to_complement = true_indices_in_col(matrix, index)
     ground_set = Polymake.Set(range(1, Polymake.nrows(matrix), step=1))
     set_complement(ground_set, to_complement)
end

end

max_stars_repo_name: FGlazov/ChowRings

max_issues_repo_name: FGlazov/ChowRings

max_forks_repo_name: FGlazov/ChowRings


text:
# 2D Advection-Diffusion problem
# LBM
#
# Hongyang Zhou, hyzhou@umich.edu 05/27/2019

# Let us consider flow in a square porous domain.
# Initially the domain is at zero temperature. The upper boundary of the domain
# is kept at zero temperature, while the bottom and vertical right boundaries
# are assumed to be adiabatic. The left vertical wall is subjected to unity
# temperature at time = 0. Find the temperature distribution in the domain at
# time = 200, for material of diffusivity of 1.0. The velocity field of u = 0.1
# and v = 0.2 is kept constant in the domain.

using PyPlot
include("Parameters.jl")
include("lattice.jl")
include("BC.jl")

using .Parameters, .BC, .lattice
using LinearAlgebra

function solveAdvectDiffuse2D!(param::Param,lat::Lattice)

   # Set aliases
   nI,nJ,ω,Δt,TotalTime = param.nI,param.nJ,param.ω,param.Δt,param.TotalTime
   u = param.u

   nw,w,cssq,c = lat.nw, lat.w, lat.cssq, lat.c
   T     = zeros(1,nI,nJ)
   feq   = Array{Float64,3}(undef,nw,nI,nJ)
   ftemp = Array{Float64,3}(undef,nw,nI,nJ)
   t     = 0.0     # Current simulation time

   # Advance
   while t < TotalTime

      # Calculate equilibrium distribution
      for iw = 1:nw
         feq[iw,:,:] .= w[iw].*T[1,:,:].*(1 + dot(c[:,iw],u)/cssq)
      end

      # Collision
      for iw = 1:nw
         lat.f[iw,:,:] .= (1.0-ω).*lat.f[iw,:,:] .+ ω.*feq[iw,:,:]
      end

      # Streaming
      for j=1:nJ, i=1:nI, iw = 1:nw
         iNew = i + c[1,iw]
         jNew = j + c[2,iw]
         (iNew > nI || iNew < 1 || jNew > nJ || jNew < 1) && continue
         ftemp[iw,iNew,jNew] = lat.f[iw,i,j]
      end
      copy!(lat.f,ftemp)

      # Remaining BCs
      setBCs!(lat,param)

      # Calculate macroscopic T
      T .= sum(lat.f, dims=1)

      t += Δt
   end

end

function init_lattice()

   param = setParameters()

   if param.DoEcho printParam(param) end

   Δx,Δt,nI,nJ,LT = param.Δx,param.Δt,param.nI,param.nJ,param.LatticeType

   if LT == 1
      lat = LatticeD2Q4(Δx,Δt,nI,nJ)
   elseif LT == 2
      lat = LatticeD2Q5(Δx,Δt,nI,nJ)
   elseif LT == 3
      lat = LatticeD2Q9(Δx,Δt,nI,nJ)
   end

   return param,lat

end

function printParam(param::Param)
   println("Input parameters for LBM advection-diffusion:")
   println("Lattice type: ",param.LatticeType)
   println("nI, nJ = $(param.nI), $(param.nJ)")
   println("BC type = $(param.BCtype)")
   # may be more to display?
end

function plotSol(lat::Lattice)
   # Visualization
   x = range(0,stop=lat.nI-1,length=lat.nI)
   y = range(0,stop=lat.nJ-1,length=lat.nJ)
   T = sum(lat.f, dims=1)

   # Python follows row major order
   contourf(x,y,T[1,:,:]'); colorbar()
   title("2D Heat Advection-Diffusion in a Plain")
   xlabel("x")
   ylabel("y")
   grid(true)

   # 1D cut plot
   plot(0:lat.nI-1,T[1,:,floor(Int,lat.nJ/2)],
      linestyle="-",marker="o",label="LBM")
end

function main()
   param,lat = init_lattice()

   solveAdvectDiffuse2D!(param,lat)

   plotSol(lat)

end


text:
function inverse_vectorofvariables_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(model, MOI.VectorOfVariables, CP.Inverse)

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(2))
    c2 = MOI.add_constraint(model, x2, MOI.EqualTo(1))

    c3 = MOI.add_constraint(model, MOI.VectorOfVariables([x1, x2, x3, x4]), CP.Inverse(2))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 1
    end
end

function inverse_vectoraffinefunction_test(
    model::MOI.ModelLike,
    config::MOIT.Config,
)
    MOI.empty!(model)

    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.EqualTo{Int})
    @test MOI.supports_constraint(
        model,
        MOI.VectorAffineFunction{Int},
        CP.Inverse,
    )

    x1, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x2, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x3, _ = MOI.add_constrained_variable(model, MOI.Integer())
    x4, _ = MOI.add_constrained_variable(model, MOI.Integer())

    c1 = MOI.add_constraint(model, x1, MOI.EqualTo(2))
    c2 = MOI.add_constraint(model, x2, MOI.Interval(1, 2))

    c3 = MOI.add_constraint(model, MOIU.vectorize([x1, x2, x3, x4]), CP.Inverse(2))

    @test MOI.is_valid(model, x1)
    @test MOI.is_valid(model, x2)
    @test MOI.is_valid(model, x3)
    @test MOI.is_valid(model, x4)
    @test MOI.is_valid(model, c1)
    @test MOI.is_valid(model, c2)
    @test MOI.is_valid(model, c3)

    if config.solve
        MOI.optimize!(model)
        @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
        @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT

        @test MOI.get(model, MOI.ResultCount()) >= 1
        @test MOI.get(model, MOI.VariablePrimal(), x1) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x2) == 1
        @test MOI.get(model, MOI.VariablePrimal(), x3) == 2
        @test MOI.get(model, MOI.VariablePrimal(), x4) == 1
    end
end

const inversetests = Dict(
    "inverse_vectorofvariables" => inverse_vectorofvariables_test,
    "inverse_vectoraffinefunction" => inverse_vectoraffinefunction_test,
)

MOIT.@moitestset inverse


text:
module DirectDetections
using ComponentArrays # TODO: remove last uses of component arrays
using Distributions
using Bijectors
using AbstractMCMC
using AdvancedHMC
using NamedTupleTools
using ForwardDiff
using Logging
using Statistics
using StatsBase
using NamedTupleTools
using DirectImages

using DirectOrbits
# Re-export the mjd helper. It's useful!
export mjd, KeplerianElements, KeplerianElementsDeg

using Base.Threads: @threads
using StaticArrays
using MCMCChains: MCMCChains, Chains
using Random
using DataDeps
using RecipesBase
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

const mjup2msol = 0.0009543

include("types.jl")
include("distributions.jl")
include("models.jl")
include("sampling.jl")
include("analysis.jl")
include("macros.jl")
include("hgca.jl")
include("sonora.jl")

function __init__()
    init_plots()

    # List DataDeps here.
    # These are optional, automatic downloads required for certain
    # functionality.
    register(DataDep("HGCA_eDR3",
        """
        Dataset: Hipparcos Gaia Catalog of Accelerations
        Author: Brandt et al
        License: arXiv.org Non-exclusive license to distribute
        Website: https://arxiv.org/abs/2105.11662

        A catalog by Brandt et al containing cross calibrated proper motions
        between Hipparcos and GAIA eDR3.

        File size: 19MiB
        """,
        "http://physics.ucsb.edu/~tbrandt/HGCA_vEDR3.fits",
        "23684d583baaa236775108b360c650e79770a695e16914b1201f290c1826065c"
    ))

    register(DataDep("SonoraBobcatEvoPhot",
        """
        Dataset: Sonora Bobcat: cloud-free, substellar atmosphere models, spectra, photometry, evolution, and chemistry
        Author: Marley et al
        License: Creative Commons Attribution 4.0 International
        Website: https://zenodo.org/record/5063476

        ``Presented here are models for non-irradiated, substellar mass objects belonging to the Sonora model series, described in Marley et al. (2021).''

        This download contains just the thermal evolution and photometry tables.

        File size: 1MiB
        """,
        "https://zenodo.org/record/5063476/files/evolution_and_photometery.tar.gz?download=1",
        "2198426d1ca0e410fda7b63c3b7f45f3890a8d9f2fcf0a3a1e36e14185283ca5",
        post_fetch_method=unpack
    ))

    return
end

include("precompile.jl")

end


text:
module OperatorClients
using gRPCClient
using ProtoBuf:google

include("dapr.jl")
using .dapr.proto.operator.v1

import Base: show

# begin service: dapr.proto.operator.v1.Operator

export OperatorBlockingClient, OperatorClient

struct OperatorBlockingClient
    controller::gRPCController
    channel::gRPCChannel
    stub::OperatorBlockingStub

    function OperatorBlockingClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = OperatorBlockingStub(channel)
        new(controller, channel, stub)
    end
end

struct OperatorClient
    controller::gRPCController
    channel::gRPCChannel
    stub::OperatorStub

    function OperatorClient(api_base_url::String; kwargs...)
        controller = gRPCController(; kwargs...)
        channel = gRPCChannel(api_base_url)
        stub = OperatorStub(channel)
        new(controller, channel, stub)
    end
end

show(io::IO, client::OperatorBlockingClient) = print(io, "OperatorBlockingClient(", client.channel.baseurl, ")")
show(io::IO, client::OperatorClient) = print(io, "OperatorClient(", client.channel.baseurl, ")")

import .dapr.proto.operator.v1: ComponentUpdate
"""
    ComponentUpdate

- input: dapr.proto.operator.v1.ComponentUpdateRequest
- output: Channel{dapr.proto.operator.v1.ComponentUpdateEvent}
"""
ComponentUpdate(client::OperatorBlockingClient, inp::dapr.proto.operator.v1.ComponentUpdateRequest) = ComponentUpdate(client.stub, client.controller, inp)
ComponentUpdate(client::OperatorClient, inp::dapr.proto.operator.v1.ComponentUpdateRequest, done::Function) = ComponentUpdate(client.stub, client.controller, inp, done)

import .dapr.proto.operator.v1: ListComponents
"""
    ListComponents

- input: dapr.proto.operator.v1.ListComponentsRequest
- output: dapr.proto.operator.v1.ListComponentResponse
"""
ListComponents(client::OperatorBlockingClient, inp::dapr.proto.operator.v1.ListComponentsRequest) = ListComponents(client.stub, client.controller, inp)
ListComponents(client::OperatorClient, inp::dapr.proto.operator.v1.ListComponentsRequest, done::Function) = ListComponents(client.stub, client.controller, inp, done)

import .dapr.proto.operator.v1: GetConfiguration
"""
    GetConfiguration

- input: dapr.proto.operator.v1.GetConfigurationRequest
- output: dapr.proto.operator.v1.GetConfigurationResponse
"""
GetConfiguration(client::OperatorBlockingClient, inp::dapr.proto.operator.v1.GetConfigurationRequest) = GetConfiguration(client.stub, client.controller, inp)
GetConfiguration(client::OperatorClient, inp::dapr.proto.operator.v1.GetConfigurationRequest, done::Function) = GetConfiguration(client.stub, client.controller, inp, done)

import .dapr.proto.operator.v1: ListSubscriptions
"""
    ListSubscriptions

- input: google.protobuf.Empty
- output: dapr.proto.operator.v1.ListSubscriptionsResponse
"""
ListSubscriptions(client::OperatorBlockingClient, inp::google.protobuf.Empty) = ListSubscriptions(client.stub, client.controller, inp)
ListSubscriptions(client::OperatorClient, inp::google.protobuf.Empty, done::Function) = ListSubscriptions(client.stub, client.controller, inp, done)

# end service: dapr.proto.operator.v1.Operator

end # module DaprClients


text:
"""
Authors: Yaroslav Salii, 2020+
         Kara Ignatenko, 2021

Submodule responsible for preprocessing tract data, i.e.
linking tracts with corresponding flight data.

tracts.jl
2021-08-13 v.0.1: First modular version
2021-08-17 v.0.2: Removed dependency on io.jl by getting rid of default args
"""
module Tracts
using DataFrames
using Statistics: mean
using Distances: haversine

export assignDsgAPs,
       assignPsgShares,
       mkAP_pop_dict

#===ASSIGNING DESIGNATED AIRPORTS===#

#= Distance
Starting “haversine”/orthodromic/big-circle distance
Should get the almost true (ball vs. geoid)
when coupled with Earth's radius, or rather, AP's and loc's altitudes
CAVEAT: getting both altitudes is almost as bad
as scrubbing “road distances” from the net
TODO: route through JuliaGeo/Geodesy
=#

"""Earth's radius in km, IAU 2015"""
const R_Earth = 6371

"""
Calculate big-circle distance between two points on Earth.
`x` and `y` are `Union{AbstractVector{T}, NTuple{2, T}}` where T<:Real
"""
myDist(x,y) = haversine(x,y,R_Earth)

"""
Get the designated AP for a node-row `[:ID(:Ste,:Cty,:Tra),:Pop,:LAT,:LNG],
return the *nearest* AP's :IATA_Code and the distance :dst to it in km
"""
#NB: worked slow-ish in Jupyter Notebook; might rework without sorting all this thing.
#UPD: this slow-down was probably due to lack of caching with arguments by default
function getDsgAP(node::DataFrameRow,APs::DataFrame)
    alle=[(dst = myDist((node.LAT,node.LNG)
            ,(ap.LAT,ap.LNG))
            ,IATA_Code = ap.IATA_Code
            ,psg= ap.IN+ap.OUT)  for ap ∈ eachrow(APs)]
    sort!(alle, by = first) # sort by distance to APs
    return (dsg=alle[1],all=alle)
end

"""
Assign a designated AP to each node in `nodes`, return a DF with info on both nodes and
their designated APs, `[:ID,:Pop,:LAT,:LNG,:IATA_Code]`; `:IATA_Code` for the designated AP
dbg columns: `:dst` distance to the dsg AP, `:psg` = `:IN` + `:OUT` passengers through `dsg_AP`
"""
function assignDsgAPs(nodes::DataFrame,APs::DataFrame)
    dsgAPs = map(n -> getDsgAP(n,APs).dsg,eachrow(nodes)) |> DataFrame
    hcat(nodes,dsgAPs)
end

#------------AUX-----------#
"""
Calculate how many people fly through each AP; default to per-state aggregation
`nodes` must have designated APs (:IATA_Code) and populations (:Pop)
"""
function mkClusterPops(nodes::DataFrame)
    #group by the same designated AP: each group is this AP's __catchment area__
    gd = groupby(nodes,[:IATA_Code])
    #= the prime here is Σpopulation by designated AP,
    but it's also a clusterization so here's `mean` :LAT and :LNG as dumb centroids =#
    out = combine(gd,:Pop => sum, :LAT => mean, :LNG => mean,renamecols = false)
end

#TODO: consider leaving it *inside* `assignDsgAPs`, it has to be called anyway
"""Make a `Dict` mapping each AP's ID (:IATA_Code) to the :Pop of its *catchment area*"""
function mkAP_pop_dict(nodes::DataFrame)
    cache=mkClusterPops(nodes)
    return Dict(cache.IATA_Code .=> cache.Pop)
end

"""
GLEaM-like aggregation: lump together all nodes in AP's catchment area (start with *census tracts*)
difference with GLEaM: some airports may get lost, whereas GLEaM starts with APs (gotta re-check that!)
"""
function aggByAPs(tracts::DataFrame, APs::DataFrame)
    return assignDsgAPs(tracts, APs) |> mkClusterPops
end

"""
For a node `n`, find the fraction of pop in `dsg_n`'s catchment area.
#the node MUST have an `:IATA_Code` column (its *designated AP*)
"""
function nodePsgShare(node::DataFrameRow,dAP_pop::Dict)
    return node.Pop / dAP_pop[node.IATA_Code]
end

"""
Apply `nodePsgShare` to all nodes in `nodes` and store them in the `shr` column of the output.
each node MUST have an `:IATA_Code` column (its *designated AP*)
"""
function assignPsgShares(nodes::DataFrame,dAP_pop::Dict)
    psgShares = map(n -> nodePsgShare(n,dAP_pop), eachrow(nodes)) #compute the shares
    out = hcat(nodes,DataFrame("shr" => psgShares)) #add them as a :shr column
end

end


text:
let
 tot2 = 0
 tot3 = 0
 open("d2data.txt") do file
   lines = readlines(file)

   for l1 in lines, l2 in lines
     if l1 == l2
       continue
     end
     diffCount = 0
     for i = 1:length(l1)
       if l1[i] != l2[i]
         diffCount += 1
       end
     end
     if diffCount == 1
       println(l1)
       println(l2)
     else
       #println(diffCount)
       #println(l1)
       #println(l2)
     end
   end
end

end