text:
struct Marginalize{T} end
Base.@pure Marginalize(marg...) = Marginalize{marg}()

struct Marginal{T} end
Base.@pure Marginal(marg...) = Marginal{marg}()

struct SparseGrid{T,N,L} <: AbstractArray{T,N}
    nodes::Vector{T}
end
@generated grid_length(::Val{N}, ::Val{L}) where {N,L} = _grid_length(Val{N}(), L)
@generated function _grid_length(::Val{N}, L::Int) where {N}
    quote
        ilim_0 = L
        s_0 = 0
        ind = 0
        out = 0
        @nloops $N i p -> begin
            0:ilim_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            ilim_{$N-p+1} = L - s_{$N-p+1}
        end begin
            @nloops $N j p -> begin
                1:rule_lengths[i_p+1]
            end begin
                out += 1
            end
        end
        out
    end
end

@generated function SparseGrid(f, data, ::Val{N}, ::Val{L}, ::Type{T} = Float64) where {T,N,L}
    quote
        ilim_0 = $L
        s_0 = 0
        ind = 0
        out = Vector{$T}(grid_length(Val{$N}(), Val{$L}()))
        @nloops $N i p -> begin
            0:ilim_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            ilim_{$N-p+1} = $L - s_{$N-p+1}
        end begin
            #Maxmimum length left over to the marg-outs
            # max_l = $L - sum( $tupj_marg_course )
            # l_remaining = max_l - sum( $tupj_margout_course )
            l_remaining = $L - sum( ( @ntuple $N i ) )
            @nloops $N j p -> begin
                HermiteQuadratureRules.gk_nodes[i_p+1]
            end begin
                ind += 1
                out[ind] = f(data, (@ntuple $N j)...)
            end
        end
        SparseGrid{$T,$N,$L}(out)
    end
end

const HasTypedLength{N,T} = Union{NTuple{N,T}, SVector{N,T}}
function grid_adjust(node::HasTypedLength{N,T}, center::HasTypedLength{N,T}, scale::HasTypedLength{N,T}) where {N,T}
    @. node * scale + center
end

@generated function SparseGrid(f, data, ::Val{N}, ::Val{L}, neg_max::T, center::HasTypedLength{N,T}, scale::HasTypedLength{N,T}) where {T,N,L}
    quote
        ilim_0 = $L
        s_0 = 0
        ind = 0
        out = Vector{$T}(undef,grid_length(Val{$N}(), Val{$L}()))
        @nloops $N i p -> begin
            0:ilim_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            ilim_{$N-p+1} = $L - s_{$N-p+1}
        end begin
            #Maxmimum length left over to the marg-outs
            # max_l = $L - sum( $tupj_marg_course )
            # l_remaining = max_l - sum( $tupj_margout_course )
            l_remaining = $L - sum( ( @ntuple $N i ) )
            @nloops $N j p -> begin
                HermiteQuadratureRules.gk_nodes_v[i_p+1]
            end begin
                ind += 1
                out[ind] = exp(neg_max + f(data, grid_adjust(SVector(@ntuple $N j), center, scale) ))
            end
        end

        SparseGrid{$T,$N,$L}(out)
    end
end

@generated function SparseGrid(::Type{T}, ::Val{N}, ::Val{L}) where {T <: Real, N, L}
    SparseGrid{T,N,L}(sparse_grid_mat(T, Val{N}(), L))
end
@generated function SparseGrid(::Type{A}, ::Val{N}, ::Val{L}) where {T, A <: AbstractArray{T}, N, L}
    SparseGrid{T,N,L}(A(sparse_grid_mat(T, Val{N}(), L)))
end

@generated function SparseGrid(sg::SparseGrid{T,Nold,L}, ::Marginal{marg}) where {T,Nold,L,marg}
    Nnew = length(marg)
    Ndrop = Nold - Nnew
    quote
        increment = 0
        @nexprs $Nd i -> begin
            while (i + increment) in retain
                increment += 1
            end
            drop_i = i + increment
        end
        drop = @ntuple $Nd i -> drop_i

        out = Matrix{T}($Nn, polycount(Val{Nn}(), L))
    end
end
@generated function SparseGrid(sg::SparseGrid{T,No,L}, retain::NTuple{Nn,I}) where {T,No,L,Nn,I<:Integer}
    Nd = No - Nn
    quote
        increment = 0
        @nexprs $Nd i -> begin
            while (i + increment) in retain
                increment += 1
            end
            drop_i = i + increment
        end
        drop = @ntuple $Nd i -> drop_i
        out = Matrix{T}($Nn, polycount(Val{Nn}(), L))
    end
end

sparse_grid_mat(::Type{T}, ::Val{N}, l::Int) where {T,N} = sparse_grid_mat!(Matrix{T}(N, polycount(Val{N}(), l)), Val{N}(), l)

@generated function sparse_grid_mat!(out::Matrix{T}, ::Val{N}, l::Int) where {T,N}
    quote
        j_0 = l
        l_0 = 1
        s_0 = 0
        ind = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N j p -> begin
                p == $N ? $hermite_nodes[i_p] : $hermite_nodes_diffs[i_p]
            end begin
                ind += 1
                @nexprs $N k -> begin
                    out[k,ind] = j_k
                end
            end
        end
        out
    end
end

function gen_sub(marg::NTuple{N,I}) where {N,I<:Integer}
    ntuple(i -> Symbol(), Val{N}())
end
function symbolify(x::NTuple{N,T}, sym::Symbol = :j_) where {N,T}
    ntuple( p -> Symbol(sym, x[p]) , Val{N}())
end

@generated function missing_inds(x::NTuple{N,T}, ::Val{Ntotal}) where {N,Ntotal,T <: Integer}
    Ndiff = Ntotal - N
    quote
        nfound = 0
        @nexprs $Ndiff i -> begin
            local ind_i
        end
        # for i ∈ 1:Ntotal
        i = 0
        while nfound < $Ndiff
            i += 1
            if i ∉ x
                nfound += 1
                @nif $Ndiff j -> j == nfound j -> ind_j = i
                # if nfound == $Ndiff
                #     break
                # end
            end
        end
       @ntuple $Ndiff ind
    end
end
@generated function MargDiff(::Marginal{marg}, ::Val{N}) where {marg,N}
    Marginalize{missing_inds(marg, Val{N}())}()
end
@generated function MargDiff(::Marginalize{margout}, ::Val{N}) where {margout,N}
    Marginal{missing_inds(margout, Val{N}())}()
end

@generated function marginalize(sg::SparseGrid{T,No,L}, ::Marginal{marg}) where {T,No,L, marg}
    _marginalize(sg, Marginal{marg}(), MargDiff( Marginal{marg}(), Val{No}() ) )
end
@generated function marginalize(sg::SparseGrid{T,No,L}, ::Marginalize{margout}) where {T,No,L, margout}
    _marginalize(sg, MargDiff( Marginalize{margout}(), Val{No}() ), Marginalize{margout}())
end

#@generated function _marginalize(sg::SparseGrid{T,No,L}, marg::NTuple{Nin}, margout::NTuple{Nout}) where {T,No,L, Nin, Nout}
#end
@generated function marginalize(sg::SparseGrid{T,No,L}, ::Marginal{marg}, ::Marginalize{margout}) where {T,No,L, marg, margout}
#    @assert length(marg) + length(margout) == No
    _marginalize(sg, Marginal{marg}(), Marginalize{margout}())
end

function marginal_symbol(course_inds::NTuple{N}, fine_inds::NTuple{N}, symbolstring = "m") where N
    for (i,j) in zip(course_inds, fine_inds)
        symbolstring *= "_" * string(i) * "_" * string(j)
    end
    Symbol(symbolstring)
end

function set_symbol_to_zero(course_inds::NTuple{N}, fine_inds::NTuple{N}, ::Type{T} = Float64, symbolstring = "m") where {N,T}
    symbol = marginal_symbol(course_inds, fine_inds, symbolstring)
    z = zero(T)
    :($symbol = $z)
end

## Need to replace Marginal{marg} and Marginalize{margout} with marg::NTuple and margout::NTuple
## add out and in_interp
# @generated function _marginalize(sg::SparseGrid{T,N,L}, marg::NTuple{Nin,Int}, margout::NTuple{Nout}) where {T,N,L,Nin,Nout}
@generated function _marginalize(::Type{SparseGrid{T,N,L}}, ::Marginal{marg}, ::Marginalize{margout}) where {T,N,L,marg,margout}
    Nin = length(marg)
    Nout = length(margout)
    # calc_l = :(sum($tupj_margout_course) - $(length(margout)))
    tupj_marg = Expr(:tuple, symbolify(marg, :j_)... ) #This tuple is for the fine indices we keep.
    tupj_margout = Expr(:tuple, symbolify(margout, :j_)... ) #This tuple gives the fine indices we're marginalizing out.
    tupj_marg_course = Expr(:tuple, symbolify(marg, :i_)... ) #This tuple gives the course indices we keep.
    tupj_margout_course = Expr(:tuple, symbolify(margout, :i_)...) #This tuple gives the course indices we're marginalizing out.
    quote
        ex = quote @fastmath begin in_interp = sg.nodes end end
        exa = ex.args[2].args[3].args# 0.7
        # exa = ex.args[2].args[2].args# 0.6

        ilim_0 = $L
        s_0 = 0

        marg_length = 0
        #Here we set all the marginal values to zero, so that later we can "+=" them without worry.
        @nloops $Nin i p -> begin
            0:ilim_{$Nin-p}
        end p -> begin
            s_{$Nin-p+1} = s_{$Nin-p} + i_p
            ilim_{$Nin-p+1} = $L - s_{$Nin-p+1}
        end begin
            @nloops $Nin j p -> begin
                1:rule_lengths[i_p+1]
            end begin
                marg_length += 1
                push!(exa,  set_symbol_to_zero( ( @ntuple $Nin i ) , ( @ntuple $Nin j ), T ) )
            end
        end

        ilim_0 = $L
        # s_0 = 0
        ind = 0

        @nloops $N i p -> begin
            0:ilim_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            ilim_{$N-p+1} = $L - s_{$N-p+1}
        end begin

            #Maxmimum length left over to the marg-outs
            # max_l = $L - sum( $tupj_marg_course )
            # l_remaining = max_l - sum( $tupj_margout_course )
            l_remaining = $L - sum( ( @ntuple $N i ) )

            @nloops $N j p -> begin
                1:rule_lengths[i_p+1]
            end begin
                ind += 1
                marg_symbol = marginal_symbol( $tupj_marg_course, $tupj_marg )
                weight_val = calc_weights( l_remaining, $tupj_margout_course, $tupj_margout )
                expr = :( $marg_symbol += in_interp[$ind] * $weight_val )
                push!(exa, expr)
            end
        end

        push!(exa, :(marginal_interp = Vector{$T}(undef, $marg_length)))
        marg_length = 0
        #Here we fill in the marginal vector.
        @nloops $Nin i p -> begin
            0:ilim_{$Nin-p}
        end p -> begin
            s_{$Nin-p+1} = s_{$Nin-p} + i_p
            ilim_{$Nin-p+1} = $L - s_{$Nin-p+1}
        end begin
            @nloops $Nin j p -> begin
                1:rule_lengths[i_p+1]
            end begin
                marg_length += 1
                marg_symbol = marginal_symbol( ( @ntuple $Nin i ), ( @ntuple $Nin j ) )
                push!(exa,  :(marginal_interp[$marg_length] = $marg_symbol ) )
            end
        end
        Nin = $Nin
        output = :(SparseGrid{$T,$Nin,$L}(marginal_interp))
        push!(exa, output )
        ex
    end
