text:
################################################################################
### Lattice Iteration Cache
################################################################################

#=
Many of the iterators implemented here do some sorting or filtering based on
distance vectors. For them to run efficiently they may need various maps
connecting distance vector indices, source site indices and target site indices.

In order to avoid duplicating data across different iterators we construct a
cache here. Each iterator will fetch the data it needs from the cache or
generate the data if it is not available yet.
=#

abstract type LICacheKeys end
struct Dir2SrcTrg <: LICacheKeys end
struct Src2DirTrg <: LICacheKeys end
struct SrcTrg2Dir <: LICacheKeys end
# remaps dirs: filter -> 1:length(filter) to avoid empty data in output
struct FilteredSrc2IdxTrg{T} <: LICacheKeys
    filter::T
end

struct LatticeIteratorCache
    cache::Dict{Any, Any}
end
LatticeIteratorCache() = LatticeIteratorCache(Dict{Any, Any}())

# for simplicity
function Base.getindex(mc::MonteCarloFlavor, key::LICacheKeys)
    get!(mc.lattice_iterator_cache, key, lattice(mc))
end

Base.getindex(cache::LatticeIteratorCache, key::LICacheKeys) = cache.cache[key]
function Base.get!(cache::LatticeIteratorCache, key::LICacheKeys, lattice::AbstractLattice)
    !haskey(cache, key) && push!(cache, key, lattice)
    cache[key]
end
Base.empty!(cache::LatticeIteratorCache) = empty!(cache.cache)
Base.haskey(cache::LatticeIteratorCache, key::LICacheKeys) = haskey(cache.cache, key)

# For shifting sites across periodic bounds
function generate_combinations(vs::Vector{<: Vector})
    out = [zeros(length(vs[1]))]
    for v in vs
        out = vcat([e.-v for e in out], out, [e.+v for e in out])
    end
    out
end

# norm + ϵ * angle(v, e_x)
function directed_norm(v, ϵ)
    l = norm(v)
    if length(v) == 2 && l > ϵ
        angle = acos(dot([1, 0], v) / l)
        v[2] < 0 && (angle = 2pi - angle)
        return l + ϵ * angle
    else
        return l
    end
end

function Base.push!(cache::LatticeIteratorCache, key::Dir2SrcTrg, lattice::AbstractLattice, ϵ=1e-6)
    if !haskey(cache.cache, key)
        _positions = positions(lattice)
        wrap = generate_combinations(lattice_vectors(lattice))
        directions = Vector{Float64}[]
        # (src, trg), first index is dir, second index irrelevant
        bonds = [Tuple{Int64, Int64}[] for _ in 1:length(lattice)]

        for origin in 1:length(lattice)
            for (trg, p) in enumerate(_positions)
                d = _positions[origin] .- p .+ wrap[1]
                for v in wrap[2:end]
                    new_d = _positions[origin] .- p .+ v
                    if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                        d .= new_d
                    end
                end
                # The rounding will allow us to use == here
                idx = findfirst(dir -> isapprox(dir, d, atol=ϵ), directions)
                if idx === nothing
                    push!(directions, d)
                    if length(bonds) < length(directions)
                        push!(bonds, Tuple{Int64, Int64}[])
                    end
                    push!(bonds[length(directions)], (origin, trg))
                else
                    push!(bonds[idx], (origin, trg))
                end
            end
        end

        temp = sortperm(directions, by = v -> directed_norm(v, ϵ))
        cache.cache[key] = bonds[temp]
    end
    nothing
end

function Base.push!(cache::LatticeIteratorCache, key::SrcTrg2Dir, lattice::AbstractLattice)
    if !haskey(cache.cache, key)
        push!(cache, Dir2SrcTrg(), lattice)
        dir2srctrg = cache[Dir2SrcTrg()]
        srctrg2dir = [-1 for _ in 1:length(lattice), __ in 1:length(lattice)]
        for dir in eachindex(dir2srctrg)
            for (src, trg) in dir2srctrg[dir]
                srctrg2dir[src, trg] = dir
            end
        end
        cache.cache[key] = srctrg2dir
    end
end
function Base.push!(cache::LatticeIteratorCache, key::Src2DirTrg, lattice::AbstractLattice)
    if !haskey(cache.cache, key)
        push!(cache, Dir2SrcTrg(), lattice)
        dir2srctrg = cache[Dir2SrcTrg()]
        trg_from_src = [Tuple{Int64, Int64}[] for _ in 1:length(lattice)]
        for dir in eachindex(dir2srctrg)
            for (src, trg) in dir2srctrg[dir]
                push!(trg_from_src[src], (dir, trg))
            end
        end
        cache.cache[key] = trg_from_src
    end
    nothing
end

function Base.push!(cache::LatticeIteratorCache, key::FilteredSrc2IdxTrg, lattice::AbstractLattice)
    if !haskey(cache.cache, key)
        push!(cache, Src2DirTrg(), lattice)

        dir2idx = [-1 for _ in 1:length(lattice)]
        for (idx, dir) in enumerate(key.filter)
            dir2idx[dir] = idx
        end

        cache.cache[key] = map(cache[Src2DirTrg()]) do dir_trg_list
            [(dir2idx[dir], trg) for (dir, trg) in dir_trg_list if dir in key.filter]
        end
    end
    nothing
end

################################################################################
### Abstract Iterator Types
################################################################################

abstract type AbstractLatticeIterator end
# All indices are sites
abstract type DirectLatticeIterator <: AbstractLatticeIterator end
# first index is a meta index (e.g. direction), rest for sites
abstract type DeferredLatticeIterator <: AbstractLatticeIterator end
# Wraps a lattice iterator to do change what happens with the output
abstract type LatticeIterationWrapper{LI <: AbstractLatticeIterator} <: AbstractLatticeIterator end

abstract type AbstractLatticeIteratorTemplate end
abstract type DirectLatticeIteratorTemplate <: AbstractLatticeIteratorTemplate end
abstract type DeferredLatticeIteratorTemplate <: AbstractLatticeIteratorTemplate end

# is this valid?
function (iter::AbstractLatticeIteratorTemplate)(mc::MonteCarloFlavor, ::Model)
    iter(mc.lattice_iterator_cache, lattice(mc))
end

################################################################################
### RangeIterator (1 index)
################################################################################

"""
    EachSiteAndFlavor()

Creates an iterator template which iterates through the diagonal of the Greensfunction.
"""
struct EachSiteAndFlavor <: DirectLatticeIteratorTemplate end
function (::EachSiteAndFlavor)(mc::MonteCarloFlavor, model::Model)
    RangeIterator(1 : length(lattice(mc)) * nflavors(field(mc)))
end

"""
    EachSite()

Creates an iterator template which iterates through every site of a given lattice.
"""
struct EachSite <: DirectLatticeIteratorTemplate end
(::EachSite)(_, l::AbstractLattice) = RangeIterator(1 : length(l))

struct RangeIterator <: DirectLatticeIterator
    range::UnitRange{Int64}
end
Base.iterate(iter::RangeIterator) = iterate(iter.range)
Base.iterate(iter::RangeIterator, state) = iterate(iter.range, state)
Base.length(iter::RangeIterator) = length(iter.range)
Base.eltype(iter::RangeIterator) = eltype(iter.range)

################################################################################
### OnSite (src, trg)
################################################################################

"""
    OnSite()

Creates an iterator template which iterates through every site of a given lattice,
returning the linear index equivalent to (site, site) at every step.
"""
struct OnSite <: DirectLatticeIteratorTemplate end
(::OnSite)(_, l::AbstractLattice) = _OnSite(length(l))

struct _OnSite <: DirectLatticeIterator
    N::Int64
end

Base.iterate(iter::_OnSite, i=1) = i ≤ iter.N ? ((i, i), i+1) : nothing
Base.length(iter::_OnSite) = iter.N
Base.eltype(::_OnSite) = Tuple{Int64, Int64}

################################################################################
### EachSitePair (src, trg)
################################################################################

"""
    EachSitePair()

Creates an iterator template which returns every pair of sites `(s1, s2)` with
`s1, s2 ∈ 1:Nsites`.
"""
struct EachSitePair <: DirectLatticeIteratorTemplate end
(::EachSitePair)(_, l::AbstractLattice) = _EachSitePair(length(l))

struct _EachSitePair <: DirectLatticeIterator
    N::Int64
end

function Base.iterate(iter::_EachSitePair, i=1)
    i ≤ iter.N^2 ? ((div(i-1, iter.N)+1, mod1(i, iter.N)), i+1) : nothing
end
Base.length(iter::_EachSitePair) = iter.N^2
Base.eltype(::_EachSitePair) = NTuple{2, Int64}

################################################################################
### EachSitePairByDistance (dir, src, trg)
################################################################################

"""
    EachSitePairByDistance()

Creates an iterator template which returns triplets
`(direction index, source, target)` sorted by distance. The `direction index`
identifies each unique direction `position(target) - position(source)`.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachSitePairByDistance <: DeferredLatticeIteratorTemplate end
function (::EachSitePairByDistance)(cache::LatticeIteratorCache, l::AbstractLattice)
    push!(cache, Dir2SrcTrg(), l)
    # this should only keep a reference I think?
    dir2srctrg = cache[Dir2SrcTrg()]
    _EachSitePairByDistance(length(l)^2, dir2srctrg)
end
function ndirections(mc, ::EachSitePairByDistance)
    push!(mc.lattice_iterator_cache, Dir2SrcTrg(), lattice(mc))
    Ndir = length(mc.lattice_iterator_cache[Dir2SrcTrg()])
    (Ndir, )
end

struct _EachSitePairByDistance <: DeferredLatticeIterator
    N::Int64
    dir2srctrg::Vector{Vector{Tuple{Int64, Int64}}}
end

function Base.iterate(iter::_EachSitePairByDistance, state = (1, 1))
    dir, i = state
    # next_dir = dir + div(i, length(iter.pairs[dir]))
    # next_i = mod1(i+1, length(iter.pairs[dir]))
    # if next_dir > length(iter.pairs)
    #     return nothing
    # else
    #     src, trg = iter.pairs[dir][i]
    #     return (dir, src, trg), (next_dir, next_i)
    # end

    if dir ≤ length(iter.dir2srctrg)
        if i ≤ length(iter.dir2srctrg[dir])
            src, trg = iter.dir2srctrg[dir][i]
            return (dir, src, trg), (dir, i+1)
        else
            return iterate(iter, (dir+1, 1))
        end
    else
        return nothing
    end
end
Base.length(iter::_EachSitePairByDistance) = iter.N
Base.eltype(::_EachSitePairByDistance) = NTuple{3, Int64}

################################################################################
### EachLocalQuadByDistance
################################################################################

"""
    EachLocalQuadByDistance([directions])()

Creates an iterator template which yields 4-tuples of sites with 3 directions
inbetween. The output is given as `(combined_dir, src1, trg1, src2, trg2)`
where `combined_dir` is a linear index of `(dir12, dir1, dir2)` corresponding to
the distance vectors `(src2 - src1, trg1 - src1, trg2 - src2)`.

The target-source directions can be limited with `directions`. If an integer is
passed it will include `1:directions` as dir1 and dir2. If a range or vector is
passed it will use that range. Note that if that range or vector is
discontinuous, the output will be remapped to avoid storing unnecessary data.
I.e. `directions = [2, 5, 6]` will pick target sites matching directions 2, 5
and 6, but return 1, 2 and 3 as directions.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachLocalQuadByDistance{T} <: DeferredLatticeIteratorTemplate
    directions::T
    EachLocalQuadByDistance(N::Integer) = EachLocalQuadByDistance(1:N)
    EachLocalQuadByDistance(dirs::T) where T = new{T}(dirs)
end

function EachLocalQuadByDistance{K}(args...) where {K}
    EachLocalQuadByDistance(K)(args...)
end
function Base.:(==)(l::EachLocalQuadByDistance, r::EachLocalQuadByDistance)
    l.directions == r.directions
end
function ndirections(mc, iter::EachLocalQuadByDistance)
    push!(mc.lattice_iterator_cache, Dir2SrcTrg(), lattice(mc))
    Ndir = length(mc.lattice_iterator_cache[Dir2SrcTrg()])
    K = length(iter.directions)
    (Ndir, K, K)
end

struct _EachLocalQuadByDistance <: DeferredLatticeIterator
    mult::Tuple{Int64, Int64}
    N::Int64
    src_mask::Vector{Int64}
    srctrg2dir::Matrix{Int64}
    filtered_src2dirtrg::Vector{Vector{Tuple{Int64, Int64}}}
