text:
abstract type AbstractAbundanceTable{T <: Real,
                                     F <: AbstractFeature,
                                     S <: AbstractSample} <: EcoBase.AbstractAssemblage{T, F, S}
end

"""
    CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}

An `AbstractAssemblage` from [EcoBase.jl](https://github.com/EcoJulia/EcoBase.jl)
that uses a `SparseMatrixCSC` under the hood.

`CommunityProfile`s are tables with `AbstractFeature`-indexed rows and
`AbstractSample`-indexed columns.
Note - we can use the `name` of samples and features to index.
"""
mutable struct CommunityProfile{T, F, S} <: AbstractAbundanceTable{T, F, S}
    abundances::SparseMatrixCSC{T}
    features::AbstractVector{F}
    samples::AbstractVector{S}
    fidx::Dictionary{String, Int}
    sidx::Dictionary{String, Int}

    function CommunityProfile(abunds::AbstractSparseMatrix,
                              feats::AbstractVector{<:AbstractFeature},
                              smpls::AbstractVector{<:AbstractSample})

        length(feats) == size(abunds, 1) || throw(DimensionMismatch("Number of features must equal number of rows in matrix"))
        length(smpls) == size(abunds, 2) || throw(DimensionMismatch("Number of samples must equal number columns in matrix"))
        fidx = Dictionary(String.(feats), eachindex(feats))
        sidx = Dictionary(String.(smpls), eachindex(smpls))

        T = eltype(abunds)
        F = eltype(feats)
        S = eltype(smpls)
        return new{T, F, S}(abunds, feats, smpls, fidx, sidx)
    end
end

function CommunityProfile(tab::AbstractVecOrMat,
                          feats::AbstractVector{<:AbstractFeature},
                          smpls::AbstractVector{<:AbstractSample})
    return CommunityProfile(sparse(tab), feats, smpls)
end

# single-column CommunityProfile
function CommunityProfile(tab::AbstractVecOrMat,
                          feats::AbstractVector{<:AbstractFeature},
                          smpl::AbstractSample)
    return CommunityProfile(sparse(reshape(tab, size(tab,1), size(tab,2))), feats, [smpl])
end

# single-row CommunityProfile
function CommunityProfile(tab::AbstractVecOrMat,
                          feat::AbstractFeature,
                          smpls::AbstractVector{<:AbstractSample})
    return CommunityProfile(sparse(reshape(tab, size(tab,1), size(tab,2))), [feat], smpls)
end

## -- Convienience functions -- ##

function ==(p1::CommunityProfile, p2::CommunityProfile)
    return abundances(p1) == abundances(p2) &&
           samples(p1)    == samples(p2) &&
           features(p1)   == features(p2)
end

"""
    taxonomicprofile(mat, features, samples)
"""
function taxonomicprofile(mat, features::AbstractVector{<:AbstractString}, samples::AbstractVector{<:AbstractString})
    CommunityProfile(mat, Taxon.(features), MicrobiomeSample.(samples))
end

"""
    functionalprofile(mat, features, samples)
"""
function functionalprofile(mat, features::AbstractVector{<:AbstractString}, samples::AbstractVector{<:AbstractString})
    CommunityProfile(mat, GeneFunction.(features), MicrobiomeSample.(samples))
end

"""
    metabolicprofile(mat, features, samples)
"""
function metabolicprofile(mat, features::AbstractVector{<:AbstractString}, samples::AbstractVector{<:AbstractString})
    CommunityProfile(mat, Metabolite.(features), MicrobiomeSample.(samples))
end

@testset "String Constructors" begin
    tp = taxonomicprofile([1 0; 0 1], ["feature1", "feature2"], ["sample1", "sample2"])
    @test tp isa CommunityProfile
    @test all(f-> f isa Taxon, features(tp))
    fp = functionalprofile([1 0; 0 1], ["feature1", "feature2"], ["sample1", "sample2"])
    @test fp isa CommunityProfile
    @test all(f-> f isa GeneFunction, features(fp))
    mp = metabolicprofile([1 0; 0 1], ["feature1", "feature2"], ["sample1", "sample2"])
    @test mp isa CommunityProfile
    @test all(f-> f isa Metabolite, features(mp))
end

"""
    features(at::AbstractAbundanceTable)

Returns features in `at`. To get featurenames instead, use [`featurenames`](@ref).
"""
features(at::AbstractAbundanceTable) = at.features

"""
    samples(at::AbstractAbundanceTable)

Returns samples in `at`. To get samplenames instead, use [`samplenames`](@ref).
"""
samples(at::AbstractAbundanceTable) = at.samples

"""
    samples(at::AbstractAbundanceTable, name::AbstractString)

Returns sample in `at` with name `name`.
"""
function samples(at::AbstractAbundanceTable, name::AbstractString)
    idx = findall(==(name), samplenames(at))
    length(idx) == 0 && throw(IndexError("No samples called $name"))
    length(idx) > 1 && throw(IndexError("More than one sample matches name $name"))
    return samples(at)[axes(at.abundances, 2)][first(idx)]
end

profiletype(at::AbstractAbundanceTable) = eltype(features(at))
ranks(at::AbstractAbundanceTable) = taxrank.(features(at))

Base.size(at::AbstractAbundanceTable, dims...) = size(at.abundances, dims...)

Base.copy(at::AbstractAbundanceTable) = CommunityProfile(copy(abundances(at)), copy(features(at)), deepcopy(samples(at)))

# -- Indexing -- #

function Base.getindex(at::AbstractAbundanceTable, rowind, colind)
    rows = _toind(at.fidx, rowind)
    cols = _toind(at.sidx, colind)

    mat = copy(abundances(at)[rows, cols])

    isempty(size(mat)) && return mat

    feat = copy(features(at))[rows]
    smpl = deepcopy(samples(at))[cols]

    feat isa AbstractFeature && (mat = reshape(mat, 1, length(mat)))
    return CommunityProfile(mat, feat, smpl)
end

# For integers, or vectors of integers, just return them
_toind(_, ind) = ind
_toind(_, inds::AbstractVector) = inds

# for strings and regex, look for matches
_toind(d, ind::AbstractString) = only((d[i] for i in findall(key-> key == ind, keys(d))))
_toind(d, ind::Regex)          = [d[i] for i in findall(key-> contains(key, ind), keys(d))]

_toind(d, inds::AbstractVector{<:AbstractString}) = [d[i] for i in findall(key-> any(ind-> key == ind, inds), keys(d))]
_toind(d, inds::AbstractVector{<:Regex})          = [d[i] for i in findall(key-> any(ind-> contains(key, ind), inds), keys(d))]

# For samples and features, look for string representation matches
_toind(d, ind::Union{AbstractSample, AbstractFeature}) = only((d[i] for i in findall(key-> key == String(ind), keys(d))))
_toind(d, inds::AbstractVector{<:Union{AbstractSample, AbstractFeature}}) = [d[i] for i in findall(key-> any(ind-> key == String(ind), inds), keys(d))]

## -- EcoBase Translations -- ##
# see src/ecobase.jl for Microbiome function names
# thing => feature
# place => sample
# occurrences => abundances (or totals)

EcoBase.thingnames(at::AbstractAbundanceTable) = name.(features(at))
EcoBase.placenames(at::AbstractAbundanceTable) = name.(samples(at))
EcoBase.occurrences(at::AbstractAbundanceTable) = at.abundances
EcoBase.nthings(at::AbstractAbundanceTable) = size(at, 1)
EcoBase.nplaces(at::AbstractAbundanceTable) = size(at, 2)
## todo
# EcoBase.thingoccurrences(at::AbstractAbundanceTable, things) = nothing
# EcoBase.placeoccurrences(at::AbstractAbundanceTable, places) = nothing

# for custom printing
EcoBase.thingkind(asm::AbstractAbundanceTable) = "feature"
EcoBase.placekind(asm::AbstractAbundanceTable) = "sample"
## not needed for now
# EcoBase.thingkindplural(asm::AbstractAbundanceTable) = "$(thingkind(asm))s"
# EcoBase.placekindplural(asm::AbstractAbundanceTable) = "$(placekind(asm))s"

"""
    featuretotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a nfeatures x 1 `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(featuretotals(at))`.
"""
featuretotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=2)

"""
    sampletotals(at::AbstractAbundanceTable)

Returns sum of each row (feature) in `at`.
Note, return value is a 1 x nsamples `Matrix`, not a `Vector`.
If you need 1D `Vector`, use `vec(sampletotals(at))`.
"""
sampletotals(at::AbstractAbundanceTable) = sum(abundances(at), dims=1)

## -- Tables Interface -- ##

Tables.istable(::AbstractAbundanceTable) = true
Tables.columnaccess(::AbstractAbundanceTable) = true
Tables.rowaccess(::AbstractAbundanceTable) = true

Tables.getcolumn(at::AbstractAbundanceTable, i::Int) = i == 1 ? featurenames(at) : abundances(at[:, i-1])
Tables.getcolumn(at::AbstractAbundanceTable, i::AbstractString) = i == "features" ? features(at) : abundances(at[:, i])
Tables.getcolumn(at::AbstractAbundanceTable, i::Symbol) = Tables.getcolumn(at, string(i))

Tables.columnnames(at::AbstractAbundanceTable) = [:features, Symbol.(samplenames(at))...]

function Tables.schema(at::AbstractAbundanceTable)
    elt = eltype(abundances(at))
    coltypes = [eltype(features(at)), (elt for _ in 1:nsamples(at))...]
    return Tables.Schema(Tables.columnnames(at), coltypes)
end

Tables.columns(at::AbstractAbundanceTable) = (; (col => Tables.getcolumn(at, col) for col in Tables.columnnames(at))...)

function _makerow(row::AbstractAbundanceTable)
    size(row, 1) == 1 || error("Can't make row from table of size $(size(row))")
    NamedTuple{(:features, Symbol.(samplenames(row))...)}((first(features(row)), abundances(row)...))
end

Tables.rows(at::AbstractAbundanceTable) = (_makerow(at[i, :]) for i in 1:nfeatures(at))

## -- Methods for absolute and relative abundances -- ##

# """
#     filterabund(abun::AbstractComMatrix, n::Int=minimum(10, nfeatures(abun)))

# Filter an abundance table to the top `n` features accross all samples

# This function also adds a row for "other", which sums the abundances of the
# remaining features.
# """
# function filterabund(abun::AbstractAbundanceTable, n::Int=minimum(10, nfeatures(abun)))
#     # TODO: add prevalence filter

#     totals = featuretotals(abun)

#     srt = sortperm(totals, rev=true)

#     newabun = getfeature(abun, srt[1:n])

#     remainder = [sum(occurrences(abun)[srt[n+1:end], i]) for i in 1:size(abun, 2)]'
#     newabun = vcat(newabun, remainder)
#     newrows = cat(featurenames(abun)[srt[1:n]], ["other"], dims=1)

#     return abundancetable(newabun, samplenames(abun), newrows)
# end

"""
    relativeabundance!(a::AbstractAbundanceTable; kind::Symbol=:fraction)

Normalize each sample in AbstractAbundanceTable to the sum of the sample.

By default, columns sum to 1.0.
Use `kind=:percent` for columns to sum to 100.
"""
function relativeabundance!(at::AbstractAbundanceTable; kind::Symbol=:fraction)
    in(kind, [:percent, :fraction]) || throw(ArgumentError("Invalid kind: $kind"))
    eltype(abundances(at)) <: AbstractFloat || throw(ArgumentError("relativeabundance! requires profile to have AbstractFloat eltype. Try relativeabundance instead"))
    abund = abundances(at)
    abund ./= sampletotals(at)
    kind == :percent && (abund .*= 100)
    dropzeros!(abund) # shouldn't need dropzeros - https://github.com/JuliaLang/julia/issues/39018
    return at
end

"""
    relativeabundance(at::AbstractAbundanceTable, kind::Symbol=:fraction)

Like [`relativeabundance!`](@ref), but does not mutate original.
"""
function relativeabundance(at::AbstractAbundanceTable, kind::Symbol=:fraction)
    comm = CommunityProfile(float.(abundances(at)), deepcopy(features(at)), deepcopy(samples(at)))
    relativeabundance!(comm)
end

"""
    present(t::Union{Real, Missing}, minabundance::Real=0.0)
    present(at::AbstractAbundanceTable, minabundance::Real=0.0)

Check if a given (non-zero) value is greater than or equal to a minimum value.
If the minimum abundance is 0, just checks if value is non-zero.

If used on an `AbstractAbundanceTable`, returns a sparse boolean matrix of the same size.
"""
function present(t::Real, minabundance::Real=0.0)
    (minabundance >= 0 && t >= 0) || throw(DomainError("Only defined for positive values"))
    t == 0 ? false : t >= minabundance