end

## Need to replace Marginal{marg} and Marginalize{margout} with marg::NTuple and margout::NTuple
## add out and in_interp
@generated function _marginalize_deriv(sg::SparseGrid{T,No,L}, ::Marginal{marg}, ::Marginalize{margout}) where {T,No,L, marg, margout}
    Nin = length(marg) #@ntuple $Nin p -> j_{ $marg[p] }
    Nout = length(margout)
    tupj_marg = Expr(:tuple, symbolify(marg)... ) #For sub2ind-ing
    tupj_margout = Expr(:tuple, symbolify(margout)... ) #For marginalizing out
    tupj_margout_course = Expr(:tuple, symbolify(margout, :i_)...)
    calc_l = :(sum($tupj_margout_course) - $(length(margout)))
    quote
        ex = quote @fastmath begin end end
        exa = ex.args[2].args[2].args

        j_0 = $L
        l_0 = 1
        s_0 = 0
        ind = 0
        @nloops $No i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = $L - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*$L - s_{$N-p+1} )
        end begin
            @nloops $N j p -> begin
#                max_l = p == $N ? length(hermite_nodes[i_p]) : length(hermite_nodes_diffs[i_p])
#                1:max_l

                max_l = length(hermite_nodes[i_p])
                if $No == 1
                    min_l = i_1 == 1 ? 1 : length(hermite_nodes[i_p-1])+1
                elseif p == $N
                    min_l = 1
                elseif i_p == 1
                    min_l = 1
                else
                    min_l = length(hermite_nodes[i_p-1])+1
                end
                min_l:max_l
            end begin
                ind += 1
           #     @nexprs $N k -> begin
                l = $L + $No - sum( ( @ntuple $No i ) )

                marg_ind = sg_sub2ind( Val{L}(), $tupj_marg )
                weight_val = calc_weight( Val{L}(), l, $tupj_margout, $tupj_margout_course)
                expr = :( temp_out = in_interp[1,$ind] * $weight_val )
                push!(exa, expr)
                for j ∈ 1:$Nout #jth derivative
                    deriv_ind = $margout[j]+1
                    weight = calc_weight( Val{L}(), l, $tupj_margout, $tupj_margout_course, j)
                    expr = :( temp_out += in_interp[$deriv_ind,$ind] * $weight )
                    push!(exa, expr)
                end
                push!(exa, :(out[1,$marg_ind] += temp_out))
                for j ∈ 2:$Nin+1
                    deriv_ind = $marg[j-1] + 1
                    push!(exa, :(out[$j,$marg_ind] += in_interp[$deriv_ind,$ind] * $weight_val ))
                end

         #           out[k,ind] = j_k
         #       end

            end
        end
        ex

    end
end

macro nnexprs(N, itersym, rangeexpr, args...)
    _nnexprs(N, itersym, rangeexpr, args...)
end
import Base.Cartesian: inlineanonymous
function _nnexprs(N, itersym, rangeexpr, args...)

#    body = args[end]
#    ex = Expr(:escape, body)
    ex = args
    for dim ∈ 1:N
        itervar = inlineanonymous(itersym, dim) #for :i, dim produces i_dim
        rng = inlineanonymous(rangeexpr, dim) #subs the dim appropriately into expression and simplifies
        ex = quote
            @nexprs $(esc(rng)) $(esc(itervar)) -> $ex
        end

    end
    ex
end
nnt = quote
        @nexprs 4 i_4 -> begin
            @nexprs 5 i_3 -> begin
                @nexprs 6 i_2 -> begin
                    @nexprs 7 i_1 -> begin

                        @ntuple 4 i
                    end
                end
            end
        end
    end

function Base.start(sg::SparseGrid{T,N}) where {T,N}

end
function Base.next(sg::SparseGrid{T,N}, state) where {T,N}

end
function Base.done(sg::SparseGrid{T,N}, state) where {T,N}

end
Base.IteratorSize(sg::SparseGrid) = HaseLength()
Base.IteratorEltype(sg::SparseGrid) = HasEltype()
Base.eltype(sg::SparseGrid{T,N}) where {T,N} = Tuple{Tuple{NTuple{N,T},T},SparseGrid{T,N}}

@generated function Base.size(sg::SparseGrid{T,N,L}) where {T,N,L}
    ntuple( i -> HermiteQuadratureRules.rules_total[L+1] , Val{N}())
end
# @generated function Base.length(sg::SparseGrid{T,N,L}) where {T,N,L}
#     polycount(Val{N}(), L)
# end

Base.length(sg::SparseGrid) = length(sg.nodes)

function Base.show(io::IO, sg::SparseGrid)
    println(io, "Sparse grid of size:", size(sg))
    FillDensity = length(sg.nodes) / prod(size(sg))
    println("The fill density: ", FillDensity)
    nothing
end

@inline function Base.getindex(sg::SparseGrid{T,N}, i::Int) where {T,N}
    ntuple( j -> sg.nodes[j,i], Val{N}())
end
@inline function Base.getindex(sg::SparseGrid{T,N}, sub::Vararg{I, N}) where {T,N,I<:Integer}

end

function Base.setindex!(sg::SparseGrid, v, i)

end

function Base.setindex!(sg::SparseGrid{T,N}, v, i::Vararg{I, N}) where {T,N,I<:Integer}

end

function Base.firstindex(sg::SparseGrid)

end

function Base.lastindex(sg::SparseGrid)

end

Base.IndexStyle(::SparseGrid) = IndexLinear()

@generated function SparseIterator(f, sg::SparseGrid{T,N}, args) where {T,N}

    quote
        j_0 = l
        s_0 = 0
        @nloops $N i p -> begin
            0:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p
            j_{$N-p+1} = l - s_{$N-p+1}
        end begin
            f(args, (@ntuple $N j -> T(i_{j}))... )
        end
        out

    end
end

@generated function polycount(::Val{N}, l::Int) where N
    quote
        l == 1 && return 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                p == $N ? $hermite_nodes[i_p] : $hermite_nodes_diffs[i_p]
            end begin
                out += 1
            end
        end
        out
    end
end
@generated function pc2(::Val{N}, l::Int) where N
    quote
        l == 1 && return 1
        j_0 = l
        l_0 = 1
        s_0 = 0
        out = 0
        @nloops $N i p -> begin
            l_{$N-p}:j_{$N-p}
        end p -> begin
            s_{$N-p+1} = s_{$N-p} + i_p - 1
            j_{$N-p+1} = l - s_{$N-p+1}
            l_{$N-p+1} = max(1, (3-p)*l - s_{$N-p+1} )
        end begin
            @nloops $N k p -> begin
                max_l = length(hermite_nodes[i_p])
                if $N == 1
                    min_l = i_1 == 1 ? 1 : length(hermite_nodes[i_p-1])+1
                elseif p == $N
                    min_l = 1
                elseif i_p == 1
                    min_l = 1
                else
                    min_l = length(hermite_nodes[i_p-1])+1
                end
                min_l:max_l
            end begin
                location = @ntuple $N d -> k_d
                @show  location
                out += 1
            end
        end
        out
    end
end


text:
# the fast radixsort for internedstrings
function radixsort(svec::AbstractVector{InternedString})
    l = length(InternedStrings.pool)

    # get the list of pointers
    pointer_pool = UInt.(pointer.(collect(keys(InternedStrings.pool))))
    # sorting the pointers will make it easy to look up
    sort!(pointer_pool)

    # get the pointer to the InternedStrings and sort them as well
    pointer_svec = UInt.(pointer.(svec))
    sort!(pointer_svec)

    # now both pointer to svec and pointer to pointer to pool arr
    cnts = zeros(UInt32, l)
    last_ps = pointer_svec[1]
    lo = 1
    @inbounds for i = 2:length(pointer_svec)
        next_ps = pointer_svec[i]
        if  next_ps != last_ps
            idx = searchsortedfirst(pointer_pool, last_ps)
            cnts[idx] += i - lo
            last_ps = next_ps
            lo = i
        end
    end

    idx = searchsortedfirst(pointer_pool, pointer_svec[end])
    cnts[idx] = length(pointer_svec) - lo + 1
    # cnts

    # pointer_pool, cnts

    # # # # load the strings based on pointers
    pstr = Base.unsafe_pointer_to_objref.(Ptr{UInt8}.(pointer_pool) .- 8)
    # # pstr, cnts

    # this index here will tell you how to sort
    # sort the strings and return the order
    # once you look up the position of pointer then can use this to look up
    # the right index
    indexes = fsortperm(pstr)
    pstr .= Base.unsafe_pointer_to_objref.(Ptr{UInt8}.(pointer_pool[indexes]) .- 8)
    # pstr, cnts[indexes]
    StatsBase.inverse_rle(pstr, cnts[indexes])

    # cnts .= cnts[indexes]
    # pstr .= pstr[indexes]
    # pstr, cnts

    # bad combination
    # res = StatsBase.inverse_rle(@view(pointer_pool[indexes]), @view(cnts[indexes]))
    # Base.unsafe_pointer_to_objref.(Ptr{UInt8}.(res) .- 8)
end

# this is slow due to string searchedsortedfirst not being very fast for strings
# but would be good at some point in the future to test again once it is faster
function radixsort1(svec::AbstractVector{InternedString})
    p = InternedStrings.pool

    pstr = collect(keys(p))
    radixsort!(pstr)

    cnts = zeros(UInt32, length(p))

    @inbounds for ss in svec
        idx = searchsortedfirst(pstr, ss.value)
        cnts[idx] += 1
    end
    cnts
end

ext: jl


text:
export CircularCircular, PriorCircular, PackedCircularCircular, PackedPriorCircular