end
function (config::EachLocalQuadByDistance)(cache::LatticeIteratorCache, lattice::AbstractLattice)
    key = FilteredSrc2IdxTrg(config.directions)
    push!(cache, key, lattice)
    push!(cache, SrcTrg2Dir(), lattice)

    Ndir = length(cache[Dir2SrcTrg()])
    srctrg2dir = cache[SrcTrg2Dir()]
    filtered_src2dirtrg = cache[key]

    src_mask = [i for i in eachindex(filtered_src2dirtrg) if !isempty(filtered_src2dirtrg[i])]
    N = mapreduce(length, +, filtered_src2dirtrg)

    _EachLocalQuadByDistance(
        (Ndir, Ndir * length(config.directions)),
        N^2, src_mask, srctrg2dir, filtered_src2dirtrg
    )
end

function Base.iterate(iter::_EachLocalQuadByDistance)
    src = iter.src_mask[1]
    dir12 = iter.srctrg2dir[src, src]
    dir, trg = iter.filtered_src2dirtrg[src][1]
    combined_dir = dir12 + iter.mult[1] * (dir-1) + iter.mult[2] * (dir-1)
    # state = (src1 mask index, src2 mask index, filter1 index, filter2 index)
    return ((combined_dir, src, trg, src, trg), (1, 1, 1, 1))
end

function Base.iterate(iter::_EachLocalQuadByDistance, state)
    midx1, midx2, fidx1, fidx2 = state
    src1 = iter.src_mask[midx1]
    src2 = iter.src_mask[midx2]

    # for reference:
    # if fidx2 == length(iter.filtered_src2dirtrg[src2])
    #     if fidx1 == length(iter.filtered_src2dirtrg[src1])
    #         if midx2 == length(iter.src_mask)
    #             if midx1 == length(iter.src_mask)
    #                 return nothing
    #             else
    #                 fidx1 = fidx2 = midx2 = 1
    #                 midx1 += 1
    #             end
    #         else
    #             fidx1 = fidx2 = 1
    #             midx2 += 1
    #         end
    #     else
    #         fidx2 = 1
    #         fidx1 += 1
    #     end
    # else
    #     fidx2 += 1
    # end

    # Branchless V2
    b1 = fidx2 == length(iter.filtered_src2dirtrg[src2])
    b2 = b1 && (fidx1 == length(iter.filtered_src2dirtrg[src1]))
    b3 = b2 && (midx2 == length(iter.src_mask))
    b4 = b3 && (midx1 == length(iter.src_mask))
    b4 && return nothing
    fidx2 = Int64(b1 || (fidx2 + 1))
    fidx1 = Int64(b2 || (fidx1 + b1))
    midx2 = Int64(b3 || (midx2 + b2))
    midx1 = Int64(midx1 + b3)

    src1 = iter.src_mask[midx1]
    src2 = iter.src_mask[midx2]
    dir12 = iter.srctrg2dir[src1, src2]
    dir1, trg1 = iter.filtered_src2dirtrg[src1][fidx1]
    dir2, trg2 = iter.filtered_src2dirtrg[src2][fidx2]
    combined_dir = dir12 + iter.mult[1] * (dir1-1) + iter.mult[2] * (dir2-1)
    # state = (src1 mask index, src2 mask index, filter1 index, filter2 index)
    return ((combined_dir, src1, trg1, src2, trg2), (midx1, midx2, fidx1, fidx2))
end

Base.length(iter::_EachLocalQuadByDistance) = iter.N
Base.eltype(::_EachLocalQuadByDistance) = NTuple{5, Int64}

################################################################################
### EachLocalQuadBySyncedDistance
################################################################################

"""
    EachLocalQuadBySyncedDistance(directions)

Creates an iterator template which yields 5-tuples of one combined direction and
4 sites `(combined_dir, src1, trg1, src2, trg2)`. The former is the linear
index associated with `(dir12, dir)` where `dir12` refers to the distance
vector between sources and `dir` to the distance vector between both target
source pairs.

The target-source directions can be limited with `directions`. If an integer is
passed it will include `1:directions` as dir1 and dir2. If a range or vector is
passed it will use that range. Note that if that range or vector is
discontinuous, the output will be remapped to avoid storing unnecessary data.
I.e. `directions = [2, 5, 6]` will pick target sites matching directions 2, 5
and 6, but return 1, 2 and 3 as directions.

Requires `lattice` to implement `positions` and `lattice_vectors`.
"""
struct EachLocalQuadBySyncedDistance{T} <: DeferredLatticeIteratorTemplate
    directions::T
    EachLocalQuadBySyncedDistance(N::Integer) = EachLocalQuadBySyncedDistance(1:N)
    EachLocalQuadBySyncedDistance(dirs::T) where T = new{T}(dirs)
end

function EachLocalQuadBySyncedDistance{K}(args...) where {K}
    EachLocalQuadBySyncedDistance(K)(args...)
end
function Base.:(==)(l::EachLocalQuadBySyncedDistance, r::EachLocalQuadBySyncedDistance)
    l.directions == r.directions
end
function ndirections(mc, iter::EachLocalQuadBySyncedDistance)
    push!(mc.lattice_iterator_cache, Dir2SrcTrg(), lattice(mc))
    Ndir = length(mc.lattice_iterator_cache[Dir2SrcTrg()])
    K = length(iter.directions)
    (Ndir, K)
end

struct _EachLocalQuadBySyncedDistance <: DeferredLatticeIterator
    N::Int64
    Ndir::Int64
    directions::Vector{Int64}
    dir2srctrg::Vector{Vector{Tuple{Int64, Int64}}}
    srctrg2dir::Matrix{Int64}
end

function (config::EachLocalQuadBySyncedDistance)(
        cache::LatticeIteratorCache, lattice::AbstractLattice
    )
    push!(cache, Dir2SrcTrg(), lattice)
    push!(cache, SrcTrg2Dir(), lattice)

    dir2srctrg = cache[Dir2SrcTrg()]
    srctrg2dir = cache[SrcTrg2Dir()]

    # Loop structure
    # for sync_dir in directions
    #     for (src1, trg1) in dir2srctrg[sync_dir]
    # 	    for (src2, trg2) in dir2srctrg[sync_dir]
    # 		    dir12 = srctrg2dir[src1, src2]
    N = mapreduce(dir -> length(dir2srctrg[dir])^2, +, config.directions)

    _EachLocalQuadBySyncedDistance(
        N, length(dir2srctrg), collect(config.directions), dir2srctrg, srctrg2dir
    )
end

function Base.iterate(iter::_EachLocalQuadBySyncedDistance, state = (1,1,0))
    sync_idx, idx1, idx2 = state

    # Branchless increments
    sync_dir = iter.directions[sync_idx]
    N = length(iter.dir2srctrg[sync_dir])

    b1 = idx2 == N
    b2 = b1 && (idx1 == N)
    b3 = b2 && (sync_idx == length(iter.directions))

    b3 && return nothing

    idx2 = Int64(b1 || (idx2 + 1))
    idx1 = Int64(b2 || (idx1 + b1))
    sync_idx = Int64(sync_idx + b2)

    sync_dir = iter.directions[sync_idx]
    src1, trg1 = iter.dir2srctrg[sync_dir][idx1]
    src2, trg2 = iter.dir2srctrg[sync_dir][idx2]
    dir12 = iter.srctrg2dir[src1, src2]
    combined_dir = dir12 + iter.Ndir * (sync_dir - 1)

    # state = (src1 mask index, src2 mask index, filter1 index, filter2 index)
    return ((combined_dir, src1, trg1, src2, trg2), (sync_idx, idx1, idx2))
end

Base.length(iter::_EachLocalQuadBySyncedDistance) = iter.N
Base.eltype(::_EachLocalQuadBySyncedDistance) = NTuple{5, Int64}

################################################################################
### Sum Wrapper
################################################################################

"""
    Sum(iteration_template)

Sum is a wrapper for iteration templates that indicates summation. It generates
`_Sum(iter)` when constructed (constructing whatever template it holds) so that
it can be dispatched on.
"""
struct Sum{T} <: AbstractLatticeIteratorTemplate
    template::T
end
# compat
Sum{T}(args...) where {T} = _Sum(T()(args...))
struct _Sum{LI} <: LatticeIterationWrapper{LI}
    iter::LI
end
function (config::Sum)(cache::LatticeIteratorCache, lattice::AbstractLattice)
    if config.template isa EachSitePairByDistance
        _Sum(EachSitePair()(cache, lattice))
    else
        _Sum(config.template(cache, lattice))
    end
end

Base.iterate(s::_Sum) = iterate(s.iter)
Base.iterate(s::_Sum, state) = iterate(s.iter, state)
Base.length(s::_Sum) = length(s.iter)
Base.eltype(s::_Sum) = eltype(s.iter)

################################################################################
### Symmetry Wrapper
################################################################################

# This is a weird monster :(
# The idea is that we construct a thin wrapper
# li = ApplySymmetries{EachLocalQuadByDistance{5}}(sym1, sym2, ...)
# which contains weights for different neighbors in the given symmetries
# The backend then bundles all of these and constructs one thick wrapper by
# calling
# li(dqmc, model)
# This then actually contains the udnerlying lattice iterator

"""
    ApplySymmetries{lattice_iterator_type}(symmetries...)

`ApplySymmetries` is a wrapper for a `DeferredLatticeIterator`. It is meant to
specify how results from different directions are to be added up.

For example `T = EachLocalQuadByDistance{5}` specifies 4 site tuples where two
sites have set distances between them - one of the first 5 smallest ones. In a
square lattice these would be on-site (1) and the four nearest neighbors (2-5).
We may use `iter = ApplySymmetries{}([1], [0, 1, 1, 1, 1])` to specify how
results in these directions should be added up. The first rule would be s-wave,
the second extended s-wave.
These rules will be applied for DQMCMeasurements during the simulation. I.e.
first, the normal iteration and summation from `EachLocalQuadByDistance` is
performed. After that we have 5 values for each direction. These are then
weighted by each "symmetry" in `ApplySymmetries` to give the final result, saved
in the DQMCMeasurement.
"""
struct ApplySymmetries{IT, N, T} <: AbstractLatticeIteratorTemplate
    template::IT
    symmetries::NTuple{N, Vector{T}}
end
function ApplySymmetries{LI}(symmetries::Vector...) where {LI}
    ApplySymmetries(LI(), symmetries)
end
function ApplySymmetries(template::DeferredLatticeIteratorTemplate, symmetries::Vector...)
    ApplySymmetries(template, symmetries)
end

struct _ApplySymmetries{LI <: DeferredLatticeIterator, N, T} <: LatticeIterationWrapper{LI}
    iter::LI
    symmetries::NTuple{N, T}
end
function (x::ApplySymmetries)(cache::LatticeIteratorCache, lattice::AbstractLattice)
    iter = x.template(cache, lattice)
    _ApplySymmetries(iter, x.symmetries)
end

Base.iterate(s::_ApplySymmetries) = iterate(s.iter)
Base.iterate(s::_ApplySymmetries, state) = iterate(s.iter, state)
Base.length(s::_ApplySymmetries) = length(s.iter)
Base.eltype(s::_ApplySymmetries) = eltype(s.iter)

################################################################################
### Utility (directions)
################################################################################

# function directions(::EachSitePair, lattice::AbstractLattice)
#     pos = positions(lattice)
#     [p2 .- p1 for p2 in pos for p1 in pos]
# end

function directions(iter::_EachSitePairByDistance, lattice::AbstractLattice, ϵ=1e-6)
    pos = MonteCarlo.positions(lattice)
    wrap = generate_combinations(lattice_vectors(lattice))

    map(iter.dir2srctrg) do pairs
        src, trg = pairs[1]
        _d = pos[src] - pos[trg]
        # Find lowest distance w/ periodic bounds
        d = _d .+ wrap[1]
        for v in wrap[2:end]
            new_d = _d .+ v
            if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                d .= new_d
            end
        end
        d
    end
end

directions(dqmc::MonteCarloFlavor, ϵ=1e-6) = directions(lattice(dqmc), ϵ)
directions(model::Model, ϵ=1e-6) = directions(lattice(model), ϵ)
function directions(lattice::AbstractLattice, ϵ = 1e-6)
    _positions = positions(lattice)
    wrap = generate_combinations(lattice_vectors(lattice))
    directions = Vector{Float64}[]
    for origin in 1:length(lattice)
        for (trg, p) in enumerate(_positions)
            d = _positions[origin] .- p .+ wrap[1]
            for v in wrap[2:end]
                new_d = _positions[origin] .- p .+ v
                if directed_norm(new_d, ϵ) + ϵ < directed_norm(d, ϵ)
                    d .= new_d
                end
            end
            idx = findfirst(dir -> isapprox(dir, d, atol=ϵ), directions)
            if idx === nothing
                push!(directions, d)
            end
        end
    end
    # temp = sortperm(directions, by=norm)
    # directions[temp]
    sort!(directions, by = v -> directed_norm(v, ϵ))
end


text:
# ported from my MATLAB utils
struct BSpline{T}
    k               # degree of each piecewise polynomial
    t::Vector{T}    # knots; breakpoints
    bdknots         # number of knots at boundaries
end

Base.length(s::BSpline) = length(s.t)