end

present(::Missing, m::Real=0.0) = missing

function present(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mat = spzeros(Bool, size(at)...)
    for i in eachindex(mat)
        mat[i] = present(at[Tuple(i)...], minabundance)
    end
    return mat
end

"""
    prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0)
    prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)

Return the fraction of values that are greater than or equal to a minimum.
If the minimum abundance is 0, returns the fraction of non-zero values.

If used on an `AbstractAbundanceTable`,
returns a prevalence value for each `feature` accross the `sample`s.
"""
prevalence(a::AbstractArray{<:Real}, minabundance::Real=0.0) = mean(x-> present(x, minabundance), a)

# makes it work for any iterable
prevalence(a, minabundance::Real=0.0) = mean(x-> present(x, minabundance), (y for y in a))

function prevalence(at::AbstractAbundanceTable, minabundance::Real=0.0)
    mean(x-> present(x, minabundance), abundances(at), dims=2)
end

"""
    prevalence_filter(comm::AbstractAbundanceTable; minabundance=0.0; minprevalence=0.05, renorm=false)

Return a filtered `CommunityProfile` where features with prevalence lower than `minprevalence` are removed.
By default, a feature is considered "present" if > 0, but this can be changed by setting `minabundance`.

Optionally, set `renorm = true` to calculate relative abundances after low prevalence features are removed.
"""
function prevalence_filter(comm::AbstractAbundanceTable; minabundance=0.0, minprevalence=0.05, renorm=false)
    comm = comm[vec(prevalence(comm, minabundance) .>= minprevalence), :]
    return renorm ? relativeabundance(comm) : comm
end

"""
    rankfilter(comm::AbstractAbundanceTable, cl::Union{Symbol, Int}; keepempty=false)

Return a copy of `comm`, where only rows that have `taxrank(feature) == cl` are kept.
Use `keepempty = true` to also keep features that don't have a `rank` (eg "UNIDENTIFIED").
"""
function rankfilter(comm::AbstractAbundanceTable, cl::Symbol; keepempty=false)
    in(cl, keys(_ranks)) ||  error("Invalid rank $cl, must be one of $(keys(_ranks))")
    if keepempty
        return filter(f-> !hasrank(f) || taxrank(f) == cl, comm)
    else
        return filter(f-> hasrank(f) && taxrank(f) == cl, comm)
    end
end

function rankfilter(comm::AbstractAbundanceTable, rank::Int; keepempty=false)
    0 <= rank <= 9 ||  error("Invalid rank $rank, must be one of $_ranks")
    return rankfilter(comm, keys(_ranks)[rank+1]; keepempty)
end

## Metadata

"""
    metadata(commp::CommunityProfile)

Returns iterator of `NamedTuple` per sample, where keys are `:sample`
and each metadata key found in `commp`.
Samples without given metadata are filled with `missing`.

Returned values can be passed to any Tables.rowtable - compliant type,
eg `DataFrame`.
"""
function metadata(commp::CommunityProfile)
    ss = samples(commp)
    cols = unique(reduce(vcat, collect.(keys.(metadata.(samples(commp))))))
    return Tables.rowtable(merge((; sample=name(s)),
                     NamedTuple(c => get(s, c, missing) for c in cols)
                    ) for s in ss)
end

"""
    set!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    set!(commp::CommunityProfile, sample::AbstractString, md::Union{AbstractDict, NamedTuple})

Update or insert a value `val` to the metadata of `sample` in the CommunityProfile `commp` using a Symbol `prop`.
If you want an error to be thrown if the value already exists, use [`insert!`](@ref).

Can also pass a Dictionary or NamedTuple containing key=> value pairs,
all of which will be `set!`.
"""
function set!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot set! $prop for $(typeof(sample)).")
    set!(sample.metadata, prop, val)
    return sample
end

function set!(commp::CommunityProfile, sample::AbstractString, md::Dictionary)
    for (key, value) in pairs(md)
        set!(commp, sample, key, value)
    end
    return nothing
end

function set!(commp::CommunityProfile, sample::AbstractString, md::Union{AbstractDict, NamedTuple})
    md = Dictionary(md)
    set!(commp, sample, md)
    return nothing
end

"""
    set!(cp::CommunityProfile, md; namecol=:sample)

Add metadata (in the form of a `Tables.jl` table) a `CommunityProfile`.
One column (`namecol`) should contain sample names that exist in `commp`,
and other columns should contain metadata that will be added to the metadata of each sample.
"""
function set!(commp::CommunityProfile, md; namecol=:sample)
    Tables.istable(md) || throw(ArgumentError("Metadata must be a Tables.table"))
    sns = Set(samplenames(commp))
    md = filter(row-> row[namecol] in sns, md)
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            set!(sample, k, row[k])
        end
    end
    return nothing
end

"""
    insert!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)

Insert a value `val` to the metadata of `sample` in a CommunityProfile `commp` using a Symbol `prop`,
and it will throw an error if `prop` exists.
If you don't want an error to be thrown if the value exists, use [`set!`](@ref).
"""
function insert!(commp::CommunityProfile, sample::AbstractString, prop::Symbol, val)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot insert! $prop for $(typeof(sample)).")
    insert!(sample.metadata, prop, val)
    return sample
end

function insert!(commp::CommunityProfile, sample::AbstractString, md::Dictionary)
    isempty(Set(keys(commp, sample)) ∩ Set(keys(md))) || throw(IndexError("Duplicate keys found. Use `set!` to overwrite"))
    for (key, value) in pairs(md)
        insert!(commp, sample, key, value)
    end
    return nothing
end

function insert!(commp::CommunityProfile, sample::AbstractString, md::Union{<:AbstractDict, NamedTuple})
    md = Dictionary(md)
    insert!(commp, sample, md)
    return nothing
end

"""
    insert!(cp::CommunityProfile, md; namecol=:sample)

Add metadata (in the form of a `Tables.jl` table) a `CommunityProfile`.
One column (`namecol`) should contain sample names that exist in `commp`,
and other columns should contain metadata that will be added to the metadata of each sample.

Before starting, this will check that every value in every row is `insert!`able,
and will throw an error if not.
This requires iterating over the metadata table twice, which may be slow.
If performance matters, you can use `set!` instead,
though this will overwrite existing data.
"""
function insert!(commp::CommunityProfile, md; namecol=:sample, careful=true)
    Tables.istable(md) || throw(ArgumentError("Metadata must be a Tables.table"))
    sns = Set(samplenames(commp))
    md = filter(row-> row[namecol] in sns, md)
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            haskey(sample, k) && throw(IndexError("Duplicate metadata detected. Use `set!` to force overwrite."))
        end
    end
    for row in Tables.rows(md)
        sample = samples(commp, row[namecol])
        ks = filter(!=(namecol), keys(first(md)))
        for k in ks
            set!(sample, k, row[k])
        end
    end
    return nothing
end

"""
unset!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)

Delete a metadata entry in `sample` from CommunityProfile `commp` using the Symbol `prop`.
If you want an error to be thrown if the value does not exist, use [`delete!`](@ref).
"""
function unset!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot unset! $prop for $(typeof(sample)).")
    unset!(sample.metadata, prop)
    return sample
end

"""
    delete!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)

Delete a metadata entry in `sample` from CommunityProfile `commp` using the Symbol `prop` if it exists, or throw an error otherwise.
If you don't want an error to be thrown if the value does not exist, use [`unset!`](@ref).
"""
function delete!(commp::CommunityProfile, sample::AbstractString, prop::Symbol)
    sample = samples(commp, sample)
    prop in _restricted_fields(sample) && error("Cannot delete! $prop for $(typeof(sample)).")
    delete!(sample.metadata, prop)
    return sample
end

"""
    keys(commp::CommunityProfile, sample::AbstractString)

Return an iterator over all keys of the metadata attached to `sample` in a CommunityProfile `commp`.
`collect(keys(commp, sample))` returns an array of keys.
"""
Base.keys(commp::CommunityProfile, sample::AbstractString) = keys(metadata(samples(commp, sample)))

"""
    haskey(commp::CommunityProfile, sample::AbstractString, key::Symbol)

Determine whether the metadata of `sample` in a CommunityProfile `commp` has a mapping for a given `key`.
Use `!haskey` to determine whether a `sample` in a CommunityProfile doesn't have a mapping for a given `key`
"""
Base.haskey(commp::CommunityProfile, sample::AbstractString, key::Symbol) = in(key, keys(samples(commp, sample)))

"""
    get(commp::CommunityProfile, sample::AbstractString, key::Symbol, default)

Return the value of the metadata in a `sample` stored for the given `key`, or the given `default` value if no mapping for the key is present.
"""
Base.get(commp::CommunityProfile, sample::AbstractString, key::Symbol, default=missing) = get(metadata(samples(commp, sample)), key, default)

"""
    get(commp::CommunityProfile, key::Symbol, default)

Return the value of the metadata in a `sample` stored for the given `key`, or the given `default` value if no mapping for the key is present.
"""
Base.get(commp::CommunityProfile, key::Symbol, default=missing) = [get(commp, sample, key, default) for sample in samplenames(commp)]

"""
    filter(f, comm::CommunityProfile)

Apply `f` to the features of `comm`,
and return a copy where `f(feature)` is `true`.
"""
function Base.filter(f::Function, commp::CommunityProfile)
    ridx = findall(f, features(commp))
    isempty(ridx) && error("Can't return empty profile")
    return copy(commp[ridx, :])
end


text:
ApproxFun.treecount(::Union{ConcreteHilbert,OffHilbert})=1
ApproxFun.treecount(M::HilbertWrapper)=1+ApproxFun.opcount(M.op)

ApproxFun.texname(D::ConcreteHilbert)=(labels[nd]=string(nd)*":"*(D.order==1?"\${\\cal H}":"\${\\cal H}\^"*string(D.order))*"\$")
ApproxFun.texname(D::OffHilbert)=(labels[nd]=string(nd)*":\${\\cal OH}\$")
ApproxFun.texname(D::HilbertWrapper)=(D.order==1?"\$({\\cal H}":"\$({\\cal H}\^"*string(D.order))*")\$"

@eval ApproxFun.add_edges!(A::HilbertWrapper,nd,M,labels)=treeadd_edges!(string(nd)*":"*texname(A),[A.op],nd,M,labels)

ext: jl


text:
export by, aggregate

# TODO: remove definitions in 2.0 release
by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))
aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))

# TODO: remove deprecation in 2.0 release
import Base.delete!
@deprecate delete!(df::DataFrame, inds) deleteat!(df::DataFrame, inds)


text:
export insert_bulk

## Bulk insert
function insert_bulk(collection::Mongo.MongoCollection,
documents::Vector{Mongo.BSONObject},
flags::Int = Mongo.MongoInsertFlags.None)
bulk = ccall((:mongoc_collection_create_bulk_operation, Mongo.libmongoc),
 Ptr{Void}, (Ptr{Void}, Bool, Ptr{Void}),
 collection._wrap_, true, C_NULL)

for doc in documents
    ccall((:mongoc_bulk_operation_insert, Mongo.libmongoc),
    Void, (Ptr{Void}, Ptr{Void}),
    bulk, doc._wrap_)
end

reply = LibBSON.BSONObject()

ret = ccall((:mongoc_bulk_operation_execute, Mongo.libmongoc),
Bool, (Ptr{Void}, Ptr{Void}, Ptr{Void}),
bulk,
      reply._wrap_,
      C_NULL
            )

ccall((:mongoc_bulk_operation_destroy, Mongo.libmongoc),
      Void, (Ptr{Void},), bulk)

return ret
end


text:
function update!(   pssrb::PotentialSimulationSetupRB{T}; use_nthreads::Int = Base.Threads.nthreads(),
        depletion_handling::Val{depletion_handling_enabled} = Val{false}(), only2d::Val{only_2d} = Val{false}(),
        is_weighting_potential::Val{_is_weighting_potential} = Val{false}())::Nothing where {T, depletion_handling_enabled, only_2d, _is_weighting_potential}
    outerloop!(pssrb, use_nthreads, Val{true}(), depletion_handling, is_weighting_potential, only2d)
    apply_boundary_conditions!(pssrb, Val{true}(), only2d)
    outerloop!(pssrb, use_nthreads, Val{false}(), depletion_handling, is_weighting_potential, only2d)
    apply_boundary_conditions!(pssrb, Val{false}(), only2d)
    nothing
end