"""
$(TYPEDEF)

Factor between two Sphere1 variables.

Related

[`Sphere1`](@ref), [`PriorSphere1`](@ref), [`Polar`](@ref), [`ContinuousEuclid`](@ref)
"""
mutable struct CircularCircular{T<: SamplableBelief} <: AbstractManifoldMinimize
  Z::T
  # Sphere1Sphere1(z::T=Normal()) where {T <: SamplableBelief} = new{T}(z)
end

const Sphere1Sphere1 = CircularCircular

CircularCircular(::UniformScaling) = CircularCircular(Normal())

DFG.getManifold(::CircularCircular) = RealCircleGroup()

function getSample(cf::CalcFactor{<:CircularCircular})
  # return (sampleTangent(getManifold(cf.factor), cf.factor.Z), )
  # FIXME workaround for issue #TBD with manifolds CircularGroup
  return sampleTangent(getManifold(cf.factor), cf.factor.Z, [0.0])

end

function (cf::CalcFactor{<:CircularCircular})(X, p, q)
  #
  M = getManifold(cf.factor)
  return distanceTangent2Point(M, X, p, q)
end

Base.convert(::Type{<:MB.AbstractManifold}, ::InstanceType{CircularCircular}) = Manifolds.RealCircleGroup()

"""
$(TYPEDEF)

Introduce direct observations on all dimensions of a Circular variable:

Example:
--------
```julia
PriorCircular( MvNormal([10; 10; pi/6.0], diagm([0.1;0.1;0.05].^2)) )
```

Related

[`Circular`](@ref), [`Prior`](@ref), [`PartialPrior`](@ref)
"""
mutable struct PriorCircular{T<: SamplableBelief} <: AbstractPrior
  Z::T
end

PriorCircular(::UniformScaling) = PriorCircular(Normal())

DFG.getManifold(::PriorCircular) = RealCircleGroup()

function getSample(cf::CalcFactor{<:PriorCircular})
  # FIXME workaround for issue #TBD with manifolds CircularGroup
  return samplePoint(getManifold(cf.factor), cf.factor.Z, [0.0])
end

Base.convert(::Type{<:MB.AbstractManifold}, ::InstanceType{PriorCircular}) = Manifolds.RealCircleGroup()

"""
$(TYPEDEF)

Serialized object for storing PriorCircular.
"""
mutable struct PackedPriorCircular  <: IncrementalInference.PackedInferenceType
  datastr::String
  # PackedPriorCircular() = new()
  # PackedPriorCircular(x::String) = new(x)
end
function convert(::Type{PackedPriorCircular}, d::PriorCircular)
  return PackedPriorCircular(convert(PackedSamplableBelief, d.Z))
end
function convert(::Type{PriorCircular}, d::PackedPriorCircular)
  distr = convert(SamplableBelief, d.datastr)
  return PriorCircular{typeof(distr)}(distr)
end

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

"""
$(TYPEDEF)

Serialized object for storing CircularCircular.
"""
mutable struct PackedCircularCircular  <: IncrementalInference.PackedInferenceType
  datastr::String
  # PackedCircularCircular() = new()
  # PackedCircularCircular(x::String) = new(x)
end
function convert(::Type{CircularCircular}, d::PackedCircularCircular)
  return CircularCircular(convert(SamplableBelief, d.datastr))
end
function convert(::Type{PackedCircularCircular}, d::CircularCircular)
  return PackedCircularCircular(convert(PackedSamplableBelief, d.Z))
end

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


text:
@json_read mutable struct TextDocumentIdentifier
    uri::DocumentUri
end

@json_read mutable struct TextDocumentItem
    uri::DocumentUri
    languageId::String
    version::Int
    text::String
end

@json_read mutable struct VersionedTextDocumentIdentifier
    uri::DocumentUri
    version::Int
end

@json_read mutable struct TextDocumentPositionParams
    textDocument::TextDocumentIdentifier
    position::Position
end

@json_read mutable struct DocumentFilter
    language::Union{Nothing,String}
    scheme::Union{Nothing,String}
    pattern::Union{Nothing,String}
end

const DocumentSelector = Vector{DocumentFilter}

mutable struct TextDocumentRegistrationOptions
    documentSelector::DocumentSelector
end

mutable struct TextDocumentChangeRegistrationOptions
    documentSelector::DocumentSelector
    syncKind::Int
end

mutable struct TextDocumentEdit
    textDocument::VersionedTextDocumentIdentifier
    edits::Vector{TextEdit}
end

mutable struct WorkspaceEdit
    changes
    documentChanges::Union{Nothing,Vector{TextDocumentEdit}}
end

@json_read mutable struct DidOpenTextDocumentParams
    textDocument::TextDocumentItem
end

@json_read mutable struct TextDocumentContentChangeEvent
    range::Union{Nothing,Range}
    rangeLength::Union{Nothing,Int}
    text::String
end

@json_read mutable struct DidChangeTextDocumentParams
    textDocument::VersionedTextDocumentIdentifier
    contentChanges::Vector{TextDocumentContentChangeEvent}
end

@json_read mutable struct DidSaveTextDocumentParams
    textDocument::TextDocumentIdentifier
    text::Union{Nothing,String}
end

@json_read mutable struct DidCloseTextDocumentParams
    textDocument::TextDocumentIdentifier
end

mutable struct FileEvent
    uri::String
    _type::Int
end
FileEvent(d::Dict) = FileEvent(d["uri"], d["type"])

@json_read mutable struct DidChangeWatchedFilesParams
    changes::Vector{FileEvent}
end

mutable struct FileSystemWatcher
    globPattern::String
    kind::Union{Nothing,Int}
end

mutable struct DidChangeWatchedFilesRegistrationOptions
    watchers::Vector{FileSystemWatcher}
end

mutable struct WillSaveTextDocumentParams
    textDocument::TextDocumentIdentifier
    reason::Int
end


text:
# First millisecond of 2015.
const DISCORD_EPOCH = 1420070400000

# Discord's form of ID.
const Snowflake = UInt64

snowflake(s::Integer) = Snowflake(s)
snowflake(s::AbstractString) = parse(Snowflake, s)

# TODO: Put these in helpers?
snowflake2datetime(s::Snowflake) = unix2datetime(((s >> 22) + DISCORD_EPOCH) / 1000)
worker_id(s::Snowflake) = (s & 0x3E0000) >> 17
process_id(s::Snowflake) = (s & 0x1F000) >> 12
increment(s::Snowflake) = s & 0xFFF

# Discord sends both Unix and ISO timestamps.
datetime(s::Int) = unix2datetime(s / 1000)
datetime(s::AbstractString) = DateTime(replace(s, "+" => ".000+")[1:23], ISODateTimeFormat)

lowered(x::DateTime) = round(Int, datetime2unix(x))
lowered(x::Union{Integer, Bool}) = x
lowered(x::Vector) = lowered.(x)
lowered(x::Nothing) = nothing
lowered(x) = JSON.lower(x)

macro lower(T)
    if supertype(eval(T)) <: Enum{<:Integer}
        quote
            JSON.lower(x::$T) = Int(x)
        end
    else
        quote
            function JSON.lower(x::$T)
                d = Dict{String, Any}()

                for f in fieldnames($T)
                    v = getfield(x, f)
                    if !ismissing(v)
                        d[string(f)] = lowered(v)
                    end
                end

                return d
            end
        end
    end
end

macro merge(T)
    quote
        function Base.merge(a::$T, b::$T)
            vals = []

            for f in fieldnames($T)
                va = getfield(a, f)
                vb = getfield(b, f)
                push!(vals, ismissing(vb) ? va : vb)
            end

            return $T(vals...)
        end
        Base.merge(::Missing, x::$T) = x
        Base.merge(x::$T, ::Missing) = x
    end
end

field(k::String, ::Type{Any}) = :(d[$k])
field(k::String, ::Type{Snowflake}) = :(snowflake(d[$k]))
field(k::String, ::Type{DateTime}) = :(datetime(d[$k]))
field(k::String, ::Type{T}) where T = :($T(d[$k]))
field(k::String, ::Type{Vector{Snowflake}}) = :(snowflake.(d[$k]))
field(k::String, ::Type{Vector{DateTime}}) = :(datetime.(d[$k]))
field(k::String, ::Type{Vector{T}}) where T = :($T.(d[$k]))
function field(k::String, ::Type{Union{T, Missing}}) where T
    return :(haskey(d, $k) ? $(field(k, T)) : missing)
end
function field(k::String, ::Type{Union{T, Nothing}}) where T
    return :(d[$k] === nothing ? nothing : $(field(k, T)))
end
function field(k::String, ::Type{Union{T, Nothing, Missing}}) where T
    return :(haskey(d, $k) ? $(field(k, Union{T, Nothing})) : missing)
end

macro dict(T)
    TT = eval(T)
    args = map(f -> field(string(f), fieldtype(TT, f)), fieldnames(TT))

    quote
        function $(esc(T))(d::Dict{String, Any})
            $(esc(T))($(args...))
        end
    end
end

function doctype(s::String)
    s = replace(s, "UInt64" => "Snowflake")
    s = replace(s, "Int64" => "Int")
    s = replace(s, "Discord." => "")
    s = replace(s, "Dates." => "")
    m = match(r"Array{([^{}]+),1}", s)
    m === nothing || (s = replace(s, m.match => "Vector{$(m.captures[1])}"))
    m = match(r"Union{Missing, Nothing, (.+)}", s)
    m === nothing || return replace(s, m.match => "Union{$(m.captures[1]), Missing, Nothing}")
    m = match(r"Union{Missing, (.+)}", s)
    m === nothing || return replace(s, m.match => "Union{$(m.captures[1]), Missing}")
    m = match(r"Union{Nothing, (.+)}", s)
    m === nothing || return replace(s, m.match => "Union{$(m.captures[1]), Nothing}")
    return s
end

macro fielddoc(T)
    TT = eval(T)
    ns = collect(string.(fieldnames(TT)))
    width = maximum(length, ns)
    map!(n -> rpad(n, width), ns, ns)
    ts = collect(map(f -> string(fieldtype(TT, f)), fieldnames(TT)))
    map!(doctype, ts, ts)
    docs = join(map(t -> "$(t[1]) :: $(t[2])", zip(ns, ts)), "\n")

    quote
        doc = string(@doc $T)
        docstring = doc * "\n# Fields\n\n```\n" * $docs * "\n```\n"

        Base.CoreLogging.with_logger(Base.CoreLogging.NullLogger()) do
            @doc docstring $T
        end
    end
end

macro boilerplate(T, exs...)
    macros = map(e -> e.value, exs)

    quote
        @static if :docs in $macros
            @fielddoc $T
        end
        @static if :dict in $macros
            @dict $T
        end
        @static if :lower in $macros
            @lower $T
        end
        @static if :merge in $macros
            @merge $T
        end
    end
