text:
# choleskyinv.jl
# Cholesky factor and its inverse (complex-conjugate) transposed in one pass.
# For small matrices this is faster than computing the Cholesky factor
# using the `cholesky` function in LinearAlgebra and then invert it.

"""
Cholesky factorization and inverse Cholesky factorization, accessible
in fields `.c` ans `.ci`, respectively.
"""
struct CholeskyInv{T, F<:Factorization{T}} <: Factorization{T}
     c::F
     ci::F
end

"""
    choleskyinv(P::AbstractMatrix{T};
		check::Bool=true, tol::Real = √eps(T)) where T<:Union{Real, Complex}

 Compute the Cholesky factorization of a dense positive definite
 matrix P and return a `Choleskyinv` object, holding in field `.c`
 the Cholesky factorization and in field `ci` the inverse of the Cholesky
 factorization.

 The two factorizations are obtained in one pass and this is faster
 then calling Julia's [chelosky](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.cholesky)
 function and inverting the lower factor for small matrices.

 Input matrix `P` may be of type `Matrix` or `Hermitian`. Since only the
 lower triangle is used, `P` may also be a `LowerTriangular` view of a
 positive definite matrix.
 If `P` is real, it can also be of the `Symmetric` type.

 The algorithm is a *multiplicative Gaussian elimination*.

 **Notes:**
 The inverse Cholesky factor ``L^{-1}``, obtained as `.ci.L`,
 is an inverse square root (whitening matrix) of `P`, since
 ``L^{-1}PL^{-H}=I``. It therefore yields the inversion of ``P`` as
 ``P^{-1}=L^{-H}L^{-1}``. This is the fastest whitening matrix to be computed,
 however it yields poor numerical precision, especially for large matrices.

 The following relations holds:
 - ``L=PL^{-H}``,
 - ``L^{H}=L^{-1}P``,
 - ``L^{-H}=P^{-1}L``
 - ``L^{-1}=L^{H}P^{-1}``.

 We also have ``L^{H}L=L^{-1}P^{2}L^{-H}=UPU^H``, with ``U`` orthogonal
 (see below) and ``L^{-1}L^{-H}=L^{H}P^{-2}L=UP^{-1}U^H``.
 ``LL^{H}`` and ``L^{H}L`` are unitarily similar, that is,

 ``ULL^{H}U^H=L^{H}L``,

 where ``U=L^{-1}P^{1/2}``, with ``P^{1/2}=H`` the *principal* (unique symmetric)
 square root of ``P``. This is seen writing
 ``PP^{-1}=HHL^{-H}L^{-1}``; multiplying both sides on the left by ``L^{-1}``
 and on the right by ``L`` we obtain
 ``L^{-1}PP^{-1}L=L^{-1}HHL^{-H}=I=(L^{-1}H)(L^{-1}H)^H`` and since
 ``L^{-1}H`` is square it must be unitary.

 From these expressions we have
 - ``H=LU=U^HL^H,
 - ``L=HU^H``,
 - ``H^{-1}=U^HL^{-1}``
 - ``L^{-1}=UH^{-1}``.

 ``U`` is the *polar factor* of ``L^{H}``, *i.e.*, ``L^{H}=UH``,
 since ``LL^{H}=HU^HUH^H=H^2=P``.
 From ``L^{H}L=UCU^H`` we have ``L^{H}LU=UC=ULL^{H}`` and from
 ``U=L^{-1}H`` we have ``L=HU^H``.

 ## Examples
	using PosDefManifold, Test
	n = 40
	etol = 1e-9
	Y=randP(n)
	Yi=inv(Y)
	a=cholesky(Y)

	C=choleskyinv!(copy(Matrix(Y)))
	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)

	# repeat the test for complex matrices
	Y=randP(ComplexF64, n)
	Yi=inv(Y)
	C=choleskyinv!(copy(Matrix(Y)))
	@test(norm(C.c.L*C.c.U-Y)/√n < etol)
	@test(norm(C.ci.U*C.ci.L-Yi)/√n < etol)

	# Benchmark
	using BenchmarkTools

	# computing the Cholesky factor and its inverse using LinearAlgebra
	function linearAlgebraWay(P)
		C=cholesky(P)
		Q=inv(C.L)
	end

	Y=randP(n)
	@benchmark(choleskyinv(Y))
	@benchmark(linearAlgebraWay(Y))

"""
choleskyinv(P::Union{Hermitian, Symmetric, Matrix, LowerTriangular};
	   		check::Bool = true,
	   		tol::Real = √eps(real(eltype(P)))) =
    choleskyinv!(copy(Matrix(P)); check=check, tol=tol)