"""
  bspline(k, t)
create bspline object of degree k, and knots -- possibly coin-cident.
Constructor will assume discontinuity at boundaries.

  bspline(k, t, 'boundaries', bb)
as above, but constructor will adapt continuity at the
boundaries:

* bb = k     => discontinuity
* bb = k - a => discontinuity in a'th derivative.

with methods
  * basis_eval(bspline, x)

The resulting basis can then be used with your favourite loss function as reqd.
Unfortunately I have not taken advantage of any sparsity properties (e.g. the
bounded support of the spline basis), because I have better things to do :D.
The Interpolations.jl package does indeed do this, but I can't find an easy way
to extract the basis from this for use in e.g. penalised fits.

follows de Boor - A Practical Guide to Splines (1978 / 2001).
"""
function BSpline(k::Int, t::AbstractVector{T}; bdknots=k) where T <: Real
    @argcheck k > 0
    @argcheck bdknots <= k
    @assert issorted(t)
    @argcheck length(unique(t)) > 1

    if bdknots > 0
        # add additional boundary knots to beginning of vector (if reqd)
        first_not_bd = findfirst(t .!= t[1]) - 1
        numaddl = max(0, bdknots - first_not_bd);
        t       = vcat(repeat(t[1:1], numaddl), t);

        # add additional boundary knots to end of vector (if reqd)
        last_not_bd = findlast(t .!= t[end]) + 1
        numaddl = max(0, bdknots - (length(t) - last_not_bd + 1));
        t       = vcat(t, repeat(t[end:end], numaddl));
    end
    return BSpline{eltype(t)}(k, t, bdknots)
end

function basis_eval(obj::BSpline, x::AbstractVector{T}) where T <: Real
    """
        basis_eval(obj::BSpline, x [,nderiv])
    evaluate the given basis at positions x = [x1, x2, ..., xm].
    Uses recurrence B_{jk} = w_{jk} B_{j,degk-1} + (1-w_{jk}) B_{j+1,degk-1}.

    Derivatives can be calculated instead by specifying the order in the third
    argument (i.e. 0 = fn, 1 = 1st deriv, 2 = 2nd deriv, ...).

    This function owes a lot to James Ramsay's implementation in the FDA
    toolbox.
    """
    N               = length(x);
    # if do_slow === nothing
    #     do_slow = N < 100 ? true : false
    # end

    # recurrence assumes nondecreasing input sequence. Ensure this
    # is so. (cf James Ramsay)
    if minimum(diff(x)) < 0
        reorder_perm = sortperm(x)
        x = x[reorder_perm]
        reordered = true;
    else
        reordered = false;
    end
    nt          = length(obj);

    nbasis      = nt - obj.bdknots;
    degk        = obj.k;
    knotslower  = obj.t[1:nbasis];
    index       = sortperm(vcat(knotslower, x));
    pointer     = findall(index .> nbasis) - (1:N);
    left        = max.(pointer, degk);
    b           = repeat(hcat(1, zeros(1, degk)), N, 1);

    # recursion directly from de Boor (and James Ramsay).
    for jj = 1:degk
        saved = zeros(N, 1);
        for r = 1:jj
            leftpr    = left .+ r;
            tr        = obj.t[leftpr] - x;
            tl        = x - obj.t[leftpr .- jj];

            term      = b[1:N, r]./(tr+tl);
            b[1:N,r]  = saved + tr.*term;
            saved     = tl.*term;
        end
        b[1:N, jj+1]    = saved;
    end

    # now each end column in b corresponds to B_{i-degk+1,degk}, ..., B_{ik}
    # for each i, where i corresponds to the knot to the left of
    # the interval in which each x falls. So to place them all in
    # the same (comparable matrix, we need to shift the columns
    # about.

    # James Ramsay has a cute MATALAB hack to restructure the matrix quickly
    # (benchmarking in MATLAB demonstrates substantial improvement). However the
    # naive (and simple) version here appears to perform just as quickly in julia.
    # if do_slow
    B = zeros(N, nt - obj.bdknots+1);
    for nn = 1:N
        B[nn, (left[nn]-degk+1):left[nn]+1] = b[nn,:];
    end

    # if reordered elements if they were not monotonically increasing, put back in original
    reordered && (B[reorder_perm,:] = B);
    return B
end