end

include(joinpath("types", "overwrite.jl"))
include(joinpath("types", "role.jl"))
include(joinpath("types", "guild_embed.jl"))
include(joinpath("types", "attachment.jl"))
include(joinpath("types", "voice_region.jl"))
include(joinpath("types", "activity.jl"))
include(joinpath("types", "embed.jl"))
include(joinpath("types", "user.jl"))
include(joinpath("types", "ban.jl"))
include(joinpath("types", "integration.jl"))
include(joinpath("types", "connection.jl"))
include(joinpath("types", "emoji.jl"))
include(joinpath("types", "reaction.jl"))
include(joinpath("types", "presence.jl"))
include(joinpath("types", "channel.jl"))
include(joinpath("types", "webhook.jl"))
include(joinpath("types", "invite_metadata.jl"))
include(joinpath("types", "member.jl"))
include(joinpath("types", "voice_state.jl"))
include(joinpath("types", "message.jl"))
include(joinpath("types", "guild.jl"))
include(joinpath("types", "invite.jl"))
include(joinpath("types", "audit_log.jl"))


text:
using ParallelAccelerator

@acc function fdtd(tmax::Int, nx::Int, ny::Int, nxm1::Int, nym1::Int,
                   ex::Matrix{Float64}, ey::Matrix{Float64},
                   hz::Matrix{Float64}, _fict_::Vector{Float64})

     for t = 1:tmax
         for j = 1:ny
             ey[1, j] = _fict_[j]
         end
         for i = 2:nx
             for j = 1:ny
                 ey[i, j] -= 0.5 * (hz[i, j] - hz[i - 1, j])
             end
         end
         for i = 1:nx
             for j = 1:ny
                 ex[i, j] -= 0.5 * (hz[i, j] - hz[i, j - 1])
             end
         end
         for i = 1:nxm1
             for j = 1:nym1
                 hz[i, j] -= 0.7 * (ex[i, j + 1] - ex[i, j] + ex[i + 1, j] - ey[i, j])
             end
         end
     end
end

tmax = 50
nx = 50
ny = nx

nxm1 = nx - 1
nym1 = ny - 1

_fict_ = rand(tmax)
ex1 = rand(nx, ny)
ex2 = copy(ex1)
ey1 = rand(nx, ny)
ey2 = copy(ey1)
hz1 = rand(nx, ny)
hz2 = copy(hz1)

fdtd(tmax, nx, ny, nxm1, nym1, ex1, ey1, hz1, _fict_)

max_stars_repo_name: malekbr/julia-to-tiramisu

max_issues_repo_name: malekbr/julia-to-tiramisu

max_forks_repo_name: malekbr/julia-to-tiramisu


text:
"""

  given a path p and a flattened scattering transform flat, return the transform of that row.
"""
function findp(flat::Vector{T}, p::Vector{S}, layers::stParallel) where {T<:Number,S<:Integer}

end

"""
  p = computePath(layers::stParallel, layerM::Int, λ::Int; stType::String="full"))

compute the path p that corresponds to location λ in the stType of transform
"""
function computePath(layers::stParallel, layered1D::ScatteredFull{T,1}, layerM::Int, λ::Int; stType::String = "full") where {T<:Number}
  nScalesLayers = [numScales(layers.shears[i], layered1D.data[i]) for i = 1:layerM]
  p = zeros(layerM)
  if stType == "full"
    # start at the outermost level and work in
    for i = 1:layerM
      p[i] = div(tmpλ, prod(nScalesLayers[i+1:end])) + 1
      tmpλ -= (p[i] - 1) * prod(nScalesLayers[i+1:end])
    end
  elseif stType == "decreasing"
    @error "not implemented"
  end
end

@doc """
    numChildren = numChildren(keeper::Array{Int}, layers::stParallel, nScalesLayers::Array{Int})

given a keeper, determine the number of children it has.
"""
function numChildren(keeperOriginal::Array{Int}, scalingFactors::Array{Float64}, nScalesLayers::Array{Int})
  keeperOriginal = [keeperOriginal; 1]
  keeper = copy(keeperOriginal)
  m = length(keeper)
  if m == 0
    return 1 # not sure why you bothered
  elseif m == 1
    numThisLayer = 1
  else
    numThisLayer = 0
  end
  while true
    if keeperOriginal == [12, 12, 1]
    end
    # we're done if either if we've hit the maximum number of scales, or if at least one of the old indices is larger than it used to be. In the first case, count it, in the second case, don't
    if keeper[m] == nScalesLayers[m]
      return numThisLayer += 1
    elseif reduce((a, b) -> a | b, [keeper[i] > keeperOriginal[i] for i = 1:m-1])
      if keeperOriginal == [12, 12, 1]
      end
      return numThisLayer
    end
    numThisLayer += 1
    isLast, keeper = incrementKeeper(keeper, m, scalingFactors, nScalesLayers)
    if isLast
      break
    end
  end
  numThisLayer
end
numChildren(keeperOriginal::Array{Int}, layers::stParallel, nScalesLayers::Array{Int}) = numChildren(keeperOriginal::Array{Int}, [shear.scalingFactor for shear in layers.shears], nScalesLayers::Array{Int})

@doc """
numInLayer(m::Int, layers::stParallel, nScalesLayers::Array{Int})

A stupidly straightforward way of counting the number of decreasing paths in a given layer.
"""
function numInLayer(m::Int, scalingFactors::Array{Float64}, nScalesLayers::Array{Int})
  keeper = [1 for i = 1:m]
  if m == 0
    return 1 # not sure why you bothered
  else
    numThisLayer = 1
  end
  while true
    isLast, keeper = incrementKeeper(keeper, m, scalingFactors, nScalesLayers)
    if isLast
      break
    end
    numThisLayer += 1
  end
  numThisLayer
end

numInLayer(m::Int, layers::stParallel, nScalesLayers::Array{Int}) = numInLayer(m, [shear.scalingFactor for shear in layers.shears], nScalesLayers)


text:
# Test isbit-ness of VariableIndex and terms.
# It is important that these struct remain isbits as otherwise,
# the performance of function will be deteriored.
# These tests explicit this design choice and makes sure that it remains the case.
@testset "isbits" begin
    x = @inferred MOI.VariableIndex(1)
    @test isbits(x)
    at = @inferred MOI.ScalarAffineTerm(1.0, x)
    @test isbits(at)
    @test isbits(@inferred MOI.VectorAffineTerm(1, at))
    qt = @inferred MOI.ScalarQuadraticTerm(1.0, x, x)
    @test isbits(qt)
    @test isbits(@inferred MOI.VectorQuadraticTerm(1, qt))
end


text:
cd(@__DIR__)
D = KitBase.read_dict("config.txt")
for key in keys(D)
    s = key
    @eval $s = $(D[key])
end

#--- settings ---#
Setup(
    "gas",
    "cylinder",
    "2d2f",
    "kfvs",
    "bgk",
    1, # species
    2, # order of accuracy
    "vanleer", # limiter
    "maxwell",
    0.5, # cfl
    10.0, # time
) |> show
KitBase.Scalar(1.0, 1e-3)
KitBase.Radiation(1e-2, 1.0, 0.0, 1e-3, 1000)
gas = KitBase.Gas(knudsen, mach, prandtl, inK, 3.0, omega, alphaRef, omegaRef, 0.01)
show(gas)
KitBase.Gas(
    knudsen,
    mach,
    prandtl,
    inK,
    3.0,
    omega,
    alphaRef,
    omegaRef,
    0.01,
    1e-4,
    10000,
) |> show
KitBase.Mixture([0.1, 0.5], mach, prandtl, inK, 3.0, 1.0, 0.5, 0.5, 0.5) |> show
KitBase.Plasma1D(
    [0.1, 0.5],
    mach,
    prandtl,
    inK,
    3.0,
    1.0,
    0.5,
    0.5,
    0.5,
    0.01,
    0.01,
    100.0,
    1.0,
    1.0,
) |> show
KitBase.Plasma2D(
    [0.1, 0.5],
    mach,
    prandtl,
    inK,
    3.0,
    1.0,
    0.5,
    0.5,
    0.5,
    0.01,
    0.01,
    100.0,
    1.0,
    1.0,
) |> show
KitBase.DiatomicGas(
    knudsen,
    mach,
    prandtl,
    inK,
    inK,
    1.4,
    0.81,
    1.0,
    0.5,
    1e-3,
    89.1 / 273,
    18.1,
    1 / 1.55,
    0.2354,
    0.3049,
) |> show

begin
    prim = [1.0, 0.0, 1.0]
    w = KitBase.prim_conserve(prim, 1.4)
    u = Float64.(collect(umin:nu:umax)) # u should be float
    h = KitBase.maxwellian(u, prim)
    b = h .* inK ./ (2.0 * prim[end])
    x = Float64(x0) # x & dx should be of same type
    dx = x0 / nx
end

fw = (args...) -> [1.0, 0.0, rand()]
KitBase.IB(fw, gas) |> show

#--- control volume ---#
KitBase.ControlVolume1D(w, prim) |> show
KitBase.ControlVolume1D1F(w, prim, h) |> show
KitBase.ControlVolume1D2F(w, prim, h, b) |> show
KitBase.ControlVolume1D3F(
    hcat(w, w),
    hcat(prim, prim),
    zeros(nu, nu, 2),
    zeros(nu, nu, 2),
    zeros(nu, nu, 2),
    zeros(3),
    zeros(3),
    zeros(3, 2),
) |> show
KitBase.ControlVolume1D3F(
    zeros(5, 7, 2), # indexed with [flow entry, uq, species]
    zeros(5, 7, 2),
    zeros(nu, nu, 7, 2),
    zeros(nu, nu, 7, 2),
    zeros(nu, nu, 7, 2),
    zeros(3, 7),
    zeros(3, 7),
    zeros(3, 7, 2),
) |> show
# Rykov
KitBase.ControlVolume1D3F(rand(4), rand(5), rand(nu), rand(nu), rand(nu)) |> show
KitBase.ControlVolume1D4F(
    hcat(w, w),
    hcat(prim, prim),
    hcat(h, h),
    hcat(h, h),
    hcat(h, h),
    hcat(h, h),
    zeros(3),
    zeros(3),
    zeros(3, 2),
) |> show
KitBase.ControlVolume1D4F(
    zeros(5, 7, 2), # indexed with [flow entry, uq, species]
    zeros(5, 7, 2),
    zeros(nu, 7, 2),
    zeros(nu, 7, 2),
    zeros(nu, 7, 2),
    zeros(nu, 7, 2),
    zeros(3, 7),
    zeros(3, 7),
    zeros(3, 7, 2),
) |> show
KitBase.ControlVolume2D(w, prim) |> show
KitBase.ControlVolume2D1F(w, prim, h) |> show
KitBase.ControlVolume2D2F(w, prim, h, b) |> show
KitBase.ControlVolume2D3F(w, prim, h, b, b, zeros(3), zeros(3), zeros(3, 2)) |> show
# Rykov
KitBase.ControlVolume2D3F(rand(4), rand(4), rand(8, 8), rand(8, 8), rand(8, 8)) |> show
# unstructured
KitBase.ControlVolumeUS([1, 0], x, dx, w, prim)
KitBase.ControlVolumeUS1F([1, 0], x, dx, w, prim, h)
KitBase.ControlVolumeUS2F([1, 0], x, dx, w, prim, h, b)