@fastmath function outerloop!( pssrb::PotentialSimulationSetupRB{T}, use_nthreads::Int,
    update_even_points::Val{even_points},
    depletion_handling::Val{depletion_handling_enabled},
    is_weighting_potential::Val{_is_weighting_potential},
    only2d::Val{only_2d})::Nothing where {T, even_points, depletion_handling_enabled, _is_weighting_potential, only_2d}
    @inbounds begin
        rb_tar_idx::Int, rb_src_idx::Int = even_points ? (rb_even::Int, rb_odd::Int) : (rb_odd::Int, rb_even::Int)

        @onthreads 1:use_nthreads for idx3 in workpart(2:1:(size(pssrb.potential, 3) - 1), 1:use_nthreads, Base.Threads.threadid())
            middleloop!( idx3, rb_tar_idx, rb_src_idx, pssrb,
                         update_even_points, depletion_handling,
                         is_weighting_potential, only2d, Val(iseven(idx3))::Union{Val{true}, Val{false}})
        end
    end
    nothing
end

function _guess_optimal_number_of_threads_for_SOR(gs::NTuple{3, Integer}, max_nthreads::Integer, S::Union{Type{Cylindrical}, Type{Cartesian}})::Int
    max_nthreads = min(Base.Threads.nthreads(), max_nthreads)
    n = S == Cylindrical ? gs[2] * gs[3] : gs[1] * gs[2] # Number of grid points to be updated in each iteration of the outer loop
    #=
        Due to threading overhead calculating small grids with too many threads results in a worse performance.
        Thus, we limit them depending on the grid size.
        These ranges may depend on the hardware and are set by educated guesses.
        # nt = if n < 100
        #     4
        # elseif n < 200
        #     8
        # elseif n < 400
        #     16
        # elseif n < 800
        #     32
        # elseif n < 1600
        #     64
        # else
        #     max_nthreads
        # end
    =#
    return min(nextpow(2, max(cld(n+1, 25), 4)), max_nthreads)
end


text:
tests = [
    "Memory_Tests",
    "Utility_Tests",
]

for test in tests
    include(test * ".jl")
end


text:
module HPCManager

using Reexport
@reexport using Distributed
export init_cluster, close_workers

function init_cluster(method = :worker_per_task; kwargs...)
    host_name = gethostname()
    node_list = get_node_list()
    tpn = get_task_per_node()
    tpn = [parse(Int, i) for i in tpn]
    if (length(node_list) == 1) && (method == :worker_per_node)
        @info "Only one nodes exists, skip cluster initialization."
        nothing
    elseif (length(node_list) == 1) && (method == :worker_per_task)
        @info "Task is assigned to a single node. Starting with local manager..."
        addprocs(tpn[1]; kwargs...)
    elseif (length(node_list) > 1) && (method == :worker_per_node)
        machine = []
        for i in zip(node_list, tpn)
            push!(machine, (i[1], 1))
        end
        @info "Task is assigned to multiple nodes. Starting with SSH manager..."
        @info "Machine configuration: " machine
        addprocs(machine; kwargs...)
    elseif (length(node_list) > 1) && (method == :worker_per_task)
        if length(tpn) == 1
            tpn = fill(tpn[1], length(node_list))
        end
        machine = []
        for i in zip(node_list, tpn)
            push!(machine, i)
        end
        @info "Task is assigned to multiple nodes. Starting with SSH manager..."
        @info "Machine configuration: " machine
        addprocs(machine; kwargs...)
    end
end

function close_workers()
    @info "Closing all workers..."
    for i in workers()
        if i != 1
            rmprocs(i)
        end
    end
end

function get_node_list()
    slurm_node = ENV["SLURM_JOB_NODELIST"]
    @debug "Slurm node list: " slurm_node
    re = match(r"\[.*\]", slurm_node)
    if re === nothing
        return [slurm_node]
    else
        header = slurm_node[1:re.offset-1]
        raw_str = strip(re.match, ['[', ']'])
        raw_str = split(raw_str, ',')
        res = []
        for str_iter in raw_str
            if occursin("-", str_iter)
                sp = split(str_iter, '-')
                st = parse(Int, sp[1])
                en = parse(Int, sp[2])
                for n = st:en
                    push!(res, header * lpad(n, 2, '0'))
                end
            else
                push!(res, header * str_iter)
            end
        end
        return res
    end
end

function get_task_per_node()
    tpn = ENV["SLURM_TASKS_PER_NODE"]
    @debug "Slurm tasks per node: " tpn
    res = []
    raw_str = split(tpn, ',')
    for str_iter in raw_str
        if occursin("x", str_iter)
            x_idx = findfirst(isequal('x'), str_iter)
            str = str_iter[1:x_idx-2]
            num = parse(Int, str_iter[x_idx+1:end-1])
            for i = 1:num
                push!(res, str)
            end
        else
            push!(res, str_iter)
        end
    end
    res
end

end # module


text:
#test2
yb=1974
tmax,jmax=30,16
w=zeros(tmax,jmax)
x=zeros(tmax,jmax)
y=zeros(tmax,jmax)
w=readdlm("w.txt")
y=readdlm("y.txt")
x=readdlm("x.txt")
#println(size(x))
#println(tmax," ",jmax)

#println(csmtx(yb,tmax,jmax,x))
#println(csmtx2csvs(csmtx(yb,tmax,jmax,x)))
x==csvs2csmtx(csmtx2csvs(csmtx(yb,tmax,jmax,x)))

max_stars_repo_name: picestcode/Fishmetica

max_issues_repo_name: picestcode/Fishmetica

max_forks_repo_name: picestcode/Fishmetica


text:
# Constraints
const UncertainConstraintRef = JuMP.ConstraintRef{RobustModel,ConstraintIndex}
function JuMP.add_constraint(
    model::RobustModel,
    c::JuMP.AbstractConstraint,
    name::String = "",
)
    model.nextconidx += 1
    index = ConstraintIndex(model.nextconidx)
    #print("INDEX: ", index)
    cref = JuMP.ConstraintRef(model, index, JuMP.shape(c))
    model.uncertainConstraints[index] = c
    JuMP.set_name(cref, name)
    return cref
end
function JuMP.delete(model::RobustModel, constraint_ref::UncertainConstraintRef)
    @assert JuMP.is_valid(model, constraint_ref)
    delete!(model.uncertainConstraints, constraint_ref.index)
    return delete!(model.con_to_name, constraint_ref.index)
end
function JuMP.delete(model::RobustModel, con_refs::Vector{<:UncertainConstraintRef})
    return JuMP.delete.(model, con_refs)
end
function JuMP.is_valid(model::RobustModel, constraint_ref::UncertainConstraintRef)
    return (
        model === constraint_ref.model &&
        constraint_ref.index in keys(model.uncertainConstraints)
    )
end
function JuMP.constraint_object(cref::UncertainConstraintRef)
    return cref.model.uncertainConstraints[cref.index]
end
function JuMP.num_constraints(
    model::RobustModel,
    F::Type{<:JuMP.AbstractJuMPScalar},
    S::Type{<:MOI.AbstractSet},
)
    return count(
        con -> con isa JuMP.ScalarConstraint{F,S},
        values(model.uncertainConstraints),
    )
end
function JuMP.num_constraints(
    model::RobustModel,
    ::Type{<:Vector{F}},
    S::Type{<:MOI.AbstractSet},
) where {F<:JuMP.AbstractJuMPScalar}
    return count(
        con -> con isa JuMP.VectorConstraint{F,S},
        values(model.uncertainConstraints),
    )
end


text:
using FieldDistributionNonuniformMedium
using Test

const C = 299792458

all_eq(x::AbstractArray) = all(xᵢ -> xᵢ == x[1], x)

@testset "FieldDistributionNonuniformMedium.jl" begin
    include("simulation.jl")
    include("plot.jl")
end


text:
function train!(agent::AbstractAgent, game::SnakeAI.Game)
    # Get the current step
    old_state = SnakeAI.get_state(game)

    # Get the predicted move for the state
    move = get_action(agent, old_state)
    SnakeAI.send_inputs!(game, move)

    # Play the step
    reward, done, score = SnakeAI.play_step!(game)
    new_state = SnakeAI.get_state(game)

    # Train the short memory
    train_short_memory(agent, old_state, move, reward, new_state, done)

    # Remember
    remember(agent, old_state, move, reward, new_state, done)

    if done
        # Reset the game
        train_long_memory(agent)
        SnakeAI.reset!(game)
        agent.n_games += 1

        if score > agent.record
            agent.record = score
            # save_model(joinpath(MODELS_PATH, "model_$(agent.n_games).bson"), agent.model)
        end
    end

    return done
end

function remember(
    agent::AbstractAgent,
    state::S,
    action::S,
    reward::T,
    next_state::S,
    done::Bool
) where {T<:Integer,S<:AbstractArray{<:T}}
    push!(agent.memory.data, (state, action, [reward], next_state, convert.(Int, [done])))
end

function train_short_memory(
    agent::AbstractAgent,
    state::S,
    action::S,
    reward::T,
    next_state::S,
    done::Bool
) where {T<:Integer,S<:AbstractArray{<:T}}
    update!(agent, state, action, reward, next_state, done)
end

function train_long_memory(agent::AbstractAgent)
    if length(agent.memory.data) > BATCH_SIZE
        mini_sample = sample(agent.memory.data, BATCH_SIZE)
    else
        mini_sample = agent.memory.data
    end

    states, actions, rewards, next_states, dones = map(x -> getfield.(mini_sample, x), fieldnames(eltype(mini_sample)))

    update!(agent, states, actions, rewards, next_states, dones)
end

function get_action(agent::SnakeAgent, state::AbstractArray{<:Integer}; rand_range=1:200)
    agent.ϵ = 80 - agent.n_games
    final_move = zeros(Int, 3)

    if rand(rand_range) < agent.ϵ
        move = rand(1:3)
        final_move[move] = 1
    else
        pred = agent.model(state)
        final_move[Flux.onecold(pred)] = 1
    end

    return final_move
end

function update!(
    agent::SnakeAgent,
    state::Union{A,AA},
    action::Union{A,AA},
    reward::Union{T,AA},
    next_state::Union{A,AA},
    done::Union{Bool,AA};
    α::Float32=0.9f0    # Step size
) where {T<:Integer,A<:AbstractArray{<:T},AA<:AbstractArray{A}}
    # Batching the states and converting data to Float32 (done implicitly otherwise)
    state = Flux.batch(state) |> x -> convert.(Float32, x)
    next_state = Flux.batch(next_state) |> x -> convert.(Float32, x)
    action = Flux.batch(action) |> x -> convert.(Float32, x)
    reward = Flux.batch(reward) |> x -> convert.(Float32, x)
    done = Flux.batch(done)

    # Model's prediction for next state
    y = agent.model(next_state)

    # Get the model's params for back propagation
    ps = Flux.params(agent.model)

    # Calculate the gradients
    gs = Flux.gradient(ps) do
        # Forward pass
        ŷ = agent.model(state)

        # Creating buffer to allow mutability when calculating gradients
        Rₙ = Buffer(ŷ, size(ŷ))

        # Adjusting values of current state with next state's knowledge
        for idx in 1:length(done)
            # Copy preds into buffer
            Rₙ[:, idx] = ŷ[:, idx]

            Qₙ = reward[idx]
            if done[idx] == false
                Qₙ += α * maximum(y[:, idx])
            end

            # Adjusting the expected reward for selected move
            Rₙ[argmax(action[:, idx]), idx] = Qₙ
        end
        # Calculate the loss
        agent.criterion(ŷ, copy(Rₙ))
    end

    # Update model weights
    Flux.Optimise.update!(agent.opt, ps, gs)
end


text:
@testset "Construction" begin

    # basic constructor
    xs = rand(5)
    ys = rand(5)
    p = Lagrange(xs, ys)
    @test all(p.(xs) .≈ ys)

    # test weights
    @test all(p.ws .≈ SP.lagrange_barycentric_weights(xs))
    zs = rand(5)
    q = Lagrange(xs, p.ws, zs)
    @test all(q.(xs) .≈ zs)

    # fit is alternate
    q = fit(Lagrange, xs, ys)
    us = rand(5)
    @test all(p.(us) .≈ q.(us))

    # test one, zero, variable
    x = variable()
    q = fit(Lagrange, [1, 2, 3], [2, 3, 1])
    @test iszero(zero(q))
    @test convert(Polynomial, one(q)) == one(x)
    @test convert(Polynomial, variable(q)) == x
end

@testset "Conversion" begin
    xs = [1, 2, 3]
    ys = xs .^ 2
    p = fit(Lagrange, xs, ys)
    q = convert(Polynomial, p)
    @test coeffs(q) == [0, 0, 1]

    pp = fit(Lagrange, xs, q)
    @test coeffs(pp) == coeffs(p)
end