function bsplineM(x::AbstractVector{T}, breaks::AbstractVector, norder::Int=4,
        nderiv::Int=0, sparsewrd::Bool=false) where T
    """
        bsplineM(x, breaks, norder, nderiv, sparsewrd)

    Original James Ramsay function ported from MATLAB upon which the `basis_eval`
    function was based. This has greater functionality (does derivatives too).
    I've changed the nature of the sparse return value but o.w. its ~same.

    Computes values or derivatives of B-spline basis functions.

    Arguments:
    X         ... Argument values for which function values are computed
    BREAKS    ... Increasing knot sequence spanning argument range
    NORDER    ... Order of B-spline (one greater than degree) max = 19
                 Default 4.
    NDERIV    ... Order of derivative required, default 0.
    SPARSEWRD ... if true, return in (a non-standard) sparse form

    Return:
    BSPLINEMAT ... length(X) times number of basis functions matrix
                  of Bspline values

    last modified 11 February 2015

    ported from James Ramsay (et al.?)'s FDA toolbox.
    """

    n = length(x);
    nbreaks = length(breaks);

    @assert nbreaks >= 2 "Number of knots less than 2."
    @assert issorted(breaks) "breaks are not strictly increasing"
    @assert norder >= 1 "Order of basis less than one."
    @argcheck nderiv >= 0
    @assert nderiv < norder "NDERIV cannot be as large as order of B-spline."

    if minimum(diff(x)) < 0
        reorder_perm = sortperm(x)
        x = x[reorder_perm]
        reorder_reverse = sortperm(reorder_perm)
        sortwrd = true
    else
        sortwrd = false
    end

    if x[1] - breaks[1] < -1e-10 || x[n] - breaks[nbreaks] > 1e-10
        display([x[1] x[n]; breaks[1] breaks[nbreaks]])
        error("Argument values out of range.")
    end

    if norder == 1
        @argcheck !sparsewrd "norder=1 not yet implemented in sparse form"
        bsplinemat = zeros(T, n, nbreaks-1);
        for ibreaks = 2:nbreaks
            ind = (x >= breaks[ibreaks-1]) & (x <= breaks[ibreaks]);
            bsplinemat[ind, ibreaks-1] = 1;
        end
        return bsplinemat
    end

    # set some abbreviations

    k   = norder;         # order of splines
    km1 = k-1;
    nb  = length(breaks); # number of break points
    nx  = length(x);      # number of argument values
    nd  = nderiv+1;       # ND is order of derivative plus one
    ns  = nb - 2 + k;     # number of splines to compute
    if ns < 1
       @warn "There are no B-splines for the given input."
       return zeros(T, 0, 0)
    end

    #  augment break sequence to get knots by adding a K-fold knot at each end.
    knots  = vcat(repeat(breaks[1:1], km1), breaks, repeat(breaks[nb:nb], km1));
    nbasis = length(knots) - k;

    #  For each  i , determine  left(i)  so that  K <= left(i) < nbasis+1 , and,
    #    within that restriction,
    #        knots(left(i)) <= pts(i) < knots(left(i)+1) .

    knotslower      = knots[1:nbasis];
    index           = sortperm(vcat(knotslower, x));
    pointer         = findall(index .> nbasis) - (1:length(x));
    left            = max.(pointer, k);

    # compute bspline values and derivatives, if needed:
    # (run the recurrence simultaneously for all  x(i).)
    # =====================================================

    # initialize the  b  array.
    b    = repeat(hcat(1, zeros(1, km1)), nd*nx, 1)
    nxs  = nd*(1:nx);

    # First, bring it up to the intended level:
    for j=1:k-nd
       saved = zeros(T, nx, 1);
       for r=1:j
          leftpr   = left .+ r;
          tr       = knots[leftpr] - x;
          tl       = x - knots[leftpr .- j];
          term     = b[nxs, r] ./ (tr + tl);
          b[nxs,r] = saved + tr .* term;
          saved    = tl .* term;
       end
       b[nxs, j+1]  = saved;
    end

    # save the B-spline values in successive blocks in  b .
    for jj=1:nd-1
       j = k - nd + jj;
       saved = zeros(T, nx, 1);
       nxn   = nxs .- 1;
       for r=1:j
          leftpr   = left .+ r;
          tr       = knots[leftpr] - x;
          tl       = x - knots[leftpr .- j];
          term     = b[nxs, r] ./ (tr + tl);
          b[nxn,r] = saved + tr .* term;
          saved    = tl .* term;
       end
       b[nxn,j+1]  = saved;
       nxs = nxn;
    end
    nxn = copy(nxs)

    # now use the fact that derivative values can be obtained by differencing:
    for jj=nd-1:-1:1
       j = k - jj;
       nxs = (repeat(jj:nd-1, 1, nx) + repeat(nxn', nd-jj, 1))[:];
       for r=j:-1:1
          leftpr     = left .+ r;
          temp       = repeat(knots[leftpr] - knots[leftpr .- j], 1, nd-jj)' / j;
          b[nxs, r]   = -b[nxs, r] ./ temp[:];
          b[nxs, r+1] =  b[nxs, r+1] - b[nxs, r];
       end
    end

    # Finally, zero out all rows of  b  corresponding to x outside the basic
    # interval,  [breaks[1] .. breaks[nb]] .
    index = findall((x .< breaks[1]) .| (x .> breaks[nb]));
    if length(index) > 0
       temp = repeat(1-nd:0, 1, length(index)) + nd*repeat(index', nd, 1);
       b[temp[:], :] .= 0
    end

    # If sparse return value, return now.
    if sparsewrd
        elements = b[nd*(1:nx),:];
        first_columns = left .- (km1)
        sortwrd && (elements = elements[reorder_reverse,:])
        sortwrd && (first_columns = first_columns[reorder_reverse,:])

        return elements, first_columns
    end

    # set up output matrix bsplinemat
    width = max(ns,nbasis) + km1 + km1;
    cc    = zeros(nx*width, 1);
    index = repeat(1-nx:0, 1, k) + nx*(repeat(left, 1, k) +
                                       repeat((-km1:0)', nx, 1));

    cc[index] = b[nd*(1:nx),:]
    # (This uses the fact that, for a column vector  v  and a matrix  A ,
    #  v(A)(i,j)=v(A(i,j)), all i,j.)  [#AB: no idea what v(A) etc mean...]
    bsplinemat = reshape(cc[repeat(1-nx:0, 1, ns) +
                            nx*repeat((1:ns)', nx, 1)], nx, ns);

    sortwrd && (bsplinemat = bsplinemat[reorder_reverse, :])

    return bsplinemat

end

ext: jl


text:
for i = 1:100
    if mod(i, 15) == 0
        println("FizzBuzz")
    elseif mod(i, 3) == 0
        println("Fizz")
    elseif mod(i, 5) == 0
        println("Buzz")
    else
        println(i)
    end
end

max_stars_repo_name: jdan/fizzbuzz-polyglot

max_issues_repo_name: jdan/webapp-polyglot

max_forks_repo_name: jdan/webapp-polyglot


text:
## set up McMC
using Distributed
nsamples, nchains, nchainsatone = 100001, 7, 2
Tmax = 2.5
addprocs(nchains)
##init packages on workers
@info "workers are $(workers())"
@everywhere using Distributed
@everywhere using transD_GP
## run McMC
@time transD_GP.main(opt, optn, tempest, Tmax=Tmax, nsamples=nsamples, nchains=nchains, nchainsatone=nchainsatone)

## clean up
rmprocs(workers())


text:
"""
    normalform!(g::FPGroupElement)
Compute the normal form of `g`, possibly modifying `g` in-place.
"""
@inline function normalform!(g::AbstractFPGroupElement)
    isnormalform(g) && return g

    let w = one(word(g))
        w = normalform!(w, g)
        resize!(word(g), length(w))
        copyto!(word(g), w)
    end

    _setnormalform!(g, true)
    _setvalidhash!(g, false)
    @assert isnormalform(g)
    return g
end

"""
    normalform!(res::GEl, g::GEl) where GEl<:FPGroupElement
Compute the normal fom of `g`, storing it in `res`.
"""
function normalform!(res::GEl, g::GEl) where {GEl<:AbstractFPGroupElement}
    @boundscheck @assert parent(res) === parent(g)
    if isnormalform(g)
        copyto!(res, g)
    else
        resize!(word(res), 0)
        normalform!(word(res), g)
        _setnormalform!(res, true)
        _setvalidhash!(res, false)
    end
    return res
end

"""
    normalform!(res::AbstractWord, g::FPGroupElement)
Append the normal form of `g` to word `res`, modifying `res` in place.

Defaults to the rewriting in the free group.
"""
@inline function normalform!(res::AbstractWord, g::AbstractFPGroupElement)
    isone(res) && isnormalform(g) && return append!(res, word(g))
    return KnuthBendix.rewrite_from_left!(res, word(g), rewriting(parent(g)))
end


text:
 abstract type RemarkableObject end

@with_kw struct Document <: RemarkableObject
    ID::String = string(uuid4())
    Version::Int = 1
    Message::String = ""
    Success::Bool = true
    BlobURLGet::String = ""
    BlobURLGetExpires::String = string(DateTime(0))
    BlobURLPut::String = ""
    BlobURLPutExpires::String = string(DateTime(0))
    ModifiedClient::String = string(DateTime(now(UTC))) * "Z"
    VissibleName::String = "new document"
    CurrentPage::Int = 1
    Bookmarked::Bool = false
    Type::String = "DocumentType"
    Parent::String = ""
end

doc_color = crayon"blue"
pdf_color = crayon"red"
col_color = crayon"green"
reset_color = Crayon(reset=true)

@with_kw struct Collection <: RemarkableObject
    ID::String = string(uuid4())
    Version::Int = 1
    Message::String = ""
    Success::Bool = true
    BlobURLGet::String = ""
    BlobURLGetExpires::String = string(DateTime(0))
    BlobURLPut::String = ""
    BlobURLPutExpires::String = string(DateTime(0))
    ModifiedClient::String = string(DateTime(now(UTC))) * "Z"
    VissibleName::String = "new folder"
    CurrentPage::Int = 0
    Bookmarked::Bool = false
    Type::String = "CollectionType"
    Parent::String = ""
    objects::Vector{RemarkableObject} = RemarkableObject[]
end

Document(dict::Dict{String, Any}) = Document(;Dict(Symbol(key)=>value for (key, value) in dict)...)
Collection(dict::Dict{String, Any}) = Collection(;Dict(Symbol(key)=>value for (key, value) in dict)...)

Base.getindex(c::Collection, i::Int) = c.objects[i]
Base.iterate(c::Collection, state) = iterate(c.objects, state)
Base.iterate(c::Collection) = iterate(c.objects)
Base.length(c::Collection) = length(c.objects)

AbstractTrees.children(::Document) = ()
AbstractTrees.children(c::Collection) = c.objects

title(x::RemarkableObject) = x.VissibleName

ispdf(d::Document) = endswith(d.VissibleName, ".pdf")

AbstractTrees.printnode(io::IO, d::Document) = print(io, ispdf(d) ? pdf_color : doc_color, d.VissibleName, reset_color)
AbstractTrees.printnode(io::IO, c::Collection) = print(io, col_color, c.VissibleName, reset_color)

function create_tree(docs::AbstractVector{<:RemarkableObject})
    root = Collection(ID = "", VissibleName = "Root")
    push!(root.objects, Collection(ID = "Trash", VissibleName = "Trash"))
    update_obj!(root, docs) # Recursive loop on documents
    return root
end

function update_obj!(col::Collection, docs)
    for doc in docs
        if doc.Parent == col.ID
            update_obj!(doc, docs)
            push!(col.objects, doc)
        end
    end
end

update_obj!(::Document, ::Any) = nothing

function obj_to_dict(doc::Document)
    dict = type2dict(doc)
    return Dict(string(key)=>value for (key, value) in dict)
end

function obj_to_dict(col::Collection)
    dict = type2dict(col)
    delete!(dict, :objects)
    return Dict(string(key)=>value for (key, value) in dict)
end

function find_col(f, client::RemarkableClient, col::Collection=list_items(client))
    if f(col)
        return col
    end
    for c in col
        res = find_col(f, client, c)
        isnothing(res) ? nothing : return res
    end
    return nothing
end

find_col(f, client::RemarkableClient, col::Document) = nothing


text:
zero_map(::Type{LinearMap}, M::Integer, N::Integer) =
    LinearMap((y,x) -> (y .= false), M, N,
              ismutating=true, issymmetric=true, ishermitian=true)

block_eltype(::Type{LinearMap{T}}) where T = T
block_spy_block(b::LinearMaps.WrappedMap{<:Any,<:AbstractMatrix}) =
    block_spy_block(b.lmap)

isblockentryzero(b::LinearMaps.WrappedMap{<:Any,<:AbstractMatrix}, args...) =
    isblockentryzero(b.lmap, args...)

unit_vec(::Type{T}, m, i) where T = Vcat(Zeros{T}(i-1), one(T), Zeros{T}(m-i))

function block_getindex(A::LinearMap{T}, i, j) where T
    m,n = size(A)

    ej = unit_vec(T, n, j)
    tmp = Vector{T}(undef, m)

    mul!(tmp, A, ej)

    tmp[i]
end

# This is necessary type piracy since SparseArrays checks if entry
# iszero on assignment. Should be upstreamed.
Base.iszero(::LinearMap) = false


text:
"""
    weight(sequence; ignoreunknown = false)

Estimate molecular weight of `sequence` in Dalton.

If `ignoreunknown` is set to `false` (default), non-specific residues (X: any AA, Z: glutamine or glutamic acid, B: asparagine or aspartic acid, J: leucine or isoleucine) are calculated using the average of all possible residues in that category. Otherwise they are ignored. Only standard residues are used to calculate the average for X.
"""
function weight(sequence; ignoreunknown = false)
    aacomposition = BioSequences.composition(sequence)
    mW = 18.01528 # ≈ weight of one water molecule in Dalton
    for (aa, n) in aacomposition
        if aa ∉ [AA_X, AA_Z, AA_B, AA_J] || !ignoreunknown
            mW += get(aminoacidweights, aa, 0.0) * n
        end
    end
    return mW
end


text:
cell_to_vtkcell(::Type{LineCell}) = WriteVTK.VTKCellTypes.VTK_LINE
cell_to_vtkcell(::Type{RectangleCell}) = WriteVTK.VTKCellTypes.VTK_QUAD
cell_to_vtkcell(::Type{TriangleCell}) = WriteVTK.VTKCellTypes.VTK_TRIANGLE
cell_to_vtkcell(::Type{HexagonCell}) = WriteVTK.VTKCellTypes.VTK_POLYGON
cell_to_vtkcell(::Type{HexahedronCell}) = WriteVTK.VTKCellTypes.VTK_HEXAHEDRON
cell_to_vtkcell(::Type{TetrahedronCell}) = WriteVTK.VTKCellTypes.VTK_TETRA

"""
    vtk_grid(filename::AbstractString, grid::Grid)

Create a unstructured VTK grid from a `PolytopalMesh`. Return a `DatasetFile`
which data can be appended to, see `vtk_point_data` and `vtk_cell_data`.
"""
function vtk_grid(filename::AbstractString, mesh::AbstractPolytopalMesh{dim,T}; compress::Bool=true) where {dim,T}
    cls = WriteVTK.MeshCell[]
    for cell_idx in 1:getncells(mesh)
        celltype = cell_to_vtkcell(getCellType(mesh, cell_idx))
        push!(cls, WriteVTK.MeshCell(celltype, getverticesindices(mesh,getcell(mesh,cell_idx))))
    end
    coords = get_vertices_matrix(mesh)'
    return WriteVTK.vtk_grid(filename, coords, cls; compress=compress)
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
"""
        AbstractMetrics

AbstractMetrics is an abstract type that allows the user the stock relevant information
of consecutive search along learning (either on different instances or on the same ones).

SeaPearl allows the user to define his own CustomMetrics as long as it contains:
1) a constructor `metrics(model::CPModel, heuristic::ValueSelection)`.
2) a function `(::CustomMetrics)(model::CPmodel,dt::Float64)` which is a generic call to
AbstractMetrics made just after the search.
The user can add new features on its CustomMetrics such as advanced plotting.
"""
abstract type AbstractMetrics end

"""
        AbstractTakeObjective

The structure TakeObjective and DontTakeObjective are used as parametric type for the
BasicMetrics{OBJ<:AbstractTakeObjective, H<:ValueSelection} definition. It indicates if the
CP problem (knapsack, graphcoloring, tsptw ...) that the metrics is attached to deals with
an objective function or not.

If yes, after a search on a given problem instance, the metrics will also retrieve the solution scores along the search.
"""
abstract type AbstractTakeObjective end
struct TakeObjective<:AbstractTakeObjective end
struct DontTakeObjective<:AbstractTakeObjective end

include("basicmetrics.jl")


text:
# Author: Youngjun Kim, youngjun@stanford.edu
# Date: 06/05/2014

# Correlated Encounter Model for Cooperative Aircraft in the National Airspace
# System Version 1.0

module CorrAEMImpl

export
    AddObserver,

    initialize,
    update,

    generateEncounter,
    getInitialState,
    getNextCommand,
    getTrajectory,

    CorrAEM,
    CorrAEMInitialState,
    CorrAEMCommand,

    generateEncountersToFile,
    setInitialDistributions,
    validate,
    getInitialSample

import Compat.ASCIIString

using AbstractEncounterModelImpl
using AbstractEncounterModelInterfaces
using CommonInterfaces
using Util
using RLESUtils, Observers
using Base.Test

import CommonInterfaces.addObserver
import CommonInterfaces.initialize
import CommonInterfaces.update
import AbstractEncounterModelInterfaces.generateEncounter
import AbstractEncounterModelInterfaces.getInitialState
import AbstractEncounterModelInterfaces.getNextCommand
import AbstractEncounterModelInterfaces.getTrajectory

import AbstractEncounterModelInterfaces.generateEncountersToFile
import AbstractEncounterModelInterfaces.setInitialDistributions
import AbstractEncounterModelInterfaces.validate
import AbstractEncounterModelInterfaces.getInitialSample

include("corr_aem_load_params.jl")
include("corr_aem_sample.jl")
include("corr_aem_validate.jl")
include("corr_aem_save_scripts.jl")
include("corr_aem_load_scripts.jl")

type CorrAEMInitialState

    t::Float64
    x::Float64
    y::Float64
    h::Float64
    v::Float64
    psi::Float64 #degrees, x-axis is zero, increasing toward positive y
    h_d::Float64
end

type CorrAEMCommand

    t::Float64
    v_d::Float64
    h_d::Float64
    psi_d::Float64

end

type CorrAEMParameters

    labels_initial
    n_initial

    G_initial
    r_initial
    N_initial

    labels_transition
    n_transition

    G_transition
    r_transition
    N_transition

    boundaries
    resample_rates

    temporal_map
    zero_bins

    CorrAEMParameters() = new()
end

type CorrAEM <: AbstractEncounterModel

    parameter_filename::ASCIIString
    parameters::CorrAEMParameters

    initial_sample_filename::ASCIIString
    number_of_initial_samples::Int
    f_init::Union{IOStream,Void}

    transition_sample_filename::ASCIIString
    number_of_transition_samples::Int
    f_tran::Union{IOStream, Void}

    b_read_from_file::Bool
    b_write_to_file::Bool
    file_format::ASCIIString

    number_of_encounters_generated::Int
    print_header::Bool

    number_of_aircraft::Int

    initial::Vector{Float64}

    A::Int
    L::Int
    C::Vector{Int}
    v_init::Vector{Float64}

    geometry_at_TCA::Vector{Float64}
    state_tca::Array{Float64, 2}

    states::Array{Float64, 3}
    state_index::Vector{Int}

    dynamic_states::Array{Float64, 3}
    dn_state_index::Vector{Int}

    initial_distributions::Union{Vector{Union{Function, Void}}, Void}
    ISInfo::Array{Float64, 2}

    const default_number_of_initial_samples = 1
    const default_number_of_transition_samples = 50

    observer::Observer

    function CorrAEM(
                    parameter_filename::AbstractString,
                    initial_sample_filename::AbstractString,
                    number_of_initial_samples::Int,
                    transition_sample_filename::AbstractString,
                    number_of_transition_samples::Int;
                    b_read_from_file = false,
                    b_write_to_file = false,
                    file_format = "text")

        obj = new()

        obj.parameter_filename = parameter_filename
        obj.parameters = CorrAEMParameters()
        em_read(obj.parameters, obj.parameter_filename)

        obj.initial_sample_filename = initial_sample_filename
        obj.number_of_initial_samples = number_of_initial_samples
        obj.f_init = nothing

        obj.transition_sample_filename = transition_sample_filename
        obj.number_of_transition_samples = number_of_transition_samples
        obj.f_tran = nothing

        obj.b_read_from_file = b_read_from_file
        obj.b_write_to_file = b_write_to_file
        obj.file_format = file_format

        obj.number_of_encounters_generated = 0
        obj.print_header = true

        @test (b_read_from_file && b_write_to_file) == false

        if b_read_from_file || b_write_to_file
            @test initial_sample_filename != ""
            @test transition_sample_filename != ""
        end

        if b_read_from_file
            @test file_format == "text"
        end

        obj.number_of_aircraft = 2

        obj.C = zeros(obj.number_of_aircraft)
        obj.v_init = zeros(obj.number_of_aircraft)

        obj.state_tca = zeros(obj.number_of_aircraft, 5)

        obj.state_index = zeros(obj.number_of_aircraft)
        obj.states = zeros(obj.number_of_aircraft, obj.number_of_transition_samples, 4)

        obj.dn_state_index = zeros(obj.number_of_aircraft)
        obj.dynamic_states = zeros(obj.number_of_aircraft, obj.number_of_transition_samples + 1, 6)

        obj.initial_distributions = nothing

        obj.observer = Observer()

        return obj
    end

    function CorrAEM(
                    parameter_filename::AbstractString,
                    initial_sample_filename::AbstractString,
                    transition_sample_filename::AbstractString;
                    b_write_to_file::Bool = false,
                    file_format::AbstractString = "text")

        if b_write_to_file
            CorrAEM(parameter_filename,
                    initial_sample_filename,
                    default_number_of_initial_samples,
                    transition_sample_filename,
                    default_number_of_transition_samples,
                    b_write_to_file = true,
                    file_format = file_format)
        else
            CorrAEM(parameter_filename,
                    initial_sample_filename,
                    default_number_of_initial_samples,
                    transition_sample_filename,
                    default_number_of_transition_samples,
                    b_read_from_file = true)
        end
    end

    CorrAEM(
        parameter_filename::AbstractString
        ) = CorrAEM(
                parameter_filename,
                "",
                default_number_of_initial_samples,
                "",
                default_number_of_transition_samples)

    CorrAEM() = new()
end

addObserver(aem::CorrAEM, f::Function) = add_observer(aem.observer, f)
addObserver(aem::CorrAEM, tag::AbstractString, f::Function) = add_observer(aem.observer, tag, f)

function generateEncounter(aem::CorrAEM; sample_number = 0, b_simulate = true)
    params = aem.parameters

    if aem.b_write_to_file && aem.file_format == "text"
        f_initial = open(aem.initial_sample_filename, "a+")
        #f_initial_cmp = open(aem.initial_sample_filename * "_cmp", "w")
        f_transition = open(aem.transition_sample_filename, "a+")

        if aem.print_header
            # print initial headers
            @printf(f_initial, "id ")
            #@printf(f_initial_cmp, "id ")
            for i = 1:params.n_initial
                @printf(f_initial, "%s ", params.labels_initial[i])
                #@printf(f_initial_cmp, "%s ", params.labels_initial[i])
            end
            @printf(f_initial, "\n")
            #@printf(f_initial_cmp, "\n")

            # print transition headers
            @printf(f_transition, "initial_id t ")
            for i = 1:(params.n_transition - params.n_initial)
                @printf(f_transition, "%s ", params.labels_transition[params.temporal_map[i, 2]])
            end
            @printf(f_transition, "\n")

            aem.print_header = false
        end
    end

    if aem.b_read_from_file
        println("read_from_file is true")
        if sample_number > 0
            reset_sample_from_file(aem)
            for i = 1:sample_number
                states = read_sample_from_file(aem, aem.number_of_transition_samples)

                if states == nothing
                    return -1
                end
            end
        else
            println("reading directly")
            states = read_sample_from_file(aem, aem.number_of_transition_samples)
        end
    else
        if aem.initial_distributions == nothing
            states = em_sample(aem, aem.number_of_transition_samples)
        else
            states, aem.ISInfo = em_sample(aem, aem.number_of_transition_samples, initial_dist = aem.initial_distributions)
        end

        aem.number_of_encounters_generated += 1

        states = states[:, 2:end]

        if aem.b_write_to_file && aem.file_format == "text"
            # print initial sample
            @printf(f_initial, "%d ", aem.number_of_encounters_generated)
            #@printf(f_initial_cmp, "%d ", aem.number_of_encounters_generated)
            for j = 2:size(states, 2)
                @printf(f_initial, "%s ", outputGFormatString(states[1, j]))
                #@printf(f_initial_cmp, "%s ", states[1, j])
            end
            @printf(f_initial, "\n")
            #@printf(f_initial_cmp, "\n")

            # print transition samples
            for j = 1:aem.number_of_transition_samples
                @printf(f_transition, "%d %d ", aem.number_of_encounters_generated, j - 1)
                for k = params.temporal_map[:, 1]
                    @printf(f_transition, "%s ", outputGFormatString(states[j, k + 1]))
                end
                @printf(f_transition, "\n")
            end
        end
    end

    if b_simulate
        convert_units(states)

        t, A, L, chi, beta_, C1, C2, v1, v2, v1d, v2d, h1d, h2d, psi1d, psi2d, hmd, vmd = states[1, :]

        aem.initial = vec(states[1, 2:end])

        aem.A = A
        aem.L = L
        aem.C[1] = C1
        aem.C[2] = C2
        aem.v_init[1] = v1
        aem.v_init[2] = v2

        aem.geometry_at_TCA = [chi, beta_, hmd, vmd]

        aem.states[1, :, 1] = states[:, 1]  # time
        aem.states[1, :, 2] = states[:, 10] # v1d
        aem.states[1, :, 3] = states[:, 12] # h1d
        aem.states[1, :, 4] = states[:, 14] # psi1d
        aem.state_index[1] = 0

        aem.states[2, :, 1] = states[:, 1]  # time
        aem.states[2, :, 2] = states[:, 11] # v2d
        aem.states[2, :, 3] = states[:, 13] # h2d
        aem.states[2, :, 4] = states[:, 15] # psi2d
        aem.state_index[2] = 0

        aem.dn_state_index[1] = 0
        aem.dn_state_index[2] = 0

        simulate_tracks(aem)

        transform_regarding_TCA(aem, aem.L, aem.geometry_at_TCA)

        for i = 1:aem.number_of_aircraft
            aem.state_index[i] = 0
        end

        #for i = 1:(aem.number_of_transition_samples + 1)
        #    print(reshape(aem.dynamic_states[1, i, :], 6)')
        #end

        #for i = 1:aem.number_of_transition_samples
        #    print(reshape(aem.states[1, i, :], 4)')
        #end
    end

    if aem.b_write_to_file && aem.file_format == "binary"
        encounters = Array(Any, aem.number_of_aircraft)

        for i = 1:aem.number_of_aircraft
            encounters[i] = Dict()
        end

        # initial
        # airspeed, ft/s, double
        # north position, ft, double
        # east position, ft, double
        # altitude, ft, double
        # heading angle, radians, double
        # flight path angle, radians, double
        # roll angle, radians, double
        # airspeed acceleration, ft/s^2, double

        for i = 1:aem.number_of_aircraft
            t, x, y, h, v, psi = aem.dynamic_states[i, 1, :]
            t_n, x_n, y_n, h_n, v_n, psi_n, = aem.dynamic_states[i, 2, :]

            t, v_d, h_d, psi_d = aem.states[i, 1, :]

            theta = atand((h_n - h) / norm(x_n - x, y_n - y)) |> to_plusminus_180
            phi = 0.0

            encounters[i]["initial"] = Float64[v, x, y, h, deg2rad(psi), deg2rad(theta), deg2rad(phi), v_d]
        end

        # time, s, double
        # vertical rate, ft/s, double
        # turn rate, rad/s, double
        # airspeed acceleration, ft/s^2, double

        for i = 1:aem.number_of_aircraft
            update = Array(Float64, 4, aem.number_of_transition_samples)

            k = 0

            t, v_d, h_d, psi_d = aem.states[i, 1, :]

            k += 1
            update[:, k] = [t, h_d, psi_d * pi / 180, v_d]

            v_d_p = v_d
            h_d_p = h_d
            psi_d_p = psi_d

            for j = 2:aem.number_of_transition_samples
                t, v_d, h_d, psi_d = aem.states[i, j, :]

                if v_d != v_d_p || h_d != h_d_p || psi_d != psi_d_p
                    k += 1
                    update[:, k] = [t, h_d, psi_d * pi / 180, v_d]
                end

                v_d_p = v_d
                h_d_p = h_d
                psi_d_p = psi_d
            end

            encounters[i]["update"] = update[:, 1:k]
        end

        save_scripts(aem.initial_sample_filename, encounters, append = true)
    end

    if aem.b_write_to_file && aem.file_format == "text"
        close(f_initial)
        #close(f_initial_cmp)
        close(f_transition)
    end
end

function generateEncountersToFile(aem::CorrAEM; file_format = "")

    if isfile(aem.initial_sample_filename)
        println("Removing ", aem.initial_sample_filename)
        rm(aem.initial_sample_filename)
    end

    if isfile(aem.transition_sample_filename)
        println("Removing ", aem.transition_sample_filename)
        rm(aem.transition_sample_filename)
    end

    println("write_to_file is ", aem.b_write_to_file, " file_format_orig is ", aem.file_format)
    b_write_to_file_orig = aem.b_write_to_file
    file_format_orig = aem.file_format

    aem.b_write_to_file = true
    if file_format != ""
        aem.file_format = file_format
    end

    for i = 1:aem.number_of_initial_samples
        if aem.file_format == "text"
            generateEncounter(aem, b_simulate = false)
        elseif aem.file_format == "binary"
            generateEncounter(aem)
        end
    end

    aem.b_write_to_file = b_write_to_file_orig
    aem.file_format = file_format_orig
end

function getInitialState(aem::CorrAEM, aircraft_number::Int)

    state = aem.dynamic_states[aircraft_number, 1, :]   # t, x, y, h, v, psi
    command = aem.states[aircraft_number, 1, :]         # t, v_d, h_d, psi_d

    return CorrAEMInitialState(state[1], state[2], state[3], state[4], state[5], state[6], command[3])
end

getNextCommand(aem::CorrAEM, aircraft_number::Int) = get_next_state(aem, aircraft_number)

function initialize(aem::CorrAEM)

  for i = 1:aem.number_of_aircraft
    aem.state_index[i] = 0
  end

end

function update(aem::CorrAEM, aircraft_number::Int)

    command = getNextCommand(aem, aircraft_number)

    if command != nothing
        return CorrAEMCommand(command...)
    else
        return nothing
    end
end

function getTrajectory(aem::CorrAEM, aircraft_number::Int)

    return reshape(aem.dynamic_states[aircraft_number, :, 1:4], aem.dn_state_index[aircraft_number], 4)
end

function setInitialDistributions(aem::CorrAEM, n::Int, f::Function)

    if aem.initial_distributions == nothing
        aem.initial_distributions = Array(Union{Function,Void}, 16)

        for i = 1:16
            aem.initial_distributions[i] = nothing
        end
    end

    aem.initial_distributions[n] = f
end

function getInitialSample(aem::CorrAEM, index::Int)

    return aem.initial[index]
end

function getInitialSample(aem::CorrAEM, var::Symbol=:all)

    if var == :all
        return aem.initial
    elseif var == :A
        return aem.initial[1]
    elseif var == :L
        return aem.initial[2]
    elseif var == :chi
        return aem.initial[3]
    elseif var == :beta
        return aem.initial[4]
    elseif var == :C1
        return aem.initial[5]
    elseif var == :C2
        return aem.initial[6]
    elseif var == :v1
        return aem.initial[7]
    elseif var == :v2
        return aem.initial[8]
    elseif var == :v1d
        return aem.initial[9]
    elseif var == :v2d
        return aem.initial[10]
    elseif var == :h1d
        return aem.initial[11]
    elseif var == :h2d
        return aem.initial[12]
    elseif var == :psi1d
        return aem.initial[13]
    elseif var == :psi2d
        return aem.initial[14]
    elseif var == :hmd
        return aem.initial[15]
    elseif var == :vmd
        return aem.initial[16]
    else
        return nothing
    end
end

function transform_regarding_TCA(aem, L, geometry_at_TCA)

    chi, beta_, hmd, vmd = geometry_at_TCA

    x1_tca, y1_tca, h1_tca, v1_tca, psi1_tca = aem.state_tca[1, :]
    x2_tca, y2_tca, h2_tca, v2_tca, psi2_tca = aem.state_tca[2, :]

    AC1_state = reshape(aem.dynamic_states[1, 1:aem.dn_state_index[1], :], aem.dn_state_index[1], 6)
    AC2_state = reshape(aem.dynamic_states[2, 1:aem.dn_state_index[2], :], aem.dn_state_index[2], 6)

    AC1_state[:, 2] -= x1_tca
    AC1_state[:, 3] -= y1_tca

    AC1_state[:, 2:3] = transpose([cosd(-psi1_tca) -sind(-psi1_tca); sind(-psi1_tca) cosd(-psi1_tca)] * AC1_state[:, 2:3]')

    if L == 1
        h_trans = 1000 + (3000 - 1000) * rand()
    elseif L == 2
        h_trans = 3000 + (10000 - 3000) * rand()
    elseif L == 3
        h_trans = 10000 + (18000 - 10000) * rand()
    elseif L == 4
        h_trans = 18000 + (29000 - 18000) * rand()
    elseif L == 5
        h_trans = 29000 + (50000 - 29000) * rand()
    end

    AC1_state[:, 4] += (-h1_tca + h_trans)

    psi = update_psi(AC1_state[1, 2:3], AC1_state[2, 2:3])
    AC1_state[:, 6] = map(to_plusminus_180, AC1_state[:, 6] + psi) #make sure we're in [-180, 180]

    AC2_state[:, 2] -= x2_tca
    AC2_state[:, 3] -= y2_tca

    psi2_0 = beta_ - psi2_tca
    AC2_state[:, 2:3] = transpose([cosd(psi2_0) -sind(psi2_0); sind(psi2_0) cosd(psi2_0)] * AC2_state[:, 2:3]')

    v_r_tca = [v2_tca * cosd(beta_) - v1_tca, v2_tca * sind(beta_)]

    x_trans = hmd / norm(v_r_tca) * [-v_r_tca[2], v_r_tca[1]]

    if chi == 1     # front
        if x_trans[1] < 0
            x_trans = -x_trans
        end
    else
        if x_trans[1] >= 0
            x_trans = -x_trans
        end
    end

    AC2_state[:, 2] += x_trans[1]
    AC2_state[:, 3] += x_trans[2]

    AC2_state[:, 4] += (-h2_tca + h_trans - vmd)

    psi = update_psi(AC2_state[1, 2:3], AC2_state[2, 2:3])
    AC2_state[:, 6] = map(to_plusminus_180, AC2_state[:, 6] + psi) #make sure we're in [-180, 180]

    aem.dynamic_states[1, 1:aem.dn_state_index[1], :] = reshape(AC1_state, 1,
        size(AC1_state, 1), 6)
    aem.dynamic_states[2, 1:aem.dn_state_index[2], :] = reshape(AC2_state, 1,
        size(AC2_state, 1), 6)
end

function update_psi(state1, state2)

    x1, y1 = state1
    x2, y2 = state2

    x = x2 - x1
    y = y2 - y1

    if x >= 0 && y >= 0
        psi = atand(y / x)
    elseif x <= 0 && y >= 0
        psi = atand(y / x) + 180
    elseif x <= 0 && y <= 0
        psi = atand(y / x) + 180
    elseif x >= 0 && y <= 0
        psi = atand(y / x) + 360
    else
        error("None worked state1 is ", state1, " and state2 is ", state2)
    end

    return psi
end

function simulate_tracks(aem)

    for i = 1:aem.number_of_aircraft
        t_curr, v_d_curr, h_d_curr, psi_d_curr = get_next_state(aem, i)
        v_curr = aem.v_init[i]

        aem.dn_state_index[i] += 1
        aem.dynamic_states[i, aem.dn_state_index[i], :] = [t_curr, 0., 0., 0, v_curr, 0.]

        while true
            state_next = get_next_state(aem, i)

            if state_next != nothing
                t_next, v_d_next, h_d_next, psi_d_next = state_next

                sim_time = t_next
            else    # last second
                sim_time += 1
            end

            #if i == 1
            #    print([sim_time, v_d_curr, h_d_curr, psi_d_curr]')
            #end

            state = update_state(aem, i, sim_time, v_d_curr, h_d_curr, psi_d_curr)

            if abs(sim_time - 40) < 0.001
                aem.state_tca[i, :] = [state...]
            end

            if state_next == nothing
                break
            else
                v_d_curr = v_d_next
                h_d_curr = h_d_next
                psi_d_curr = psi_d_next
            end
        end
    end
end

function update_state(aem, aircraft_number, t, v_d, h_d, psi_d)

    theta_regulated = 45    # degree

    t_, x, y, h, v, psi = aem.dynamic_states[aircraft_number, aem.dn_state_index[aircraft_number], :]

    if abs(h_d) > abs(v) * sind(theta_regulated)
        h_d = sign(h_d) * abs(v) * sind(theta_regulated)
    end

    dt = t - t_

    x_n = x + sqrt(v^2 - h_d^2) * cosd(psi) * dt
    y_n = y + sqrt(v^2 - h_d^2) * sind(psi) * dt

    v_n = v + v_d * dt  # ft/s
    h_n = h + h_d * dt  # ft
    psi_n = psi + psi_d * dt    # deg

    aem.dn_state_index[aircraft_number] += 1
    aem.dynamic_states[aircraft_number, aem.dn_state_index[aircraft_number], :] = [t, x_n, y_n, h_n, v_n, psi_n]

    return x_n, y_n, h_n, v_n, psi_n
end

function convert_units(values)

    values[:, 8] *= 1.68780
    values[:, 9] *= 1.68780
    values[:, 10] *= 1.68780
    values[:, 11] *= 1.68780
    values[:, 12] /= 60
    values[:, 13] /= 60
    values[:, 16] *= 6076.12

    return values
end

function get_next_state(aem, aircraft_number)

    aem.state_index[aircraft_number] += 1

    if aem.state_index[aircraft_number] > aem.number_of_transition_samples
        return nothing
    else
        return reshape(aem.states[aircraft_number, aem.state_index[aircraft_number], :], 4)
    end
end

# TODO support to read number of samples
# returns a sample at a time, currently

function read_sample_from_file(aem, number_of_transition_samples)

    # 1       id
    # 2       A = [1:4]
    # 3       L = [1:5]
    # 4       chi = [90, 270]   : deg
    # 5       beta = [0:30:360] : deg
    # 6       C1 = [1, 2]
    # 7       C2 = [1, 2]
    # 8,9     v1, v2 = [50, 100, 200, 300, 400, 500, 600]   : kt
    # 10, 11  v1d, v2d = [-5, -2, -0.25, 0.25, 2, 5]        : kt/s
    # 12, 13  h1d, h2d = [-5000, -3000, -2000, -1000, -400, 400, 1000, 2000, 3000, 5000]    : ft/min
    # 14, 15  psi1d, psi2d = [-8, -6, -3.5, -1, -0.25, 0.25, 1, 3.5, 6, 8]  : deg/s
    # 16      hmd = [0, 0.0822896, 0.5, 1, 3]   : nm
    # 17      vmd = [0, 100, 200, 300, 400, 500, 600, 700, 800, 900, 6000]  : ft

    if aem.f_init == nothing
        aem.f_init = open(aem.initial_sample_filename, "r")
        readline(aem.f_init)
    end

    line = readline(aem.f_init)

    if line == ""
        close(aem.f_init)
        aem.f_init = nothing
        initial = nothing
    else
        initial = float(split(chomp(line)))
    end

    # 1     id
    # 2     time step
    # 3, 4  h1d, h2d        : ft/min
    # 5, 6  psi1d, psi2d    : deg/s

    if aem.f_tran == nothing
        aem.f_tran = open(aem.transition_sample_filename, "r")
        readline(aem.f_tran)
    end

    line = readline(aem.f_tran)

    if line == ""
        close(aem.f_tran)
        aem.f_tran = nothing

        transitions = nothing
    else
        values = float(split(chomp(line)))

        transitions = Array(Float64, aem.number_of_transition_samples, length(values))
        transitions[1, :] = values

        for i = 2:aem.number_of_transition_samples
            transitions[i, :] = float(split(chomp(readline(aem.f_tran))))
        end
    end

    if initial != nothing && transitions != nothing
        states = repmat(initial', size(transitions, 1))
        states[:, 1] = transitions[:, 2]
        states[:, aem.parameters.temporal_map[:, 1] + 1] = transitions[:, 3:(3 + length(aem.parameters.temporal_map[:, 1]) - 1)]

        return states
    else
        return nothing
    end
end

function reset_sample_from_file(aem)

    if aem.f_init != nothing
        close(aem.f_init)
        aem.f_init = nothing
    end

    if aem.f_tran != nothing
        close(aem.f_tran)
        aem.f_tran = nothing
    end

end

#mods x to the range [-b, b]
function to_plusminus_b(x::AbstractFloat, b::AbstractFloat)
  z = mod(x, 2 * b)
  return (z > b) ? (z - 2 * b) : z
end

to_plusminus_180(x::AbstractFloat) = to_plusminus_b(x, 180.0)

end


text:
module TestAsTableFlag

using Test
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@astable with just assignments" begin
    df = DataFrame(a = 1, b = 2)

    d = @rtransform df @astable begin
        :x = 1
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 1)

    d = @rselect df @astable begin
        :x = 1
        y = 100
        nothing
    end

    @test d == DataFrame(x = 1)

    d = @transform df @astable begin
        :x = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 5)

    d = @select df @astable begin
        :x = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(x = 5)
end

@testset "@astable with strings" begin
    df = DataFrame(a = 1, b = 2)

    x_str = "x"
    d = @rtransform df @astable begin
        $x_str = 1
        y = 100
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 1)

    d = @rselect df @astable begin
        $x_str = 1
        y = 100
        nothing
    end

    @test d == DataFrame(x = 1)

    d = @transform df @astable begin
        $x_str = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 5)

    d = @select df @astable begin
        $x_str = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(x = 5)
end

@testset "Re-using variables" begin
    df = DataFrame(a = 1, b = 2)

    d = @rtransform df @astable begin
        :x = 1
        y = 5
        :z = :x + y
    end

    @test d == DataFrame(a = 1, b = 2, x = 1, z = 6)

    d = @rselect df @astable begin
        :x = 1
        y = 5
        :z = :x + y
    end

    @test d == DataFrame(x = 1, z = 6)

    x_str = "x"
    d = @rtransform df @astable begin
        $x_str = 1
        y = 5
        :z = $x_str + y
    end

    @test d == DataFrame(a = 1, b = 2, x = 1, z = 6)

    d = @rselect df @astable begin
        $x_str = 1
        y = 5
        :z = $x_str + y
    end

    @test d == DataFrame(x = 1, z = 6)
end

@testset "@astable with mutation" begin
    df = DataFrame(a = 1, b = 2)

    df2 = copy(df)
    d = @rtransform! df2 @astable begin
        :x = 1
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 1)
    @test d === df2

    df2 = copy(df)
    d = @rselect! df2 @astable begin
        :x = 1
        y = 100
        nothing
    end

    @test d == DataFrame(x = 1)
    @test d === df2

    df2 = copy(df)
    d = @transform! df2 @astable begin
        :x = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(a = 1, b = 2, x = 5)
    @test d === df2

    df2 = copy(df)
    d = @select! df2 @astable begin
        :x = [5]
        y = 100
        nothing
    end

    @test d == DataFrame(x = 5)
    @test d === df2
end

@testset "grouping astable flag" begin
    df = DataFrame(a = [1, 1, 2, 2], b = [5, 6, 7, 8])

    gd = groupby(df, :a)

    d = @combine gd @astable begin
        ex = extrema(:b)
        :b_min = ex[1]
        :b_max = ex[2]
    end

    res_sorted = DataFrame(a = [1, 2], b_min = [5, 7], b_max = [6, 8])

    @test sort(d, :b_min) == res_sorted

    d = @combine gd @astable begin
        ex = extrema(:b)
        $"b_min" = ex[1]
        $"b_max" = ex[2]
    end

    @test sort(d, :b_min) == res_sorted

    d = @by df :a @astable begin
        ex = extrema(:b)
        :b_min = ex[1]
        :b_max = ex[2]
    end

    @test sort(d, :b_min) == res_sorted

    d = @by df :a @astable begin
        ex = extrema(:b)
        $"b_min" = ex[1]
        $"b_max" = ex[2]
    end

    @test sort(d, :b_min) == res_sorted
end

@testset "errors with passmissing" begin
    @eval df = DataFrame(y = 1)
    @test_throws LoadError @eval @transform df @passmising @byrow @astable :x = 2
    @test_throws LoadError @eval @transform df @byrow @astable @passmissing :x = 2
    @test_throws LoadError @eval @transform df @astable @passmissing @byrow :x = 2
    @test_throws LoadError @eval @rtransform df @astable @passmissing :x = 2
    @test_throws LoadError @eval @rtransform df @passmissing @astable :x = 2
end

@testset "bad assignments" begin
    @eval df = DataFrame(y = 1)
    @test_throws ArgumentError @eval @transform df @astable cols(1) = :y
    @test_throws ArgumentError @eval @transform df @astable cols(AsTable) = :y
end

end # module


text:
# Chapter 15. Structs and Objects

include("ch15module.jl")

# Does not work?
# using .ch15
# p = Point(3.0, 4.0)

import .AllnixThinkJulia
ns = AllnixThinkJulia # alias

# Composite Types
p = ns.Point(3.0, 4.0)
q = ns.Point(12, 13)

# Mutable Structs
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_mutable_structs
blank = ns.MPoint(0.0, 0.0)
origin = ns.MPoint(0.0, 0.0)

# Rectangles
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_rectangles
box = ns.Rectangle(100., 200., origin)

# Instances as Arguments
ns.printpoint(p)
distance = ns.distancebetweenpoints(p, q)
ns.movepoint!(origin, 3, 4)
# ns.movepoint!(p, 4, 5) # setfield! immutable struct of type Point cannot be changed
ns.moverectangle!(box, 1, 2)

# Instances as Return Values
center = ns.findcenter(box)

# Copying
# https://benlauwens.github.io/ThinkJulia.jl/latest/book.html#_copying
p1 = ns.MPoint(3.0, 4.0)
p2 = deepcopy(p1)
p1 ≡ p2   # false
p1 === p2 # false
p1 == p2  # false

# Exercise 15-2
p1 = ns.Point(3.0, 4.0)
p2 = deepcopy(p1)
p1 ≡ p2   # true
p1 === p2 # true
p1 == p2  # true

# Debugging
fieldnames(ns.Point)
isdefined(p, :x)


text:
using JuMP
using GLPK
using ModelGraphs

m1 = Model(with_optimizer(GLPK.Optimizer))

@variable(m1, xm[i in 1:2],Bin)
@constraint(m1, xm[1] + xm[2] <= 1)
@objective(m1, Max, 16xm[1] + 10xm[2])

m2 = Model(solver=GLPKSolverMIP())
@variable(m2, xs[i in 1:2],Bin)
@variable(m2, y[i in 1:2], Bin)
@constraint(m2, y[1] + y[2] <= 1)
@constraint(m2, 8xs[1] + 2xs[2] + y[1] + 4y[2] <= 10)
@objective(m2, Max, 4y[2])

## Model Graph
graph = ModelGraph()
heur(g) = 16
setsolver(graph, LagrangeSolver(update_method=:subgradient,max_iterations=30,lagrangeheuristic=heur))
n1 = add_node(graph)
setmodel(n1,m1)
n2 = add_node(graph)
setmodel(n2,m2)

## Linking
# m1[x] = m2[x]  ∀i ∈ {1,2}
@linkconstraint(graph, [i in 1:2], n1[:xm][i] == n2[:xs][i])

solution = solve(graph)

true


text:
export replace_sjreg, replace_sjhand, replace_sjhand2

function replace_sjreg(s::AbstractString, p::Pair{<:AbstractString,<:AbstractString}...)
    d = Dict(p...)
    buf = IOBuffer()
    frst = true
    for pair in p
        frst || print(buf, '|')
        frst = false
        escape_string(buf, first(pair), "[\\^\$.|?*+()")
    end
    r = Regex(String(take!(buf)))
    replace(s, r => (s -> d[s]))
end

function replace_sjhand(s::AbstractString)
    io = IOBuffer()
    for c in s
        if c == 'x'
            c = 'A'
        elseif c == 'b'
            c = 'B'
        elseif islower(c)
            c = uppercase(c)
        end
        write(io, c)
    end
    String(take!(io))
end

function replace_sjhand2(s::AbstractString)
    buf = IOBuffer()
    for c in s
        if c == 'c'
            print(buf, "BOOM")
        elseif !isupper(c)
            print(buf, c)
        end
    end
    String(take!(buf))
end


text:
#*************************************************************************#
#
# Predictor Error Models
#
#*************************************************************************#

"""
    narx(mdl, yterms, uterms, iddata, [init, kwargs...]) -> Θ, results, e

Compute a NARX model with the given structure for the data
provided in `iddata`.

Estimate the parameters minimizing the one-step-ahead prediction error,
between the predicted output and the measured one.

Use the given initializer `init` to compute the initial guess used
by the optimization algorithm. Aditional arguments will be passed
on the to the nonlinear least squares solver (levenberg-marquadt).

Return the estimated parameter vector Θ and a structure `results`
containing information about the optimization procedure.
"""
function narx(mdl::ParametricModel, yterms::Vector{Vector{Int}},
              uterms::Vector{Vector{Int}}, iddata::IdData;
              init::Initializer=DefaultInitializer(), kwargs...)

    pred = one_step_ahead(mdl, yterms, uterms, iddata)

    e = SimpleError(pred, iddata)

    Θ, results = estimate_parameters(e, init; kwargs...)

    return Θ, results, e
end

"""
    noe(mdl, yterms, uterms, iddata, [init, use_extended,  kwargs...]) -> Θ, results, e

Compute a NOE  model with the given structure for the data
provided in `iddata`.

Estimate the parameters minimizing the free-run simulation error,
between the predicted output and the measured one.

If `extended_parameter` is true include initial conditions
on the optimization problem. Otherwise the measured values
will be used as initial conditions. By default it is false.

Use the given Initializer `init` to compute the initial guess
used by the optimization algorithm. Aditional arguments will
be passed on the to the nonlinear least squares solver
(levenberg-marquadt).

Return the estimated parameter vector Θ and a structure `results`
containing information about the optimization procedure.
"""
function noe(mdl::ParametricModel, yterms::Vector{Vector{Int}},
             uterms::Vector{Vector{Int}}, iddata::IdData;
             init::Initializer=DefaultInitializer(),
             use_extended::Bool=false,
             kwargs...)
    pred = free_run_simulation(mdl, yterms, uterms, iddata)

    if use_extended
        e = ExtendedParamsError(pred, iddata)
    else
        e = SimpleError(pred, iddata)
    end

    Θ, results = estimate_parameters(e, init; kwargs...)

    return Θ, results, e
end

#*************************************************************************#
#
# PredictorError
#
#*************************************************************************#
abstract type PredictorError end

function estimate_parameters(pred_error::PredictorError, init::Initializer=DefaultInitializer();
                             kwargs...)
    # Get Initial Guess
    ϕ₀ = initial_guess(pred_error, init)

    # Error Lenght
    len = pred_error.pred.length*pred_error.pred.dynmdl.ny

    # Create Buffer
    e_buffer = Vector{Float64}(len)

    # Optimize
    f!(ϕ, storage) = error_function!(pred_error, storage, ϕ)
    g!(ϕ, storage) = error_function!(pred_error, e_buffer, ϕ; dϕ=storage)
    results = levenberg_marquardt(f!, g!, ϕ₀, len; kwargs...)

    Θ = get_Θ(pred_error, Optim.minimizer(results))
    return Θ, results
end

#*************************************************************************#
#
# SimpleError
#
#*************************************************************************#
struct SimpleError <: PredictorError
    pred::Predictor
    iddata::IdData
end

get_Θ{R<:Real}(e::SimpleError, Θ::Vector{R}) = Θ

function error_function!{R<:Real,
                         S<:Real,
                         T<:Real}(pred_error::SimpleError, e::AbstractVector{R},
                                  ϕ::AbstractVector{S};
                                  dϕ::AbstractMatrix{T}=Matrix{Float64}(0, 0))
    # Get info
    pred = pred_error.pred
    iddata = pred_error.iddata
    nparams = length_params(pred.dynmdl.mdl)
    ny = pred.dynmdl.ny
    n = pred.dynmdl.n
    n0 = pred.start
    nd = pred.length

    Θ = get_Θ(pred_error, ϕ)

    # Reshape buffers
    y = reshape(e, ny, nd)
    if !isempty(dϕ)
        dΘ = reshape(dϕ, ny, nd, nparams)
    else
        dΘ = Array{Float64}(0, 0, 0)
    end

    # Compute output and its jacobian
    predict!(pred, y, Θ; dΘ=dΘ)

    # Measured values
    y_measured = get_slice(pred, iddata.y)
    y .-= y_measured

    return
end

function initial_guess(e::SimpleError, init::Initializer=DefaultInitializer())
    initial_guess(e.pred.dynmdl.mdl, init)
end

#*************************************************************************#
#
# ExtendedParamsError
#
#*************************************************************************#
struct ExtendedParamsError <: PredictorError
    pred::Predictor
    iddata::IdData
end

function get_Θ{R<:Real}(e::ExtendedParamsError, ϕ::Vector{R})
    # Get parameters length
    nparams = length_params(e.pred.dynmdl.mdl)

    # Get parameters from extended parameters
    Θ = ϕ[1:nparams]

    return Θ
end

function get_y0{R<:Real}(e::ExtendedParamsError, ϕ::Vector{R})
    # Get initial condition length
    ny, n = size(e.pred.y0)

     # Get parameters length
    nparams = length_params(e.pred.dynmdl.mdl)

    # Get initial conditions from extended parameters
    y0 = reshape(ϕ[nparams+1:end], ny, n)

    return y0
end

function assemble_ϕ{R<:Real, S<:Real}(e::ExtendedParamsError,
                                      Θ::Vector{R},
                                      y0::Matrix{S})
    y0 = vec(y0)
    ϕ = [Θ; y0]
    return ϕ
end

function error_function!{R<:Real,
                         S<:Real,
                         T<:Real}(pred_error::ExtendedParamsError, e::AbstractVector{R},
                                  ϕ::AbstractVector{S};
                                  dϕ::AbstractMatrix{T}=Matrix{Float64}(0, 0))
    # Get info
    pred = pred_error.pred
    iddata = pred_error.iddata
    nparams = length_params(pred.dynmdl.mdl)
    ny = pred.dynmdl.ny
    n = pred.dynmdl.n
    n0 = pred.start
    nd = pred.length

    # Get initial conditions and parameters from extended parameters
    Θ = get_Θ(pred_error, ϕ)
    y0 = get_y0(pred_error, ϕ)

    # Reshape buffers
    y = reshape(e, ny, nd)
    if !isempty(dϕ)
        dΘ = view(dϕ, :, 1:nparams)
        dΘ = reshape(dΘ, ny, nd, nparams)

        dy0 = view(dϕ, :, nparams+1:nparams+n*ny)
        dy0 = reshape(dy0, ny, nd, ny, n)
    else
        dΘ = Array{Float64}(0, 0, 0)
        dy0 = Array{Float64}(0, 0, 0, 0)
    end

    # Compute output and its jacobian
    predict!(pred, y, Θ, y0; dΘ=dΘ, dy0=dy0)

    # Measured values
    y_measured = get_slice(pred, iddata.y)
    y .-= y_measured

    return
end

function initial_guess(e::ExtendedParamsError, init::Initializer=DefaultInitializer())
    Θ = initial_guess(e.pred.dynmdl.mdl, init)

    return assemble_ϕ(e, Θ, e.pred.y0)
end


text:
# Provides dynamic risk measures for SDEs

"""
    dynamicRM(problem::SDEProblem,
                eval = x-> sum(x),
                RM=1.0,
                dt = 0.01,
                MCs = 200,
                iter = 100,
                u0::Flux.Chain=Flux.Chain() )
implements the unique dynamic risk measure for the SDE ``problem``associated to the nonlinear term
```math
s_\\rho(t,x) \\lVert \\sigma(t,x) \\nabla u(t,x) \\rVert.
```
For more details on dynamic risk measures for SDEs, its relationship to BSDEs and parabolic PDEs as well as discrete time dynamic risk measures
see the Section [BSDE]().
"""
function dynamicRM(problem::SDEProblem,eval = x-> sum(x), RM=1.0, dt = 0.01, MCs = 200, iter = 100,
                    u0::Flux.Chain=Flux.Chain() )

    time_steps = div(problem.tspan[2]-problem.tspan[1],dt, RoundNearest)

    f(X,u,σᵀ∇u,p,t) = RM*abs.(σᵀ∇u)[1]                 # nonlinear riskMeasure part

    BSDE = TerminalPDEProblem(eval, f, problem.f, problem.g, problem.u0, problem.tspan);

    d = length(problem.u0)
    hls  = 10 + d #hide layer size
    opt = Flux.ADAM(0.001)

    if (u0 == Flux.Chain())
        u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
    end

    σᵀ∇u = [Flux.Chain(Dense(d,hls,relu),
        Dense(hls,hls,relu),
        Dense(hls,d)) for i in 1:time_steps]

    pdealg = NNPDEHan(u0,σᵀ∇u;opt=Flux.ADAM(0.1))
    ans = solve(BSDE, pdealg, verbose=true, maxiters=iter, trajectories=MCs, dt=dt)

    return(ans)
end

# Provide a jump diffusion version of dynamicRM


text:
"""
Case 2:
This case study a three bus system with 2 machines (One d- One q-: 4th order model) and an infinite source.
The fault drop the connection between buses 1 and 3, eliminating the direct connection between the infinite source
and the generator located in bus 3.
"""

##################################################
############### LOAD DATA ########################
##################################################

############### Data Network ########################

nodes_case234   = [Bus(1 , "Bus 1"  , "REF" , 0 , 1.02  , (min=0.94, max=1.06), 138),
                    Bus(2 , "Bus 2"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138),
                    Bus(3 , "Bus 3"  , "PV" , 0 , 1.00 , (min=0.94, max=1.06), 138)]

branch_case234  =  [Line("Line1", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                    Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case234[2], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

#Trip of Line 1.
branch_case234_fault = [Line("Line2", true, 0.0, 0.0, Arc(from=nodes_case234[1], to=nodes_case234[2]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04),
                        Line("Line3", true, 0.0, 0.0, Arc(from=nodes_case234[2], to=nodes_case234[3]), 0.01, 0.12, (from=0.0, to=0.0), 100, 1.04)]

loads_case234 =   [PowerLoad("Bus1", true, nodes_case234[1], PowerSystems.ConstantPower, 1.5, 0.8, 1.5, 0.8),
                   PowerLoad("Bus2", true, nodes_case234[2], PowerSystems.ConstantPower, 1.5, 0.7, 1.5, 0.8),
                   PowerLoad("Bus3", true, nodes_case234[3], PowerSystems.ConstantPower, 0.5, 0.3, 0.5, 0.3)]

############### Data devices ########################

inf_gen_case234 = StaticSource(1, #number
                :InfBus, #name
                nodes_case234[1], #bus
                1.02, #VR
                0.0, #VI
                0.000001) #Xth

######## Machine Data #########

### Case 2: 4th Order Model with AVR (3-bus case) ###
case2_machine2 =  OneDOneQMachine(0.0, #R
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  100.0)   #MVABase

case2_machine3 =  OneDOneQMachine(0.0, #R
                                  1.3125, #Xd
                                  1.2578, #Xq
                                  0.1813, #Xd_p
                                  0.25, #Xq_p
                                  5.89, #Td0_p
                                  0.6, #Tq0_p
                                  100.0)   #MVABase

######## Shaft Data #########

### Shafts for Gens 2 and 3: Cases 2, 3 and 4 ###
case234_shaft2 = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                            0.0) #D

case234_shaft3 = SingleMass(3.01, #H (M = 6.02 -> H = M/2)
                            0.0) #D

######## PSS Data #########
cases_no_pss = PSSFixed(0.0)

######## TG Data #########

### No TG for Cases 1, 2, 3, 4 ###
case1234_no_tg = TGFixed(1.0) #eff

########  AVR Data #########
### AVRs for Case 2, 3, 4 and 5 ###
case2345_avr2 = AVRTypeI(20.0, #Ka - Gain
                        0.01, #Ke
                        0.063, #Kf
                        0.2, #Ta
                        0.314, #Te
                        0.35, #Tf
                        0.001, #Tr
                        5.0, #Vrmax
                        -5.0, #Vrmin
                        0.0039, #Ae - 1st ceiling coefficient
                        1.555) #Be - 2nd ceiling coefficient

case2345_avr3 = AVRTypeI(20.0, #Ka - Gain
                        0.01, #Ke
                        0.063, #Kf
                        0.2, #Ta
                        0.314, #Te
                        0.35, #Tf
                        0.001, #Tr
                        5.0, #Vrmax
                        -5.0, #Vrmin
                        0.0039, #Ae - 1st ceiling coefficient
                        1.555) #Be - 2nd ceiling coefficient

### Case 2 Generators ###
case2_gen2 = DynGenerator(1, #Number
                         :Case2Gen2,
                         nodes_case234[2], #bus
                         1.0, # ω_ref,
                         1.0142, #V_ref
                         1.0, #P_ref
                         case2_machine2, #machine
                         case234_shaft2, #shaft
                         case2345_avr2, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

case2_gen3 = DynGenerator(2, #Number
                         :Case2Gen3,
                         nodes_case234[3], #bus
                         1.0, # ω_ref,
                         1.0059, #V_ref
                         1.0, #P_ref
                         case2_machine3, #machine
                         case234_shaft3, #shaft
                         case2345_avr3, #avr
                         case1234_no_tg, #tg
                         cases_no_pss) #pss

######################### Dynamical System ########################

case2_DynSystem = DynamicSystem(nodes_case234,
                                branch_case234,
                                [case2_gen2, case2_gen3],
                                vcat(inf_gen_case234,loads_case234),
                                100.0,
                                60.0)

##################################################
############### SOLVE PROBLEM ####################
##################################################

#Compute Y_bus after fault
Ybus_fault = PSY.Ybus(branch_case234_fault, nodes_case234)[:,:]

#Initialize variables
dx0 = zeros(LITS.get_total_rows(case2_DynSystem))
x0 = [1.02, 1.0, 1.0, 0.0, -0.01, -0.01,
      1.0, #eq_p
      0.47, #ed_p
      0.6, #δ
      1.0, #ω
      2.1, #Vf
      0.28, #Vr1
      -0.39, #Vr2,
      1.0, #Vm
      0.81, #eq_p
      0.59, #ed_p
      0.86, #δ
      1.0, #ω
      1.7, #Vf
      0.11, #Vr1
      -0.31, #Vr2,
      1.0] #Vm
tspan = (0.0, 30.0);

#Find initial condition
inif! = (out,x) -> LITS.system_model!(out, dx0 ,x, (Ybus_fault,case2_DynSystem), 0.0)
sys_solve = nlsolve(inif!, x0)
x0_init = sys_solve.zero
cb = DiffEqBase.DiscreteCallback(LITS.change_t_one, LITS.Y_change!)

sim = DynamicSimulation(case2_DynSystem, tspan, Ybus_fault, cb, x0_init)

#Solve problem in equilibrium
run_simulation!(sim, IDA());

#Obtain data for angles
series = get_state_series(sim, (:Case2Gen2, :δ));

@test sim.solution.retcode == :Success


text:
using Pkg.TOML
using Markdown

const all_config = open(TOML.parse, joinpath(@__DIR__, "pretrains.toml"))

function description(config::Dict)
  desc = config["description"]
  host = config["host"]
  link = config["link"]
  if haskey(config, "cite")
    cite = config["cite"]
    description(desc, host, link, cite)
  else
    description(desc, host, link)
  end
end
function description(description::String, host::String, link::String, cite=nothing)
  """
  $description
  Released by $(host) at $(link).
  $(isnothing(cite) ? "" : "\nCiting:\n$cite")"""
end

subword_type(config::Dict) = subword_type(Val(Symbol(config["subword"])), config)
function subword_type(::Val{:bpe}, config)
  vocab = config["vocab"]
  case = config["case"]
  "subword with $(vocab)-vocabularies-sized $(case) bpe tokenizer."
end
function subword_type(::Val{:wordpiece}, config)
  case = config["case"]
  "subword with google wordpiece $case tokenizer."
end

function model_summary(mt::String, config::Dict)
  head = config["head"]
  hidden = config["hidden"]
  layer = config["layer"]
  tk_desc = subword_type(config)
  desc = get(config, "desc", nothing)
  """
  $mt model with $(layer) layers, $(hidden) hidden units, and $(head) attention heads.
  $tk_desc$(isnothing(desc) ? "" : "\n$desc" )"""
end

function register_config(configs)
  for (model_type, model_config) in pairs(configs)
    models = filter(!isequal("description"), keys(model_config))
    mt_desc = model_config["description"]
    for model in models
      config = model_config[model]
      model_configs = config["models"]
      model_desc = description(config)

      for (model_name, model_detail) in model_configs
        summary = model_summary(model_type, model_detail)
        depname = "$(uppercase(model_type))-$(model_name)"
        depdesc = join([mt_desc, summary, model_desc], "\n\n")
        url = model_detail["url"]
        checksum = model_detail["checksum"]
        dep = DataDep(depname, depdesc, url, checksum; fetch_method=download_gdrive)
        DataDeps.register(dep)
      end
    end
  end
end

match(model, name, query) = isequal(lowercase(query), model) || isequal(query, name)

"""
  pretrains(query::String = ""; detailed::Bool = false)

Show all available models. you can also query a specific `model` or `model name`.
show more detail with `detailed = true`.
"""
function pretrains(query::String = ""; detailed::Bool = false)
  global all_config
  rows = [Any["Type", "model", "model name", "support items"]]
  detailed &&
    push!(rows[1], "detail description")

  for (model_type, model_config) in pairs(all_config)
    models = filter(!isequal("description"), keys(model_config))
    mt_desc = model_config["description"]
    for model in models
      config = model_config[model]
      model_configs = config["models"]
      model_desc = description(config)

      for (model_name, model_detail) in model_configs
        summary = model_summary(model_type, model_detail)
        name = model_detail["name"]
        items = join(model_detail["items"], ", ")
        if isempty(query) || match(model, name, query)
          row = Any[uppercasefirst(model_type), model, name, items]
        else
          continue
        end
        detailed &&
          push!(row, summary)

        push!(rows, row)
      end
    end
  end
  Markdown.MD(Any[Markdown.Table(rows, detailed ? [:l,:l,:l,:l,:l] : [:l,:l,:l,:l])])
end

function parse_model(str)
  m = Base.match(r"([a-zA-Z]+)[-:]([a-zA-Z\d\-_]+):?(\S*)", str)
  type, name, item = m.captures
  item = isempty(item) ? :all : Symbol(item)
  return lowercase(type), name, item
end


text:
using ForwardDiff
using LinearAlgebra
using Optim
using Pathfinder
using Test

include("test_utils.jl")

function lbfgs_inverse_hessian_explicit(H₀, S, Y)
    B = [H₀ * Y S]
    R = triu(S'Y)
    E = Diagonal(R)
    D = [0*I -inv(R); -inv(R)' R' \ (E + Y' * H₀ * Y)/R]
    return H₀ + B * D * B'
end

@testset "L-BFGS inverse Hessian construction" begin
    @testset "lbfgs_inverse_hessian" begin
        #! format: off
        S0 = [[1.719573, 3.294037, -2.008877, 3.901275, 0.214324, 0.400382, 0.113598, -1.804262, 0.465563, 2.465748], [-0.445476, -0.514915, 1.069617, -1.505506, -0.036535, -0.11386, -0.029737, 0.579662, -0.118694, -1.726258], [-0.013354, 0.1981, 0.081886, -0.194172, 0.010499, -0.007944, -0.001039, 0.061604, -0.002621, 0.11145], [0.00648, 0.255961, -0.011901, -0.059042, 0.013587, -0.002064, 0.000302, 0.023457, 0.00257, -0.002343], [-0.016408, 0.015005, 0.009654, 0.016735, -0.000245, -0.002825, -0.00106, 0.004272, -0.004578, 0.002275]]
        Y0 = -[[-2.357935, -3.343312, 4.659008, -7.065433, -0.228045, -0.584164, -0.156837, 2.863289, -0.631753, -7.152021], [0.610851, 0.522617, -2.480668, 2.726559, 0.038874, 0.166123, 0.041055, -0.9199, 0.161064, 5.007094], [0.018312, -0.201064, -0.189911, 0.351657, -0.011171, 0.011591, 0.001434, -0.097763, 0.003557, -0.323265], [-0.008886, -0.25979, 0.027601, 0.106929, -0.014456, 0.003011, -0.000418, -0.037225, -0.003488, 0.006795], [0.022499, -0.015229, -0.02239, -0.030308, 0.00026, 0.004122, 0.001463, -0.00678, 0.006212, -0.006598]]
        #! format: on
        S = reduce(hcat, S0)
        Y = reduce(hcat, Y0)
        N, history_length = size(S)
        α = rand(N)
        H₀ = Diagonal(α)

        @test @inferred(Pathfinder.lbfgs_inverse_hessian(H₀, S, Y, 0, 0)) ≈ H₀

        H = Pathfinder.lbfgs_inverse_hessian(H₀, S, Y, 3, 3)
        Hexp = lbfgs_inverse_hessian_explicit(H₀, S[:, 1:3], Y[:, 1:3])
        @test H ≈ Hexp

        S2 = S[:, [4:history_length; 1:3]]
        Y2 = Y[:, [4:history_length; 1:3]]
        ilast = argmax(axes(S, 2)[[4:history_length; 1:3]])
        H = Pathfinder.lbfgs_inverse_hessian(H₀, S2, Y2, ilast, history_length)
        Hexp = lbfgs_inverse_hessian_explicit(H₀, S, Y)
        @test H ≈ Hexp

        H = Pathfinder.lbfgs_inverse_hessian(H₀, S, Y, history_length, history_length)
        Hexp = lbfgs_inverse_hessian_explicit(H₀, S, Y)
        @test H ≈ Hexp
    end

    @testset "lbfgs_inverse_hessians" begin
        n = 10
        history_length = 5
        logp(x) = logp_banana(x)
        ∇logp(x) = ForwardDiff.gradient(logp, x)
        nocedal_wright_scaling(α, s, y) = fill!(similar(α), dot(y, s) / sum(abs2, y))
        θ₀ = 10 * randn(n)

        optimizer = Optim.LBFGS(;
            m=history_length, linesearch=Optim.LineSearches.MoreThuente()
        )
        θs, logpθs, ∇logpθs = Pathfinder.maximize_with_trace(logp, ∇logp, θ₀, optimizer)

        # run lbfgs_inverse_hessians with the same initialization as Optim.LBFGS
        Hs = Pathfinder.lbfgs_inverse_hessians(
            θs, ∇logpθs; history_length=history_length, Hinit=nocedal_wright_scaling
        )
        ss = diff(θs)
        ps = (Hs .* ∇logpθs)[1:(end - 1)]
        # check that next direction computed from Hessian is the same as the actual
        # direction that was taken
        @test all(≈(1), dot.(ps, ss) ./ norm.(ss) ./ norm.(ps))
    end
end