#--- interface ---#
KitBase.Interface1D(w) |> show
KitBase.Interface1D1F(w, h) |> show
KitBase.Interface1D2F(w, h) |> show
KitBase.Interface1D3F(w, h, zeros(3)) |> show
KitBase.Interface1D3F(zeros(5, 7, 2), zeros(21, 21, 7, 2), zeros(3, 7)) |> show
KitBase.Interface1D3F(w, h) |> show # Rykov
KitBase.Interface1D4F(w, h, zeros(3)) |> show
KitBase.Interface1D4F(zeros(5, 7, 2), zeros(21, 7, 2), zeros(3, 7)) |> show
cosa = 1 / √2
sina = 1 / √2
KitBase.Interface2D(dx, cosa, sina, w) |> show
KitBase.Interface2D1F(dx, cosa, sina, w, h) |> show
KitBase.Interface2D2F(dx, cosa, sina, w, h) |> show
KitBase.Interface2D3F(dx, cosa, sina, w, h, zeros(3)) |> show
KitBase.Interface2D3F(dx, cosa, sina, w, h, zeros(3, 7)) |> show
KitBase.Interface2D3F(dx, cosa, sina, w, h) |> show

#--- solution ---#
sol_w = [w for i = 1:2]
sol_prim = [prim for i = 1:2]
sol_h = [h for i = 1:2]
KitBase.Solution1D(sol_w, sol_prim)
KitBase.Solution1D1F(sol_w, sol_prim, sol_h)
KitBase.Solution1D1F(sol_w, sol_prim, similar(sol_w), sol_h, similar(sol_h))
KitBase.Solution1D2F(sol_w, sol_prim, sol_h, sol_h)
KitBase.Solution1D2F(
    sol_w,
    sol_prim,
    similar(sol_w),
    sol_h,
    sol_h,
    similar(sol_h),
    similar(sol_h),
)

sol_w = [w for i = 1:2, j = 1:2]
sol_prim = [prim for i = 1:2, j = 1:2]
sol_h = [h for i = 1:2, j = 1:2]
KitBase.Solution2D(sol_w, sol_prim)
KitBase.Solution2D1F(sol_w, sol_prim, sol_h)
KitBase.Solution2D1F(sol_w, sol_prim, similar(sol_w), sol_h, similar(sol_h))
KitBase.Solution2D2F(sol_w, sol_prim, sol_h, sol_h)
KitBase.Solution2D2F(
    sol_w,
    sol_prim,
    similar(sol_w),
    sol_h,
    sol_h,
    similar(sol_h),
    similar(sol_h),
)

#--- flux ---#
KitBase.Flux1D(w, w)
KitBase.Flux1D1F(w, w, h)
KitBase.Flux1D2F(w, w, h, h)
KitBase.Flux2D(zeros(2), w, w, zeros(2), w, w)
KitBase.Flux2D1F(zeros(2), w, w, h, zeros(2), w, w, h)
KitBase.Flux2D2F(zeros(2), w, w, h, h, zeros(2), w, w, h, h)


text:
"""
    lsnormal(A,b)

Solve a linear least squares problem by the normal equations. Returns the
minimizer of ||b-Ax||.
"""
function lsnormal(A,b)

N = A'*A;  z = A'*b;
R = cholesky(N).U
w = forwardsub(R',z)                   # solve R'z=c
x = backsub(R,w)                       # solve Rx=z

return x
end

"""
    lsqrfact(A,b)

Solve a linear least squares problem by QR factorization. Returns the
minimizer of ||b-Ax||.
"""
function lsqrfact(A,b)

Q,R = qr(A)
c = Q'*b
x = backsub(R,c)

return x
end

max_stars_repo_name: snowdj/fnc-extras

max_issues_repo_name: snowdj/fnc-extras

max_forks_repo_name: snowdj/fnc-extras


text:
const Match = Tuple{Rule, Sub, Int64}
const MatchesBuf = Vector{Match}

import ..genrhsfun
import ..options
import ..@log

using .Schedulers

function inst(pat, sub::Sub)
    # TODO interface istree (?)
    if haskey(sub, pat)
        (eclass, lit) = sub[pat]
        lit != nothing ? lit : eclass
    else
        pat
    end
end

function instantiate(p, sub::Sub; skip_assert=false)
    # remove type assertions
    if skip_assert
        p = df_walk( x -> (isexpr(x, :(::)) ? x.args[1] : x), p; skip_call=true )
    end

    df_walk(inst, p, sub; skip_call=true)
end

"""
Core algorithm of the library: the equality saturation step.
"""
function eqsat_step!(egraph::EGraph, theory::Vector{Rule};
        scheduler=SimpleScheduler(), mod=@__MODULE__,
        match_hist=MatchesBuf(), sizeout=0, stopwhen=()->false,
        matchlimit=5000 # max number of matches
        )
    matches=MatchesBuf()
    EMPTY_DICT = Sub()

    readstep!(scheduler)

    report = Report(egraph)

    search_stats = @timed for rule ∈ theory
        # don't apply banned rules
        shouldskip(scheduler, rule) && continue

        if rule.mode ∉ [:rewrite, :dynamic, :equational]
            error("unsupported mode in rule ", rule)
        end

        # outermost symbol in lhs
        sym = getfunsym(rule.left)
        if istree(rule.left)
            ids = get(egraph.symcache, sym, [])
        else
            ids = keys(egraph.M)
        end

        for id ∈ ids
            for sub in ematch(egraph, rule.left, id, EMPTY_DICT)
                # display(sub); println()
                !isempty(sub) && push!(matches, (rule, sub, id))
            end
        end

        if rule.mode == :equational
            sym = getfunsym(rule.right)
            ids = get(egraph.symcache, sym, [])
            for id ∈ ids
                for sub in ematch(egraph, rule.right, id, EMPTY_DICT)
                    # display(sub); println()
                    !isempty(sub) && push!(matches, (rule, sub, id))
                end
            end
        end
    end
    report.search_stats = report.search_stats + discard_value(search_stats)

    # mmm = unique(matches)
    # mmm = symdiff(match_hist, matches)

    mmm = setdiff(matches, match_hist)

    if length(mmm) > matchlimit
        mmm = mmm[1:matchlimit]
        # mmm = Set(collect(mmm)[1:matchlimit])
        #
        # report.reason = :matchlimit
        # @goto quit_rebuild
    end

    # println("============ WRITE PHASE ============")
    # println("\n $(length(matches)) $(length(match_hist))")
    # println(" diff length $(length(mmm))")

    i = 0

    apply_stats = @timed for match ∈ mmm
        i += 1
        (rule, sub, id) = match

        if i % 300 == 0
            if sizeout > 0 && length(egraph.U) > sizeout
                @log "E-GRAPH SIZEOUT"
                report.reason = :sizeout
                @goto quit_rebuild
            end

            if stopwhen()
                @log "Halting requirement satisfied"
                report.reason = :condition
                @goto quit_rebuild
            end
        end

        writestep!(scheduler, rule)

        if rule.mode == :rewrite || rule.mode == :equational # symbolic replacement
            l = instantiate(rule.left, sub; skip_assert=true)
            r = instantiate(rule.right, sub)
            lc = addexpr!(egraph, l)
            rc = addexpr!(egraph, r)
            merge!(egraph, lc.id, rc.id)
        elseif rule.mode == :dynamic # execute the right hand!
            l = instantiate(rule.left,sub; skip_assert=true)
            lc = addexpr!(egraph, l)

            (params, f) = rule.right_fun[mod]
            actual_params = map(params) do x
                (eclass, literal) = sub[x]
                literal != nothing ? literal : eclass
            end
            new = f(egraph, actual_params...)
            rc = addexpr!(egraph, new)
            merge!(egraph,lc.id,rc.id)
        else
            error("unsupported rule mode")
        end
    end
    report.apply_stats = report.apply_stats + discard_value(apply_stats)

    union!(match_hist, mmm)

    # display(egraph.parents); println()
    # display(egraph.M); println()
    @label quit_rebuild
    report.saturated = isempty(egraph.dirty) && report.reason == nothing
    rebuild_stats = @timed rebuild!(egraph)
    report.rebuild_stats = report.rebuild_stats + discard_value(rebuild_stats)

    total_time!(report)

    return report, egraph
end

# TODO plot how egraph shrinks and grows during saturation
"""
Given an [`EGraph`](@ref) and a collection of rewrite rules,
execute the equality saturation algorithm.
"""
function saturate!(egraph::EGraph, theory::Vector{Rule};
    mod=@__MODULE__,
    timeout=0, stopwhen=(()->false), sizeout=2^12,
    matchlimit=5000,
    scheduler::Type{<:AbstractScheduler}=BackoffScheduler)

    if timeout == 0
        timeout = length(theory)
    end

    curr_iter = 0

    # prepare the dynamic rules in this module
    for rule ∈ theory
        if rule.mode == :dynamic && !haskey(rule.right_fun, mod)
            rule.right_fun[mod] = genrhsfun(rule.left, rule.right, mod)
        end
    end

    # init the scheduler
    sched = scheduler(egraph, theory)

    match_hist = MatchesBuf()

    tot_report = Report()

    while true
        curr_iter+=1

        options[:printiter] && @info("iteration ", curr_iter)

        report, egraph = eqsat_step!(egraph, theory;
            scheduler=sched, mod=mod,
            match_hist=match_hist, sizeout=sizeout,
            matchlimit=matchlimit,
            stopwhen=stopwhen)

        tot_report = tot_report + report

        report.reason == :matchlimit && break
        cansaturate(sched) && report.saturated && (tot_report.saturated = true; break)
        curr_iter >= timeout && (tot_report.reason = :timeout; break)
        sizeout > 0 && length(egraph.U) > sizeout && (tot_report.reason = :sizeout; break)
        stopwhen() && (tot_report.reason = :condition; break)
    end
    # println(match_hist)
    tot_report.iterations = curr_iter
    @log tot_report

    return tot_report
end


text:
using Revise
using Makie, ZigZagBoomerang, SparseArrays, LinearAlgebra

function ϕ(x, i, μ)
    x[i] - μ[i]
end
κ = 1.0
n = 100
μ = zeros(n)
x0 = randn(n)
θ0 = rand((-1.0,1.0), n)
T = 1000.0
c = 10ones(n)

#@time trace0, _ = ZigZagBoomerang.spdmp(ϕ, 0.0, x0, θ0, T, c, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)

@time trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, κ, ZigZag(sparse(1.0I,n,n), zeros(n)), μ)
ts0, xs0 = splitpairs(trace0)