@testset "Evaluation" begin
    xs = rand(5)
    ys = sin.(xs)
    p = Lagrange(xs, ys)
    @test all(p.(xs) == ys)   # shortcuts  for xs
    us = rand(5)
    @test all(abs.(p.(us) - sin.(us)) .<= 1e-1)

    xs, ws = SP.lagrange_barycentric_nodes_weights(Chebyshev, 64)
    ys = sin.(xs)
    p = Lagrange(xs, ws, ys)
    us = 0.1 .+ 0.8 * rand(5) # avoid edges...
    @test all(abs.(sin.(us) - p.(us)) .<= 1e-3)
end

@testset "Arithmetic" begin
    p = fit(Lagrange, [1, 2, 3], x -> x^2)
    q = fit(Lagrange, [2, 3, 4, 5], x -> x^3)

    us = rand(5)

    @test maximum(abs, p.(us) + q.(us) - (p + q).(us)) <= 1e-5
    @test maximum(abs, p.(us) .* q.(us) - (p * q).(us)) <= 1e-2
end

@testset "Elementwise operations" begin
    p = fit(Lagrange, [1, 2, 3], x -> x^2)
    p + 1
    -p
    2p
end

@testset "roots" begin
    q = fit(Lagrange, [2, 3, 4, 5], x -> (x - 1) * (x - 1 / 2) * (x + 1 / 2))
    @test all(roots(q) .≈ [-1 / 2, 1 / 2, 1])
end


text:
@testset "Synchronization" begin
    # TODO: Remove dummy argument
    function synckern(x)
        sync_workgroup()
        nothing
    end

    iob = IOBuffer()
    AMDGPU.code_gcn(iob, synckern, Tuple{Int}; kernel=true)
    @test occursin("s_barrier", String(take!(iob)))
end


text:
using Gadfly

plot([sin, cos], color=["sin", "cos"], 0, 25)


text:
m = 21          # The number of 'space' nodes
n = 401         # The number of 'time' nodes
τ = 10/(n-1)    # Time step size
h = 1/(m-1)     # Space step size

λ = 0.3
α = 0.3+0.1*(1-cos(1))         # Integral of g(x) = 0.3 +0.1sin(x) = 0.3+0.1(1-cos1)
x = vec(linspace(0,1,m))
t = vec(linspace(0,10,n))

α = 0.3+0.1*(1-cos(1))         # Integral of g(x) = 0.3 +0.1sin(x) = 0.3+0.1(1-cos1)
μ = 1./(1-x).^2; μ[1] = 0; μ[m] = 0;
β = zeros(m,m)
γ = (1-λ*τ+h*τ/2)

optv = 0.85 + 0.05cos(2π*t)   # The optimal distribution we are aiming at

## Calcuate A
A = zeros(m,m)

for i=2:m-1
    A[i,i-1] = 1 - τ*μ[i] + τ/h
    A[i,i+1] = 1 - τ*μ[i] - τ/h
end

## Calculate g
g = (0.3+0.1*sin(x))

## Quadrature
function mu(x)
    return 1/(1-x)^2
end

function f(x)
    ## Use just the x part
    ## Remaining cos(πt) should be multiplied later
    return 0.1*sin(1-x)*sin(1-x)
end

function fmu(x)
    return f(x)mu(x)
end

weights = [(18+√30)/36,(18-√30)/36]
point1 = √(3-2(√6)/(√5))/√7
point2 = √(3+2(√6)/(√5))/√7

## Integrate wrt x
fx = weights[1]*fmu((point1+1)/2)+weights[1]*fmu((-point1+1)/2)+weights[2]*fmu((point2+1)/2)+weights[2]*fmu((-point2+1)/2)
fx = fx/2
## Multiply by time vector
fxt = fx*cos(π*t)

## Calculate β
β = zeros(n,n)
β[1,1] = α
for j=2:n
    for k = 1:j-1
       a = (A^(j-1-k))'*μ
       β[j,k] = dot(μ,g)[1]
    end
       β[j,j] = α
end

B = zeros(n,n)
B[1,1] = α
for j=2:n
    for k=1:j
        B[j,:] += γ^(j-k)*β[k,:]
    end
end
B = τ*B
d = zeros(n,1)
for i=1:n
  d[i] = γ^(i-1)*0.9
end

estb = pinv(B)*vec(d-optv)
estv = -B*estb + d
errv = norm(estv-optv,2);
c = estb + fxt - 0.3optv

w = zeros(m,n)
w[1,:] = λ*estv
w[1,1] = 0
for i=2:n
    w[:,i] = A*w[:,i]-τ*estb[i]*g
end
M = ones(m)
p1hat = M'*w*h

#Pkg.add("Gadfly")
#Pkg.add("Cairo")
using Gadfly
using Cairo
draw(PDF("estv.pdf",6inch, 3inch), plot(x=t, y=estv, Guide.XLabel("Time"), Guide.YLabel("Estimated p0"),Guide.XTicks(ticks=[0:1:10])))
draw(PDF("error.pdf",6inch, 3inch), plot(x=t, y=abs(estv-optv), Guide.XLabel("Time"), Guide.YLabel("Absolute Error"),Guide.XTicks(ticks=[0:1:10])))
draw(PDF("estb.pdf",6inch, 3inch), plot(x=t,y=estb, Guide.XLabel("Time"), Guide.YLabel("b(t)"),Guide.XTicks(ticks=[0:1:10])))
draw(PDF("p1hat.pdf",6inch, 3inch), plot(x=t,y=p1hat, Guide.XLabel("Time"), Guide.YLabel("p1hat(t)"), Guide.YTicks(ticks=[0:0.0025:0.0225]), Guide.XTicks(ticks=[0:1:10])))


text:
function Simulation_PlayPause(p_state, method_type, solver_type, n_iterations=-1, n_iterations_log=-1, idx_image=-1, idx_chain=-1)
  ccall((:Simulation_PlayPause, spiritlib),  Void, (Ptr{Void}, Cstring, Cstring, Cint, Cint, Cint, Cint), p_state, method_type, solver_type, n_iterations, n_iterations_log, idx_image, idx_chain)
end

max_stars_repo_name: ddkn/spirit

max_issues_repo_name: ddkn/spirit

max_forks_repo_name: ddkn/spirit


text:
using StochasticDelayDiffEq
using Random
using SparseArrays