"""
    choleskyinv!(P::AbstractMatrix{T};
		kind::Symbol = :LLt, tol::Real = √eps(T)) where T<:RealOrComplex
 The same thing as [`choleskyinv`](@ref), but destroys the input matrix.
"""
function choleskyinv!(P::Matrix{T};
			check::Bool = true,
			tol::Real = √eps(real(T))) where T<:Union{Real, Complex}
	LinearAlgebra.require_one_based_indexing(P)
	n = LinearAlgebra.checksquare(P)

	@inbounds for j=1:n-1
		check && abs2(P[j, j])<tol && throw(LinearAlgebra.PosDefException(1))
		for i=j+1:n
			θ = conj(P[i, j] / -P[j, j])
			for k=i:n P[k, i] += θ * P[k, j] end # update P and write D
			P[i, j] = conj(-θ) # write Cholesky factor
			for k=1:j-1 P[k, i] += θ * P[k, j] end # write inv Cholesky factor
			P[j, i] = θ # write inv Cholesky factor
		end
	end

	LT, ULT, UUT = LowerTriangular, UnitLowerTriangular, UnitUpperTriangular
	D=sqrt.(Diagonal(P))
	return CholeskyInv(Cholesky(LT(ULT(P)*D), :L, 0), Cholesky(LT(Matrix((UUT(P)*inv(D))')) , :L, 0 ))
end

function show(io::IO, ::MIME{Symbol("text/plain")}, C::CholeskyInv)
    println(io, "Cholesky (.c) and inverse Cholesky (.ci) factorizations")
	#show(io, C.c)
	#show(io, C.ci)
end


text:
# Routines for SPN evaluation and sampling

export
    logpdf,
    rand

"""
    (spn::SumProductNetwork)(x::AbstractVector{<:Real})
    (spn::SumProductNetwork)(x...)

Evaluates the sum-product network at a given instantiation `x` of the variables.
Summed-out variables are represented as `NaN`s.

# Parameters

- `x`: vector of values of variables (integers or reals).

# Examples

To compute the probability of ``P(b=1)`` using spn `S`, use
```julia
julia> S = SumProductNetwork(IOBuffer("1 + 2 0.2 3 0.5 4 0.3\n2 * 5 7\n3 * 5 8\n4 * 6 8\n5 categorical 1 0.6 0.4\n6 categorical 1 0.1 0.9\n7 categorical 2 0.3 0.7\n8 categorical 2 0.8 0.2"));

julia> S([NaN, 2])
0.3

julia> S(NaN, 2)
0.3
```
"""
function (spn::SumProductNetwork)(x::AbstractVector{<:Real})
    return exp(logpdf(spn,x))
end
function (spn::SumProductNetwork)(x...)
    return exp(logpdf(spn,[x...]))
end

"""
    logpdf(spn, X)

Returns the sums of the log-probabilities of instances `x` in `X`. Uses multithreading to speed up computations if `JULIA_NUM_THREADS > 1`.

# Parameters

- `spn`: Sum-Product Network
- `X`: matrix of values of variables (integers or reals). Summed-out variables are represented as `NaN`s.
"""
function logpdf(spn::SumProductNetwork, X::AbstractMatrix{<:Real})::Float64
    # single-threaded version
    if Threads.nthreads() == 1
        values = Array{Float64}(undef,length(spn))
        return sum(logpdf!(values,spn,view(X,i,:)) for i=1:size(X,1))
    end
    # multi-threaded version
    values = Array{Float64}(undef,length(spn),Threads.nthreads())
    s = Threads.Atomic{Float64}(0.0)
    Threads.@threads for i=1:size(X,1)
        Threads.atomic_add!(s, logpdf!(view(values,:,Threads.threadid()),spn,view(X,i,:)))
    end
    s[]
end

"""
    logpdf!(values, spn, X)

Computes log-probabilities of instances `x` in `X` and stores the results in a given vector. Uses multithreading to speed up computations.

# Parameters

- `results`: vector to store results (log-probabilities)
- `spn`: Sum-Product Network
- `X`: matrix of values of variables (integers or reals). Summed-out variables are represented as `NaN`s.
"""
function logpdf!(results::AbstractVector{<:Real}, spn::SumProductNetwork, X::AbstractMatrix{<:Real})
    @assert length(results) == size(X,1)
    # multi-threaded version
    values = Array{Float64}(undef,length(spn),Threads.nthreads())
    Threads.@threads for i=1:size(X,1)
        @inbounds results[i] = logpdf!(view(values,:,Threads.threadid()), spn, view(X,i,:))
    end
    Nothing
end

"""
    logpdf(spn, x)

Evaluates the sum-product network `spn` in log domain at configuration `x`.

# Parameters

- `x`: vector of values of variables (integers or reals). Summed-out variables are represented as `NaN`s

# Examples

To compute the probability of ``P(b=1)`` using spn `S`, use
```julia
julia> S = SumProductNetwork(IOBuffer("1 + 2 0.2 3 0.5 4 0.3\n2 * 5 7\n3 * 5 8\n4 * 6 8\n5 categorical 1 0.6 0.4\n6 categorical 1 0.1 0.9\n7 categorical 2 0.3 0.7\n8 categorical 2 0.8 0.2"));
julia> logpdf(S,[NaN, 2])
-1.2039728043259361

```
"""
function logpdf(spn::SumProductNetwork, x::AbstractVector{<:Real})::Float64
    values = Array{Float64}(undef,length(spn))
    return logpdf!(values,spn,x)
end

"""
    logpdf!(values,spn,x)

Evaluates the sum-product network `spn` in log domain at configuration `x` and stores values of each node in the vector `values`.
"""
function logpdf!(values::AbstractVector{Float64}, spn::SumProductNetwork, x::AbstractVector{<:Real})::Float64
    # @assert length(values) == length(spn)
    # traverse nodes in reverse topological order (bottom-up)
    @inbounds for i in length(spn):-1:1
        node = spn[i]
        if isprod(node)
            lval = 0.0
            for j in node.children
                lval += values[j]
            end
            values[i] = isfinite(lval) ? lval : -Inf
        elseif issum(node)
            # log-sum-exp trick to improve numerical stability
            m = -Inf
            # get maximum incoming value
            for j in node.children
                m = values[j] > m ? values[j] : m
            end
            lval = 0.0
            for (k,j) in enumerate(node.children)
                # ensure exp in only computed on nonpositive arguments (avoid overflow)
                lval += exp(values[j]-m)*node.weights[k]
            end
            # if something went wrong (e.g. incoming value is NaN or Inf) return -Inf
            values[i] = isfinite(lval) ? m + log(lval) : -Inf
        else # is a leaf node
            values[i] = logpdf(node, x[node.scope])
        end
    end
    @inbounds return values[1]
end

"""
    plogpdf!(values,spn,layers,x)

Evaluates the sum-product network `spn` in log domain at configuration `x` using the scheduling in `nlayers` as obtained by the method `layers(spn)`. Stores values of each node in the vector `values`.

# Parameters

- `values`: vector to cache node values
- `spn`: the sum product network
- `nlayers`: Vector of vector of node indices determinig the layers of the `spn`; each node in a layer is computed based on values of nodes in smaller layers.
- `x`: Vector containing assignment

"""
function plogpdf!(values::AbstractVector{Float64}, spn::SumProductNetwork, nlayers, x::AbstractVector{<:Real})::Float64
    # visit layers from last (leaves) to first (root)
    @inbounds for l in length(nlayers):-1:1
        # parallelize computations within layer
        Threads.@threads for i in nlayers[l]
            node = spn[i]
            if isprod(node)
                lval = 0.0
                for j in node.children
                    lval += values[j]
                end
                values[i] = isfinite(lval) ? lval : -Inf
            elseif issum(node)
                # log-sum-exp trick to improve numerical stability (assumes weights are normalized)
                m = -Inf
                # get maximum incoming value
                for j in node.children
                    m = values[j] > m ? values[j] : m
                end
                lval = 0.0
                for (k,j) in enumerate(node.children)
                    # ensure exp in only computed on nonpositive arguments (avoid overflow)
                    lval += exp(values[j]-m)*node.weights[k]
                end
                # if something went wrong (e.g. incoming value is NaN or Inf) return -Inf, otherwise
                # return maximum plus lval (adding m ensures signficant digits are numerically precise)
                values[i] = isfinite(lval) ? m+log(lval) : -Inf
            else # is a leaf node
                values[i] = logpdf(node, x[node.scope])
            end
        end
    end
    @inbounds return values[1]
end

"""
    plogpdf(spn, x)

Parallelized version of `logpdf(spn, x)`. Set `JULIA_NUM_THREADS` > 1 to make this effective.
"""
function plogpdf(spn::SumProductNetwork, x::AbstractVector{<:Real})::Float64
    values = Array{Float64}(undef,length(spn))
    return plogpdf!(values,spn,layers(spn),x)
end

"""
    sample(weights)::UInt

Sample integer with probability proportional to given `weights`.

# Parameters

- `weights`: vector of nonnegative reals.
"""
function sample(weights)::UInt
    Z = sum(weights)
    u = rand()
    cum = 0.0
    for i = 1:length(weights)
        @inbounds cum += weights[i]/Z
        if u < cum
            return i
        end
    end
end

"""
    rand(n::IndicatorFunction)
    rand(n::CategoricalDistribution)
    rand(n::GaussianDistribution)

Sample values from sum-product network leaves.
"""
@inline Base.rand(n::IndicatorFunction) = n.value
@inline Base.rand(n::CategoricalDistribution) = sample(n.values)
@inline Base.rand(n::GaussianDistribution) = n.mean + sqrt(n.variance)*randn()

"""
    rand(spn)

Returns a sample of values of the variables generated according
to the probability defined by the network `spn`. Stores the sample
as a vector of values

# Example

```julia
julia> rand(spn)
[2, 1]
```
"""
function Base.rand(spn::SumProductNetwork)
    if length(scope(spn)) > 0
        numvars = length(scope(spn))
    else
        numvars = length(union(n.scope for n in leaves(spn)))
    end
    a = Vector{Float64}(undef, numvars)
    # sample induced tree
    queue = [1]
    while !isempty(queue)
        n = popfirst!(queue)
        if issum(spn[n])
            # sample one child to visit
            c = sample(spn[n].weights)
            push!(queue, spn[n].children[c])
        elseif isprod(spn[n])
            # visit every child
            append!(queue, children(spn,n))
        else
            # draw value from node distribution
            a[ spn[n].scope ] = rand(spn[n])
        end
    end
    return a
end

"""
  rand(spn::SumProductNetwork, N::Integer)

Returns a matrix of samples generated according to the probability
defined by the network `spn`.
"""
function Base.rand(spn::SumProductNetwork, N::Integer)
    if length(scope(spn)) > 0
        numvars = length(scope(spn))
    else
        numvars = length(union(n.scope for n in leaves(spn)))
    end
    Sample = Array{Float64}(undef,N,numvars)
    # get sample
    for i=1:N
        queue = [1]
        while length(queue) > 0
            n = popfirst!(queue)
            if issum(spn[n])
                # sample one child to visit
                c = sample(spn[n].weights) # sparse array to vector inserts 0 at first coordinate
                push!(queue, spn[n].children[c])
            elseif isprod(spn[n])
                # visit every child
                append!(queue, children(spn,n))
            else
                # draw value from distribution
                Sample[ i, spn[n].scope ] = rand(spn[n])
            end
        end
    end
    return Sample
end

"""
    ncircuits(spn::SumProductNetwork)

Counts the number of induced circuits of the sum-product network `spn`.
"""
ncircuits(spn::SumProductNetwork) = ncircuits!(Array{Int}(undef,length(spn)),spn)
"""
    ncircuits!(values,spn)

Counts the number of induced circuits of the sum-product network `spn`, caching intermediate values.
"""
function ncircuits!(values::AbstractVector{Int}, spn::SumProductNetwork)
    @assert length(values) == length(spn)
    # traverse nodes in reverse topological order (bottom-up)
    for i in length(spn):-1:1
        @inbounds node = spn[i]
        if isa(node,ProductNode)
            @inbounds values[i] = mapreduce(j -> values[j], *, node.children)
        elseif isa(node,SumNode)
            @inbounds values[i] = sum(values[node.children])
        else # is a leaf node
            @inbounds values[i] = 1
        end
    end
    @inbounds values[1]
end

"""
    NLL(spn::SumProductNetwork,data::AbstractMatrix{<:Real})

Computes the average negative loglikelihood of a dataset `data` assigned by spn.
"""
NLL(spn::SumProductNetwork,data::AbstractMatrix{<:Real}) = -logpdf(spn,data)/size(data,1)

"""
    MAE(S1::SumProductNetwork,S2::SumProductNetwork,data::AbstractMatrix{<:Real})

Computes the Mean Absolute Error of sum-product networks `S1` and `S2` on given `data`.
"""
MAE(S1::SumProductNetwork,S2::SumProductNetwork,data::AbstractMatrix{<:Real}) = sum( abs(logpdf(S1,view(data,i,:))-logpdf(S2,view(data,i,:))) for i=1:size(data,1))/size(data,1)

ext: jl


text:
using DataFrames
using Grids
using Crayons

action_space_size(g::MuZeroGame)                        = action_space_size(g.game)
to_play(g::MuZeroGame)                                  = gettoplay(g.game)
to_play_at(g::MuZeroGame, idx::Integer)                 = to_play_at(g.game, idx)
next_to_play(g::MuZeroGame)                             = getnexttoplay(g.game)
discount(g::MuZeroGame)                                 = discount(g.game)
reward_matrix(g::MuZeroGame)                            = reward_matrix(g.game)
value_matrix(g::MuZeroGame, idx::Integer; kwargs...)    = value_matrix(g.game, idx; kwargs...)
uniform_policy(g::MuZeroGame)                           = uniform_policy(g.game)
history(g::MuZeroGame)                                  = gethistory(g.game)
is_terminal(g::MuZeroGame)                              = isterminal(g.game)
legal_actions(g::MuZeroGame)                            = getlegalactions(g.game)
action_history(g::MuZeroGame, idx)                      = action_history(g.game, idx)
action_history(g::MuZeroGame)                           = action_history(g.game)
apply!(g::MuZeroGame, action::Integer)                  = step!(g.game, action)
make_image(g::MuZeroGame, state_idx::Integer)           = make_image(g.game, state_idx)
#getstate(g::MuZeroGame, state_idx::Integer)             = get_state(g, state_idx)
get_state(g::MuZeroGame, state_idx::Integer)            = getstate(g.game, state_idx)
make_image_label(g::MuZeroGame, state_idx::Integer)     = make_image_label(g.game, state_idx)
get_state(g::MuZeroGame)                                = getstate(g.game)

struct MuZeroGameSVG
    game::MuZeroGame
end

getgame(mzg::MuZeroGame) = mzg.game
alltrue(itr) = length(itr)==0 || all(identity, itr)
is_terminal(games::Vector) = is_terminal.(games) |> alltrue
movesmaxedout(game::MuZeroGame, max_moves) = length(history(game)) >= max_moves
movesmaxedout(games::Vector, max_moves) = movesmaxedout.(games, max_moves) |> alltrue

function svg(g::MuZeroGame)
    MuZeroGameSVG(g)
end

function Base.show(io::IO, m::MIME"image/svg+xml", g::MuZeroGameSVG)
    show(io, m, getmoves(g.game.game))
end

Base.getindex(x::MuZeroGame, idx::Integer) = getindex(x.game, idx)

function df(x::MuZeroGame; range=1:length(x.game.states), limit=3)
    g = x.game
    state_col = "State #"
    action_col = "Action"
    reward_col = "Reward"
    d = Dict(
        state_col=>[],
        action_col=>[],
        reward_col=>[],
    )
    for s_idx in range
        action = s_idx > length(g.actions) ? nothing : g.actions[s_idx]
        reward = s_idx > length(g.rewards) ? nothing : g.rewards[s_idx]
        push!(d[state_col], s_idx)
        push!(d[action_col], isnothing(action) ? nothing : action.idx)
        push!(d[reward_col], reward)
    end
    DataFrame(
        state_col=>d[state_col],
        action_col=>d[action_col],
        reward_col=>d[reward_col],
    )
end

function Base.show(io::IO, m::MIME"text/plain", x::MuZeroGame)
    child_visit_strings = []
    for visits in x.child_visits
        from = 1
        limit = 5
        to = min(limit, length(visits))
        actions = reverse(sortperm(visits))[from:to]
        prs = visits[actions]
        action_prs = map(
            x->begin
                action, pr = x
                pr_s = @sprintf "%.2f" pr
                "$action ($(pr_s))"
            end,
            zip(actions, prs)
        )
        push!(child_visit_strings, join(action_prs, ", "))
    end
    game_df = df(x)
    root_values_str = map(root_value->begin
        @sprintf "%.2f" root_value
    end, x.root_values)
    game_df[!, "V(S)"] = [root_values_str..., nothing]
    game_df[!, "Sampled Actions (Top 5)"] = [child_visit_strings..., nothing]
    game_df = game_df[!, [1, 4, 5, 2, 3]]
    println()
    pretty_game(
        io,
        game_df;
        title="MuZero Game History"
    )
end

function pretty_game(io::IO, df::DataFrame; title)
    reward_hl = Highlighter((data, i, j)->(j==5) && !isnothing(data[i,j]) && data[i,j] > 0 , crayon"green");
    pretty_table(
        io,
        df;
        show_row_number=false,
        nosubheader=true,
        title=title,
        highlighters=(reward_hl),
    )
end

# @todo make this more readable
# @todo use exclamation point to indicate it is mutating
function store_search_statistics(game::MuZeroGame, root::Node)
    sum_visits = sum([child.visit_count for child in values(root.children)])
    action_space = [idx for idx in 1:action_space_size(game)]
    push!(game.child_visits,
        [in(a, keys(root.children)) ? (root.children[a].visit_count / sum_visits) : 0 for a in action_space]
    )
    push!(game.root_values, calc_value(root))
end

# @todo This function needs extensive unit testing for edge cases, especially for non-endgame rewards.
function make_target(game::MuZeroGame, state_index::Integer, num_unroll_steps::Integer, td_steps::Integer)
    # The value target is the discounted root value of the search tree N steps
    # into the future, plus the discounted sum of all rewards until then.
    targets = []
    rewards = reward_matrix(game)
    num_rewards = size(rewards, 1)

    for (target_id, current_index) in enumerate(state_index:(state_index + num_unroll_steps))
        if current_index <= num_rewards
            player = to_play_at(game, current_index-1)
            next_player = to_play_at(game, current_index)
            bootstrap_index = current_index + td_steps
            value = 0
            if bootstrap_index <= length(game.root_values)
                @assert false "This should not be executing for board games"
                value = game.root_values[bootstrap_index] * discount(game)^td_steps
            end
            # For simplicity the network always predicts the most recently received
            # reward, even for the initial representation network where we already
            # know this reward.
            reward = reward_matrix(game)[current_index, player]
            value += value_matrix(game, current_index, td_steps=td_steps)[1][next_player]
            try
                push!(targets, Target(value, reward, policy_matrix(game)[current_index]))
            catch e
                @show game.uuid, current_index, policy_matrix(game)
                throw(e)
            end
        else
            # States past the end of games are treated as absorbing states.
            push!(targets, Target(0, 0, uniform_policy(game)))
        end
    end
    targets
end

function policy_matrix(g::MuZeroGame)
    [g.child_visits..., uniform_policy(g)]
end

function Base.show(io::IO, m::MIME"text/plain", g::MuZeroGame)
    show(io, "text/plain", g.game)
    println(io, "Root Values: $(join(g.root_values, ", "))")
    action_visits = zip(g.game.history, g.child_visits)
    child_visits_str = join(["[$(join(child_visits, ", "))] -> A($action)" for (action, child_visits) in action_visits], "\n")
    println(io, "Child Visits:\n$child_visits_str")
end

# Base.copy(g::AbstractBoardGame) = copy_at(g, length(g.states))

# function copy_at(g::AbstractBoardGame, state_idx::Integer)
#     typeof(g)(
#         copy(g.history[1:(state_idx-1)]),
#         copy(g.states[1:state_idx]),
#         copy(g.rewards[1:(state_idx-1)]),
#         # @todo child_visits and root_values should not be
#         # copied for AbstractEnv
#         copy(g.child_visits[1:(state_idx-1)]),
#         copy(g.root_values[1:(state_idx-1)]),
#     )
# end


text:
using DifferentialEquations
using SimpleDiffEq
using Tables
using DataFrames
using StatsPlots
using BenchmarkTools

function sir_ode!(du,u,p,t)
    (S,I,R) = u
    (β,c,γ) = p
    N = S+I+R
    @inbounds begin
        du[1] = -β*c*I/N*S
        du[2] = β*c*I/N*S - γ*I
        du[3] = γ*I
    end
    nothing
end;

δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)

u0 = [990.0,10.0,0.0]; # S,I,R

p = [0.05,10.0,0.25]; # β,c,γ

prob_ode = ODEProblem(sir_ode!, u0, tspan, p);

sol_ode = solve(prob_ode, dt = δt);

df_ode = DataFrame(Tables.table(sol_ode'))
rename!(df_ode,["S","I","R"])
df_ode[!,:t] = sol_ode.t;

@df df_ode plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number")

@benchmark solve(prob_ode, dt = δt);


text:
# Run all the code up until the variable definitions to load functions to memory

using TestImages #You don't need this line when using local files
using Images

##### Define the fitness function

function fitness(a,b)
    score = sqrt(512*512-sum((a - b).^2))
    return score
end

# gets the difference between the value (0 is black, 1 is white) of all pixels in one of the population arrays and the pixels in the target array

function all_fit(pop, t)
    fitness_vec =  zeros(N)
    for i in 1:N
        fitness_vec[i] = fitness(pop[i,:,:], t)
        #println(fitness_vec[i])
    end
    return fitness_vec
end

# Generates a child using parameters set below

function genfilho(cdf,pop,crossover, prob_m, hillclimbing, sigma_hc, dim)
    if prob_m > rand()
        filho = rand(dim[1],dim[2])
    else
        min_val = rand(1)
        cdf[cdf .> min_val]
        index_pai = sum((cdf .< min_val))+1
        pai = pop[index_pai]
        mae = zeros(dim[1],dim[2])
        if crossover
            index_mae = index_pai
            while index_mae == index_pai
                min_val = rand(1)
                cdf[cdf .> min_val]
                index_mae = sum((cdf .< min_val))+1
            end
            mae = pop[index_mae]
        end
        filho = 0.5*pai+0.5*mae
        if hillclimbing
            filho+= sigma_hc*randn(dim[1],dim[2])
        end
    end
    return filho
end

# get the best members of the population
function get_best_members(fitness_vec_prob,pop,N,percentage_best)
    ord_vec = reverse(sortperm(fitness_vec_prob))
    filhos = pop[ord_vec[1:Int(floor(N*percentage_best))],:,:]
    return filhos
end

#Create the next generation

function epoch(pop, dim, target, N, percentage_best, sigma, crossover, prob_m, hillclimbing, sigma_hc)

    top_idx = Int(floor(N*percentage_best))
    fitness_vec = all_fit(pop, target_float)
    println(mean(fitness_vec))
    fitness_vec_prob=(fitness_vec)/sum(fitness_vec)
    cdf = cumsum(fitness_vec_prob)

    new_pop = zeros(N,dim[1], dim[2])
    new_pop[1:top_idx,:,:] = get_best_members(fitness_vec_prob, pop,N, percentage_best)

    num_best_members = Int(floor(N*percentage_best))

    while num_best_members < N
        new_pop[num_best_members,:,:] = (genfilho(cdf,pop,crossover, prob_m, hillclimbing, sigma_hc, dim))
        num_best_members +=1
    end
    return new_pop
end

############################# STOP FIRST RUN HERE #######################

#Show a sample start image
sample_img = rand(512,512)
colorview(Gray, sample_img)

#Define your variables
target = testimage("livingroom")
target_float = float(target) #convert img to array
dim = size(target_float) #get size of array
N = 100
percentage_best = 0.15
sigma = 0.05
pop = rand(N,dim[1], dim[2]) #empty pop vector
nr_epochs = 300
crossover = true
prob_m = 0.05
hillclimbing = false
sigma_hc = 0.05

#loop all the code "nr_epochs" times

for i in 1:nr_epochs
    colorview(Gray, pop[1,:,:])
    pop = epoch(pop, dim, target_float, N, percentage_best, sigma, crossover, prob_m, hillclimbing, sigma_hc)
end

max_stars_repo_name: quendera/EvolveImage

max_issues_repo_name: quendera/EvolveImage

max_forks_repo_name: quendera/EvolveImage


text:
module Lp_glpk
  using Lp
  # Waiting for https://github.com/JuliaLang/julia/pull/6884
  importall JuMP
  importall GLPKMathProgInterface

  function install()
    Pkg.add("JuMP")
    Pkg.add("GLPK")
    Pkg.add("GLPKMathProgInterface")
  end

  type Dat{T}
    prob::Lp.Cf0_problem{T}

    Dat() = new()
  end

  function construct_dat(T::DataType) return Dat{T}() end

  function fill_dat{T}(prob::Lp.Cf0_problem, dat::Dat{T})
    dat.prob = prob
  end

  function solve_dat{T}(dat::Dat{T}, sol::Lp.Solution{T})
    m = Model(solver = GLPKSolverLP(method=:Exact))

    lp_prob = dat.prob
    @defVar(m, lp_x[1:lp_prob.n] >= 0 )
    for ri = 1:lp_prob.m
      aff = AffExpr(lp_x[1:lp_prob.n], reshape(lp_prob.A[ri,1:lp_prob.n], lp_prob.n), 0.0)
      @addConstraint(m, aff <= lp_prob.b[ri])
    end
    obj = AffExpr(lp_x[1:lp_prob.n], lp_prob.c, 0.0)
    setObjective(m, :Min, obj)

    @time status = JuMP.solve(m)

    status_dict = { :Optimal => :Optimal, :Unbounded => :Unbounded, :Infeasible => :Infeasible, :UserLimit => :Maxit, :Error => :Error, :NotSolved => :Created }
    sol.status = status_dict[status]
    sol.solved = (sol.status == :Optimal)
    sol.z = getObjectiveValue(m)
    if (sol.solved)
      sol.x = Array(lp_prob.conv.t, lp_prob.n)
      for i=1:lp_prob.n
        sol.x[i] = JuMP.getValue(lp_x[i])
      end
    end
  end

end

max_stars_repo_name: jadnohra/mlcp-lib

max_issues_repo_name: jadnohra/mlcp-lib

max_forks_repo_name: jadnohra/mlcp-lib


text:
module aoc_23
    using CircularArrays
    using CircularList
    using Test

    const input = "685974213"
    const input_array = CircularArray( parse.(UInt8, Vector{Char}(input) ) )
    const demo = "389125467"
    const demo_array = CircularArray( parse.(UInt8, Vector{Char}(demo) ) )

    const demo_destinations = [2, 7, 3, 7, 3, 9, 8, 1, 5, 3]

    # CircularArray delete / insert enhancement
    include("circulararray.jl")

    # CircularArrayList actually subtypes AbstractArray
    # We also added a dict into it
    include("CircularListArray.jl")

    function play( input = input_array , n = 100)
        select = 1
        #A = deepcopy(input)
        A = input
        # @info "Finding min and max"
        min_cup = minimum(A)
        # @info "Minimum", min_cup
        max_cup = maximum(A)
        # @info "Maximum", max_cup

        # println(A)
        # println(typeof(A))

        for i=1:n
            if i % 1_000_000 == 0
                @info i
            end
            selected = select+1:select+3
            current = A[ select ]
            destination = current - 1
            # For CircularListArray the following will return a CircularList
            snip = A[ selected ]
            if !isa(A, CircularListArray)
                # We will does this while we insert
                deleteat!(A, selected)
            end
            if destination < min_cup
                destination = max_cup
            end
            while destination ∈ snip
                destination -= 1
                #println("Destination", destination)
                if destination < min_cup
                    destination = max_cup
                end
            end
            #println("Final destination: ", destination )
            #println(A)
            insertaftervalue!(A, destination, snip)
            #println(A)
            if isa( A, CircularListArray )
                jump!( A.list, next(A.list) )
                select = 1
            else
                select = findfirst(isequal(current), A) + 1
            end
        end
        A
    end

    function insertaftervalue!(A, value, inserter)
        dest_ind = findfirst(isequal(value), A) + 1
        insert!(A, dest_ind, inserter)
    end

    function insertaftervalue!(A::CircularListArray, value, inserter)
        # c = current(A.list)
        node = findfirstnode(isequal(value), A)
        # node = A.dict[value]
        # jump!(A.list, node.next)
        insert!(A, node, inserter)
        # jump!(A.list, c)
    end

    function report( A )
        cup_1 = findfirst(isequal(1), A)
        print.(A[cup_1+1:cup_1+8])
        nothing
    end

    part1() = report( play() ) # 82635947

    function test_part1( input = input_array )
        part1_real_answer = [8,2,6,3,5,9,4,7]

        A = play( deepcopy( input ) )
        cup_1 = findfirst(isequal(1), A)
        @test all( .==( A[cup_1+1:cup_1+8], part1_real_answer ) )

        A = play( CircularListArray( deepcopy( input ) ) )
        cup_1 = findfirst(isequal(1), A)
        @test all( .==( A[cup_1+1:cup_1+8], part1_real_answer ) )
    end

    function part2( input = input_array, n = 10_000_000 )
        input = deepcopy( input )
        A = CircularListArray([input; maximum(input)+1:1_000_000])
        @info "One million array created"
        play( A , n )
        cup_1 = findfirst(isequal(1), A)
        println( A[cup_1+1:cup_1+2] )
        A
        # 157047826689
    end
end


text:
import LowLevelJPEG: @code_str

@testset "DCSymbol" begin
    @testset "Constructors" begin
        @test DCSymbol(0x12).bits === UInt8(0x12)
    end

    @testset "Traits" begin
        @test LowLevelJPEG.amplitude(DCSymbol(0x0A)) == 10
    end

    @testset "Display" begin
        @test shown_string(DCSymbol(1)) == "DCSymbol(0x01)"
        @test shown_string([DCSymbol(2), DCSymbol(10)]) == "DCSymbol[0x02, 0x0a]"
    end

    @testset "Conversions" begin
        @test DCSymbol[0x00, 0x01] == [DCSymbol(0x00), DCSymbol(0x01)]
    end
end

@testset "ACSymbol" begin
    @testset "Constructors" begin
        @test ACSymbol(0x12).bits === UInt8(0x12)
        @test ACSymbol(10, 15).bits === UInt8(0xAF)
        @test ACSymbol(10, 0) === LowLevelJPEG.EOB10
    end

    @testset "Traits" begin
        @test LowLevelJPEG.runlength(ACSymbol(10, 15)) == 10
        @test LowLevelJPEG.amplitude(ACSymbol(10, 15)) == 15

        @test LowLevelJPEG.iszrl(LowLevelJPEG.ZRL)
        @test !LowLevelJPEG.iszrl(ACSymbol(10,15))

        @test LowLevelJPEG.iseob(LowLevelJPEG.EOB)
        @test !LowLevelJPEG.iseob(LowLevelJPEG.EOB1)

        @test LowLevelJPEG.iseobn(LowLevelJPEG.EOB)
        @test !LowLevelJPEG.iseobn(LowLevelJPEG.ZRL)
        @test LowLevelJPEG.iseobn(ACSymbol(10,0))
        @test !LowLevelJPEG.iseobn(ACSymbol(10,1))
    end

    @testset "Display" begin
        @test shown_string(ACSymbol(1, 2)) == "ACSymbol(0x12)"
        @test shown_string([ACSymbol(3, 4), ACSymbol(5, 6)]) == "ACSymbol[0x34, 0x56]"
    end

    @testset "Conversions" begin
        @test ACSymbol[0x00, 0x01] == [ACSymbol(0x00), ACSymbol(0x01)]
    end
end

@testset "DCCoeff" begin
    @testset "Constructors" begin
        @test DCCoeff( 0) === DCCoeff(DCSymbol(0x00), code"")
        @test DCCoeff(-1) === DCCoeff(DCSymbol(0x01), code"0")
        @test DCCoeff(+1) === DCCoeff(DCSymbol(0x01), code"1")
        @test DCCoeff(-3) === DCCoeff(DCSymbol(0x02), code"00")
        @test DCCoeff(-2) === DCCoeff(DCSymbol(0x02), code"01")
        @test DCCoeff(+2) === DCCoeff(DCSymbol(0x02), code"10")
        @test DCCoeff(+3) === DCCoeff(DCSymbol(0x02), code"11")
        @test DCCoeff(-2047) === DCCoeff(DCSymbol(0x0B), code"00000000000")
        @test DCCoeff(+2047) === DCCoeff(DCSymbol(0x0B), code"11111111111")

        # for 12-bit sample precision
        @test DCCoeff(-32767) === DCCoeff(DCSymbol(0x0F), code"000000000000000")
        @test DCCoeff(+32767) === DCCoeff(DCSymbol(0x0F), code"111111111111111")
    end

    @testset "Display" begin
        @test shown_string(DCCoeff(-10)) == "DCCoeff(-10)"
        @test shown_string(DCTCoeff[DCCoeff(0), DCCoeff(1)]) ==
            "DCTCoeff[DCCoeff(0), DCCoeff(1)]"
        @test shown_string([DCCoeff(-1), DCCoeff(2)]) ==
            "DCTCoeff{DCSymbol}[-1, 2]"
        @test shown_plaintext([DCCoeff(-1), DCCoeff(2)]) ==
            """
            2-element Array{DCTCoeff{DCSymbol},1}:
             -1
             2"""
    end

    @testset "Conversions" begin
        @test all(v -> Int(DCCoeff(v)) == v, -2047:2047)
    end
end

@testset "ACCoeff" begin
    @testset "Constructors" begin
        @test ACCoeff(0) === ACCoeff(ACSymbol(0x00), code"") # EOB

        @test ACCoeff(-1) === ACCoeff(ACSymbol(0x01), code"0")
        @test ACCoeff(+1) === ACCoeff(ACSymbol(0x01), code"1")
        @test ACCoeff(-3) === ACCoeff(ACSymbol(0x02), code"00")
        @test ACCoeff(-2) === ACCoeff(ACSymbol(0x02), code"01")
        @test ACCoeff(+2) === ACCoeff(ACSymbol(0x02), code"10")
        @test ACCoeff(+3) === ACCoeff(ACSymbol(0x02), code"11")
        @test ACCoeff(-1023) === ACCoeff(ACSymbol(0x0A), code"0000000000")
        @test ACCoeff(+1023) === ACCoeff(ACSymbol(0x0A), code"1111111111")

        @test ACCoeff(-4, 5) === ACCoeff(ACSymbol(0x53), code"011")

        # fixes amplitude
        @test ACCoeff(4, ACSymbol(0xAA)) === ACCoeff(ACSymbol(0xA3), code"100")

        # EOBn
        @test ACCoeff(5, LowLevelJPEG.EOB2) === ACCoeff(ACSymbol(0x20), code"01")
        @test ACCoeff(16384, LowLevelJPEG.EOB14) === ACCoeff(ACSymbol(0xE0), code"00000000000000")
        @test ACCoeff(32767, LowLevelJPEG.EOB14) === ACCoeff(ACSymbol(0xE0), code"11111111111111")
        # the result symbol will be fixed to the proper EOBn.
        @test ACCoeff(7, LowLevelJPEG.EOB) === ACCoeff(LowLevelJPEG.EOB2, code"11")
    end

    @testset "Traits" begin
        @test length(ACCoeff(5, 15)) == 16
        @test length(ACCoeff(5, LowLevelJPEG.EOB4)) == 0
    end

    @testset "Display" begin
        @test shown_string(ACCoeff(-10)) == "ACCoeff(-10)"
        @test shown_string(ACCoeff(-10, 2)) == "ACCoeff(-10, 2)"
        @test shown_plaintext(ACCoeff(-10, 2)) ==   """
                                                    3-element DCTCoeff{ACSymbol}:
                                                     0
                                                     0
                                                     -10"""
        @test shown_string(ACCoeff(0)) == "ACCoeff(1, EOB)"
        @test shown_plaintext(ACCoeff(10, LowLevelJPEG.EOB)) ==
            "DCTCoeff{ACSymbol}(10, EOB3)"
        @test shown_string([ACCoeff(10, LowLevelJPEG.EOB), ACCoeff(1,5)]) ===
            "DCTCoeff{ACSymbol}[ACCoeff(10, EOB3), ACCoeff(1, 5)]"
        @test shown_plaintext([ACCoeff(10, LowLevelJPEG.EOB) ACCoeff(1,5)]) ===
            """
            1×2 Array{DCTCoeff{ACSymbol},2}:
             ACCoeff(10, EOB3)  ACCoeff(1, 5)"""
    end

    @testset "Conversions" begin
        @test all(v -> Int(ACCoeff(v, 1)) == v || v == 0, -1023:1023)
        @test all(v -> Int(ACCoeff(v, LowLevelJPEG.EOB)) == v, 1:1024)
    end
end


text:
thisfile = joinpath(datadir, "GermanMotorPremiums", "data1Small.csv")
@test isfile(thisfile)
@time df_tmp = CSV.read(thisfile, allowmissing=:none, types=eltypesData1, categorical=false, rows_for_type_detect=10000);
df_tmp_orig = deepcopy(df_tmp)

# example with only char columns
selected_explanatory_vars_char_only = ["ART_DES_WOHNEIGENTUM","FAMILIENSTAND","GESCHLECHT","FINANZIERUNGSART","STADT","KENNZEICHEN"]
dtmtableMini, sett, df_prepped = prepare_dataframe_for_dtm!(df_tmp[1:100,:], weightcol="EXPOSURE", numcol="LOSS20HALF", denomcol="PREMIUM66", independent_vars=selected_explanatory_vars_char_only);
sett.minWeight = -.2
@info("DTM:Testing build tree (dtm mini; only character independent variables)")
strs, resm = dtm(dtmtableMini, sett)

features = dtmtableMini.features

featuresSubArr = view(features, collect(5:8), :)

i = 2

f = featuresSubArr[i]
featuresSubArr[i].parent.pool

fieldnames(typeof(f))

f.indices
f.indices

szTMP = size(f, 1)
scores = rand(szTMP);
numeratorEst = rand(szTMP);
numerator = rand(szTMP);
denominator = rand(szTMP);
weight = rand(szTMP);


text:
#=
Test the wrapper function.

Reference
---------
1. https://github.com/stevengogogo/FindSteadyStates.jl/issues/29
=#

include("model.jl")

function fixedpoint_gen(func, u0, p, ranges)

    # Define a problem
    de = DEsteady(func=func, p=p, u0= u0, method=SSRootfind())

    j_gen = jacobian(de) # jacobian generator

    # Searching method and domain
    param_gen = ParameterGrid(ranges)

    # Solve
    sols = solve(de, param_gen)

    # Remove similar solutions
    steadies = unique(sols)

    # Jacobian
    jac_ms = j_gen.(steadies)

    # Stability
    stab_modes = StabilityType.(jac_ms)

    return steadies, jac_ms, jac_ms

end

func_ = bistable_model.func
p_ = bistable_model.p
u0_ = bistable_model.u0

fixedpoint_gen(func_, u0_, p_, [1:0.1:2, 1:0.1:2])


text:
using Documenter, ExcelReaders

makedocs(modules = [ExcelReaders],
	sitename = "ExcelReaders.jl",
	analytics = "UA-132838790-1",
	pages = [
        "Introduction" => "index.md"
    ])

deploydocs(repo = "github.com/queryverse/ExcelReaders.jl.git")


text:
"""
    reflectinmodel(x::AbstractMillNode, fm=d -> Dense(d, 10), fa=BagCount ∘ SegmentedMeanMax;
        fsm=Dict(), fsa=Dict(), single_key_identity=true, single_scalar_identity=true)

Build a `Mill.jl` model capable of processing `x`.

All inner `Dense` layers are constructed using `fm`, a function accepting input dimension `d` and
returning suitable model. All aggregation operators are constructed using `fa` in a similar manner.

More fine-grained control can be achieved with `fsm` and `fsa` keyword arguments, which should be
`Dict`s of `c => f` pairs, where `c` is a `String` traversal code from [HierarchicalUtils.jl](@ref) and
`f` is a function. These definitions override `fm` and `fa`.

If a [`ProductNode`](@ref) with only a single child (subtree) is encountered, its final `m` model
is instantiated as `identity` instead of using `fm` and `fsm`. This can be controlled with `single_key_identity`.

Similarly, if an [`ArrayNode`](@ref) contains data `X` where `size(X, 1)` is `1`, the corresponding
model is instantiated as `identity` unless `single_scalar_identity` is `false`.

# Examples
```jldoctest
julia> n1 = ProductNode((; a=ArrayNode(NGramMatrix(["a", "b"]))))
ProductNode with 2 obs
  └── a: ArrayNode(2053×2 NGramMatrix with Int64 elements)

julia> n2 = ProductNode((ArrayNode([0 1]), BagNode(ArrayNode([0 1; 2 3]), bags([1:1, 2:2]))))
ProductNode with 2 obs
  ├── ArrayNode(1×2 Array with Int64 elements)
  └── BagNode with 2 obs
        └── ArrayNode(2×2 Array with Int64 elements)

julia> n = ProductNode((n1, n2))
ProductNode with 2 obs
  ├── ProductNode with 2 obs
  │     └── a: ArrayNode(2053×2 NGramMatrix with Int64 elements)
  └── ProductNode with 2 obs
        ├── ArrayNode(1×2 Array with Int64 elements)
        └── BagNode with 2 obs
              ⋮

julia> printtree(n; trav=true)
ProductNode with 2 obs [""]
  ├── ProductNode with 2 obs ["E"]
  │     └── a: ArrayNode(2053×2 NGramMatrix with Int64 elements) ["M"]
  └── ProductNode with 2 obs ["U"]
        ├── ArrayNode(1×2 Array with Int64 elements) ["Y"]
        └── BagNode with 2 obs ["c"]
              └── ArrayNode(2×2 Array with Int64 elements) ["e"]

julia> reflectinmodel(n) |> printtree
ProductModel … ↦ ArrayModel(Dense(20, 10))
  ├── ProductModel … ↦ ArrayModel(identity)
  │     └── a: ArrayModel(Dense(2053, 10))
  └── ProductModel … ↦ ArrayModel(Dense(11, 10))
        ├── ArrayModel(identity)
        └── BagModel … ↦ BagCount([SegmentedMean(10); SegmentedMax(10)]) ↦ ArrayModel(Dense(21, 10))
              └── ArrayModel(Dense(2, 10))

julia> reflectinmodel(n, d -> Dense(d, 3), SegmentedMean) |> printtree
ProductModel … ↦ ArrayModel(Dense(6, 3))
  ├── ProductModel … ↦ ArrayModel(identity)
  │     └── a: ArrayModel(Dense(2053, 3))
  └── ProductModel … ↦ ArrayModel(Dense(4, 3))
        ├── ArrayModel(identity)
        └── BagModel … ↦ SegmentedMean(3) ↦ ArrayModel(Dense(3, 3))
              └── ArrayModel(Dense(2, 3))

julia> reflectinmodel(n, d -> Dense(d, 3), SegmentedMean;
                        fsm=Dict("e" => d -> Chain(Dense(d, 2), Dense(2, 2))),
                        fsa=Dict("c" => SegmentedLSE),
                        single_key_identity=false,
                        single_scalar_identity=false) |> printtree
ProductModel … ↦ ArrayModel(Dense(6, 3))
  ├── ProductModel … ↦ ArrayModel(Dense(3, 3))
  │     └── a: ArrayModel(Dense(2053, 3))
  └── ProductModel … ↦ ArrayModel(Dense(6, 3))
        ├── ArrayModel(Dense(1, 3))
        └── BagModel … ↦ SegmentedLSE(2) ↦ ArrayModel(Dense(2, 3))
              └── ArrayModel(Chain(Dense(2, 2), Dense(2, 2)))
```

See also: [`AbstractMillNode`](@ref), [`AbstractMillModel`](@ref), [`ProductNode`](@ref), [`ArrayNode`](@ref).
"""
function reflectinmodel(x, fm=d -> Dense(d, 10), fa=BagCount ∘ SegmentedMeanMax; fsm=Dict(),
        fsa=Dict(), single_key_identity=true, single_scalar_identity=true, all_imputing=false)
    _reflectinmodel(x, fm, fa, fsm, fsa, "", single_key_identity, single_scalar_identity, all_imputing)[1]
end

function _reflectinmodel(x::AbstractBagNode, fm, fa, fsm, fsa, s, args...)
    c = stringify(s)
    im, d = _reflectinmodel(x.data, fm, fa, fsm, fsa, s * encode(1, 1), args...)
    agg = haskey(fsa, c) ? fsa[c](d) : fa(d)
    d = size(BagModel(im, agg)(x).data, 1)
    bm = haskey(fsm, c) ? fsm[c](d) : fm(d)
    m = BagModel(im, agg, bm)
    m, size(m(x).data, 1)
end

_remap(data::NamedTuple, ms) = (; zip(keys(data), ms)...)
_remap(::Tuple, ms) = tuple(ms...)

function _reflectinmodel(x::AbstractProductNode, fm, fa, fsm, fsa, s, ski, args...)
    c = stringify(s)
    n = length(x.data)
    ks = keys(x.data)
    ms, d = zip([_reflectinmodel(x.data[k], fm, fa, fsm, fsa, s * encode(i, n), ski, args...)
                  for (i, k) in enumerate(ks)]...) |> collect
    ms = _remap(x.data, ms)
    m = if haskey(fsm, c)
        ArrayModel(fsm[c](sum(d)))
    elseif ski && n == 1
        identity_model()
    else
        _reflectinmodel(ProductModel(ms)(x), fm, fa, fsm, fsa, s, ski, args...)[1]
    end
    m = ProductModel(ms, m)
    m, size(m(x).data, 1)
end

function _reflectinmodel(x::ArrayNode, fm, fa, fsm, fsa, s, ski, ssi, ai)
    c = stringify(s)
    r = size(x.data, 1)
    m = if haskey(fsm, c)
        fsm[c](r)
    elseif ssi && r == 1
        identity
    else
        fm(r)
    end |> ArrayModel
    m = _make_imputing(x.data, m, ai)
    m, size(m(x).data, 1)
end

function _reflectinmodel(ds::LazyNode{Name}, fm, fa, fsm, fsa, s, args...) where Name
    pm, d = _reflectinmodel(unpack2mill(ds), fm, fa, fsm, fsa, s * encode(1, 1), args...)
    LazyModel{Name}(pm), d
end

_make_imputing(x, t::ArrayModel, ai) = _make_imputing(x, t.m, ai) |> ArrayModel
_make_imputing(x, t::Chain, ai) = Chain(_make_imputing(x, t[1], ai), t[2:end]...)

_make_imputing(x::AbstractArray, t::Dense, ai) = ai ? preimputing_dense(t) : t
_make_imputing(x::AbstractArray{Maybe{T}}, t::Dense, ai) where T <: Number = preimputing_dense(t)
_make_imputing(x::AbstractArray, ::typeof(identity), ai) = ai ? preimputing_dense(_identity_dense(x)) : identity
_make_imputing(x::AbstractArray{Maybe{T}}, ::typeof(identity), ai) where T <: Number = preimputing_dense(_identity_dense(x))

_make_imputing(x::NGramMatrix, t::Dense, ai) = ai ? postimputing_dense(t) : t
_make_imputing(x::NGramMatrix{Maybe{T}}, t::Dense, ai) where T <: Sequence = postimputing_dense(t)
_make_imputing(x::NGramMatrix, ::typeof(identity), ai) = ai ? postimputing_dense(_identity_dense(x)) : identity
_make_imputing(x::NGramMatrix{Maybe{T}}, ::typeof(identity), ai) where T <: Sequence = postimputing_dense(_identity_dense(x))

_make_imputing(x::MaybeHotArray, t::Dense, ai) = postimputing_dense(t)
_make_imputing(x::MaybeHotArray, ::typeof(identity), ai) = postimputing_dense(_identity_dense(x))

_identity_dense(x) = Dense(Matrix{Float32}(I, size(x, 1), size(x, 1)), zeros(Float32, size(x, 1)))


text:
################################################################################
# Rectangular
################################################################################
abstract type AbstractRectangular <: Style end

#####
##### General rectangular methods
#####
boxx(x, width) = [x, x + width, x + width, x, x]
boxy(y, height) = [y, y, y + height, y + height, y]

cartesian(::AbstractRectangular, x, y) = (1.5*x, 1.5*y)
polygon(::AbstractRectangular, x, y) = (boxx(x, 1), boxy(y, 1))

function iorules(style::AbstractRectangular)

    offset_skeletons = (
        ( ( 0, 1), "east", "west"),
        ( ( 0,-1), "west", "east"),
    )

    input_rules = squash([
        Offset(a, "out[$i]", "$(c)_in[$i]")
        for (a,b,c) in offset_skeletons, i in 0:style.io_links-1
    ])

    output_rules = squash([
        Offset(a, "$(b)_out[$i]", "in[$i]")
        for (a,b,c) in offset_skeletons, i in 0:style.io_links-1
    ])

    return ConnectionRule(
        vcat(input_rules, output_rules),
        source_filter = filter_io,
        dest_filter = filter_io,
    )
end

function memory_request_rules(::AbstractRectangular)

    # 2 port memories
    proc_to_mem2 = ConnectionRule([
            Offset((1,0), "memory_out", "in[0]"),
            Offset((1,-1), "memory_out", "in[1]")
        ],
        source_filter = filter_memproc,
        dest_filter = filter_memory(2),
    )

    # 1 port memories
    proc_to_mem1 = ConnectionRule(
        [Offset((1,0), "memory_out", "in[0]")],
        source_filter = filter_memproc,
        dest_filter = filter_memory(1),

    )
    return (proc_to_mem1, proc_to_mem2)
end

function memory_return_rules(::AbstractRectangular)

    # 2 port memories
    mem2_to_proc = ConnectionRule([
            Offset((-1,0), "out[0]", "memory_in"),
            Offset((-1,1), "out[1]", "memory_in"),
        ],
        source_filter = filter_memory(2),
        dest_filter = filter_memproc,
    )

    # 1 port memories
    mem1_to_proc = ConnectionRule(
        [Offset((-1,0), "out[0]", "memory_in")],
        source_filter = filter_memory(1),
        dest_filter = filter_memproc,
    )

    return (mem2_to_proc, mem1_to_proc)
end

#####
##### Standard Rectangular
#####
struct Rectangular <: AbstractRectangular
    # Number of interprocessor links
    links :: Int
    io_links :: Int
end

directions(::Rectangular) = ("east", "north", "south", "west")
function procrules(style::Rectangular)

    offset_skeleton = (
        ( (-1, 0), "north", "south"),
        ( ( 1, 0), "south", "north"),
        ( ( 0, 1), "east",  "west"),
        ( ( 0,-1), "west",  "east"),
    )

    offsets = squash([
        Offset(a, "$(b)_out[$i]", "$(c)_in[$i]")
        for (a,b,c) in offset_skeleton, i in 0:style.links
    ])

    return ConnectionRule(
        offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
    )
end

# Port annotations
function port_boundaries(::Rectangular, orientation)
    coords = Dict(
        "east" => ((1,0), (1,1)),
        "south" => ((1,1), (0,1)),
        "west" => ((0,0), (0,1)),
        "north" => ((1,0), (0,0)),
    )

    return coords[orientation]
end

function initial_offset(::Rectangular, orientation, direction)
    # Select where to start between the start and stop coordinates depending
    # on the direction selected
    starts = Dict(
        ( "north", Output ) => 0.55,
        ( "east" , Output ) => 0.55,
        ( "south", Input ) => 0.55,
        ( "west" , Input ) => 0.55,

        ( "north", Input ) => 0.05,
        ( "east" , Input ) => 0.05,
        ( "south", Output ) => 0.05,
        ( "west" , Output ) => 0.05,
    )

    return starts[(orientation, direction)]
end

#####
##### 8-4 rectangular
#####

# Neighbors are still in the four cardinal directions, but now there are long distance
# links also that hop over neighbors
struct Rectangular84 <: AbstractRectangular
    # Number of interprocessor links
    links :: Int
    io_links :: Int
end

# Have two types of links in each direction - one for 1 hop links and another for 2 hop links
directions(::Rectangular84) = [
    "east",
    "east_far",
    "north",
    "north_far",
    "south",
    "south_far",
    "west",
    "west_far",
]
function procrules(style::Rectangular84)

    offset_skeleton = (
        ( (-1, 0), "north", "south"),
        ( ( 1, 0), "south", "north"),
        ( ( 0, 1), "east",  "west"),
        ( ( 0,-1), "west",  "east"),
        ( (-2, 0), "north_far", "south_far"),
        ( ( 2, 0), "south_far", "north_far"),
        ( ( 0, 2), "east_far",  "west_far"),
        ( ( 0,-2), "west_far",  "east_far"),
    )

    offsets = squash([
        Offset(a, "$(b)_out[$i]", "$(c)_in[$i]")
        for (a,b,c) in offset_skeleton, i in 0:style.links
    ])

    return ConnectionRule(
        offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
    )
end

# Port annotations
function port_boundaries(::Rectangular84, orientation)
    # Split the near are far links halfway across the border
    coords = Dict(
        "east"      =>  ((1,0),    (1,0.5)),
        "east_far"  =>  ((1, 0.5), (1, 1)),
        "south"     =>  ((1,1),    (0.5, 1)),
        "south_far" =>  ((0.5, 1), (0, 1)),
        "west"      =>  ((0,0),    (0, 0.5)),
        "west_far"  =>  ((0,0.5),  (0,1)),
        "north"     =>  ((1,0), (0.5,0)),
        "north_far" =>  ((0.5,0), (0,0)),
    )

    return coords[orientation]
end

function initial_offset(::Rectangular84, orientation, direction)
    # Select where to start between the start and stop coordinates depending
    # on the direction selected
    starts = Dict(
        ( "north", Output ) => 0.55,
        ( "east" , Output ) => 0.55,
        ( "north_far", Output ) => 0.55,
        ( "east_far" , Output ) => 0.55,
        ( "south", Input ) => 0.55,
        ( "west" , Input ) => 0.55,
        ( "south_far", Input ) => 0.55,
        ( "west_far" , Input ) => 0.55,

        ( "north", Input ) => 0.05,
        ( "east" , Input ) => 0.05,
        ( "north_far", Input ) => 0.05,
        ( "east_far" , Input ) => 0.05,
        ( "south", Output ) => 0.05,
        ( "west" , Output ) => 0.05,
        ( "south_far", Output ) => 0.05,
        ( "west_far" , Output ) => 0.05,
    )

    return starts[(orientation, direction)]
end

#####
##### 8-8 rectangular
#####

# 8 nearest neighbors
struct Rectangular88 <: AbstractRectangular
    # Number of interprocessor links
    links :: Int
    io_links :: Int
end

# Have two types of links in each direction - one for 1 hop links and another for 2 hop links
directions(::Rectangular88) = [
    "east",
    "north_east",
    "north",
    "north_west",
    "west",
    "south_west",
    "south",
    "south_east",
]
function procrules(style::Rectangular88)

    offset_skeleton = (
        ( (-1, 0), "north", "south"),
        ( ( 1, 0), "south", "north"),
        ( ( 0, 1), "east",  "west"),
        ( ( 0,-1), "west",  "east"),
        ( (-1, 1), "north_east", "south_west"),
        ( (-1,-1), "north_west", "south_east"),
        ( ( 1, 1), "south_east",  "north_west"),
        ( ( 1,-1), "south_west",  "north_east"),
    )

    offsets = squash([
        Offset(a, "$(b)_out[$i]", "$(c)_in[$i]")
        for (a,b,c) in offset_skeleton, i in 0:style.links
    ])

    return ConnectionRule(
        offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
    )
end

# Port annotations
function port_boundaries(::Rectangular88, orientation)
    # Split the near are far links halfway across the border
    coords = Dict(
        "east"        =>  ((1,0),    (1,0.5)),
        "south_east"  =>  ((1, 0.5), (1, 1)),
        "south"       =>  ((1,1),    (0.5, 1)),
        "south_west"  =>  ((0.5, 1), (0, 1)),
        "north_west"  =>  ((0,0),    (0, 0.5)),
        "west"        =>  ((0,0.5),  (0,1)),
        "north_east"  =>  ((1,0),    (0.5,0)),
        "north"       =>  ((0.5,0),  (0,0)),
    )

    return coords[orientation]
end

function initial_offset(::Rectangular88, orientation, direction)
    # Select where to start between the start and stop coordinates depending
    # on the direction selected
    starts = Dict(
        ( "north", Output ) => 0.55,
        ( "east" , Output ) => 0.55,
        ( "north_east", Output ) => 0.55,
        ( "south_east" , Output ) => 0.55,
        ( "south", Input ) => 0.55,
        ( "west" , Input ) => 0.55,
        ( "south_west", Input ) => 0.55,
        ( "north_west" , Input ) => 0.55,

        ( "north", Input ) => 0.05,
        ( "east" , Input ) => 0.05,
        ( "north_east", Input ) => 0.05,
        ( "south_east" , Input ) => 0.05,
        ( "south", Output ) => 0.05,
        ( "west" , Output ) => 0.05,
        ( "south_west", Output ) => 0.05,
        ( "north_west" , Output ) => 0.05,
    )

    return starts[(orientation, direction)]
end


text:
using BeeEncoder
using Suppressor
using Test

@testset "BeeEncoder.jl" begin
    @testset "simple" begin
        BeeEncoder.reset()

        example="""new_int(w, 0, 10)
        new_int(x, 0, 5)
        new_int(z, -5, 10)
        new_int(y, -4, 9)
        new_bool(x1)
        new_bool(x4)
        new_bool(x2)
        new_bool(x3)
        int_plus(x, y, z)
        bool_eq(x1, -x2)
        bool_eq(x2, true)
        bool_array_sum_eq([-x1, x2, -x3, x4], w)
        solve satisfy
        """
        ret = @capture_out include("../example/simple-example.jl")
        @test ret == example

        @test "new_int(w, 0, 10)\n" == capture_render(w)

        c = x+y == z
        @test "int_plus(x, y, z)\n" == capture_render(c)

        c = xl[1] == -xl[2]
        @test "bool_eq(x1, -x2)\n" == capture_render(c)

        BeeEncoder.reset()
        @constrain xl[1] == -xl[2]
        @test "bool_eq(x1, -x2)\nsolve satisfy\n" == capture_render()

        BeeEncoder.reset()
        constrain(xl[1] == -xl[2])
        @test "bool_eq(x1, -x2)\nsolve satisfy\n" == capture_render()
    end

    @testset "Declaring Variable" begin
        BeeEncoder.reset()

        @beebool x1
        @test "new_bool(x1)\n" == capture_render(x1)

        x2 = beebool("x2")
        @test "new_bool(x2)\n" == capture_render(x2)

        x3 = beebool("x3")
        @test "new_bool(x3)\n" == capture_render(x3)

        @test hasbool("x3")
        @test !hasbool("x4")

        x5 = fetchbool("x5")
        @test isa(x5, BeeBool)
        @test x5.name == "x5"
        @test x5 === fetchbool("x5")

        @test x5 === getbool("x5")
        @test_throws KeyError getbool("x6")

        xx1 = fetchbool("xx")
        xx2 = fetchbool("xx")

        @test isequal(xx1, xx2)

        yl = @beebool y1 y2 y3
        for i in 1:3
            @test "new_bool(y$i)\n" == capture_render(yl[i])
        end

        zl = @beebool z[1:10]
        for i in 1:10
            @test "new_bool(z$i)\n" == capture_render(zl[i])
        end

        z1 = fetchbool("z1")
        @test z1 in zl

        (a, b, cl, d) = @beebool a b c[1:10] d
        for i in 1:10
            @test "new_bool(c$i)\n" == capture_render(cl[i])
        end

        BeeEncoder.reset()

        @beeint xx 3 55

        @test isequal(xx, xx)

        @test "new_int(xx, 3, 55)\n" == capture_render(xx)

        il = @beeint i[1:10] 3 7
        for i in 1:10
            @test "new_int(i$i, 3, 7)\n" == capture_render(il[i])
        end

        @test "[i1, i2, i3]" == capture_render(il[1:3])

        i1 = getint("i1")
        @test i1 in il

        @test hasint("i4")

        @test !hasint("i12")

        @test xx === getint("xx")
        @test_throws KeyError getint("x6")
    end

    @testset "Boolean statements" begin
        BeeEncoder.reset()

        @beebool x1
        @beebool x2
        x3 = BeeBool("x3")

        c = x1 == x2
        @test "bool_eq(x1, x2)\n" == capture_render(c)

        c = -x1 == x2
        @test "bool_eq(-x1, x2)\n" == capture_render(c)

        c = true == and([x1, -x2, x3])
        @test "bool_array_and([x1, -x2, x3])\n" == capture_render(c)

        c = and([x1, x2]) == -x3
        @test "bool_array_and_reif([x1, x2], -x3)\n" == capture_render(c)

        c = -x3 == BeeEncoder.xor(-x1, x2)
        @test "bool_xor_reif(-x1, x2, -x3)\n" == capture_render(c)
    end

    @testset "Integer statements" begin
        BeeEncoder.reset()

        @beeint x1 3 7
        @beeint x2 4 6
        x3 = beeint("x3", 10, 15)

        c = (x1 < x2) == true
        @test "int_lt(x1, x2)\n" == capture_render(c)

        c = true == (x1 < x2)
        @test "int_lt(x1, x2)\n" == capture_render(c)

        c = false == (x1 < x2)
        @test "int_lt_reif(x1, x2, false)\n" == capture_render(c)

        c = (x1 < x2) == false
        @test "int_lt_reif(x1, x2, false)\n" == capture_render(c)

        c = alldiff([x1, x2, x3]) == true
        @test "int_array_allDiff([x1, x2, x3])\n" == capture_render(c)

        c = true == alldiff([x1, x2, x3])
        @test "int_array_allDiff([x1, x2, x3])\n" == capture_render(c)

        c = alldiff([x1, x2]) == x3
        @test "int_array_allDiff_reif([x1, x2], x3)\n" == capture_render(c)

        c = x3 == alldiff([x1, x2])
        @test "int_array_allDiff_reif([x1, x2], x3)\n" == capture_render(c)

        c = x1 * x2 == x3
        @test "int_times(x1, x2, x3)\n" == capture_render(c)

        c = x3 == x1 * x2
        @test "int_times(x1, x2, x3)\n" == capture_render(c)

        c = x1 + x2 == 33
        @test "int_plus(x1, x2, 33)\n" == capture_render(c)

        c = 33 == x1 + x2
        @test "int_plus(x1, x2, 33)\n" == capture_render(c)

        c = max([x1, x2]) == x3
        @test "int_array_max([x1, x2], x3)\n" == capture_render(c)

        c = x3 == max([x1, x2])
        @test "int_array_max([x1, x2], x3)\n" == capture_render(c)
    end

    @testset "Cardinality" begin
        BeeEncoder.reset()

        xl = [beebool("x$i") for i in 1:3]

        il = [beeint("i$i", 3, 5) for i in 1:3]

        c = sum(il) == il[1]
        @test "int_array_sum_eq([i1, i2, i3], i1)\n" == capture_render(c)

        c = sum(xl) >= il[1]
        @test "bool_array_sum_geq([x1, x2, x3], i1)\n" == capture_render(c)

        c = sum(xl) > il[1]
        @test "bool_array_sum_gt([x1, x2, x3], i1)\n" == capture_render(c)
    end

    @testset "Boolean arrays relation" begin
        BeeEncoder.reset()

        xl = [beebool("x$i") for i in 1:3]

        il = [beebool("i$i") for i in 1:3]

        @beebool b

        c = (xl == il) == true
        @test "bool_arrays_eq([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl != il) == true
        @test "bool_arrays_neq([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl <= il) == true
        @test "bool_arrays_lex([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl < il) == true
        @test "bool_arrays_lexLt([x1, x2, x3], [i1, i2, i3])\n" == capture_render(c)

        c = (xl <= il) == b
        @test "bool_arrays_lex_reif([x1, x2, x3], [i1, i2, i3], b)\n" == capture_render(c)

        c = (xl < il) == b
        @test "bool_arrays_lexLt_reif([x1, x2, x3], [i1, i2, i3], b)\n" == capture_render(c)
    end
end

max_stars_repo_name: newptcai/BeeEncoder

max_issues_repo_name: newptcai/BeeEncoder

max_forks_repo_name: newptcai/BeeEncoder


text:
using Bio

macro output_timings(t,name)
    quote
        @printf "BioJulia,%s,%f\n" $name minimum($t)
    end
end

macro timeit(ex,name)
    quote
        t = Float64[]
        for i in 1:6
            e = 1000*(@elapsed $(esc(ex)));
            if i > 1
                # warm up on first iteration
                push!(t, e)
            end
        end
        @output_timings t $name
    end
end

max_stars_repo_name: diegozea/mitos-benchmarks

max_issues_repo_name: diegozea/mitos-benchmarks

max_forks_repo_name: diegozea/mitos-benchmarks


text:
x = randn(16,5)
wt = wavelet(WT.haar)
xw = cat([wpd(x[:,i], wt) for i in axes(x,2)]..., dims=3)
xsw = cat([swpd(x[:,i], wt) for i in axes(x,2)]..., dims=3)
xacw = cat([acwpd(x[:,i], wt) for i in axes(x,2)]..., dims=3)

# bb
@test isvalidtree(x[:,1], bestbasistree(xw[:,:,1], BB()))
@test isvalidtree(x[:,1], bestbasistree(xw, method=BB())[:,1])
@test isvalidtree(x[:,1],
    bestbasistree(xw, BB(LogEnergyEntropyCost(), false))[:,1])
@test isvalidtree(x[:,1], bestbasistree(xsw, BB(redundant=true))[:,1])
@test isvalidtree(x[:,1], bestbasistree(xacw, BB(redundant=true))[:,1])

# jbb
@test isvalidtree(x[:,1], bestbasistree(xw))
@test isvalidtree(x[:,1], bestbasistree(xw, JBB(NormCost(), false)))
@test isvalidtree(x[:,1], bestbasistree(xsw, JBB(redundant=true)))

# lsdb
@test isvalidtree(x[:,1], bestbasistree(xw, method=LSDB()))
@test isvalidtree(x[:,1], bestbasistree(xw, LSDB()))
@test isvalidtree(x[:,1], bestbasistree(xsw, LSDB(redundant=true)))

# siwpd_bb
x = randn(16)
xw0 = siwpd(x, wt)
xw4 = siwpd(circshift(x,4), wt)
bt0 = map(node -> sum(node), bestbasistree(xw0, 4, SIBB()))
bt4 = map(node -> sum(node), bestbasistree(xw4, 4, SIBB()))
@test bt0 == bt4

# misc
@test_throws ArgumentError bestbasis_treeselection(randn(15), 8, :fail)
@test_throws AssertionError bestbasis_treeselection(randn(7), 3, :fail) # true n=4


text:
function plot_example_track_results(trackresults; f_if=0, M=10,
                                    save_to_file=missing,
                                    code_length=1023,
                                    top=0.9,
                                    bottom=0.1,
                                    left=0.1,
                                    right=0.9,
                                    wspace=0.1,
                                    hspace=0.1,
                                    marker_size=1,
                                    line_width=1,
                                    use_ylims=true,
                                    truth_snr=missing,
                                    truth_doppler=missing,
                                    truth_doppler_rate=missing,
                                    truth_t=missing,
                                    truth_phi=missing,
                                    truth_code_phase=missing)
    t = trackresults.t
    ZP = trackresults.ZP
    SNR = trackresults.SNR
    x = trackresults.x
    P = trackresults.P
    ϕ = x[1,:] ./ 2π #.* (180/π)
    σ_ϕ = sqrt.(P[1,:]) ./ 2π # .* (180/π)
    f_d = x[2,:] ./ (2π) .- f_if
    σ_f_d = sqrt.(P[2,:]) ./ (2π)
    f_dd = x[3,:] ./ (2π)
    σ_f_dd = sqrt.(P[3,:]) ./ (2π)
    dphi_meas = trackresults.dphi_meas
    dphi_filt = trackresults.dphi_filt
    code_err_meas = trackresults.code_err_meas
    code_err_filt = trackresults.code_err_filt
    code_phase_meas = trackresults.code_phase_meas
    code_phase_filt = trackresults.code_phase_filt
    N_over_2 = floor(Int, length(t)/2)
    fig = figure(figsize=(8, 7.5))
    ax1 = fig.add_subplot(4,2,1)
    ax1.plot(t, real.(ZP), label="real(ZP)", linewidth=line_width)
    ax1.plot(t, imag.(ZP), label="imag(ZP)", linewidth=line_width)
    xlabel("Time (seconds)")
    ylabel("ZP")
    legend()
    title("(a)")
    ax2 = fig.add_subplot(4,2,2)
    if ~ismissing(truth_snr)
        ax2.plot(t, SNR, "k.", markersize=marker_size, label="Measured SNR")
        ax2.axhline(y=truth_snr, c="blue", linestyle="-", label="Truth SNR",
                    linewidth=line_width)
        legend()
    else
        ax2.plot(t, SNR, "k.", markersize=marker_size)
    end
    xlabel("Time (seconds)")
    ylabel("SNR (dB)")
    title("(b)")
    ax3 = fig.add_subplot(4,2,3)
    if ~ismissing(truth_phi)
        phi_diff = truth_phi./2π .- ϕ
        phi_diff = phi_diff .- mean(phi_diff)
        ax3.plot(t,  phi_diff, "k-", label=L"ϕ - \hat{\phi}", linewidth=line_width)
        ax3.plot(t, -σ_ϕ, label=L"\pm\sigma_{\phi}",
                 color="grey", linestyle=":", linewidth=line_width)
        ax3.plot(t, σ_ϕ, color="grey", linestyle=":", linewidth=line_width)
        ylabel(string(L"ϕ - \hat{\phi}", " (cycles)"))
    else
        ax3.plot(t, ϕ, "k-", label=L"\hat{\phi}", linewidth=line_width)
        ax3.plot(t, (ϕ .- σ_ϕ), label=L"\pm\sigma_{\phi}",
                 color="grey", linestyle=":", linewidth=line_width)
        ax3.plot(t, (ϕ .+ σ_ϕ), color="grey", linestyle=":", linewidth=line_width)
        ylabel(string(L"\hat{\phi}", " (cycles)"))
    end
    xlabel("Time (seconds)")
    legend()
    title("(c)")
    ax4 = fig.add_subplot(4,2,4)
    ax4.plot(t, dphi_meas .* (180/π), "k.", label=L"\Delta\phi_{measured}", markersize=marker_size)
    ax4.plot(t, dphi_filt .* (180/π), "b-", label=L"\Delta\phi_{filtered}", linewidth=line_width)
    xlabel("Time (seconds)")
    ylabel(string(L"\Delta\phi", " (degrees)"))
    legend()
    title("(d)")
    ax5 = fig.add_subplot(4,2,5)
    ax5.plot(t, f_d, "k-", label=L"\hat{f}_d", linewidth=line_width)
    ax5.plot(t, f_d .- σ_f_d, label=L"\pm\sigma_{f_d}",
             color="grey", linestyle=":", linewidth=line_width)
    ax5.plot(t, f_d .+ σ_f_d, color="grey", linestyle=":", linewidth=line_width)
    if ~ismissing(truth_doppler)
        ax5.plot(truth_t, truth_doppler, "b-", linewidth=line_width,
                 label=string("Truth ", L"f_d"))
    end
    if use_ylims
        ylim([median(f_d .- σ_f_d.*8), median(f_d .+ σ_f_d.*8)])
    end
    xlabel("Time (seconds)")
    ylabel(string(L"\hat{f}_d", " (Hz)"))
    legend()
    title("(e)")
    ax6 = fig.add_subplot(4,2,6)
    ax6.plot(t, code_err_meas, "k.", label=L"\Delta n_{measured}", markersize=marker_size)
    ax6.plot(t, code_err_filt, "b-", label=L"\Delta n_{filtered}", linewidth=line_width)
    xlabel("Time (seconds)")
    ylabel("Code Phase Error (chips)")
    legend()
    title("(f)")
    ax7 = fig.add_subplot(4,2,7)
    ax7.plot(t, f_dd, "k-", label=L"\hat{\dot{f}}_d", linewidth=line_width)
    ax7.plot(t, f_dd .- σ_f_dd, label=L"\pm\sigma_{\dot{f}_d}",
             color="grey", linestyle=":", linewidth=line_width)
    ax7.plot(t, f_dd .+ σ_f_dd, color="grey", linestyle=":", linewidth=line_width)
    if ~ismissing(truth_doppler_rate)
        ax7.plot(truth_t, truth_doppler_rate, "b-", linewidth=line_width,
                 label=string("Truth ", L"\dot{f}_d"))
    end
    if use_ylims
        ylim([median(f_dd .- σ_f_dd.*4), median(f_dd .+ σ_f_dd.*4)])
    end
    # ylim([-408, 88])
    xlabel("Time (seconds)")
    ylabel(string(L"\hat{\dot{f}}_d", " (Hz/s)"))
    legend()
    title("(g)")
    ax8 = fig.add_subplot(4,2,8)
    ax8.plot(t, code_phase_filt .% code_length, "k-", label=L"\hat{n}", linewidth=line_width)
    intervals = floor(Int, length(t)/M)
    t_sigma_n = Array{Float64}(undef, intervals+1)
    sigma_n = Array{Float64}(undef, intervals+1)
    n = Array{Float64}(undef, intervals+1)
    for i in 1:intervals
        sigma_n[i+1] = std(code_err_meas[(i-1)*M+1:i*M])
        # n[i] = mean(code_phase_filt[(i-1)*M+1:i*M])
        # n[i] = code_phase_filt[(i-1)*M+floor(Int, M/2)]
        n[i+1] = code_phase_filt[i*M]
        t_sigma_n[i+1] = t[i*M]
        if i == 1
            t_sigma_n[1] = t[1]
            n[1] = code_phase_filt[1]
            sigma_n[1] = sigma_n[i+1]
        end
    end
    ax8.plot(t_sigma_n, (n .- sigma_n) .% code_length, label=L"\pm\sigma_n",
             color="grey", linestyle=":", linewidth=line_width)
    ax8.plot(t_sigma_n, (n .+ sigma_n) .% code_length, color="grey", linestyle=":",
             linewidth=line_width)
    if ~ismissing(truth_code_phase)
        ax8.plot(truth_t, truth_code_phase .% code_length, "b-", linewidth=line_width,
                 label=string("Truth ", L"n"))
    end
    xlabel("Time (seconds)")
    ylabel("Code Phase (chips)")
    legend()
    title("(h)")
    subplots_adjust(hspace=hspace, wspace=wspace,  top=top, left=left, right=right,
                    bottom=bottom)
    if ~ismissing(save_to_file)
        savefig(save_to_file, dpi=300)
    end
end


text:
using LinearAlgebra: I, norm
# using CairoMakie

export euclideandistance, inductiveVR, edgelist, H0Barcode, prunedataset

"""
euclideandistance(dataset)

Computes the Euclidean distance matrix for a set of points.
This corresponds to the weighted neighborhood graph when
the maximum allowable distance is greater than the maximum
distance between points.

### Arguments
 - `dataset` -- An N×M array containing N points in Rᴹ

### Returns
An N×N symmetric matrix of pairwise distances between points.
"""
function euclideandistance(dataset)
    N = size(dataset)[1]
    distances = zeros(N,N)
    for i in 1:N
        for j in i+1:N
            distances[i,j] = norm(dataset[i,:] - dataset[j,:])
            distances[j,i] = distances[i,j]
        end
    end
    return distances
end

"""
lowerneighbors(k, distances)

Returns a list of the neighboring nodes which are lower in ordering.

### Arguments
 - `k`         -- Chosen node
 - `distances` -- Matrix with zeros when nodes are not connected

### Returns
An array of node neighbors with ``index<k``
"""
function lowerneighbors(k, distances)
    return findall(!iszero, distances[k,1:k-1])
end

"""
    edgelist(distances)

Generates a list of graph edges from a weighted adjacency matrix.

### Arguments
 - `distances` --   Adjacency matrix of a graph where connections are
                    represented by the distance between the nodes.

### Returns
`(simplices, weights)`
 - `simplices`  -- 2×M matrix of M edges in the graph
 - `weights`    -- M dimensional vector of distances
"""
function edgelist(distances)
    N = size(distances)[1]
    simplices = zeros(Int,2,N^2)
    weights = zeros(N^2)
    ind = 1
    for i in 1:N-1
        for j in i+1:N
            if !iszero(distances[i,j])
                simplices[1, ind] = i
                simplices[2, ind] = j
                weights[ind] = distances[i,j]
                ind += 1
            end
        end
    end
    return simplices[:,1:ind-1], weights[1:ind-1]
end

"""
    neighborintersect(simplex,adjacency)

Computes the interesection of the neighboring nodes in a simplex.

### Arguments
 - `simplex`   -- An array of nodes
 - `adjacency` -- Adjacency matrix of a graph

### Returns
A list of nodes adjacent to all nodes

### Notes
Any non-zero elements in `adjacency` are interpretted as 1's
"""
function neighborintersect(simplex,adjacency)
    N = size(adjacency)[1]
    neighbors = [l for l in 1:N]
    for l in 1:size(simplex)[1]
        intersect!(neighbors,lowerneighbors(simplex[l],adjacency))
    end
    return neighbors
end

"""
    inductiveVR(distances, k)

Computes the Vietoris-Rips complex using the inductive algorithm.
Builds the simplices starting from 1-simplices and forming higher
dimensional faces.

### Arguments
 - `distances` -- Matrix with zeros when nodes are not connected
 - `k`         -- Maximum simplice dimension

### Returns
`(simplices, weights)`
 - `simplices`  -- An array of simplicial matrices up to dimension `k`
 - `weights`    -- Weights associated with `simplices`

### Notes
This is the most basic construction mentioned in the paper.

### Reference

Afra Zomorodian,
Fast construction of the Vietoris-Rips complex,
Computers & Graphics,
Volume 34, Issue 3,
2010,
Pages 263-271,
ISSN 0097-8493,
https://doi.org/10.1016/j.cag.2010.03.007.
"""
function inductiveVR(distances, k)
    N = size(distances)[1]

    # 1-Simplices
    edges, weights = edgelist(distances)
    V = [edges]
    W = [weights]

    # Iteratively expand the complex
    for i in 1:k
        simplices = zeros(Int,i+2,0)
        weights = zeros(0);

        simplices_t = [zeros(Int,i+2,0) for j in 1:Threads.nthreads()]
        weights_t = [zeros(0) for j in 1:Threads.nthreads()]

        # For each i-simplex
        Threads.@threads for j in 1:size(V[end])[2]
            tnum = Threads.threadid()
            simplex = V[end][:,j]

            # Compute interaction of neighbors
            neighbors = neighborintersect(simplex,distances)

            # Append each element of the intersection into the simplices
            for v in neighbors
                w = W[end][j]
                for node in simplex
                    w = max(w,distances[v,node])
                end
                push!(weights_t[tnum],w)
                new_simplice = sort([v, simplex...])
                simplices_t[tnum] = hcat(simplices_t[tnum], new_simplice)
            end
        end

        @inbounds for j in 1:Threads.nthreads()
            simplices = hcat(simplices,simplices_t[j])
            append!(weights,weights_t[j])
        end

        # Append the next level of simplices.
        push!(V,simplices)
        push!(W,weights)
    end
    return V, W
end

"""
    H0Barcode(distances)

Compute the H0 persistent homology from a euclidean distance matrix

### Arguments
 - `distances` -- A Euclidean Distance Matrix

### Returns
Sequence of H0 homology lifetimes
"""
function H0Barcode(distances)

    # Initialize the arrays
    n = size(distances)[1]
    persistence = zeros(2,n)
    dist_sort = zeros(n*(n-1) ÷ 2)
    indices = zeros(Int,2,n*(n-1) ÷ 2)
    rips_complex = Array(1:n)

    # Sort the distances for Birth-Death Process
    p = sortperm(distances[:])[n+2:2:length(distances)]
    i = 1
    @inbounds for ind in p
        a, b = divrem(ind-1,n)
        indices[1,i] = a+1
        indices[2,i] = b+1
        dist_sort[i] = distances[ind]
        i += 1
    end

    # Step through the sorted distances to record death process
    @inbounds for i in 1:length(dist_sort)
        # remove later homology class, copy into first
        hom1 = rips_complex[indices[1,i]]
        hom2 = rips_complex[indices[2,i]]

        if hom1 != hom2
            persistence[2,hom2] = dist_sort[i]
            for j in 1:n
                if rips_complex[j] == hom2
                    rips_complex[j] = hom1
                end
            end
        end
    end

    return persistence[:,1:end-1]
end

"""
    prunedataset(distances, ϵ)

Prunes the dataset based on a Euclidean distance matrix.
Given a point, by the triangle inequality, any two points within
distance ``ϵ`` of the first point are within ``2ϵ`` of each other.

Beginning with the first element of the dataset, iterate up removing
all later elements within ``ϵ``.
Should only change the persistent homology by up to ``ϵ``

### Arguments
 - `distances`  -- Euclidean Distance Matrix
 - `ϵ`          -- Pruning Distance
### Returns
Two values, `pruned_dist, remaining_indices`

 - `pruned_dist`        -- Distance matrix of remaining indices
 - `remaining_indices`  -- Remaining indices
"""
function prunedataset(distances, ϵ)
    N = size(distances)[1]

    # Find initial connections
    removals = 1 .+ findall(distances[1,2:end] .< ϵ)

    # For each of the later nodes
    for i in 2:N

        # If the node hasn't been removed yet, remove all close nodes
        if !(i in removals)
            union!(removals, i .+ findall(distances[i,i+1:end] .< ϵ))
        end
    end

    remaining_indices = setdiff(1:N,removals)
    pruned_dist = distances[remaining_indices,remaining_indices]

    return pruned_dist, remaining_indices
end


text:
reload("PyTest")
include("different_module.jl")

module tst

using PyTest

using Base.Test

@testset "PyTest tests" begin
  # FIXME: avoid this
  using ResumableFunctions

  include("base_fixtures.jl")
  include("exceptions.jl")
  include("request.jl")
  include("parametrized.jl")
  include("basetest_integration.jl")
  include("builtin/tempdir.jl")
  include("runner.jl")
end

end


text:
module PixelArts

"""
    display_js(s::String)
# Details
Sends to display in IJulia or (to be implemented) to blink electron window in Atom
# Arguments
* `s::String` the javascript to be executed
* `require` an array of strings specifying the js libraries to use, assuming already loaded, defaults to d3
# Example
```julia
display_js("console.log('hello world');")
```
"""
function display_js(s::String, require = ["d3"])
    display("text/html",
    """
    <script>
        require(['$(join(require, "'], ['"))'],  function($(join(require, ","))) {
        $s
        });
    </script>
    """)
end
export display_js

"""
    create_canvas(id::String, vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)
    create_canvas(vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)

# Details
To start a canvas called "id" where the `html svg` image will be launched use

    create_canvas(vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)

If `id` is not provided, then the function creates and returns a random id of the form `id = "canvas" * randstring(3)` using

    create_canvas(vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)

# Arguments
* `id::String` `html` id of the `svg` tag for the canvas.
* `vunits::Int` number of vertical measurement units in the canvas
* `hunits::Int` number of horizontal measurement units in the canvas
* `height::Int` the height in display pixels of the image, the `svg` image is scalated to this height
* `width::Int` the width in display pixels of the image, the `svg` image is scalated to this width
# Examples
```julia
create_canvas("canvas", 10, 10)
canvasid = create_canvas(10, 10)
```
"""
function create_canvas(id::String, vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)
    display("text/html", """
    <div class="usersvg">
    <svg id='$id' width='$(width)' height='$(height)' viewbox='0 0 $hunits $vunits'></svg>
    </div>
    <script>
    requirejs.config({paths: {d3: ["https://d3js.org/d3.v3.min.js?noext", "?$(@__DIR__)/../assets/d3.min"] }});
    </script>
    """)
end
function create_canvas(vunits::Int, hunits::Int, height::Int = 250, width::Int = 250)
    id = "canvas" * randstring(3)
    create_canvas(id, vunits, hunits, height, width)
    return id
end
export create_canvas

"""
    remove_element_by_id(id::String)
# Details
All elements created using PixelArts have unique `html` ids. This functions eliminates an id from the current canvas.
# Arguments
* `id::String` the id to be removed from the canvas (the pages's DOM)
# Examples
```julia
create_canvas("canvas", 10, 10)
remove_element_by_id("canvas")
canvasid = create_canvas(10, 10)
remove_element_by_id(canvasid)
```
"""
function remove_element_by_id(id::String)
    display("text/html", """
    <script>
    var x = document.querySelector(".usersvg #$id");
    if (!(x == null)) {
        x.parentNode.removeChild(x);
    }
    """)
end
export remove_element_by_id

"""
    add_pixels(id::String, i::Int, j::Int, colour::String = "black"; <keywords>)
    add_pixels(id::String, pos::Array{Array{Int, 1}, 1}, colour::Any = ["black"]; <keywords>)
    add_pixels(id::String, colour_array::Array; <keywords>)

# Details
Use when adding a single pixel

    add_pixels(id::String, i::Int, j::Int, colour::String = "black"; <keywords>)

# Arguments
* `id::String` the canvas id where the pixels will be added
* `i::Int` the row to place the pixel
* `j::Int` the col to place the pixel
* `colour::String` Any string representing a color interpretable by javascript (color name, rgb, hexadecimal)
* `; <keywords>` pending
# Examples
```julia
cv = create_canvas(10, 10)
add_pixels(cv, 1, 2) # added black pixels in position [1,2]
add_pixels(cv, 2, 2, "yellow")
```
"""
function add_pixels(id::String, i::Int, j::Int, colour::Any = "black"; disp = true)
    s = """
    d3.select(".usersvg #$id")
      .append("rect")
          .attr("width", "1")
          .attr("height", "1")
          .attr("transform", "translate($(j-1) $(i-1))")
          .style("fill", "$colour");
    """
    disp ? display_js(s) : return s
end
"""
# Details
Use when adding many pixels and passing an array of paris of coordiantes

    add_pixels(id::String, pos::Array{Array{Int, 1}, 1}, colour::Any = ["black"])

# Arguments
* `id::String` the canvas id where the pixels will be added
* `pos::Array{Array{Int, 1}, 1}`an array of entries `[[i1, j1], [i2,j2],...]`
* `<keywords>` dictionary of attr and style to be passed to svg elements
# Examples
```julia
cv = create_canvas(10, 10)
add_pixels(cv, [[1,2], [2, 4]]) # added black pixels in those positions
```
"""
function add_pixels(id::String, pos::Array{Array{Int, 1}, 1}; attr = Dict("fill" => "black"), style = Dict(), disp = true)
    s = """var svg = d3.select(".usersvg #$id");"""
    for p in pos
        i, j = p
        s *= """svg.append("rect").attr("width", "1").attr("height", "1").attr("transform", "translate($(j-1) $(i-1))")"""
        for key in keys(attr) s*= """.attr("$key", "$(attr[key])")""" end
        for key in keys(style) s*= """.attr("$key", "$(style[key])")""" end
        s *= ";"
    end
    disp ? display_js(s) : return s
end
"""
# Details
Use when passing a square array of colors representing the greed

    add_pixels(id::String, pos::Array{Array{Int, 1}, 1}, colour::Any = ["black"])

# Arguments
* `id::String` the canvas id where the pixels will be added
* `pos::Array{Array{Int, 1}, 1}`an array of entries `[[i1, j1], [i2,j2],...]`
* `colour::Any` An array strings of javascript interpretable colors matching the length of `pos` or with length one, in which case the same colous is used for every pixel.
# Examples
```julia
colour_array = ["green" "red";"yellow" "blue"]
cv = create_canvas(size(colour_array)...)
add_pixels(cv, colour_array) # added black pixels in those positions
```
"""
function add_pixels(id::String, colour_array::Array; disp = true)
    vunits, hunits = size(colour_array)
    pos = [[i,j] for i in 1:vunits for j in 1:hunits]
    s = ""
    for p in pos
        s *= add_pixels(id, p..., colour_array[p...]; disp = false)
    end
    display_js(s)
end
export add_pixels

"""
    render_bg(colour_array::Array, width::Int = 250, height::Int = 250)
    render_bg(array::Any, colour_dict::Dict, width::Int = 250, height::Int = 250)
# Details
Uses `create_canvas` and `add pixels` to convert a rectangular array of `html` interpretable colors in an image.

    render_bg(colour_array::Array, width::Int = 250, height::Int = 250)
# Arguments
* `colour_array::Array` an array of strings or elements that can be interpreted as html colors when converted to text
* ``height::Int`, width::Int` size in screen pixels of the image
# Value
* an id representing the created canvas
# Examples
```julia
render_bg(["green" "red";"yellow" "blue"])
```
"""
function render_bg(colour_array::Array, height::Int = 250, width::Int = 250)
    vunits, hunits = size(colour_array)
    pos = [[i,j] for i in 1:vunits for j in 1:hunits]
    colour = [colour_array[x...] for x in pos]
    new_canvas = create_canvas(vunits, hunits, height, width)
    add_pixels(new_canvas, pos, colour)
    return new_canvas
end
"""
# Details
An array of any type can be passed with a dictionary to convert the array to colors

    render_bg(array::Array, colour_dict::Dict, width::Int = 250, height::Int = 250)
# Arguments
* `array::Array` any rectangular array
* `colour_dict::Dict` a dictionary mapping the elements of render_bg to `html` interpretable colours
* ``height::Int`, width::Int` size in screen pixels of the canvas
# Value
* an id representing the created canvas
# Examples
```julia
render_bg([1 "#";3 "*"], Dict(1 => "green", "#" => "red", 3 => "yellow", "*" => "blue"))
```
"""
function render_bg(array::Array, colour_dict::Dict, width::Int = 250, height::Int = 250)
    vunits, hunits = size(array)
    pos = [[i,j] for i in 1:vunits for j in 1:hunits]
    colour = [colour_dict[array[x...]] for x in pos]
    new_canvas = create_canvas(vunits, hunits, height, width)
    add_pixels(new_canvas, pos, colour)
    return new_canvas
end
export render_bg

#function render_bg(bg_name::String, width::Int = 250, height::Int = 250)
#    if !(bg_name in ["right_turn", "circuit"])
#        error("Invalid bg_name")
#    end
#    array = readcsv(joinpath(dirname(@__FILE__), "files", "$(bg_name).csv"))
#    colour_dict = readcsv(joinpath(dirname(@__FILE__), "files", "$(bg_name)_colours.csv"))
#    colour_dict = Dict(zip(colour_dict[:,1], colour_dict[:,2]))
#    return render_bg(array, colour_dict, width, height)
#end

#function bg_colour_array(bg_name::String)
#    if !(bg_name in ["right_turn", "circuit"])
#        error("Invalid bg_name")
#    end
#    array = readcsv(joinpath(dirname(@__FILE__), "files", "$(bg_name).csv"))
#    colour_dict = readcsv(joinpath(dirname(@__FILE__), "files", "$(bg_name)_colours.csv"))
#    colour_dict = Dict(zip(colour_dict[:,1], colour_dict[:,2]))
#    vunits, hunits = size(array)
#    pos = [[i,j] for i in 1:vunits, j in 1:hunits]
#    colour_array = map(x -> colour_dict[array[x...]], pos)
#    return colour_array, pos, array
#end
#export bg_colour_array

function add_svg_element(element_id::String, canvas_id::String, element_type::String, i::Real = 1, j::Real = 1; attr = Dict(), style = Dict(), disp = true)
    # if !(element_type in ["circle", "rect", "path", "line", "image"]) error("Element type $(element_type) not supported") end
    s = """d3.select(".usersvg #$canvas_id").append("$element_type")"""
    s *= """.attr("id", "$element_id").attr("transform", "translate($(j-1) $(i-1))")"""
    for key in keys(attr) s*= """.attr("$key", "$(attr[key])")""" end
    for key in keys(style) s*= """.attr("$key", "$(style[key])")""" end
    s *= ";"
    disp ? display_js(s) : return s
end
function add_svg_element(canvas_id::String, element_type::String, i::Real = 0, j::Real = 0; attr = Dict(), style = Dict(), disp = true)
    element_id = "element" * randstring(3)
    return element_id, add_svg_element(element_id, canvas_id, element_type, i, j, attr = attr, style = style, disp = disp)
end
export add_svg_element

function set_attr(element_id::String, attr)
    s = """d3.select(".usersvg #$element_id")"""
    for key in keys(attr) s*= """.attr("$key", "$(attr[key])")""" end
    s *= ";"
    display_js(s)
end
function set_style(element_id::String, style)
    s = """d3.select(".usersvg #$element_id")"""
    for key in keys(style) s*= """.attr("$key", "$(attr[style])")""" end
    s *= ";"
    display_js(s)
end
function translate_element(element_id::String, vmove::Int, hmove::Int)
    s = """d3.select(".usersvg #$element_id").attr("transform", "translate($(hmove-1) $(vmove-1))");"""
    display_js(s)
end`
export set_attr, set_style, translate_element

function add_pixel_cross(element_id::String, canvas_id::String, i::Int, j::Int, colour::Any = "black"; disp = true)``
    attr = Dict("d" => "M0.5 0 L0.5 1 M0 0.5 L1 0.5", "transform" => "translate($(j-1) $(i-1))")
    style = Dict("stroke" => "$colour", "stroke-width" => 0.2)
    add_svg_element(element_id, canvas_id, "path", i, j, attr = attr, style = style, disp = disp)
end
function add_pixel_cross(canvas_id::String, i::Int, j::Int, colour::Any; disp = true)
    element_id = "pixel_cross" * randstring(3)
    add_pixel_cross(element_id, canvas_id, i, j, colour, disp = disp)
    return element_id
end
export add_pixel_cross

function add_pixel_arrow(element_id::String, canvas_id::String, i::Int, j::Int, rotate::Real = 0, colour::Any = "black"; disp = true)
    attr = Dict("d" => "M0 0 L0 0.5 M -0.15 0.3 L0 0.5 L0.15 0.3", "transform" => "translate($(j-0.5) $(i-0.5)) rotate($(-rotate-90))")
    style = Dict("stroke" => "$colour", "stroke-width" => 0.1)
    add_svg_element(element_id, canvas_id, "path", i, j, attr = attr, style = style, disp = disp)
end
function add_pixel_arrow(canvas_id::String, i::Int, j::Int, rotate::Real = 0, colour::Any = "black"; disp = true)
    element_id = "pixel_arrow" * randstring(3)
    add_pixel_arrow(element_id, canvas_id, i, j, rotate, colour, disp = disp)
    return element_id
end
export add_pixel_arrow

"""
    add_pixel_image(element_id::String, canvas_id::String, i::Int, j::Int, path::String; <keywords>)
    add_pixel_image(canvas_id::String, i::Int, j::Int, path::String; <keywords>)
# Details
Creates an image of unit size in a canvas.
# Arguments
* `element_id::String` the `html` id of the new image
* `canvas_id::String` the canvas in which the image will be placed
* `i::Int, j::Int` row and column to place the image of unit size
* `path` path to the image
* `<keywords>` modifiers to be passed to `add_svg_element`
# Value
If `element_id` is not passed to the function, then a random id is returned.
# Example
```julia
cv = render_bg(["green" "red";"yellow" "blue"])
add_pixel_image(cv, 1, 2, "<path>")
```
"""
function add_pixel_image(element_id::String, canvas_id::String, i::Int, j::Int, path::String; attr = Dict(), style = Dict(), disp = true)
    attr = merge(Dict(
        "xlink:href" => path,
        "x" => "0",
        "y" => "0",
        "width" => "1px",
        "height" => "1px"
    ), attr)
    add_svg_element(element_id, canvas_id, "image", i, j, attr = attr, style = style, disp = disp)
end
function add_pixel_image(canvas_id::String, i::Int, j::Int, path::String; attr = Dict(), style = Dict(), disp = true)
    element_id = "pixel_image" * randstring(3)
    add_pixel_image(element_id, canvas_id, i, j, path, attr = attr, style = style, disp = disp)
    return element_id
end
export add_pixel_image

end


text:
#
# # Example with VAR as in Sims (1980)
#

using VARmodels, TimeSeries

data = readdlm("/Users/tomaskrehlik/.julia/v0.3/VARmodels/src/simsData.txt", ';')

d = [Date(1959,1,1):Month(3):Date(2014,7,1)]
names = convert(Vector{String},vec(data[1,2:end]))

data = TimeArray(d, convert(Matrix{Float64},data[46:end,2:end]), names)

s = varEstimate(data["wages","imports","cpi"], 3, "Const")

Series: ....
Lags: ....
Type: ....

No restrictions imposed. Estimated by OLS.


text:
module TableTraits

using IteratorInterfaceExtensions

export isiterabletable
export supports_get_columns_copy, get_columns_copy
export supports_get_columns_copy_using_missing, get_columns_copy_using_missing
export supports_get_columns_view, get_columns_view

# Iterable table trait

function isiterabletable(x::T) where {T}
    isiterable(x) || return false

    if Base.IteratorEltype(x)==Base.HasEltype()
        et = Base.eltype(x)
        if et <: NamedTuple
            return true
        elseif et===Any
            return missing
        else
            return false
        end
    else
        return missing
    end
end

# Column copy trait

supports_get_columns_copy(x::T) where {T} = false

function get_columns_copy end

# Column copy trait using Missing

supports_get_columns_copy_using_missing(x::T) where {T} = false

function get_columns_copy_using_missing end

# Column view trait

supports_get_columns_view(x::T) where {T} = false

function get_columns_view end

end # module


text:
export SuffixArray

struct SuffixArray
    data::Vector
    index::Vector{Int}
    lcps::Vector{Int}
end

function SuffixArray(data::String)
    index = sais(data)
    lcps = lcparray(data, index)
end

function Base.findall(sa::SuffixArray, query::String)
    chars = Vector{Char}(query)
    i = 1
    j = length(sa)
    while i <= j
        k = (i+j) ÷ 2
        ii = sa.index[k]
        while sa.data[ii+m] == chars[m]
            m == length(chars) && return ii
            m += 1
        end

    end

end


text:
module MaternQ
using SparseArrays: SparseMatrixCSC
using LinearAlgebra: Symmetric
using OffsetArrays
export Qmatrix_matern2d

include("matern_values.jl")

""" Directly construct a SparseMatrixCSC with a 2d nearest neighbour structure used by
the sparse Matern precision matrix.
"""
function matern_grid2_CSC(nrow, ncol, nu::Int, f)
    colptr = fill(0, nrow*ncol+1)
    # Upper bound
    nnz = cld((nu*2+3)^2,2)*nrow*ncol
    ii = 1
    col = 1 # result matrix column index
    rowval = Vector{Int}(undef, nnz)
    nzval = Vector{Float64}(undef, nnz)
    li = LinearIndices((1:nrow, 1:ncol))
    for j = 1:ncol
        for i = 1:nrow
            colptr[col] = ii
            col += 1
            for l = -(nu+1):(nu+1)
                for k = -(nu+1):(nu+1)
                    if checkindex(Bool, 1:nrow, i+k) && checkindex(Bool, 1:ncol, j+l) && abs(k) + abs(l) <= nu + 1
                        @assert ii <= nnz
                        I = CartesianIndex(i, j)
                        J = CartesianIndex(i+k, j+l)
                        @inbounds nzval[ii] = f(I, J)
                        @inbounds rowval[ii] = li[J]
                        ii += 1
                    end
                end
            end
        end
    end
    colptr[end] = ii
    nnz = ii-1
    resize!(nzval, nnz)
    resize!(rowval, nnz)
    SparseMatrixCSC(ncol*nrow, ncol*nrow, colptr, rowval, nzval)
end
import GaussianRandomFields
dense_cache_matrix_size = (0,0)
dense_cache_valid = falses(1000)
dense_cache = Vector{Matrix{Float64}}(undef, 1000)
""" Constuct a dense precision matrix by inverting a dense covariance matrix
"""
function Qmatrix_matern2d_dense_cached(size, prec, range, nu::Int)
    global dense_cache_matrix_size, dense_cache, dense_cache_valid
    rangei = round(Int, range*100)
    range = rangei/100.0
    if dense_cache_matrix_size != size
        fill!(dense_cache_valid, false)
        dense_cache_matrix_size = size
    end
    if dense_cache_valid[rangei]
        Q = dense_cache[rangei]
    else
        covf = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(range/4, nu, σ=1.0))
        C = GaussianRandomFields.apply(covf, 1:size[1], 1:size[2])
        Q = inv(C)
        dense_cache[rangei] = Q
        dense_cache_valid[rangei] = true
    end
    Symmetric(Q * prec)
end
function Qmatrix_matern2d_dense(size, prec, range, nu::Int)
    covf = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(range/4, nu))
    C = GaussianRandomFields.apply(covf, 1:size[1], 1:size[2])
    Symmetric(inv(C) * prec)
end

dense = Ref(false)
function Qmatrix_matern2d(size, prec, range, nu::Int)
    dense[] && return Qmatrix_matern2d_dense_cached(size, prec, range, nu)
    nrow, ncol = size
    Q = matern_grid2_CSC(nrow, ncol, nu, Q_matern2d_fun(prec, range, nu))
    Symmetric(Q)
end

end # module


text:
using Oceananigans

# We'll want the make_vertical_slice_movie function from utils.jl to turn our
# rising thermal bubble output into a movie.
include("utils.jl")

# We'll set up a 2D model with an xz-slice so there's only 1 grid point in y.
Nx, Ny, Nz = 256, 1, 256    # Number of grid points in each dimension.
Lx, Ly, Lz = 2000, 1, 2000  # Domain size (meters).
Nt, Δt = 5000, 10           # Number of time steps, time step size (seconds).

# Set up the model and use an artificially high viscosity ν and diffusivity κ.
model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

# Get location of the cell centers in x and z and reshape them to easily
# broadcast over them when calculating hot_bubble_perturbation.
xC, zC = model.grid.xC, model.grid.zC
xC, zC = reshape(xC, (Nx, 1, 1)), reshape(zC, (1, 1, Nz))

# Set a temperature perturbation with a Gaussian profile located at the center
# of the vertical slice. It roughly corresponds to a background temperature of
# T = 282.99 K and a bubble temperature of T = 283.01 K.
hot_bubble_perturbation = @. 0.01 * exp(-100 * ((xC - Lx/2)^2 + (zC + Lz/2)^2) / (Lx^2 + Lz^2))
model.tracers.T.data .= 282.99 .+ 2 .* reshape(hot_bubble_perturbation, (Nx, Ny, Nz))

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "thermal_bubble_2d_" every 10 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="thermal_bubble_2d_", frequency=10)
push!(model.output_writers, nc_writer)

time_step!(model, Nt, Δt)

make_vertical_slice_movie(model, nc_writer, "T", Nt, Δt, model.eos.T₀)