#lines(ts0, getindex.(xs0,1))
#ts0, xs0 = splitpairs(discretize(trace0, 0.01))
#p1 = scatter(getindex.(xs0,1), getindex.(xs0,2), color=(:black, 0.4), markersize=0.05)
#save(joinpath("figures", "spikeandslab.png"), title(p1, "Spike and Slab"))
function prob0(ts0, xs0)
    p = 0.0
    t⁻ = 0.0
    first = true
    for i in eachindex(ts0)
        if xs0[i] == 0.0
            if first
                t⁻ = ts0[i]
                first = false
                continue
            else
                continue
            end
        elseif first == false && xs0[i] != 0.0
            p += ts0[i-1] - t⁻
            first = true
        else
            continue
        end
    end
    p/ ts0[end]
end
prob0(ts0, getindex.(xs0,1))

tp(κ, μ) = 1/(1 + κ*sqrt(2pi)*exp(0.5*(-μ)^2))

using Random
Random.seed!(0)
μ1 = 0.0:0.5:2.0
κ = 0.1:0.1:1.0
A = zeros(length(μ1), length(κ))
B =  zeros(length(μ1), length(κ))
for i in eachindex(μ1)
    for j in eachindex(κ)
        global κ
        global μ1
        μ = fill(μ1[i], n)
        trace0, _ = ZigZagBoomerang.sspdmp(ϕ, 0.0, x0, θ0, T, c, κ[j], ZigZag(sparse(1.0I,n,n), zeros(n)), μ)
        ts0, xs0 = splitpairs(trace0)
        A[i,j] = sum([prob0(ts0, getindex.(xs0, k)) for k in 1:100])/100
    end
end

using CairoMakie, ColorSchemes
mygradcol = ColorSchemes.deep
length(mygradcol)
p1 = Scene()
for i in eachindex(μ1)
    global μ1, κ
    p1 = plot!(κ, A[i,:], color = mygradcol[Int(floor(i*256/5))], label = "mu = 0.0")
    Makie.lines!(κ, tp.(κ, μ1[i]), color = mygradcol[Int(floor(i*256/5))])
end
save("figures/zz_test.png", title(p1, "sticky ZZ"))


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.test.behaviour.concept.thing.entity;
#
# import typedb.client.api.concept.thing.Attribute;
# import typedb.client.api.concept.thing.Entity;
# import typedb.client.common.Label;
# import io.cucumber.java.en.Then;
# import io.cucumber.java.en.When;
#
# import java.time.LocalDateTime;
#
# import static typedb.client.test.behaviour.concept.thing.ThingSteps.get;
# import static typedb.client.test.behaviour.concept.thing.ThingSteps.put;
# import static typedb.client.test.behaviour.connection.ConnectionStepsBase.tx;
# import static typedb.client.test.behaviour.util.Util.assertThrows;
# import static org.junit.Assert.assertEquals;
# import static org.junit.Assert.assertTrue;
#
# @SuppressWarnings("CheckReturnValue")
# public class EntitySteps {
#
#     @When("{var} = entity\\( ?{type_label} ?) create new instance")
#     public void entity_type_create_new_instance(String var, String typeLabel) {
#         put(var, tx().concepts().getEntityType(typeLabel).asRemote(tx()).create());
#     }
#
#     @When("entity\\( ?{type_label} ?) create new instance; throws exception")
#     public void entity_type_create_new_instance_throws_exception(String typeLabel) {
#         assertThrows(() -> tx().concepts().getEntityType(typeLabel).asRemote(tx()).create());
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {int}")
#     public void entity_type_create_new_instance_with_key(String var, String type, String keyType, int keyValue) {
#         Attribute.Long key = tx().concepts().getAttributeType(keyType).asLong().asRemote(tx()).put(keyValue);
#         Entity entity = tx().concepts().getEntityType(type).asRemote(tx()).create();
#         entity.asRemote(tx()).setHas(key);
#         put(var, entity);
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {word}")
#     public void entity_type_create_new_instance_with_key(String var, String type, String keyType, String keyValue) {
#         Attribute.String key = tx().concepts().getAttributeType(keyType).asString().asRemote(tx()).put(keyValue);
#         Entity entity = tx().concepts().getEntityType(type).asRemote(tx()).create();
#         entity.asRemote(tx()).setHas(key);
#         put(var, entity);
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) create new instance with key\\( ?{type_label} ?): {datetime}")
#     public void entity_type_create_new_instance_with_key(String var, String type, String keyType, LocalDateTime keyValue) {
#         Attribute.DateTime key = tx().concepts().getAttributeType(keyType).asDateTime().asRemote(tx()).put(keyValue);
#         Entity entity = tx().concepts().getEntityType(type).asRemote(tx()).create();
#         entity.asRemote(tx()).setHas(key);
#         put(var, entity);
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {long}")
#     public void entity_type_get_instance_with_key(String var1, String type, String keyType, long keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asLong().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {word}")
#     public void entity_type_get_instance_with_key(String var1, String type, String keyType, String keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asString().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @When("{var} = entity\\( ?{type_label} ?) get instance with key\\( ?{type_label} ?): {datetime}")
#     public void entity_type_get_instance_with_key(String var1, String type, String keyType, LocalDateTime keyValue) {
#         put(var1, tx().concepts().getAttributeType(keyType).asDateTime().asRemote(tx()).get(keyValue).asRemote(tx()).getOwners()
#                 .filter(owner -> owner.getType().getLabel().equals(Label.of(type))).findFirst().orElse(null));
#     }
#
#     @Then("entity\\( ?{type_label} ?) get instances contain: {var}")
#     public void entity_type_get_instances_contain(String typeLabel, String var) {
#         assertTrue(tx().concepts().getEntityType(typeLabel).asRemote(tx()).getInstances().anyMatch(i -> i.equals(get(var))));
#     }
#
#     @Then("entity\\( ?{type_label} ?) get instances is empty")
#     public void entity_type_get_instances_is_empty(String typeLabel) {
#         assertEquals(0, tx().concepts().getEntityType(typeLabel).asRemote(tx()).getInstances().count());
#     }
# }


text:
"""
    create_ols_approximation(y::Vector{<:Real}, x::Vector{<:Real}, base_x::Real = 0.0, degree::Integer = 1, intercept::Bool = true)
    create_ols_approximation(y::Vector{<:Real}, x::Union{Vector{DateTime},Vector{Date},Vector{Union{Date,DateTime}}}, base_x::Union{Date,DateTime} = global_base_date, degree::Integer = 1, intercept::Bool = true)

An approximation function calculated via OLS.

### Inputs
* `y` - A `Vector` with the y coordinates
* `x` - A `Vector` with the x coordinates
* `base_x` - A real that offsets the x. So a coordinate with x value of 2.0 will be converted to 1.8 if base_x is 0.2.
* `degree` - What the highest power of x should be. So if this is 3 then the equation will have x, x^2, x^3 as predictors.
* `intercept` - Should there be an x intercept.
### Returns
* A `Sum_Of_Functions` containing the approximation function.
"""
function create_ols_approximation(y::Vector{<:Real}, x::Vector{<:Real}, base_x::Real = 0.0, degree::Integer = 1, intercept::Bool = true)
    obs = length(y)
    if degree < 0
        error("Cannot approximate with OLS with a degree that is negative")
    end
    x = x .- base_x
    if intercept
        X = ones(obs)
        for i in 1:degree
            X = hcat(X, (x .^ i))
        end
    else
        X = x
        for i in 2:degree
            X = hcat(X, (x .^ i))
        end
    end

    lm1 = fit(LinearModel,  hcat(X), y)
    beta = lm1.pp.beta0
    func_array = Array{PE_Function,1}(undef,convert(Int, intercept) + degree)
    if intercept
        func_array[1] = PE_Function(beta[1], 0.0, base_x, 0)
    end
    for d in 1:degree
        func_array[d+Integer(intercept)] = PE_Function(beta[d+Integer(intercept)], 0.0, base_x, d)
    end
    return Sum_Of_Functions(func_array)
end

function create_ols_approximation(y::Vector{<:Real}, x::Union{Vector{DateTime},Vector{Date},Vector{Union{Date,DateTime}}}, base_x::Union{Date,DateTime} = global_base_date, degree::Integer = 1, intercept::Bool = true)
    base = years_from_global_base.(base_x)
    xx   = years_from_global_base.(x)
    return create_ols_approximation(y, xx, base, degree, intercept)
end

function get_cholesky_coefficients(chebyshev::Sum_Of_Functions, y::Vector{<:Real}, unnormalised_nodes::Vector{<:Real})
    chebyshev_on_nodes = evaluate.(Ref(chebyshev), unnormalised_nodes)
    a = sum(y .* chebyshev_on_nodes) / sum(chebyshev_on_nodes .^ 2)
    return a
end

"""
    create_chebyshev_approximation(func::Function, nodes::Integer, degree::Integer, left::Real, right::Real)

An function that will approximate another function via Chebyshev polynomials.

### Inputs
* `func` - A function that you want to approximation
* `nodes` - The number of approximation nodes
* `degree` - The degree of the Chebyshev polynomials.
* `left` - The left limit of the approximation
* `right` - The right limit of the approximation.
### Returns
* A `Sum_Of_Functions` containing the approximation function.
"""
function create_chebyshev_approximation(func::Function, nodes::Integer, degree::Integer, left::Real, right::Real)
    # This is all after Algorithm 6.2 from Judd (1998) Numerical Methods in Economics.
    if nodes <= degree
        error("Need to have more nodes than degree to use a chebyshev approximation")
    end
    k = 1:nodes
    unnormalised_nodes = -cos.( (((2 .* k) .- 1) ./ (2 * nodes)) .* pi    )
    normalised_nodes = ((unnormalised_nodes .+ 1) .* ((right-left)/2)) .+ left
    y = func.(normalised_nodes)
    chebyshevs = get_chevyshevs_up_to(degree, true)
    a = get_cholesky_coefficients.(chebyshevs, Ref(y), Ref(unnormalised_nodes))
    transformed_chebyshevs = convert_to_linearly_rescale_inputs.(chebyshevs, (right-left)/2, +(right+left)/2)
    # Note that these alpha and beta parameters in the convert_to_linearly_rescale_inputs function differ from those in Judd because those did not work in this context.
    all_terms = a .* transformed_chebyshevs
    final_func = Sum_Of_Functions(all_terms)
    return final_func