function sir_dde!(du,u,h,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    infection = β*c*I/N*S
    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
    Nd = Sd+Id+Rd
    recovery = β*c*Id/Nd*Sd
    @inbounds begin
        du[1] = -infection
        du[2] = infection - recovery
        du[3] = recovery
    end
    nothing
end;

# Define a sparse matrix by making a dense matrix and setting some values as not zero
A = zeros(3,2)
A[1,1] = 1
A[2,1] = 1
A[2,2] = 1
A[3,2] = 1
A = SparseArrays.sparse(A);

# Make `g` write the sparse matrix values
function sir_delayed_noise!(du,u,h,p,t)
    (S,I,R) = u
    (β,c,τ) = p
    N = S+I+R
    infection = β*c*I/N*S
    (Sd,Id,Rd) = h(p, t-τ) # Time delayed variables
    Nd = Sd+Id+Rd
    recovery = β*c*Id/Nd*Sd
    du[1,1] = -sqrt(infection)
    du[2,1] = sqrt(infection)
    du[2,2] = -sqrt(recovery)
    du[3,2] = sqrt(recovery)
end;

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[2]
end;
function affect!(integrator)
  integrator.u[2] = 0.0
end;
cb = ContinuousCallback(condition,affect!);

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
t = 0.0:δt:tmax;
u0 = [990.0,10.0,0.0]; # S,I,R

function sir_history(p, t)
    [1000.0, 0.0, 0.0]
end;

p = [0.05,10.0,4.0]; # β,c,τ
Random.seed!(1234);

prob_sdde = SDDEProblem(sir_dde!,sir_delayed_noise!,u0,sir_history,tspan,p;noise_rate_prototype=A);
sol_sdde = solve(prob_sdde,LambaEM(),callback=cb);


text:
"""
"""
struct VCycle <: AbstractCycle end

""
function _core_cycle!(problem, level, ::VCycle)
    process_fine!(problem, level)
    _descend!(problem, level, VCycle())
    process_coarse!(problem, level)
    return nothing
end


text:
function optimize(diag::Union{Diagram,Tuple,AbstractVector}, optlevel = 1; verbose = 0, kwargs...)
    single = false
    diag = deepcopy(diag)
    if diag isa Diagram
        single = true
        diag = [diag,]
    end
    removeOneChildParent!(diag, verbose = verbose)
    removeDuplicatedLeaves!(diag, verbose = verbose)
    return single ? diag[1] : diag
end

"""
    removeOneChildParent!(diags::AbstractVector; verbose = 0)

    remove duplicated nodes such as:  ---> ver4 ---> InteractionId. Leaf will not be touched!
"""
function removeOneChildParent!(diags::AbstractVector; verbose = 0)
    verbose > 0 && println("remove nodes with only one child.")
    for diag in diags
        #deep first search, remove one-child parent from the leaf level first
        removeOneChildParent!(diag.subdiagram)
        #then remove the one-child subdiagram of the current diagram
        for (si, subdiag) in enumerate(diag.subdiagram)
            if length(subdiag.subdiagram) == 1
                subdiag.subdiagram[1].factor *= subdiag.factor
                diag.subdiagram[si] = subdiag.subdiagram[1]
            end
        end
    end
    return diags
end

"""
    removeDuplicatedLeaves!(diags::AbstractVector; verbose = 0)

    remove duplicated nodes such as:  ---> ver4 ---> InteractionId. Leaf will not be touched!
"""
function removeDuplicatedLeaves!(diags::AbstractVector; verbose = 0)
    verbose > 0 && println("remove duplicated leaves.")
    leaves = []
    for diag in diags
        #leaves must be the propagators
        append!(leaves, collect(Leaves(diag)))
    end
    # println([d.hash for d in leaves])
    sort!(leaves, by = x -> x.hash) #sort the hash of the leaves in an asscend order
    unique!(x -> x.hash, leaves) #filter out the leaves with the same hash number

    for l in leaves
        #make sure all leaves are either Green's functions or interactions
        @assert l.id isa PropagatorId
    end

    function uniqueLeaves(_diags::AbstractVector)
        ############### find the unique Leaves #####################
        uniqueDiag = []
        mapping = Dict{Int,Any}()
        for diag in _diags
            flag = true
            for (ei, e) in enumerate(uniqueDiag)
                if e.factor ≈ diag.factor && e.id == diag.id
                    mapping[diag.hash] = e
                    flag = false
                    break
                end
            end
            if flag
                push!(uniqueDiag, diag)
                # push!(mapping, length(uniqueDiag))
                mapping[diag.hash] = diag
            end
        end
        return uniqueDiag, mapping
    end

    # println(leaves)
    green = [l for l in leaves if l.id isa BareGreenId]
    interaction = [l for l in leaves if l.id isa BareInteractionId]
    greenN = [l for l in leaves if l.id isa BareGreenNId]
    hopping = [l for l in leaves if l.id isa BareHoppingId]
    # println(green)
    # println(interaction)

    uniqueGreen, greenMap = uniqueLeaves(green)
    uniqueGreenN, greenNMap = uniqueLeaves(greenN)
    uniqueInteraction, interactionMap = uniqueLeaves(interaction)
    uniqueHopping, hoppingMap = uniqueLeaves(hopping)
    # println(uniqueInteraction)
    # display(greenMap)

    verbose > 0 && length(green) > 0 && println("Number of independent Greens $(length(green)) → $(length(uniqueGreen))")
    verbose > 0 && length(greenN) > 0 && println("Number of independent GreenNs $(length(greenN)) → $(length(uniqueGreenN))")
    verbose > 0 && length(interaction) > 0 && println("Number of independent Interactions $(length(interaction)) → $(length(uniqueInteraction))")
    verbose > 0 && length(hopping) > 0 && println("Number of independent Hopping $(length(hopping)) → $(length(uniqueHopping))")

    for diag in diags
        for n in PreOrderDFS(diag)
            for (si, subdiag) in enumerate(n.subdiagram)
                @assert (n.id isa PropagatorId) == false "the diagram $n with subdiagrams cannot be a proapgator!"

                if subdiag.id isa PropagatorId
                    if subdiag.id isa BareGreenId
                        n.subdiagram[si] = greenMap[subdiag.hash]
                    elseif subdiag.id isa BareInteractionId
                        n.subdiagram[si] = interactionMap[subdiag.hash]
                    elseif subdiag.id isa BareGreenNId
                        n.subdiagram[si] = greenNMap[subdiag.hash]
                    elseif subdiag.id isa BareHoppingId
                        n.subdiagram[si] = hoppingMap[subdiag.hash]
                    else
                        error("not implemented!")
                    end
                end
            end
        end
    end

    return uniqueGreen, uniqueInteraction
    # return diags
end


text:
using PooksoftOptionsKit
using PyPlot
using Dates
using Distributions

# build -
function _setup(path_to_config_file::String)

    # initialize -
    lattice_model = nothing
    asset_set = nothing

    # load the lattice model -
    result = build_simulation_lattice_data_structure(path_to_config_file)
    if (typeof(result.value) == PSError)
        return result
    else
        lattice_model = result.value
    end

    # load the asset set -
    result = build_simulation_asset_set(path_to_config_file)
    if (typeof(result.value) == PSError)
        return result
    else
        asset_set = result.value
    end

    # Setup tmp structure -
    named_tuple = (model=lattice_model,assets=asset_set)

    # return -
    return PSResult(named_tuple)
end

# setup -
path_to_design_file = "$(pwd())/tests/config/Call.json"
result = _setup(path_to_design_file)
problem_setup = nothing
if (isa(result.value,Exception) == false)
    problem_setup = result.value
end

# call to the greeks -
underlyingAssetPrice = 54.17
delta_value = 0.0
result = delta(problem_setup.assets,problem_setup.model,underlyingAssetPrice)
if (isa(result.value, Exception) == false)
    delta_value = result.value
end

# calculate theta -
theta_value = 0.0
result = theta(problem_setup.assets, problem_setup.model, underlyingAssetPrice)
if (isa(result.value, Exception) == false)
    theta_value = result.value
end

# calculate vega -
vega_value = 0.0
result = vega(problem_setup.assets, problem_setup.model, underlyingAssetPrice)
if (isa(result.value, Exception) == false)
    vega_value = result.value
end

# calculate rho -
rho_value = 0.0
result = rho(problem_setup.assets, problem_setup.model, underlyingAssetPrice)
if (isa(result.value, Exception) == false)
    rho_value = result.value
end


text:
function get_pkg_artifact_dir(pkg::Module, target::String)
    afts = first(values(load_artifacts_toml(find_artifacts_toml(pathof(pkg)))))
    target_arch, target_os, target_libc = get_arch_os_libc(target)
    candidates = Dict[]
    for info in afts
        if info isa Dict
            arch = get(info, "arch", "")
            os = get(info, "os", "")
            libc = get(info, "libc", "")
            if arch == target_arch && os == target_os && libc == target_libc
                push!(candidates, info)
            end
        else
            # this could be an "Any"-platform JLL package
            push!(candidates, afts)
            break
        end
    end
    isempty(candidates) && return ""
    length(candidates) > 1 && @warn "found more than one candidate artifacts, only use the first one: $(first(candidates))"
    info = first(candidates)
    download_info = info["download"][]
    id, url, chk = info["git-tree-sha1"], download_info["url"], download_info["sha256"]
    download_artifact(Base.SHA1(id), url, chk)
    return normpath(artifact_path(Base.SHA1(id)))
end

function get_pkg_include_dir(pkg::Module, target::String)
    artifact_dir = get_pkg_artifact_dir(pkg, target)
    incdir = isempty(artifact_dir) ? "" : joinpath(artifact_dir, "include")
    return incdir
end


text:
using OrdinaryDiffEq, DiffEqCallbacks, Test

f2(u,p,t) = 1.01*u
u0=1/2
tspan = (0.0,5.0)
prob = ODEProblem(f2,u0,tspan)

sol = solve(prob,Tsit5())
@test maximum(diff(sol.t)) > 0.46

dtFE(u,p,t) = 1/20
cb = StepsizeLimiter(dtFE)
sol = solve(prob,Tsit5(),callback=cb)
@test maximum(diff(sol.t)) < 0.046001

cb = StepsizeLimiter(dtFE;safety_factor=1)
sol = solve(prob,Tsit5(),callback=cb)
@test maximum(diff(sol.t)) < 0.050001

cb = StepsizeLimiter(dtFE;safety_factor=1)
sol = solve(prob,Tsit5(),callback=cb,adaptive=false,dt=1/10)
@test maximum(diff(sol.t)) < 0.050001

dtFE(u,p,t) = t < 1 ? 1/20 : 1/10
cb = StepsizeLimiter(dtFE;safety_factor=1)
sol = solve(prob,Tsit5(),callback=cb,adaptive=false,dt=1/10)
@test maximum(diff(sol.t)) > 0.050001
@test maximum(diff(sol.t)) < 0.100001

sol = solve(prob,Tsit5(),callback=cb,adaptive=false,dt=1/20)
@test maximum(diff(sol.t)) < 0.050001

dtFE(u,p,t) = t + 0.2
cb = StepsizeLimiter(dtFE;safety_factor=1,max_step=true)
sol = solve(prob,Tsit5(),callback=cb,adaptive=false,dt=1/20)
@test diff(sol.t) ≈ [0.2,0.4,0.8,1.6,2.0]


text:
module RBNF
using MLStyle
using DataStructures
using PrettyPrint

include("Token.jl")
include("Lexer.jl")
include("Lens.jl")
include("NaiveParserc.jl")
include("ParserGen.jl")

PrettyPrint.pp_impl(io, tk::Token{T}, indent::Int) where T = begin
    s = "Token{$T}(str=$(tk.str), lineno=$(tk.lineno), colno=$(tk.lineno))"
    print(io, s)
    return length(s) + indent
end
end # module


text:
"""
Table Schema validators
https://github.com/frictionlessdata/tableschema-jl#validators
"""
function checkrow(f::Field, val, column::Array=[])
    c::Constraints = f.constraints

    c.required && (val == "" || val == nothing) &&
        throw(ConstraintError("required", f, val))

    c.unique && length(findall(in([val]), column)) > 1 &&
        throw(ConstraintError("unique", f, val))

    if typeof(val) == String

        c.minLength !== nothing && c.minLength > -1 &&
            length(val) < c.minLength &&
                throw(ConstraintError("minLength", f, val, c.minLength))

        c.maxLength !== nothing && c.maxLength > -1 &&
            length(val) > c.maxLength &&
                throw(ConstraintError("maxLength", f, val, c.maxLength))

    end

    # c.minimum > -1 && val < c.minimum &&
    #     throw(ConstraintError("minimum", f, val, c.minimum))
    #
    # c.maximum > -1 && val > c.maximum &&
    #     throw(ConstraintError("maximum", f, val, c.maximum))

    return true
end


text:
@theme EmacsTheme Dict(
    :name => "Emacs",
    :description => "A theme inspired by Emacs 22.",
    :comments => "Based on the Emacs theme from Pygments.",
    :style => S"bg: f8f8f8",
    :tokens => Dict(
        COMMENT            => S"italic; fg: 008800",
        COMMENT_PREPROC    => S"",
        COMMENT_SPECIAL    => S"bold",

        ERROR              => S"bg: ff0000; fg: ff3300",

        GENERIC_DELETED    => S"fg: a00000",
        GENERIC_EMPH       => S"italic",
        GENERIC_ERROR      => S"fg: ff0000",
        GENERIC_HEADING    => S"bold; fg: 000080",
        GENERIC_INSERTED   => S"fg: 00a000",
        GENERIC_OUTPUT     => S"fg: 888",
        GENERIC_PROMPT     => S"bold; fg: 000080",
        GENERIC_STRONG     => S"bold",
        GENERIC_SUBHEADING => S"bold; fg: 800080",
        GENERIC_TRACEBACK  => S"fg: 04d",

        KEYWORD            => S"bold; fg: aa22ff",
        KEYWORD_PSEUDO     => S"",
        KEYWORD_TYPE       => S"bold; fg: 00bb00",

        NAME               => S"",
        NAME_ATTRIBUTE     => S"fg: bb4444",
        NAME_BUILTIN       => S"fg: aa22ff",
        NAME_CLASS         => S"fg: 0000ff",
        NAME_CONSTANT      => S"fg: 880000",
        NAME_DECORATOR     => S"fg: aa22ff",
        NAME_ENTITY        => S"bold; fg: 999999",
        NAME_EXCEPTION     => S"bold; fg: d2413a",
        NAME_FUNCTION      => S"fg: 00a000",
        NAME_LABEL         => S"fg: a0a000",
        NAME_NAMESPACE     => S"bold; fg: 0000ff",
        NAME_TAG           => S"bold; fg: 008000",
        NAME_VARIABLE      => S"fg: b8860b",

        NUMBER             => S"fg: 666666",

        OPERATOR           => S"fg: 666666",
        OPERATOR_WORD      => S"bold; fg: aa22ff",

        STRING             => S"fg: bb4444",
        STRING_DOC         => S"italic",
        STRING_ESCAPE      => S"bold; fg: bb6622",
        STRING_INTERPOL    => S"bold; fg: bb6688",
        STRING_OTHER       => S"fg: 008000",
        STRING_REGEX       => S"fg: bb6688",
        STRING_SYMBOL      => S"fg: b8860b",

        TEXT               => S"",

        WHITESPACE         => S"fg: bbbbbb",
    )
)


text:
using Pkg; Pkg.activate(pwd())
using Revise
using Test, LinearAlgebra, Random, StaticArrays
dir = "I:\\My Documents\\00 AFIT\\Research\\Julia Projects\\StrategyGames\\src"
include(joinpath(dir, "ops_build.jl"))
include(joinpath(dir, "cfrops_solve.jl"))

Random.seed!(579843)
g = AirDefenseGame()
gs = GameSet(g)
ni, ni_stage, na_stage, players = gs.ni, gs.ni_stage, gs.na_stage, gs.players
r = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # regret
s = [zeros(numactions(i, ni_stage, na_stage)) for i in 1:ni] # cumulative strategies
u = 0.0
T = 100
u1 = Vector{Float64}(undef, T)

# Benchmark to find slow functions
hs = SVector(11, 14, 5, 4, 6, 110)
h = [11, 14, 5, 4, 6]
h2 = [11, 14, 5, 4, 6]
@btime vcat([1,2,3], 4) # 2.4 μs # replace with constant size h vectors and keep track of depth? (drops to ~50 ns)
@btime hs2 = setindex(hs, 10, 4) # 30 ns
@btime leafutility(hs, 1, 0.9, 0.7, $g, $gs) # 800 ns
@btime terminal($h) # 2 ns
@btime getplayer($h) # 2 ns
@btime actions($h, $gs.A) # 2 ns
@btime chance([1, 1, 3], 4, $g, $gs) # 500 ns
@btime sample(pweights([0.1, 0.2, 0.3, 0.4])) # 50 ns
@btime infoset($h, $g, $gs) # 200 ns
@btime matchregret(1, $r, $gs) # 220 ns
@btime Vector{Float64}(undef, 100) # 40 ns
@btime U = gs.U # 35 ns

# Benchmark full CFR algorithm
# with max performance lefutility (rand) we have 26 ms => 28 sec per 1k iterations on default AirDefenseGame
# with leafutility = leafutilty_sub, we have 0.08 sec => 65 sec per 1k iterations on default AirDefenseGame
# with static h node and leafutility_sub, we have 11 ms => 11 sec per 1k iterations on default AirDefenseGame
Random.seed!(579843)
g = AirDefenseGame()
gs = GameSet(g)
T = 10_000
cfr(1, g, gs)
@btime cfr(1, g, gs)
@time u, r, s, σ = cfr(T, g, gs)
sum(u) == -83.9854058521706 # check correctness

# push h: checksum passes, 38 ms per iter, 4.4 sec per 100 iter
# static h: checksum passes exactly, 11 ms per iter, 1.1 sec per 100 iter,

# Benchmark fixes to leafutility
function leafutility_rand(h::Array{Int64}, player::Int64, πi::Float64, πo::Float64, gs::GameSet)
    # U, z = gs.U, gs.z
    # ind = findfirst(x -> x == h, z)
    # u = player == 1 ? U[ind] * πo : -1.0 * U[ind] * πo
    val = rand()
    u = player == 1 ? -val : val
    return u
end

function leafutility_lookup(h::Array{Int64}, player::Int64, πi::Float64, πo::Float64, gs::GameSet)
    U, z = gs.U, gs.z
    ind = findfirst(x -> x == h, z)
    u = player == 1 ? U[ind] * πo : -1.0 * U[ind] * πo
    return u
end

function leafutility_rawcalc(h::Array{Int64}, player::Int64, πi::Float64, πo::Float64, g::AirDefenseGame, gs::GameSet)
    c1, c2, ac, c4, dc, ap = h
    attacked_cities = @view gs.A[6][h[6]][:]
    defended_iads = getindex(g.iads, getindex(gs.A[1][h[1]], gs.A[5][h[5]]))
    attacked_iads = getindex(g.iads, getindex(gs.A[2][h[2]], gs.A[3][h[3]]))
    leaf_value = 0.0
    for city in attacked_cities
        q = 1.0
        for defender in g.iads
            if gs.C[defender, city] == 0
                q *= (1 - 0.0)
            elseif defender ∈ attacked_iads && defender ∉ defended_iads
                q *= (1 - 0.0)
            elseif defender ∉ attacked_iads
                q *= (1 - g.pp)
            elseif defender ∈ attacked_iads && defender ∈ defended_iads
                q *= (1 - g.pp * (g.pdc * g.pac)) # if a defender is cyber attacked, it has pdc * pac of being cyver-defended
            else
                @warn "unexpected defense layout"
            end
        end
        leaf_value += -1.0 * g.v[city] * q # defender's expected utility for a city
    end
    u = player == 1 ? leaf_value * πo : -1.0 * leaf_value * πo
    return u
end

# this is the version used in the actual leafutility function
function leafutility_sub(h::Array{Int64}, player::Int64, πi::Float64, πo::Float64, g::AirDefenseGame, gs::GameSet)
    c1, c2, ac, c4, dc, ap = h
    attacked_cities = @view gs.A[6][h[6]][:]
    defended_iads = getindex(g.iads, getindex(gs.A[1][h[1]], gs.A[5][h[5]]))
    attacked_iads = getindex(g.iads, getindex(gs.A[2][h[2]], gs.A[3][h[3]]))
    leaf_value = utility_defender(g, gs, attacked_cities, defended_iads, attacked_iads)
    u = player == 1 ? leaf_value * πo : -1.0 * leaf_value * πo
    return u
end

@btime leafutility_rand([11, 14, 5, 4, 6, 110], 1, 0.9, 0.7, $gs) # 40 ns
@btime leafutility_lookup([11, 14, 5, 4, 6, 9], 1, 0.9, 0.7, $gs) # 66 ms
@btime leafutility_rawcalc([11, 14, 5, 4, 6, 9], 1, 0.9, 0.7, $g, $gs) # 820 ns
@btime leafutility_sub([11, 14, 5, 4, 6, 9], 1, 0.9, 0.7, $g, $gs) # 850 ns

max_stars_repo_name: ajkeith/Cyber-Air-Defense

max_issues_repo_name: ajkeith/Cyber-Air-Defense


text:
module MatrixUtils

using LinearAlgebra, Statistics, DataFrames

export cov2corr,
    make_symmetric!,
    maxind,
    eye,
    eye_diag,
    matrix_inds,
    make_upper_triangular!,
    pivot_cov,
    pivot_mean,
    findnans,
    replace_nans!,
    permutation_matrix,
    standardize_data!,
    standardize!,
    log!,
    issquare,
    randpd,
    missing_var,
    missing_cov,
    missing_mean

function cov2corr(Σ::Matrix{Float64})
    n, p = size(Σ)
    @assert n == p
    P = zeros(n, n)
    for i = 1:n
        for j = i:n
            P[i, j] = Σ[i, j] / sqrt(Σ[i, i] * Σ[j, j])
            P[j, i] = P[i, j]
        end
    end
    return P
end

function cov2corr(Σ::Diagonal{Float64})
    p = size(Σ, 1)
    return Diagonal(ones(p))
end

function make_symmetric!(X::Matrix{T}; use_upper::Bool = true) where T <: Any
    n, m = size(X)
    @assert n == m
    for i = 1:(n - 1)
        for j = (i + 1):n
            if use_upper
                X[j, i] = X[i, j]
            else
                X[i, j] = X[j, i]
            end
        end
    end
end

function maxind(X::AbstractArray{<:Real, 1})
    n = length(X)
    m = 1
    Xm = X[m]
    for i = 2:n
        if X[i] > Xm
            m = i
            Xm = X[m]
        end
    end
    return m
end

function maxind(X::AbstractArray{<:Real, 2})
    ind = maxind(vec(X))
    m, n = size(X)
    col, r = divrem(ind - 1, m)
    col += 1
    row = mod1(r, m)
    return (row, col)
end

function eye(n::Int)
    return diagm(0 => ones(n))
end

function eye_diag(n::Int)
    return Diagonal(ones(n))
end

function matrix_inds(start::Int, stop::Int; diagonals::Bool = true)

    n = stop - start + 1
    p = Int(round(sqrt(n)))
    @assert p^2 == n
    q = div(p * (p + 1), 2)
    if !diagonals
        q = q - p
    end
    inds = zeros(Int, q)
    counter = 0
    for i = 1:p
        r = i:p
        if !diagonals
            r = (i + 1):p
        end
        for j in r
            counter += 1
            v_dim = (i - 1) * p + j
            inds[counter] = v_dim + start - 1
        end
    end
    return inds
end

function matrix_inds(n::Int; diagonals::Bool = true)
    return matrix_inds(1, n^2, diagonals = diagonals)
end

function make_upper_triangular!(X::Matrix{Float64};
        positive_diagonals::Bool = false)
    n, p = size(X)
    for i = 2:n
        for j = 1:(i - 1)
            X[i, j] = 0.0
        end
    end
    if positive_diagonals
        for i = 1:min(n, p)
            X[i, i] = abs(X[i, i])
        end
    end
end

function pivot_cov(M::Hermitian{Float64}, dims::Vector{Int};
        obs_dims::Vector{Int} = setdiff(1:size(M, 1), dims))

    Mmm = M[dims, dims]
    Mmo = M[dims, obs_dims]
    Moo = M[obs_dims, obs_dims]
    return Mmm - Mmo * inv(Moo) * Mmo'
end

function pivot_mean(Y::AbstractArray{Float64, 1}, M::Hermitian{Float64},
        μ::Vector{Float64}, dims::Vector{Int};
        obs_dims::Vector{Int} = setdiff(1:size(M, 1), dims))

    μm = μ[dims]
    Mmo = M[dims, obs_dims]
    Moo = M[obs_dims, obs_dims]
    deviation = Y[obs_dims] - μ[obs_dims]
    return μm + Mmo * (inv(Moo) * deviation)
end

function findnans(x::AbstractArray{Float64}; reverse::Bool = false)
    if reverse
        return findall(y -> !isnan(y), x)
    else
        return findall(y -> isnan(y), x)
    end
end

function permutation_matrix(x::Vector{Int})
    n = length(x)
    Q = zeros(n, n)
    for i = 1:n
        Q[i, x[i]] = 1.0
    end
    return Q
end

function replace_nans!(x::AbstractArray{Float64}; replacement::Float64 = 0.0)
    n = length(x)
    for i = 1:n
        if isnan(x[i])
            x[i] = replacement
        end
    end
end

function standardize_data!(data::Matrix{Float64})
    n, p = size(data)
    ns = zeros(Int, p)
    sums = zeros(p)
    sum_squares = zeros(p)

    for i = 1:n
        for j = 1:p
            x = data[i, j]
            if !isnan(x)
                ns[j] += 1
                sums[j] += x
                sum_squares[j] += x^2
            end
        end
    end

    means = sums ./ ns
    sds = sqrt.(sum_squares ./ ns - means.^2)

    for j = 1:p
        data[:, j] .-= means[j]
        data[:, j] ./= sds[j]
    end
end

function standardize!(x::Vector{Float64})
    μ = 0.0
    σ2 = 0.0
    n_obs = 0
    for i = 1:length(x)
        if !isnan(x[i])
            μ += x[i]
            σ2 += x[i] * x[i]
            n_obs += 1
        end
    end
    if n_obs < 2
        error("Unable to standardize. The number of non-missing elements is $n_obs.")
    end
    μ /= n_obs
    σ2 /= n_obs
    sd = sqrt(σ2 - μ * μ)
    for i = 1:length(x)
        x[i] = (x[i] - μ) / sd
    end
    return x
end

function log!(x::T) where T <: AbstractArray{Float64}
    for i = 1:length(x)
        x[i] = log(x[i])
    end
    return x
end

function log_standardize!(x::T) where T <: AbstractArray{Float64}
    return standardize!(log!(x))
end

function issquare(x::AbstractArray{T, 2}) where T <: Any
    n, m = size(x)
    return n == m
end

function missing_mean(x::AbstractMatrix{Union{Missing, Float64}})
    missing_mean(convert_to_float(x))
end

function missing_mean(x::AbstractMatrix{Float64})
    n, p = size(x)
    μ = zeros(p)
    ns = zeros(Int, p)
    for j = 1:p
        for i = 1:n
            if !isnan(x[i, j])
                μ[j] += x[i, j]
                ns[j] += 1
            end
        end
    end

    return μ ./ ns
end

function missing_var(x::AbstractArray{Union{Float64, Missing}})
    return missing_var(convert_to_float(x))
end

function missing_var(x::AbstractArray{Float64})
    μ = 0.0
    σ2 = 0.0

    n = length(x)
    for i = 1:n
        if !isnan(x[i])
            μ += x[i]
            σ2 += x[i] * x[i]
        end
    end

    μ /= n
    σ2 /= n
    σ2 -= μ * μ

    return σ2
end

function missing_cov(x::AbstractArray{Union{Float64, Missing}})
    return missing_cov(convert_to_float(x))
end

function missing_cov(X::Matrix{Float64})
    n, p = size(X)
    counts = zeros(Int, p, p)
    Σ = zeros(p, p)
    μ = zeros(p, p)
    V = zeros(p, p)
    # Σ_row = zeros(p, p)
    # Σ_col = zeros(p, p)
    # μ = zeros(p)
    # μ_row = zeros(p, p)
    # μ_col = zeros(p, p)
    for i = 1:n
        for j = 1:p
            xj = X[i, j]
            if !isnan(xj)
                Σ[j, j] += xj * xj
                μ[j, j] += xj
                counts[j, j] += 1

                for k = (j + 1):p
                    xk = X[i, k]
                    if !isnan(xk)
                        Σ[j, k] += xj * xk
                        V[j, k] += xj * xj
                        V[k, j] += xk * xk
                        μ[j, k] += xj
                        μ[k, j] += xk
                        counts[j, k] += 1
                    end
                end
            end
        end
    end

    make_symmetric!(counts)
    μ ./= counts
    Σ ./= counts
    V ./= counts

    for i = 1:p
        Σ[i, i] -= μ[i, i] * μ[i, i]
        for j = (i + 1):p
            Σ[i, j] -= μ[i, j] * μ[j, i]
            Σ[j, i] = Σ[i, j]
            V[i, j] -= μ[i, j] * μ[i, j]
            V[j, i] -= μ[j, i] * μ[j, i]
        end
    end

    return Σ, V
end

function missing_cor(X::Matrix{Float64})
    Σ, V = missing_cov(X)
    p = size(Σ, 1)
    for i = 1:p
        Σ[i, i] = 1.0
        for j = (i + 1):p
            Σ[i, j] /= sqrt(V[i, j] * V[j, i])
            Σ[j, i] = Σ[i, j]
        end
    end
    return Σ
end

function convert_to_float(X::Matrix{Float64})
    return X
end

function convert_to_float(X::Matrix{Union{Missing, Float64}})
    s = size(X)
    Y = zeros(s)
    for i = 1:length(X)
        if ismissing(X[i])
            Y[i] = NaN
        else
            Y[i] = X[i]
        end
    end
    return Y
end

function missing_cor(df::DataFrame)
    M = Matrix(df)
    M = convert_to_float(M)
    return missing_cor(M)
end

function randpd(n::Int; rescale::Bool = true)
    X = randn(n, n)
    M = X' * X
    if rescale
        diags = diag(M)
        v = var(diags, corrected=false)
        sd = sqrt(v)
        M ./= sd
    end
    return Symmetric(M)
end

end


text:
module Jive

export Mock, @mockup
include("mockup.jl")

export @skip
include("skip.jl")

export @onlyonce
include("onlyonce.jl")

export @If
include("If.jl")

export @useinside
include("useinside.jl")

export @__END__
include("__END__.jl")

export runtests
include("runtests.jl")

export watch
include("watch.jl")

end # module Jive


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for parse

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

@testset "GraphIO" begin
   ###
   # loadgraph
   ###
   g = loadgraph("testgraph.txt")
   @test nv(g) == 10
   @test ne(g) == 28

   @test encode(g) == collect(1 : 10)

   @test names(g.vdata) == [:name, :age]
   @test eltypes(g.vdata) == [String, Int]
   @test size(g.vdata) == (10, 2)

   @test names(g.edata) == [:relationship]
   @test eltypes(g.edata) == [String]
   @test size(g.edata) == (28,1)

   ###
   # storegraph
   ###
   g = completegraph(10)
   # Vertex Properties
   setvprop!(g, :, rand(Int, 10), :p1)
   setvprop!(g, :, rand(10), :p2)
   setvprop!(g, :, [randstring() for i in 1:10], :p3)
   setvprop!(g, :, rand(Bool, 10), :p4)
   setvprop!(g, :, rand(Char, 10), :p5)

   # Edge properties
   seteprop!(g, :, rand(Int,90), :p1)
   seteprop!(g, :, rand(90), :p2)
   seteprop!(g, :, [randstring() for i in 1:90], :p3)
   seteprop!(g, :, rand(Bool, 90), :p4)
   seteprop!(g, :, rand(Char, 90), :p5)

   storegraph(g, "tmp.txt")

   h = loadgraph("tmp.txt")

   @test isequal(g, h)

   rm("tmp.txt")
end


text:
using DelayEmbeddings, Statistics, TimeseriesSurrogates

"""
    construct_candidate_variables(
        source::Vector{AbstractVector},
        target::Vector{AbstractVector},
        [cond::Vector{AbstractVector}];
        k::Int = 1, include_instantaneous = true,
        τexclude::Union{Int, Nothing} = nothing,
        maxlag::Union{Int, Float64} = 0.05
    ) → ([τs_source, τs_target, τs_cond, ks_targetfuture], [js_source, js_target, js_cond, js_targetfuture])

Construct candidate variables from input time series. `source` is a vector of equal-length time series
assumed to represent the putative source process. `target` and `cond` are the same, but contains time series
of the target process and of the conditional processes, respectively. `k` is the desired prediction lag.

If `include_instantaneous == true`, then the analysis will also consider instantaneous interactions between
the variables.

If `maxlag` is an integer, `maxlag` is taken as the maximum allowed embedding lag. If `maxlag` is a float,
then the maximum embedding lag is taken as `maximum([length.(source); length.(target); length.(cond)])*maxlag`.

If `τexclude` is an integer, all variables whose embedding lag has absolute value equal to `exclude` will be
excluded.
"""
function construct_candidate_variables(source, target, cond;
        k::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "ac_min",
        maxlag::Union{Int, Float64} = 0.05)

    # Ensure all time series are of the same length.
    Ls = [length.(source); length.(target); length.(cond)]
    @assert all(Ls .== maximum(Ls))

    if maxlag isa Int
        τs = 1:maxlag
    else
        τs = 1:ceil(Int, maximum(Ls)*maxlag)
    end

    # Find the maximum allowed embedding lag for each of the candidates.
    τsmax_source = [estimate_delay(s, method_delay, τs) for s in source]
    τsmax_target = [estimate_delay(t, method_delay, τs) for t in target]
    τsmax_cond = [estimate_delay(c, method_delay, τs) for c in cond]

    # The set of candidate variables generated from the target
    # time series must not overlap with the prediction variables,
    # so the k-th lag variable is never included in the candidate set.
    τs_target = [[0:-1:-τ...,] for τ in τsmax_target]

    # Source and conditional variables may have instantaneous
    # interactions with target, so include lag `k` (if desired)
    # in the set of candidate variables.
    if include_instantaneous
        τs_source = [[k; [0:-1:-τ...,]] for τ in τsmax_source]
        τs_cond = [[k; [0:-1:-τ...,]] for τ in τsmax_cond]
    else
        τs_source = [[0:-1:-τ...,] for τ in τsmax_source]
        τs_cond = [[0:-1:-τ...,] for τ in τsmax_cond]
    end

    τs = [τs_source..., τs_target..., τs_cond...]

    # Embedding variables
    js = [[i for x in 1:length(τs[i])] for i = 1:length(τs)]
    js_targetfuture = [i for i in length(τs_source)+1:length(τs_source)+length(τs_target)]

    # Prediction variables
    ks_targetfuture = [k for i in 1:length(target)]

    # Variable filtering, if desired
    if τexclude isa Int
        τs = [filtered_τs(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
        js = [filtered_js(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
    end
    return [τs..., ks_targetfuture], [js..., js_targetfuture]
end

# Usaully, we use all lags from startlag:-\tau_max to construct variables. In some situations,
# we may want to exclude som of those variables.
function filtered_τs(τs::AbstractVector{Int}, js::AbstractVector{Int}, τexclude::Int)
    [τ for τ in τs if abs(τ) != abs.(τexclude)]
end

function filtered_js(τs::AbstractVector{Int}, js::AbstractVector{Int}, τexclude::Int)
    [j for (τ, j) in zip(τs, js) if abs(τ) != abs.(τexclude)]
end

# source & target variant
function construct_candidate_variables(source, target;
        k::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "ac_min",
        maxlag::Union{Int, Float64} = 0.05)

    # Ensure all time series are of the same length.
    Ls = [length.(source); length.(target)]
    @assert all(Ls .== maximum(Ls))

    if maxlag isa Int
        τs = 1:maxlag
    else
        τs = 1:ceil(Int, maximum(Ls)*maxlag)
    end

    # Find the maximum allowed embedding lag for each of the candidates.
    τsmax_source = [estimate_delay(s, method_delay, τs) for s in source]
    τsmax_target = [estimate_delay(t, method_delay, τs) for t in target]

    # The set of candidate variables generated from the target
    # time series must not overlap with the prediction variables,
    # so the k-th lag variable is never included in the candidate set.
    τs_target = [[0:-1:-τ...,] for τ in τsmax_target]

    # Source variables may have instantaneous
    # interactions with target, so include lag `k` (if desired)
    # in the set of candidate variables.
    if include_instantaneous
        τs_source = [[k; [0:-1:-τ...,]] for τ in τsmax_source]
    else
        τs_source = [[0:-1:-τ...,] for τ in τsmax_source]
    end

    τs = [τs_source..., τs_target...]

    ks_targetfuture = [k for i in 1:length(target)]
    js_targetfuture = [i for i in length(τs_source)+1:length(τs_source)+length(τs_target)]
    τs = [τs_source..., τs_target...,]
    js = [[i for x in 1:length(τs[i])] for i = 1:length(τs)]

    # Variable filtering, if desired
    if τexclude isa Int
        τs = [filtered_τs(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
        js = [filtered_js(τsᵢ, jsᵢ, τexclude) for (τsᵢ, jsᵢ) in zip(τs, js)]
    end

    return [τs..., ks_targetfuture], [js..., js_targetfuture]
end

# source, target & cond variant
function candidate_embedding(source, target, cond;
        η::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "mi_min",
        maxlag::Union{Int, Float64} = 0.05)

    τs, js = construct_candidate_variables(
        source, target, cond,
        k = η,
        τexclude = τexclude,
        method_delay = method_delay,
        maxlag = maxlag,
        include_instantaneous = include_instantaneous)

    # TODO: This is more efficient if not using datasets. Re-do manually.
    data = Dataset([source..., target..., cond...,]...,)
    ℰ = genembed(data, ((τs...)...,), ((js...)...,))

    # Get all variables except the target future (which are the last columns of ℰ)
    n_timeseries = size(ℰ, 2)
    n_timeseries_target = length(target)
    Ω = [ℰ[:, i] for i = 1:n_timeseries - n_timeseries_target]
    Y⁺ = ℰ[:, n_timeseries - n_timeseries_target+1:end]

    # We need to keep track of which variables are from the source, because
    # when computing the final TE, we need a marginal which is 𝒮 \ 𝒮_source.
    # Hence, we need to know which indices in `js` correspond to the source.
    idxs_source = 1:length(source)
    idxs_target = length(source)+1:length(source)+length(target)
    idxs_cond = length(source)+length(target)+1:length(source)+length(target)+length(cond)

    return Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond
end

# source & target variant
function candidate_embedding(source, target;
        η::Int = 1,
        τexclude::Union{Int, Nothing} = nothing,
        include_instantaneous = true,
        method_delay = "mi_min",
        maxlag::Union{Int, Float64} = 0.05)

    τs, js = construct_candidate_variables(
        source, target,
        k = η,
        τexclude = τexclude,
        method_delay = method_delay,
        maxlag = maxlag,
        include_instantaneous = include_instantaneous)

    # TODO: This is more efficient if not using datasets. Re-do manually.
    data = Dataset([source..., target...,]...,)
    ℰ = genembed(data, ((τs...)...,), ((js...)...,))

    # Get all variables except the target future (which are the last columns of ℰ)
    n_timeseries = size(ℰ, 2)
    n_timeseries_target = length(target)
    Ω = [ℰ[:, i] for i = 1:n_timeseries - n_timeseries_target]

    Y⁺ = ℰ[:, n_timeseries - n_timeseries_target+1:end]
    idxs_source = 1:length(source)
    idxs_target = length(source)+1:length(source)+length(target)
    idxs_cond = Int[]

    return Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond
end

function optim_te(Ω, Y⁺, τs, js, idxs_source, idxs_target, idxs_cond, est;
        q = 1, base = 2,
        α = 0.05, nsurr = 100, surr::Surrogate = RandomShuffle())

    τs_comb = [(τs...)...,]
    js_comb = [(js...)...,]

    n_candidate_variables = length(Ω)

    𝒮 = Vector{Vector{Float64}}(undef, 0)
    𝒮_τs = Vector{Int}(undef, 0)
    𝒮_js = Vector{Int}(undef, 0)

    k = 1
    while k <= n_candidate_variables
        n_remaining_candidates = length(Ω)
        CMIs_between_Y⁺_and_candidates = zeros(n_remaining_candidates)

        # At first iteration, only loop through source variable. If no source variable is found that
        # yields significant TE, terminate.
        for i = 1:n_remaining_candidates
            if k == 1 || length(𝒮) == 0
                Cᵢ = Ω[i]
                CMI_Y⁺_Cᵢ =
                    genentropy(Dataset(Y⁺, Dataset(Cᵢ)), est, q = q, base = base) -
                    genentropy(Dataset(Cᵢ), est, q = q, base = base)
            else
                Cᵢ = [Ω[i], 𝒮...]
                CMI_Y⁺_Cᵢ =
                    genentropy(Dataset(Y⁺, Dataset(Cᵢ...,)), est, q = q, base = base) -
                    genentropy(Dataset(Cᵢ...,), est, q = q, base = base)
            end
            CMIs_between_Y⁺_and_candidates[i] = CMI_Y⁺_Cᵢ
        end

        idx = findfirst(x -> x == minimum(CMIs_between_Y⁺_and_candidates), CMIs_between_Y⁺_and_candidates)
        cₖ = Ω[idx]

        # Test the significance of this candidate by using a permutation test. The type of surrogate
        # is given by `surr`, and we will use `nsurr` surrogate realizations.
        CMI_permutations = zeros(nsurr)
        s = surrogenerator(cₖ, surr)

        # If k == 1, no candidates have been selected, so CMI reduces to MI
        if k == 1
            cmiₖ = CMIs_between_Y⁺_and_candidates[idx]

            for i = 1:nsurr
                surr_cₖ = s() # Surrogate version of cₖ
                CMI_permutations[i] = mutualinfo(Y⁺, surr_cₖ, est)
            end
        # If k > 1, at least one candidate has been selected, so we compute CMI
        else
            # Precompute terms that do not change during surrogate loop
            H_Y⁺_𝒮 = genentropy(Dataset(Y⁺, Dataset(𝒮...,)), est, q = q, base = base)
            H_𝒮 = genentropy(Dataset(𝒮...), est, q = q, base = base)

            # Original TE
            cmiₖ = H_Y⁺_𝒮 +
                    genentropy(Dataset([cₖ, 𝒮...,]...,), est, q = q, base = base) -
                    genentropy(Dataset(Y⁺, Dataset([cₖ, 𝒮...,]...,)), est, q = q, base = base) -
                    H_𝒮

            for i = 1:nsurr
                surr_cₖ = s() # Surrogate version of cₖ
                CMI_permutations[i] = H_Y⁺_𝒮 +
                    genentropy(Dataset([surr_cₖ, 𝒮...]...,), est, q = q, base = base) -
                    genentropy(Dataset(Y⁺, Dataset([surr_cₖ, 𝒮...]...,)), est, q = q, base = base) -
                    H_𝒮
            end
        end

        # If the candidate passes the significance test, add it to list of selected candidates
        # and remove it from list of remaining candidates.
        if cmiₖ > quantile(CMI_permutations, 1 - α)
            push!(𝒮, cₖ)
            push!(𝒮_τs, τs_comb[idx])
            push!(𝒮_js, js_comb[idx])
            deleteat!(Ω, idx)
            deleteat!(τs_comb, idx)
            deleteat!(js_comb, idx)
            k = k + 1
        # If the candidate does not pass significance test, terminate.
        else
            k = n_candidate_variables + 1
        end
    end

    # No variables were selected at all.
    if length(𝒮) == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    # If no variables were selected from the source process, then TE is not well-defined.
    n_source_vars_picked = count(x -> x ∈ idxs_source, 𝒮_js)
    if n_source_vars_picked == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    # No variables were selected from the target or conditional processes.
    𝒮_nonX = [ts for (ts, j) in zip(𝒮, 𝒮_js) if j ∉ idxs_source]
    if length(𝒮_nonX) == 0
        return 0.0, Int[], Int[], idxs_source, idxs_target, idxs_cond
    end

    CE2 = genentropy(Dataset(Y⁺, Dataset(𝒮...,)), est, base = base, q = q) -
        genentropy(Dataset(𝒮...,), est, base = base, q = q)

    CE1 = genentropy(Dataset(Y⁺, Dataset(𝒮_nonX...,)), est, base = base, q = q) -
        genentropy(Dataset(𝒮_nonX...,), est, base = base, q = q)

    CMI = CE1 - CE2
    return CMI, 𝒮_js, 𝒮_τs, idxs_source, idxs_target, idxs_cond
end

process_input(ts::Vector{T}) where T <: Number = [ts]
process_input(ts::AbstractVector{V}) where V <: Vector{N} where N <: Number = ts
process_input(ts::Dataset) = [columns(ts)...,]


text:
@inline alwaystrue(x) = true

# Matcher patterns with Slot and Segment

# matches one term
# syntax:  ~x
struct Slot{P}
    name::Symbol
    predicate::P
end

Slot(s) = Slot(s, alwaystrue)

Base.isequal(s1::Slot, s2::Slot) = s1.name == s2.name

Base.show(io::IO, s::Slot) = (print(io, "~"); print(io, s.name))

# matches zero or more terms
# syntax: ~~x
struct Segment{F}
    name::Symbol
    predicate::F
end

ismatch(s::Segment, t) = s.predicate(t)

Segment(s) = Segment(s, alwaystrue)

Base.show(io::IO, s::Segment) = (print(io, "~~"); print(io, s.name))

makesegment(s::Symbol, keys) = (push!(keys, s); Segment(s))

function makesegment(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a segment is ~~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Segment($(QuoteNode(name)), $(esc(s.args[2]))))
end

makeslot(s::Symbol, keys) = (push!(keys, s); Slot(s))

function makeslot(s::Expr, keys)
    if !(s.head == :(::))
        error("Syntax for specifying a slot is ~x::\$predicate, where predicate is a boolean function")
    end

    name = s.args[1]

    push!(keys, name)
    :(Slot($(QuoteNode(name)), $(esc(s.args[2]))))
end

function makepattern(expr, keys)
    if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Expr && expr.args[2].args[1] == :(~)
                    # matches ~~x::predicate
                    makesegment(expr.args[2].args[2], keys)
                else
                    # matches ~x::predicate
                    makeslot(expr.args[2], keys)
                end
            else
                :(term($(map(x->makepattern(x, keys), expr.args)...); type=Any))
            end
        else
            error("Unsupported Expr of type $(expr.head) found in pattern")
        end
    else
        # treat as a literal
        return esc(expr)
    end
end

function makeconsequent(expr)
    if expr isa Expr
        if expr.head === :call
            if expr.args[1] === :(~)
                if expr.args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2]))))
                elseif expr.args[2] isa Expr && expr.args[2].args[1] == :(~)
                    @assert expr.args[2].args[2] isa Symbol
                    return :(getindex(__MATCHES__, $(QuoteNode(expr.args[2].args[2]))))
                end
            else
                return Expr(:call, map(makeconsequent, expr.args)...)
            end
        else
            return Expr(expr.head, map(makeconsequent, expr.args)...)
        end
    else
        # treat as a literal
        return esc(expr)
    end
end

abstract type AbstractRule end # Currently doesn't really do anything. Can be removed.

#-----------------------------
#### Regular Rewriting Rules

struct Rule{L, M, R} <: AbstractRule
    expr::Expr               # rule pattern stored for pretty printing
    lhs::L                   # the pattern
    matcher::M               # matcher(lhs)
    rhs::R                   # consequent
    depth::Int               # number of levels of expr this rule touches
end

getdepth(r::Rule) = r.depth

function rule_depth(rule, d=0, maxdepth=0)
    if rule isa Term
        maxdepth = reduce(max, (rule_depth(r, d+1, maxdepth) for r in arguments(rule)), init=1)
    elseif rule isa Slot || rule isa Segment
        maxdepth = max(d, maxdepth)
    end
    return maxdepth
end

function Base.show(io::IO, r::Rule)
    Base.print(io, r.expr)
end

const EMPTY_DICT = ImmutableDict{Symbol, Any}(:____, nothing)

function (r::Rule)(term)
    rhs = r.rhs

    try
        return r.matcher((term,), EMPTY_DICT) do bindings, n
            # n == 1 means that exactly one term of the input (term,) was matched
            n === 1 ? (@timer "RHS" rhs(bindings)) : nothing
        end
    catch err
        throw(RuleRewriteError(r, term))
    end
end

"""
    @rule LHS => RHS

Creates a `Rule` object. A rule object is callable, and  takes an expression and rewrites
it if it matches the LHS pattern to the RHS pattern, returns `nothing` otherwise.
The rule language is described below.

LHS can be any possibly nested function call expression where any of the arugments can
optionally be a Slot (`~x`) or a Segment (`~~x`) (described below).

If an expression matches LHS entirely, then it is rewritten to the pattern in the RHS
Segment (`~x`) and slot variables (`~~x`) on the RHS will substitute the result of the
matches found for these variables in the LHS.

**Slot**:

A Slot variable is written as `~x` and matches a single expression. `x` is the name of the variable. If a slot appears more than once in an LHS expression then expression matched at every such location must be equal (as shown by `isequal`).

_Example:_

Simple rule to turn any `sin` into `cos`:

```julia
julia> @syms a b c
(a, b, c)

julia> r = @rule sin(~x) => cos(~x)
sin(~x) => cos(~x)

julia> r(sin(1+a))
cos((1 + a))
```

A rule with 2 segment variables

```julia
julia> r = @rule sin(~x + ~y) => sin(~x)*cos(~y) + cos(~x)*sin(~y)
sin(~x + ~y) => sin(~x) * cos(~y) + cos(~x) * sin(~y)

julia> r(sin(a + b))
cos(a)*sin(b) + sin(a)*cos(b)
```

A rule that matches two of the same expressions:

```julia
julia> r = @rule sin(~x)^2 + cos(~x)^2 => 1
sin(~x) ^ 2 + cos(~x) ^ 2 => 1

julia> r(sin(2a)^2 + cos(2a)^2)
1

julia> r(sin(2a)^2 + cos(a)^2)
# nothing
```

**Segment**:

A Segment variable is written as `~~x` and matches zero or more expressions in the
function call.

_Example:_

This implements the distributive property of multiplication: `+(~~ys)` matches expressions
like `a + b`, `a+b+c` and so on. On the RHS `~~ys` presents as any old julia array.

```julia
julia> r = @rule ~x * +((~~ys)) => sum(map(y-> ~x * y, ~~ys));

julia> r(2 * (a+b+c))
((2 * a) + (2 * b) + (2 * c))
```

**Predicates**:

Predicates can be used on both `~x` and `~~x` by using the `~x::f` or `~~x::f`.
Here `f` can be any julia function. In the case of a slot the function gets a single
matched subexpression, in the case of segment, it gets an array of matched expressions.

The predicate should return `true` if the current match is acceptable, and `false`
otherwise.

```julia
julia> two_πs(x::Number) = abs(round(x/(2π)) - x/(2π)) < 10^-9
two_πs (generic function with 1 method)

julia> two_πs(x) = false
two_πs (generic function with 2 methods)

julia> r = @rule sin(~~x + ~y::two_πs + ~~z) => sin(+(~~x..., ~~z...))
sin(~(~x) + ~(y::two_πs) + ~(~z)) => sin(+(~(~x)..., ~(~z)...))

julia> r(sin(a+3π))

julia> r(sin(a+6π))
sin(a)

julia> r(sin(a+6π+c))
sin((a + c))
```

Predicate function gets an array of values if attached to a segment variable (`~~x`).

**Context**:

_In predicates_: Contextual predicates are functions wrapped in the `Contextual` type.
The function is called with 2 arguments: the expression and a context object
passed during a call to the Rule object (maybe done by passing a context to `simplify` or
a `RuleSet` object).

The function can use the inputs however it wants, and must return a boolean indicating
whether the predicate holds or not.

_In the consequent pattern_: Use `(@ctx)` to access the context object on the right hand side
of an expression.
"""
macro rule(expr)
    @assert expr.head == :call && expr.args[1] == :(=>)
    lhs,rhs = expr.args[2], expr.args[3]
    keys = Symbol[]
    lhs_term = makepattern(lhs, keys)
    unique!(keys)
    quote
        $(__source__)
        lhs_pattern = $(lhs_term)
        Rule($(QuoteNode(expr)),
             lhs_pattern,
             matcher(lhs_pattern),
             __MATCHES__ -> $(makeconsequent(rhs)),
             rule_depth($lhs_term))
    end
end

#-----------------------------
#### Associative Commutative Rules

struct ACRule{F,R} <: AbstractRule
    sets::F
    rule::R
    arity::Int
end

Rule(acr::ACRule)   = acr.rule
getdepth(r::ACRule) = getdepth(r.rule)

macro acrule(expr)
    arity = length(expr.args[2].args[2:end])
    quote
        ACRule(permutations, $(esc(:(@rule($(expr))))), $arity)
    end
end

macro ordered_acrule(expr)
    arity = length(expr.args[2].args[2:end])
    quote
        ACRule(combinations, $(esc(:(@rule($(expr))))), $arity)
    end
end

Base.show(io::IO, acr::ACRule) = print(io, "ACRule(", acr.rule, ")")

function (acr::ACRule)(term)
    r = Rule(acr)
    if !istree(term)
        r(term)
    else
        f =  operation(term)
        # Assume that the matcher was formed by closing over a term
        if f != operation(r.lhs) # Maybe offer a fallback if m.term errors.
            return nothing
        end

        T = symtype(term)
        args = arguments(term)

        itr = acr.sets(eachindex(args), acr.arity)

        for inds in itr
            result = r(Term{T}(f, @views args[inds]))
            if !isnothing(result)
                # Assumption: inds are unique
                length(args) == length(inds) && return result
                return similarterm(term, f, [result, (args[i] for i in eachindex(args) if i ∉ inds)...])
            end
        end
    end
end

struct RuleRewriteError
    rule
    expr
end

getdepth(::Any) = typemax(Int)

@noinline function Base.showerror(io::IO, err::RuleRewriteError)
    msg = "Failed to apply rule $(err.rule) on expression "
    msg *= sprint(io->showraw(io, err.expr))
    print(io, msg)
end

function timerewrite(f)
    if !TIMER_OUTPUTS
        error("timerewrite must be called after enabling " *
              "TIMER_OUTPUTS in the main file of this package")
    end
    reset_timer!()
    being_timed[] = true
    x = f()
    being_timed[] = false
    print_timer()
    println()
    x
end

"""
    @timerewrite expr

If `expr` calls `simplify` or a `RuleSet` object, track the amount of time
it spent on applying each rule and pretty print the timing.

This uses [TimerOutputs.jl](https://github.com/KristofferC/TimerOutputs.jl).

## Example:

```julia

julia> expr = foldr(*, rand([a,b,c,d], 100))
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)

julia> @timerewrite simplify(expr)
 ────────────────────────────────────────────────────────────────────────────────────────────────
                                                         Time                   Allocations
                                                 ──────────────────────   ───────────────────────
                Tot / % measured:                     340ms / 15.3%           92.2MiB / 10.8%

 Section                                 ncalls     time   %tot     avg     alloc   %tot      avg
 ────────────────────────────────────────────────────────────────────────────────────────────────
 ACRule((~y) ^ ~n * ~y => (~y) ^ (~n ...    667   11.1ms  21.3%  16.7μs   2.66MiB  26.8%  4.08KiB
   RHS                                       92    277μs  0.53%  3.01μs   14.4KiB  0.14%     160B
 ACRule((~x) ^ ~n * (~x) ^ ~m => (~x)...    575   7.63ms  14.6%  13.3μs   1.83MiB  18.4%  3.26KiB
 (*)(~(~(x::!issortedₑ))) => sort_arg...    831   6.31ms  12.1%  7.59μs    738KiB  7.26%     910B
   RHS                                      164   3.03ms  5.81%  18.5μs    250KiB  2.46%  1.52KiB
   ...
   ...
 ────────────────────────────────────────────────────────────────────────────────────────────────
(a ^ 26) * (b ^ 30) * (c ^ 16) * (d ^ 28)
```
"""
macro timerewrite(expr)
    :(timerewrite(()->$(esc(expr))))
end

Base.@deprecate RuleSet(x) Postwalk(Chain(x))