end


text:
export fieldterm, xxzmodel, hopping, z_field

_get_interaction(J::AbstractMatrix, i, j) = J[i,j]
_get_interaction(J::Number, i, j) = J

function fieldterm(N, op)
    res = op
    for i in 2:N
        res = 𝟙 ⊗ res + op ⊗ identity_op(i-1)
    end
    res
end

function hopping(N, i, j)
    entries = 2^(N-1)
    A,B,V = zeros(Int, entries), zeros(Int, entries), ones(Int, entries)
    _hopping_coordinates!(A, B, N, i, j)
    return sparse(A,B,V,2^N,2^N)
end

xxzmodel(J::AbstractMatrix, Δ) = xxzmodel(size(J, 1), J, Δ)

function xxzmodel(N,J,Δ)
    # 2^N diagonal entries + N(N-1)/2 * 2^(N-1) offdiagonal terms
    entries = 2^N + 2^(N-2)*N*(N-1)
    A,B,V = zeros(Int, entries), zeros(Int, entries), zeros(float(promote_type(eltype(J), typeof(Δ))), entries)

    @views _zcorrelator_coo!(A[1:2^N],B[1:2^N],V[1:2^N], N, J)
    V[1:2^N] .*= Δ

    @views _hopping!(A[2^N+1:end], B[2^N+1:end], V[2^N+1:end], N, J)
    V[2^N+1:end] .*= 2

    V ./= 4 # pauli -> spin operators
    sparse(A,B,V,2^N,2^N)
    ## In principle using Hermitian and only computing only
    ## half of the off-diagonal speeds up constructing the Hamiltonian.
    ## In practice this slows down everythin else dramatically, as Hermitian
    ## (sometimes) hides the fact that this is a sparse matrix and that leads to
    ## unfortunate interactions (f.e. Hermitian(sparse) ≈ dense takes AGES)
    #Hermitian(sparse(A,B,V,2^N,2^N), :L)
end

"""
    _zcorrelator_coo!(A, B, V, J)

Generate the COO vectors of ∑ᵢⱼJᵢⱼ σzⁱσzʲ.

# Parameters
 - `A`,`B`,`V` vectors for column (`A`) and row indices (`B`) respectively values (`V`)
 - `N` number of sites in the system
 - `J` interaction matrix of size `N`x`N` or scalar. Assumed symmetric!

`A`,`B` and `V` need to have length 2^N.
"""
function _zcorrelator_coo!(A,B,V, N, J)
    A[:] .= 1:2^N
    B[:] .= 1:2^N
    V[:] .= 0
    for i in 1:N
        for j in i+1:N
            _zcorrelator_values!(V, N, _get_interaction(J, j,i), i, j)
        end
    end
end

"""
    _zcorrelator_values!(V, N, Jij, i, j)

*ADD* to `V` the diagonal of Jᵢⱼ σzⁱσzʲ.

# Parameters
 - `V` array of length 2^`N`. Will be mutated!
 - `N` number of sites in the system
 - `Jij` interaction strength
 - `i`,`j` sites to couple
"""
function _zcorrelator_values!(V, N, Jij, i, j)
    # 𝟙(sec1) ⊗ σz ⊗ 𝟙(sec2) ⊗ σz ⊗ 𝟙(sec3)

    #             | 𝟙(sec2) |     0    |   | 𝟙(sec3) |    0     |
    # = 𝟙(sec1) ⊗ | ------- | -------- | ⊗ | ------- | -------- |
    #             |    0    | -𝟙(sec2) |   |    0    | -𝟙(sec3) |

    #             |           [ 𝟙(sec3) |     0    ] |                                  |
    #             | 𝟙(sec2) ⊗ [ ------- | -------- ] |              0                   |
    #             |           [    0    | -𝟙(sec3) ] |                                  |
    # = 𝟙(sec1) ⊗ | -------------------------------- | -------------------------------- |
    #             |                                  |           [ -𝟙(sec3) |    0    ] |
    #             |              0                   | 𝟙(sec2) ⊗ [ -------- | ------- ] |
    #             |                                  |           [    0     | 𝟙(sec3) ] |
    sec1 = 2^(i-1)
    sec2 = 2^(j-i-1)
    sec3 = 2^(N-j)

    blocksize3 = 2sec3 # σz ⊗ 𝟙(sec3)
    blocksize23 = 2sec2*blocksize3 # size of: σz ⊗ 𝟙(sec2) ⊗ σz ⊗ 𝟙(sec3)

    #sec1 loop
    for offset1 in blocksize23 .* (0:sec1-1)
        # upper block
        for offset2 in offset1 .+ blocksize3 .* (0:sec2-1)
            V[1 .+ offset2 .+ (0:sec3-1)] .+= Jij
            V[1 .+ offset2 .+ (sec3:2sec3-1)] .+= -Jij
        end
        # lower block
        for offset2 in offset1 .+ blocksize3 .* (sec2:2sec2-1)
            V[1 .+ offset2 .+ (0:sec3-1)] .+= -Jij
            V[1 .+ offset2 .+ (sec3:2sec3-1)] .+= Jij
        end
    end
end

"""
    _hopping!(A,B,V, J)

Generate the COO vectors for ∑ᵢⱼ Jᵢⱼ (σ₊ⁱσ₋ʲ + σ₋ⁱσ₊ʲ)

# Parameters
 - `A`,`B`,`V` vectors for column (`A`) and row indices (`B`) respectively values (`V`)
 - `J` interaction matrix of size `N`x`N`. Assumed symmetric!

 If `J` is `N`x`N`, then `A`,`B` and `V` need to have length (N-1)N*2^(N-2)
"""
function _hopping!(A,B,V, N, J)
    fill!(V, 1)
    length = 2^(N-1)
    at = 0
    for i in 1:N
        for j in i+1:N
            range = at .+ (1:length)
            @views _hopping_coordinates!(A[range], B[range], N, i, j)
            V[range] .*= _get_interaction(J, j, i)
            at += length
        end
    end
    ## _hopping values! does the transpose
    ## no measureable performance difference
    #copyto!(A,at+1,B,1,at)
    #copyto!(B,at+1,A,1,at)
    #copyto!(V,at+1,V,1,at)
end

"""
    _hopping_coordinates!(A, B, N, i, j)

Compute the coordinates of the matrix entries of the hopping operator between site i and J
in a spin chain. The values are 1 everywhere, so there is no input/output for these in this
function.

Hopping operator: σ₊ⁱσ₋ʲ + σ₋ⁱσ₊ʲ = 𝟙⊗…⊗𝟙⊗σ₊⊗𝟙⊗…⊗𝟙⊗σ₋⊗𝟙⊗…⊗𝟙 + transpose

# Parameters
 - `A`,`B` coordinate vectors. These should have length 2^(N-2) will be overwritten!
 - `N` number of sites in the system
 - `i`,`j` sites the excitation can hop between. 1 ≤ i,j ≤ N
"""
function _hopping_coordinates!(A,B,N,i,j)
    i > j && return _hopping_coordinates!(A,B,N,j,i)
    # 𝟙(sec1) ⊗ σ+ ⊗ 𝟙(sec2) ⊗ σ- ⊗ 𝟙(sec3) + transpose
    #             |                                  |                                  |
    #             |                0                 |                 0                |
    #             |                                  |                                  |
    # = 𝟙(sec1) ⊗ | -------------------------------- | -------------------------------- |
    #             |           [    0    |  𝟙(sec3) ] |                                  |
    #             | 𝟙(sec2) ⊗ [ ------- | -------- ] |                 0                |
    #             |           [    0    |    0     ] |                                  |
    sec1 = 2^(i-1)
    sec2 = 2^(j-i-1)
    sec3 = 2^(N-j)

    blocksize3 = 2sec3 # σz ⊗ 𝟙(sec3)
    blocksize23 = 2sec2*blocksize3 # size of: σz ⊗ 𝟙(sec2) ⊗ σz ⊗ 𝟙(sec3)

    currentIndex = 0
    #sec1 loop
    for offset1 in blocksize23 .* (0:sec1-1)
        for offset2 in blocksize3 .* (0:sec2-1)
            rowOffset2 = offset1 + offset2 # left columns # offset for the row-COORDINATE
            colOffset2 = offset1 + offset2 + blocksize3*sec2 # bottom rows # offset for the col-COORDINATE

            A[currentIndex .+ (1:sec3)] .= colOffset2 .+ (1:sec3) # upper rows
            B[currentIndex .+ (1:sec3)] .= rowOffset2 .+ (sec3+1:2sec3) # right columns
            currentIndex += sec3

            ## now transpose
            B[currentIndex .+ (1:sec3)] .= colOffset2 .+ (1:sec3) # upper rows
            A[currentIndex .+ (1:sec3)] .= rowOffset2 .+ (sec3+1:2sec3) # right columns
            currentIndex += sec3
        end
    end
end

z_field(N::Int) = z_field(ones(N))

function z_field(field_values::Vector)
    N = length(field_values)
    ret = zeros(2^N)
    for (i,hz) in enumerate(field_values)
        _z_field_values!(ret, N, hz, i)
    end
    return Diagonal(ret)
end

"""
    _z_field!(V, N, hz, i)

*ADD* to `V` the diagonal of hz σzⁱ.

# Parameters
 - `V` array of length 2^`N`. Will be mutated!
 - `N` number of sites in the system
 - `hz` field strength
 - `i` site to act on
"""
function _z_field_values!(V, N, hz, i)
    # σz^(i) = 𝟙(sec1) ⊗ σz ⊗  𝟙(sec2)
    sec1 = 2^(i-1)
    sec2 = 2^(N-i)

    blocksize2 = 2sec2 # σz ⊗ 𝟙(sec2)

    for offset1 in blocksize2 .* (0:sec1-1)
        V[offset1 .+ (1:sec2)] .+= hz
        V[offset1 .+ (sec2+1:2sec2)] .-= hz
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule PicoSAT_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("PicoSAT")
JLLWrappers.@generate_main_file("PicoSAT", UUID("e78fa76d-a187-569f-aede-ad11521a2edf"))
end  # module PicoSAT_jll


text:
"""
Returns the boundingbox of a patch in terms of its center and halfsize.

    function boundingbox{U,D,C,N,T}(p::Simplex{U,D,C,N,T}) -> center, halfsize
"""
function boundingbox(p::Simplex{U,D,C,N,T}) where {U,D,C,N,T}

    # ll = minimum(p.vertices)
    # ur = maximum(p.vertices)

    ll = first(p.vertices); for v in p.vertices; ll = min.(v, ll); end
    ur = first(p.vertices); for v in p.vertices; ll = max.(v, ll); end

    center = (ll + ur) / 2
    halfsize = maximum(ur - center)

    return center, halfsize
end

function boundingbox(V)

    # ll = minimum(V)
    # ur = maximum(V)

    ll = first(V); for v ∈ V; ll = min.(v, ll); end
    ur = first(V); for v ∈ V; ur = max.(v, ur); end

    c = (ll + ur) / 2
    s = maximum(ur - c)

    return c, s
end

boundingbox(p::SVector{N,T}) where {N,T<:Number} = p, zero(eltype(p))


text:
module brainvisionloader

struct dataset
    data::Array{Float64}
    labels::Array{String}
    markers
    srate::Float64
end

function loadFile(filename::String)
    dataInfo, chanInfo  = readHDR(filename)
    trial               = readEEG(filename,dataInfo,chanInfo)
    markers             = readMRK(filename)
    EEG                 = dataset(trial,chanInfo[1,:],markers,get(dataInfo,"srate",NaN));
    return EEG
end

function readEEG(filename::String,dataInfo,chanInfo)
    open(string(filename,".eeg")) do f
        filesize    = position(seekend(f))
        framenumber = Int64(filesize/( get(dataInfo,"ByteLength",2) * get(dataInfo,"ChannelNumber",32)))
        filedim     = (get(dataInfo,"channelNumber",32),framenumber)
        seekstart(f)
        trial = reinterpret(Int16,read(f))
    end

    trial = collect(trial)

    trial = reshape(trial,size(chanInfo,2),:)
    trial = trial.*repeat(chanInfo[2,:],1,size(trial,2))

    return trial
end

function readMRK(filename::String)
    doCollect   = false
    markerInfo  = Array{Any,2}(undef,0,4)
    open(string(filename,".vmrk")) do f
        while !eof(f)
            x = readline(f)
            # Read Channel-Specific Info
            if length(x)==0 || isnothing(x) || any(occursin("\r\n",x[1:2]))
                doCollect = false;
            end
            if doCollect
                markerInfo  = vcat(markerInfo,getMarkerSegment(x))
            end
            if occursin("; Commas in type or description text are coded",x)
                # once we are here, start collecting the markers
                doCollect = true
            end
        end
    end
    return markerInfo
end

function getMarkerSegment(x)
    type,code,position,duration,chan = split(split(x,"=")[2],",")
    position   = parse(Int64,position)
    return Mrk  = hcat(type,code,position,duration)
end

function readHDR(filename::String)
    Fs            = NaN;
    dataformat    = "";
    bitformat     = NaN;
    dataOrient    = "";
    binaryformat  = "";
    NumChan       = NaN;
    chanInfo      = [];
    dataInfo      = [];
    doCollect     = false;
    open(string(filename,".vhdr")) do f
        while !eof(f)
            x = readline(f)

          # Read General Info
            if occursin("Sampling Rate [Hz]:",x);
                Fs = parse(Int32,getGeneralInfo(x));
            end
            if occursin("Number of channels:",x);
                NumChan = parse(Int32,getGeneralInfo(x));
            end
            if occursin("Dataformat",x)
                dataformat = getdataInfo(x);
            end
            if occursin("Dataorientation",x)
                dataOrient = getdataInfo(x);
            end
            if occursin("Binaryformat",x);
                binaryformat,bitformat = getencoding(x)
            end

            # Read Channel-Specific Info

            if length(x)==0 || isnothing(x) || any(occursin("\r\n",x[1:2]))
                 doCollect = false;
            end
            if doCollect; chanInfo = hcat(chanInfo,getchanInfo(x)); end
            if occursin("Ch1=",x); doCollect = true; chanInfo = getchanInfo(x); end

        end
    end
    dataInfo = Dict("srate"=>Fs,"channelNumber"=>NumChan,"encoding"=>binaryformat,"byteLength"=>bitformat,"format"=>dataformat,"orientation"=>dataOrient)
    return dataInfo, chanInfo
end

function getchanInfo(x::String)
    #match(r"^Ch([0-9]+)=[A-Za-z0-9]+,[A-Za-z0-9]*,)
    id,other = split(x,"=")

    name,ref,res,unit  = split(other,",")

    res       = parse(Float64,res);
    info      = [name;res;id]

    return info
end

function getGeneralInfo(x::String)

    return split(x,':')[2];

end

function getdataInfo(x::String)
    return split(x,'=')[2];

end

function getencoding(x::String)
    TypeEncode  = Dict("INT_16" => "Int16","UINT_16" => "UInt16","FLOAT_32" => "Float32");
    BitEncode    = Dict("INT_16" => 2,"UINT_16" => 2,"FLOAT_32" => 4);
    return TypeEncode[getdataInfo(x)], BitEncode[getdataInfo(x)];
end

end


text:
using PhotoOrganizer

dry_run = false
rm_src = false
dst_root="/mnt/zpool/tank/media/Pictures"

src_dirs = [ "/mnt/zpool/tank/syncthing/Pixel4a/Glen_DCIM" ]
organize_photos(src_dirs, dst_root; rm_src, dry_run, photo_suffix="_Glen")

# Move photos off of IncomingPhotos drive to the photo library
rm_src = true
organize_photos(["/mnt/zpool/tank/IncomingPhotos"], dst_root; rm_src, dry_run)


text:
using Test
using PlayingCards
using TexasHoldem
TH = TexasHoldem

@testset "Players" begin
    players = (
        Player(Human(), 1),
        Player(Bot5050(), 2),
    )
    @test TH.name(players[1]) == "Human[1]"
    @test TH.name(players[2]) == "Bot5050[2]"
    @test TH.cards(players[1]) == players[1].cards
    @test TH.cards(players[2]) == players[2].cards
end


text:
#Demo 7: Bode Plots
#-------------------------------------------------------------------------------
using InspectDR
using Colors
import Printf: @sprintf

#==Input
===============================================================================#

#Constants
#-------------------------------------------------------------------------------
black = RGB24(0, 0, 0)
white = RGB24(1, 1, 1)
red = RGB24(1, 0, 0)
green = RGB24(0, 1, 0)
blue = RGB24(0, 0, 1)
j = im

#Fonts, Line types, etc:
default_line = line(color=blue, width=3)
markerline = line(style=:dash, width=2.5)
markerline_light = line(style=:dash, width=2.5, color=RGB24(.4,.4,.4))

#Input data
#-------------------------------------------------------------------------------
fmax = 10e9
fstep = 5e6
f = collect(0:fstep:fmax)
fmin_disp = 10e6
G = 20.0
f_p1 = 50e6
f_p2 = .9e9
f_z = Inf #.8e9
f1 = .8e9

#==Helper functions
===============================================================================#
phase(x) = rad2deg.(angle.(x))
function roots(a,b,c)
	r = sqrt(b*b-4*a*c)
	result = (-b .+ [r, -r]) / (2a)
	return (result[1], result[2])
end
xf_o2(f, G, f_z, f_p1, f_p2) = #2nd order transfer function
	G*(1 .+ j*(f./f_z))./( (1 .+ j*(f./f_p1)).*(1 .+ j*(f./f_p2)) )
function f_unitygain_f(f, G, f_z, f_p1, f_p2) #Compute unity gain frequency
	invsq(x) = (result = 1/x; return result*result)
	_z² = invsq(f_z); G² = G*G
	_p1² = invsq(f_p1); _p2² = invsq(f_p2)
	r = roots(_p1²*_p2², _p1²+_p2²-G²*_z², 1-G²)
	return sqrt(max(r...)) #Select sqrt(positive root)
end
f_3dB(f, G, f_z, f_p1, f_p2) = #Compute 3dB frequency
	f_unitygain_f(f, sqrt(2), f_z, f_p1, f_p2)

#==Calculations
===============================================================================#
X = xf_o2(f, G, f_z, f_p1, f_p2)
f0 = f_unitygain_f(f, G, f_z, f_p1, f_p2)
fBW = f_3dB(f, G, f_z, f_p1, f_p2)
phase0 = phase(xf_o2(f0, G, f_z, f_p1, f_p2))
pmargin = 180+phase0
@show pmargin

#==Generate plot
===============================================================================#
mplot = InspectDR.Multiplot(title="Sample Bode Plot")
#Bode plot looks better with wider aspect ratio:
mplot.layout[:valloc_plot] = mplot.layout[:valloc_plot]*.6

plot = add(mplot, InspectDR.bodeplot())
	strip_mag, strip_phase = plot.strips

	#Define extents & scales:
	plot.xext_full = InspectDR.PExtents1D(min=fmin_disp) #Avoid issues with log scale
	strip_mag.yext_full = InspectDR.PExtents1D(min=-10)
	strip_phase.yext_full = InspectDR.PExtents1D(min=-180)

	#Add waveforms:
	wfrm = add(plot, f, X, strip=1)
		wfrm.line = default_line
	wfrm = add(plot, f, phase(X), strip=2)
		wfrm.line = default_line

#Add specialized annotation
#-------------------------------------------------------------------------------

#Add vertical markers to both plots:
add(plot, vmarker(f0, markerline, strip=0))
add(plot, vmarker(fBW, markerline, strip=0))

#Add annotation to Magnitude plot:
afont = plot.layout[:font_annotation]
fstr = "f3dB=" * @sprintf("%.1f MHz", fBW/1e6)
	add(plot, atext(fstr, x=fBW, xoffset=-3, yoffset_rel=.5, font=afont, angle=-90, align=:bc, strip=1))
fstr = "f0=" * @sprintf("%.1f MHz", f0/1e6)
	add(plot, atext(fstr, x=f0, xoffset=-3, yoffset_rel=.5, font=afont, angle=-90, align=:bc, strip=1))
add(plot, atext("0dB", y=0, xoffset_rel=0.5, yoffset=3, font=afont, align=:bc, strip=1))
add(plot, hmarker(0, markerline_light, strip=1))

#Add annotation to Phase plot:
fstr = "PM=" * @sprintf("%.1f°", pmargin)
	add(plot, atext(fstr, y=(phase0-180)/2, xoffset_rel=0.5, font=afont, align=:cc, strip=2))
add(plot, hmarker(phase0, markerline, strip=2))

#==Render plot
===============================================================================#
gplot = display(InspectDR.GtkDisplay(), mplot)

mplot.layout[:halloc_plot] = 800; mplot.layout[:valloc_plot] = 500
InspectDR.write_png("export_bode.png", mplot)
#InspectDR.write_svg("export_bode.svg", mplot)
#InspectDR.write_eps("export_bode.eps", mplot)
#InspectDR.write_pdf("export_bode.pdf", mplot)

:DONE