text:
export reversed

function set_capture(sequence::ReversedString, index::Int, x)
    @warn "check"
    set_capture(sequence.x,index,x)
end

## caveat!
regex_string(x::ReversedString) = regex_escape(x.representation)

export PositiveLookbehind
"""
    PositiveLookbehind(parser)

Parser that succeeds if and only if `parser` succeeds **before cursor**. Consumes no input.
The match is returned.
Useful for checks like "must be preceded by `parser`, don't consume its match".
"""
@auto_hash_equals struct PositiveLookbehind{S,T,P} <: WrappedAssertion{S,T}
    parser::P
    function PositiveLookbehind(p_,reversed_parser=true)
        p = reversed_parser ? reversed(parser(p_)) : parser(p_)
        new{Tuple{Int,state_type(p)},result_type(p),typeof(p)}(p)
    end
end
# result_type(p::Type{PositiveLookbehind{T}}) where T = T
regex_prefix(x::PositiveLookbehind) = "(?<="

export NegativeLookbehind
"""
    NegativeLookbehind(parser)

Parser that succeeds if and only if `parser` does not succeed **before cursor**.  Consumes no input.
`nothing` is returned as match.
Useful for checks like "must not be preceded by `parser`, don't consume its match".

```jldoctest
julia> la=NegativeLookbehind("keep")
re"(?<!keep)"

julia> parse("peek"*la,"peek")
("peek", re"(?<!keep)")
```
"""
@auto_hash_equals struct NegativeLookbehind{P} <: WrappedAssertion{MatchState,NegativeLookbehind{P}}
    parser::P
    function NegativeLookbehind(p_,reversed_parser=true)
        p = reversed_parser ? reversed(parser(p_)) : parser(p_)
        new{typeof(p)}(p)
    end
end
regex_prefix(x::NegativeLookbehind) = "(?<!"

export Lookbehind
"""
    Lookbehind(does_match::Bool, p)

[`PositiveLookbehind`](@ref) if `does_match==true`,
[`NegativeLookbehind`](@ref) otherwise.
"""
function Lookbehind(does_match::Bool, p)
    if does_match
        PositiveLookbehind(p)
    else
        NegativeLookbehind(p)
    end
end
@deprecate look_behind(does_match,p) Lookbehind(does_match, p)

children(x::PositiveLookbehind) =
    children(x.parser)

function iterate_state(t::NegativeLookbehind, str, till, posi, next_i, state::Nothing)
    rseq=reversed(str)
    next_i < 1 && return next_i, MatchState()
    p = reverse_index(rseq,_prevind(str,next_i))
    r = iterate_state(t.parser, rseq, till, p, p, nothing)
    if r === nothing
        next_i,MatchState()
    else
        nothing
    end
end

@inline iterate_state(t::NegativeLookbehind, str, till, posi, next_i, state::MatchState) =
    nothing

@inline iterate_state(t::PositiveLookbehind, str, till, posi, next_i, state::MatchState) =
    nothing

function iterate_state(t::PositiveLookbehind, str, till, posi, next_i, state)
    rseq=reversed(str)
    ri = reverse_index(rseq,_prevind(str,next_i))
    next_i < 1 && return nothing
    r = iterate_state(t.parser, rseq, till, ri, tuple_pos(state,ri), tuple_state(state))
    if r === nothing
        nothing
    else
        next_i, r
    end
end

"""
    Base.get(parser::PositiveLookbehind, sequence, till, after, i, state)

get result of PositiveLookbehind

!!! note
    The result is currently for a `reversed` sequence, and
    you might find it difficult to [`Base.map`](@ref) a lookbehind parser match.
    If you require this functionality please open an issue for discussion.

    Assertions do not consume input, so typically these input chars are parsed/mapped outside of the assertion.

```jldoctest
julia> p = Sequence(!re"a+b", PositiveLookbehind(!re"a+b"))
🗄 Sequence
├─ 🗄 Sequence |> !
│  ├─ a+  |> Repeat
│  └─ b
└─ (?<=🗄) Sequence |> ! |> PositiveLookbehind
   ├─ b
   └─ a+  |> Repeat
::Tuple{SubString{String}, SubString{String}}

julia> p("aaab")
("aaab", "baaa")
```
"""
function Base.get(parser::PositiveLookbehind, sequence, till, after, i, state)
    rseq = reversed(sequence)
    after_ = tuple_pos(state)
    get(parser.parser, rseq, till, after_, reverse_index(rseq,prevind(sequence, i)), tuple_state(state))
end

regex_inner(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    regex_inner(reversed(x.parser))

children(x::Union{PositiveLookbehind,NegativeLookbehind}) =
    reverse(children(x.parser))
reversed(x::CombinedParser) = deepmap_parser(_reversed, x)
_reversed(x::ConstantParser{<:AbstractString}) =
    ConstantParser(reversed(x.parser))
_reversed(x) = x
_reversed(x::AtEnd) = AtStart()
_reversed(x::AtStart) = AtEnd()

_deepmap_parser(::typeof(_reversed),mem::AbstractDict,x::Sequence) =
    Sequence(( deepmap_parser(_reversed,mem,p) for p in reverse(x.parts) )...)

for T in [PositiveLookahead,NegativeLookahead,PositiveLookbehind,NegativeLookbehind]
    eval(quote
             _deepmap_parser(f::Function,mem::AbstractDict,x_::$T,a...; kw...) =
                 $T(deepmap_parser(f,mem,x_.parser,a...; kw...),false)
         end)
end

_deepmap_parser(::typeof(_reversed),mem::AbstractDict,x::NegativeLookbehind) =
    NegativeLookahead(x.parser) ##deepmap_parser(reversed,@show x.parser))
_deepmap_parser(::typeof(_reversed),mem::AbstractDict,x::NegativeLookahead) =
    NegativeLookbehind(x.parser) ##deepmap_parser(reversed,x.parser))
_deepmap_parser(::typeof(_reversed),mem::AbstractDict,x::PositiveLookbehind) =
    PositiveLookahead(x.parser) ##deepmap_parser(reversed,x.parser))
_deepmap_parser(::typeof(_reversed),mem::AbstractDict,x::PositiveLookahead) =
    PositiveLookbehind(x.parser) ##deepmap_parser(reversed,x.parser))


text:
abstract type Component end

function get_name(ch::Component)
    Base.typename(ch |> typeof).name |> Symbol
end

abstract type Compartment <: Component end
abstract type FlowChannel{Sensors<:Tuple,Actuators<:Tuple} <: Component end
abstract type Neuron <: Compartment end

abstract type Synapse <: Component end
struct EmptyConnection <: Synapse end

abstract type Geometry end

abstract type PlasticityRule{S} end

struct PlasticisedChannel{S,S2,A,C<:FlowChannel{S,A},P<:PlasticityRule{S2}} <: FlowChannel{Tuple{S,S2},A}
    channel::C
    modification::P
end

function get_name(p::PlasticisedChannel)
    Symbol(
        get_name(p.channel),
        :_with_,
        get_name(p.mutation)
    )
end

ext: jl


text:
module MapperGraphs

using ..Mapper2.Helper
Helper.@SetupDocStringTemplates

using LightGraphs

# Additions to the light graphs API
export SparseDiGraph, source_vertices, sink_vertices, linearize, make_lightgraph

# Manually export the needed items from LightGraphs instead of
# Reexporting to:
# 1. Not pollute the namespace with random unused names
# 2. Avoid Documentor looking at all the docstrings in LightGraphs and
# throwing a million "missing docstring" errors.

# Graph types
export AbstractGraph,
    SimpleDiGraph,
    DiGraph,
    # Mutating methods
    add_edge!,
    add_vertex!,
    add_vertices!,
    # Iterating methods
    edges,
    vertices,
    # Query Methods
    nv,
    ne,
    outneighbors,
    inneighbors,
    has_edge,
    has_vertex,
    src,
    dst,
    # Analysis
    is_weakly_connected,
    has_path

################################################################################
# Subgraph for for lightgraphs
################################################################################
struct AdjacencyList{T}
    neighbors_in::Vector{T}
    neighbors_out::Vector{T}
end

AdjacencyList{T}() where {T} = AdjacencyList(T[], T[])

"""
Graph representation with arbitrary vertices of type `T`.
"""
struct SparseDiGraph{T}
    vertices::Dict{T,AdjacencyList{T}}
end

SparseDiGraph{T}() where {T} = SparseDiGraph(Dict{T,AdjacencyList{T}}())

LightGraphs.has_vertex(g::SparseDiGraph, v) = haskey(g.vertices, v)

function LightGraphs.add_vertex!(g::SparseDiGraph{T}, v) where {T}
    has_vertex(g, v) && return false
    g.vertices[v] = AdjacencyList{T}()
    return true
end

function LightGraphs.add_edge!(g::SparseDiGraph, src, snk)
    has_vertex(g, src) || throw(KeyError(src))
    has_vertex(g, snk) || throw(KeyError(snk))
    push!(g.vertices[src].neighbors_out, snk)
    push!(g.vertices[snk].neighbors_in, src)
    return nothing
end

LightGraphs.vertices(g::SparseDiGraph) = keys(g.vertices)
LightGraphs.outneighbors(g::SparseDiGraph, v) = g.vertices[v].neighbors_out
LightGraphs.inneighbors(g::SparseDiGraph, v) = g.vertices[v].neighbors_in
LightGraphs.nv(g::SparseDiGraph) = length(g.vertices)

"""
    source_vertices(graph::SparseDiGraph{T}) :: Vector{T} where T

Return the vertices of `graph` that have no incoming edges.
"""
function source_vertices(g::SparseDiGraph)
    return [v for v in vertices(g) if length(inneighbors(g, v)) == 0]
end

"""
    sink_vertices(graph::SparseDiGraph{T}) :: Vector{T} where T

Return the vertices of `graph` that have no outgoing edges.
"""
function sink_vertices(g::SparseDiGraph)
    return [v for v in vertices(g) if length(outneighbors(g, v)) == 0]
end

"""
    linearize(graph::SparseDiGraph{T}) where T

Return a Vector{T} of vertices of `graph` in linearized traversal order if
`graph` is linear. Throw error if `ggraph` is not a linear graph.
"""
function linearize(g::SparseDiGraph{T}) where {T}
    # Do a check for an empty graph.
    if nv(g) == 0
        return T[]
    end

    sv = source_vertices(g)
    length(sv) != 1 && error("Expected 1 source vertex. Found $(length(sv)).")
    vertices = T[first(sv)]
    neighbors = outneighbors(g, last(vertices))
    while length(neighbors) > 0
        if length(neighbors) > 1
            error("""
                Vertex $(last(vertices)) has $(length(neighbors)) neighbors.
                Expected 1.
                """)
        end
        push!(vertices, first(neighbors))
        neighbors = outneighbors(g, last(vertices))
    end
    return vertices
end

"""
    make_lightgraph(graph::SparseDiGraph{T}) where T :: (SimpleDiGraph, Dict{T,Int})

Return tuple `(g, d)` where `g :: SimpleDiGraph` is a `LightGraph` ismorphic
to `graph` and `d` maps vertices of `graph` to vertices of `g`.
"""
function make_lightgraph(s::SparseDiGraph)
    g = DiGraph(nv(s))
    # Create a mapping dictionary
    d = Dict(v => i for (i, v) in enumerate(vertices(s)))
    # Iterate over edges, adding each edge to the lightgraph
    for i in vertices(s), j in outneighbors(s, i)
        add_edge!(g, d[i], d[j])
    end
    return g, d
end

end


text:
include("utils.jl")

build_so("001")

x = 5.3

y = ccall(
    (:foo, "binary/lib001.so"),
    Cdouble,
    (Cdouble, ),
    x
)

@assert y == x + 3.14


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DataFrames: DataFrame
using Suppressor

include("../helper/gpu.jl")

@testset "Likelihood" begin
    # Uses a PC with 4 variables, and tests 3 of the configurations to
    # match with python. Also tests all probabilities sum up to 1.

    prob_circuit = zoo_psdd("little_4var.psdd");
    @test prob_circuit isa ProbCircuit;

    # Step 1. Check Probabilities for 3 samples
    data = DataFrame(BitArray([0 0 0 0; 0 1 1 0; 0 0 1 1]));
    true_prob = [0.07; 0.03; 0.13999999999999999]

    calc_prob = EVI(prob_circuit, data)
    calc_prob = exp.(calc_prob)

    @test true_prob ≈ calc_prob atol=1e-7;

    # Step 2. Add up all probabilities and see if they add up to one
    N = 4;
    data_all = generate_data_all(N)

    calc_prob_all = EVI(prob_circuit, data_all)
    calc_prob_all = exp.(calc_prob_all)
    sum_prob_all = sum(calc_prob_all)

    @test 1 ≈ sum_prob_all atol = 1e-7;

    cpu_gpu_agree_approx(data_all) do d
        EVI(prob_circuit, d)
    end

    # make sure log-likelihoods are -Inf when the input is not satisfied
    alltrue = multiply(pos_literals(ProbCircuit,4))
    @test all(EVI(alltrue, data) .== -Inf)

    cpu_gpu_agree(data) do d
        EVI(alltrue, d)
    end

    # Test Sturdel EVI
    samples, _ = sample(prob_circuit, 100000)
    mix, weights, _ = learn_strudel(DataFrame(convert(BitArray, samples)); num_mix = 10,
                                    init_maxiter = 20, em_maxiter = 100, verbose = false)
    mix_calc_prob = exp.(EVI(mix, data, weights))

    @test true_prob ≈ mix_calc_prob atol = 0.1
    mix_calc_prob_all = exp.(EVI(mix, data_all))
    @test 1 ≈ sum(mix_calc_prob_all) atol = 0.1

    cpu_gpu_agree_approx(data_all) do d
        EVI(mix, d, weights)
    end
end

@testset "Bagging models' likelihood" begin
    dfb = DataFrame(BitMatrix([true true; true true; true true; true true]))
    r = fully_factorized_circuit(ProbCircuit,num_features(dfb))
    # bag_dfb = bagging_dataset(dfb; num_bags = 2, frac_examples = 1.0)
    bag_dfb = Vector{DataFrame}(undef, 2)
    bag_dfb[1] = dfb[[2, 1, 3, 4], :]
    bag_dfb[2] = dfb[[4, 3, 2, 1], :]

    r = compile(SharedProbCircuit, r, 2)

    estimate_parameters(r, bag_dfb; pseudocount = 1.0)

    ll = log_likelihood_per_instance(r, dfb)
    @test ll[1] ≈ -0.2107210 atol = 1e-6
    @test ll[2] ≈ -0.2107210 atol = 1e-6
    @test ll[3] ≈ -0.2107210 atol = 1e-6
    @test ll[4] ≈ -0.2107210 atol = 1e-6

    @test log_likelihood(r, dfb) ≈ -0.8428841 atol = 1e-6

    @test log_likelihood_avg(r, dfb) ≈ -0.2107210 atol = 1e-6
end


text:
using Test, GeneticAlgorithm
using Random

Random.seed!(42)

global numberOfFitnessEvaluation = 0

"""
createGene()

Create a gene value
"""
function createGeneExample()
    return floor(Int, rand() * 10)
end

"""
    fitnessExample(ind)

Example of a fitness function that determines how good an individual is. It returns a number.
"""
function fitnessExample(ind)
    global numberOfFitnessEvaluation += 1
    solution = 1:5
    allDifferent = filter(((a,b),) -> a != b, collect(zip(solution, ind)))
    return length(allDifferent)
end

@test fitnessExample([3, 2, 1, 4, 5]) == 2
@test fitnessExample(reverse(collect(1:5))) == 4
@test fitnessExample(1:5) == 0

@test pickBestIndividual(fitnessExample, [[3, 2, 1, 4, 5], reverse(collect(1:5)), 1:5]) == 1:5
@test pickBestIndividual(fitnessExample, [[3, 2, 1, 4, 5], reverse(collect(1:5))]) == [3, 2, 1, 4, 5]

@test selectIndividual(fitnessExample, [[3, 2, 1, 4, 5], reverse(collect(1:5)), 1:5]) == 1:5
@test selectIndividual(fitnessExample, [[3, 2, 1, 4, 5], reverse(collect(1:5))]) == [3, 2, 1, 4, 5]

numberOfFitnessEvaluation = 0
@test runGA(fitnessExample, createGeneExample, 5, maxNumberOfIterations=200, logging=false)[1] == [1, 2, 3, 4, 5]
@test numberOfFitnessEvaluation == 24010

## Checking cache!!!
numberOfFitnessEvaluation = 0
fitnessCache = Dict()
function fitnessExampleWithCache(ind)
    if get(fitnessCache, ind, :notPresent) !== :notPresent
        return fitnessCache[ind]
    end
    global numberOfFitnessEvaluation += 1
    solution = 1:5
    allDifferent = filter(((a,b),) -> a != b, collect(zip(solution, ind)))
    fitnessCache[ind] = length(allDifferent)
    return fitnessCache[ind]
end
@test runGA(fitnessExampleWithCache, createGeneExample, 5, maxNumberOfIterations=200, logging=false)[1] == [1, 2, 3, 4, 5]
@test numberOfFitnessEvaluation == 187


text:
mutable struct AccumulatorObsScalarGPU{A,B,C,D} <: AbstractObservableAccumulator
    res::A

    ψ_σ::B          # Buffers alls the <σ|ψ> of the denominator
    Δ_σ::C
    mel_buf::C      # ⟨σ|Ô|σ'⟩ in the buffer

    logψ_acc::D
end

function AccumulatorObsScalar(net::NeuralNetwork, hilb, v::Union{gpuAStateBatch, gpuADoubleStateBatch},
                              batch_sz)
    bnet      = cached(net, batch_sz)
    CT        = Complex{real(out_type(net))}
    v         = state(hilb, bnet)

    acc       = AccumulatorLogPsi(bnet, v, batch_sz)

    w         = trainable_first(net)
    RT        = real(eltype(w))

    res       = zero(eltype(out_similar(bnet)))
    mel_buf   = collect(similar(w, CT, 1, batch_sz))
    Δ_σ       = similar(mel_buf)

    return AccumulatorObsScalarGPU(
        res,
        res,
        Δ_σ, mel_buf,
        acc)
end

@forward AccumulatorObsScalarGPU.logψ_acc Base.length, Base.count,
    isfull

accum(a::AccumulatorObsScalarGPU) = a.logψ_acc

function init!(c::AccumulatorObsScalarGPU, σ, ψ_σ)
    c.ψ_σ     = ψ_σ
    c.mel_buf .= 0
    c.res     = 0

    init!(accum(c), σ)
    return nothing
end

# This resets the value at the end of a process_accumulatr
function reset!(c::AccumulatorObsScalarGPU, σ)
    c.mel_buf .= 0

    init!(accum(c), σ)
    return nothing
end

function (c::AccumulatorObsScalarGPU)(mel::Number, cngs_l, cngs_r, v)
    isfull(c) && process_accumulator!(c)

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return acc

    n_cngs_l = isnothing(cngs_l) ? 0 : length(cngs_l)
    n_cngs_r = isnothing(cngs_r) ? 0 : length(cngs_r)

    # If there are no changes, just sum it
    if n_cngs_l == 0 && n_cngs_r == 0
        c.res += mel
    else
        accum(c)(cngs_l, cngs_r)
        c.mel_buf[count(c)] = mel
    end

    return c
end

function (c::AccumulatorObsScalarGPU)(mel::Number, cngs, v)
    isfull(c) && process_accumulator!(c)

    # If the matrix element is zero, don't do anything
    mel == 0.0 && return c

    n_cngs = isnothing(cngs) ? 0 : length(cngs)

    # If there are no changes, just sum it
    if cngs == 0
        c.res += mel
    else
        accum(c)(cngs)
        c.mel_buf[count(c)] = mel
    end

    return c
end

function process_accumulator!(c::AccumulatorObsScalarGPU)
    count(c) == 0 && return nothing

    ψ_σp = process_accumulator!(accum(c))

    ψ_σp .= exp.(ψ_σp .- c.ψ_σ)
    copy!(c.Δ_σ, ψ_σp)

    c.mel_buf .*= c.Δ_σ
    c.res += sum(c.mel_buf)

    return nothing
end

function finalize!(c::AccumulatorObsScalarGPU)
    process_accumulator!(c)
    return c.res
end


text:
using Statistics

#lossy
function downsample_4_2_0(x::Matrix{Float32})::Matrix{Float32}
    retimg = Matrix{Float32}(undef, Int.(size(x)./2))
    for i=1:Int(size(x)[1]/2), j=1:Int(size(x)[2]/2)
        block_lower = 2*i
        block_upper = block_lower-1
        block_right = 2*j
        block_left = block_right-1
        subimg = x[block_upper:block_lower,block_left:block_right]
        retimg[i,j] = mean(subimg)
    end
    return retimg
end

function upsample_4_2_0(x::Matrix{Float32})::Matrix{Float32}
    retimg = Matrix{Float32}(undef, (size(x).*2))
    for i=1:Int(size(x)[1]), j=1:Int(size(x)[2])
        block_lower = 2*i
        block_upper = block_lower-1
        block_right = 2*j
        block_left = block_right-1
        thisval = x[i,j]
        retimg[block_upper:block_lower,block_left:block_right] .= thisval
    end
    return retimg
end


text:
# this code was closely ported from the rlpy implementaiton https://github.com/rlpy/rlpy/blob/master/rlpy/Domains/Pinball.py
# configs/*.cfg come directly from rlpy

struct PinBallConfig{T}
    start_pos::Tuple{T,T}
    target_pos::Tuple{T,T}
    target_radius::T
    ball_radius::T
    noise::T
    drag::T
    force::T

    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(0.02), T(0.), T(0.995), T(1. / 5.))
    end
    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius, ball_radius) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(ball_radius), T(0.), T(0.995), T(1. / 5.))
    end
    function PinBallConfig(::Type{T}, start_pos::Tuple, target_pos::Tuple, target_radius, ball_radius, noise, drag, force) where {T}
        new{T}(convert.(T, start_pos), convert.(T, target_pos), T(target_radius), T(ball_radius), T(noise), T(drag), T(force))
    end
end

mutable struct BallState{T} <:Any where {T}
    x::T
    y::T
    xDot::T
    yDot::T
    radius::T

    function BallState(position::Tuple{T,T}, radius::T) where {T}
        new{Float64}(Float64(position[1]), Float64(position[2]), 0., 0., Float64(radius))
    end

    function BallState(::Type{T}, position::Tuple{T,T}, radius::T) where {T}
        new{T}(T(position[1]), T(position[2]), T(0.), T(0.), T(radius))
    end
end

struct PinballObstacle{T} <: Any where {T}
    points::Array{Tuple{T,T}, 1}
    minx::T
    miny::T
    maxx::T
    maxy::T

    function PinballObstacle(::Type{T}, points::Array{Tuple{TX, TY}, 1}) where {T,TX,TY}
        pts = [T.(p) for p in points]
        minx, miny = min.(pts...)
        maxx, maxy = max.(pts...)
        new{T}(pts, minx, miny, maxx, maxy)
    end
end

mutable struct PinBall{T,Int} <: AbstractEnvironment where {T<:Real}
    state::Array{T,1}
    reward::T
    done::Bool
    ball::BallState{T}
    obstacles::Array{PinballObstacle{T}, 1}
    config::PinBallConfig{T}
    dt::T
    t::T # time
	deterministic_start::Bool

    function PinBall(T::Type, config::String, deterministic_start::Bool=false)
        obstacles, conf = read_config(T, config)
        ball = BallState(conf.start_pos, conf.ball_radius)
        new{T,Int}(zeros(T,4), T(0.0), false, ball, obstacles, conf, T(1. / 20.), T(0.), deterministic_start)
    end

	function PinBall(config::String, deterministic_start::Bool=false)
		return PinBall(Float64, config, deterministic_start)
	end

	function PinBall(::Type{T}, config::PinBallConfig{T}, deterministic_start::Bool, obstacles::Array{PinballObstacle{T}, 1}, ball::BallState{T}) where {T}
		b = BallState(T, (ball.x, ball.y), ball.radius)
		b.xDot = ball.xDot
		b.yDot = ball.yDot
		new{T,Int}(zeros(T,4), T(0.0), false, b, deepcopy(obstacles), deepcopy(config), T(1. / 20.), T(0.), deterministic_start)
	end
end

function get_gamma(env::PinBall{T})::T where {T}
	return T(1.0)
end

function randomize_env(::Type{PinBall{T}}, rng::AbstractRNG, config::String, deterministic_start::Bool=false) where {T}
	obstacles, conf = read_config(T, config)
	target_radius = conf.target_radius + rand(rng, Uniform(-0.1, 0.1))  * conf.target_radius  # scaled target radius randomly by standard deviation of 10% of the specified radius
    ball_radius = conf.ball_radius + rand(rng, Uniform(-0.1, 0.1)) * conf.ball_radius
    noise = rand(rng, T) * 0.25    # chance for random action
    drag = 1. - exp(rand(rng, Uniform(log(0.001), log(0.1))))  # friction coefficient for ball
    force =  rand(rng, Uniform(0.1, 0.3))  # force applied to the ball
	conf = PinBallConfig(T, conf.start_pos, conf.target_pos, target_radius, ball_radius, noise, drag, force)

	ball = BallState(T, (conf.start_pos[1], conf.start_pos[2]), ball_radius)
	return PinBall(T,conf, deterministic_start, obstacles, ball)
end

function get_state_desc(env::PinBall{T})::Array{T,2} where {T<:Real}
    ranges = zeros(T, (4,2))
    ranges[1,:] .= [0., 1.]  # x range
	ranges[2,:] .= [0., 1.]  # y range
	ranges[3,:] .= [-2., 2.] # xdot range
	ranges[4,:] .= [-2., 2.] # ydot range
    return ranges
end

function get_action_desc(env::PinBall{T,TA})::Int where{T<:Real,TA<:Integer}
    return 5
end

function step!(env::PinBall, action::Int, rng::AbstractRNG)
    if action <= 0 || action > 5
        error("Action needs to be an integer in [1, 5]")
    end

	if rand(rng) < env.config.noise
		action = rand(rng, Int(1):Int(5))
	end
    # add action effect
    if action == 1
        add_impulse!(env.ball, env.config.force, 0.)  # Acc x
    elseif action ==2
        add_impulse!(env.ball, 0., -env.config.force) # Dec y
    elseif action ==3
        add_impulse!(env.ball, -env.config.force, 0.) # Dec x
    elseif action == 4
        add_impulse!(env.ball, 0., env.config.force) # Acc y
    else
        add_impulse!(env.ball, 0., 0.)  # No action
    end

	dxdy = [0., 0.]

    for i in 1:20
        stepball!(env.ball, env.dt)

        ncollision = 0
        fill!(dxdy, 0.)

        for obs in env.obstacles
            hit, double_collision, intercept = collision(obs, env.ball)
            if hit
				dxdy .= dxdy .+ collision_effect(env.ball, hit, double_collision, intercept)
                ncollision += 1
            end
		end
        if ncollision == 1
            env.ball.xDot = dxdy[1]
            env.ball.yDot = dxdy[2]
            if i == 19
                stepball!(env.ball, env.dt)
            end
        elseif ncollision > 1
            env.ball.xDot = -env.ball.xDot
            env.ball.yDot = -env.ball.yDot
        end
		env.t += env.dt
        found_goal = at_goal(env.ball, env.config)
		done = found_goal || env.t > 1000.

        env.done = done
        if done
            env.reward = (10000. * found_goal) - (1 - found_goal)
            update_state!(env)
            return env.reward
        end
    end

    add_drag!(env.ball, env.config.drag)
    checkbounds!(env.ball)

    if action == 5
        env.reward = -1.
    else
        env.reward = -5.
    end

    update_state!(env)

    return env.reward
end

function checkbounds!(ball::BallState)
    if ball.x > 1.0
        ball.x = 0.95
    end
    if ball.x < 0.0
        ball.x = 0.05
    end
    if ball.y > 1.0
        ball.y = 0.95
    end
    if ball.y < 0.0
        ball.y = 0.05
    end
end

function update_state!(env::PinBall)
    env.state .= env.ball.x, env.ball.y, env.ball.xDot, env.ball.yDot
end

function at_goal(ball::BallState, config::PinBallConfig)
	norm((ball.x, ball.y) .- config.target_pos) < config.target_radius
end

function terminalCheck(env::PinBall)
    return at_goal(env.ball, env.config) || env.t > 1000.
end

function reset_ball!(ball::BallState{T}, start_pos::Tuple{T,T}, deterministic_start::Bool, rng::AbstractRNG) where {T}
	ball.x = start_pos[1]
	ball.y = start_pos[2]
	if !deterministic_start
		ball.x += 0.02 * randn(rng, T)
		ball.y += 0.02 * randn(rng, T)
	end
	ball.xDot = 0.
	ball.yDot = 0.
end

function env_reset!(env::PinBall, rng::AbstractRNG)
	reset_ball!(env.ball, env.config.start_pos, env.deterministic_start, rng)
	update_state!(env)
	env.t = 0.
	env.done = false
end

function is_terminal(env::PinBall)::Bool
    return env.done
end

function clone(env::PinBall{T})::PinBall{T} where {T}
	env2 = PinBall(T, env.config, env.deterministic_start, env.obstacles, env.ball)
	env2.state .= env.state
	env2.reward = env.reward
	env2.done = env.done
	env2.dt = env.dt
	env2.t = env.t

	return env2
end

function add_impulse!(ball::BallState{T}, Δx::T, Δy::T) where {T}
    xDot = ball.xDot + Δx
    yDot = ball.yDot + Δy
    ball.xDot = clamp(xDot, -2., 2.)
    ball.yDot = clamp(yDot, -2., 2.)
end

function add_drag!(ball::BallState{T}, drag::T) where {T}
    ball.xDot *= drag
    ball.yDot *= drag
end

function stepball!(ball::BallState{T}, dt::T) where {T}
    ball.x += ball.xDot * ball.radius * dt
    ball.y += ball.yDot * ball.radius * dt
end

function read_config(::Type{T}, source) where {T}
	source = joinpath(@__DIR__, "configs", source)
    obstacles = Array{PinballObstacle{T}, 1}()
    target_pos = Tuple{T,T}((0.,0.))
    target_radius = T(0.04)
    ball_radius = T(0.02)
    start_pos = Tuple{T,T}((0.,0.))
    noise = T(0.)
    drag = T(0.995)
    force = T(1. / 5.)
    lines = readlines(source)
    for line in lines
        tokens = split(strip(line))
        if length(tokens) <= 0
            continue
        elseif tokens[1] == "polygon"
			nums = map(x->parse(T, x), tokens[2:end])
			points = [(x,y) for (x,y) in zip(nums[1:2:end], nums[2:2:end])]
            push!(obstacles, PinballObstacle(T, points))
        elseif tokens[1] == "target"
            target_pos = Tuple{T,T}(map(x->parse(T,x), tokens[2:3]))
            target_radius = parse(T,tokens[4])
        elseif tokens[1] == "start"
            start_pos = Tuple{T,T}(map(x->parse(T,x), tokens[2:3]))
        elseif tokens[1] == "ball"
            ball_radius = parse(T, tokens[2])
        end
    end

    conf = PinBallConfig(T, start_pos, target_pos, target_radius, ball_radius, noise, drag, force)

    return obstacles, conf
end

function collision(obs::PinballObstacle{T}, ball::BallState{T}) where {T}
    if ball.x - ball.radius > obs.maxx
        return false, nothing, nothing
	end
	if ball.x + ball.radius < obs.minx
        return false, nothing, nothing
	end
    if ball.y - ball.radius > obs.maxy
        return false, nothing, nothing
	end
    if ball.y + ball.radius < obs.miny
        return false, nothing, nothing
	end
    double_collision = false
    intercept_found = false

    i = 1
    j = 2
    intercept = nothing
    while i ≤ length(obs.points)
        p1, p2 = obs.points[i], obs.points[j]
        if intercept_edge(p1, p2, ball)
            if intercept_found
                intercept = select_edge((p1, p2), intercept, ball)
                double_collision = true
            else
                intercept = (p1, p2)
                intercept_found = true
            end
        end
        i += 1
        j += 1
        if j > length(obs.points)
            j = 1
        end
    end
    return intercept_found, double_collision, intercept
end

function collision_effect(ball::BallState{T}, intercept_found::Bool, double_collision::Bool, intercept::Tuple{Tuple{T,T},Tuple{T,T}})::Tuple{T,T} where {T}
    if double_collision
        return -ball.xDot, -ball.yDot
    end

    obstacle_vector = intercept[2] .- intercept[1]
    if obstacle_vector[1] < 0.
        obstacle_vector = intercept[1] .- intercept[2]
    end

    velocity_vector = (ball.xDot, ball.yDot)
    θ = compute_angle(velocity_vector, obstacle_vector) - π
    if θ < 0.
        θ += 2π
    end

    intercept_theta = compute_angle([-1, 0], obstacle_vector)
    θ += intercept_theta

    velocity = norm(velocity_vector)

    return velocity * cos(θ), velocity * sin(θ)
end

function compute_angle(v1, v2)
    angle_diff = atan(v1[1], v1[2]) - atan(v2[1], v2[2])
    if angle_diff < 0.
        angle_diff += 2π
    end
    return angle_diff
end

function intercept_edge(p1::Tuple{T,T}, p2::Tuple{T,T}, ball::BallState{T}) where {T}
    edge = p2 .- p1
    pball = (ball.x, ball.y)
    difference = pball .- p1

    scalar_proj = dot(difference, edge) / dot(edge, edge)
    scalar_proj = clamp(scalar_proj, 0., 1.)

    closest_pt = p1 .+ (edge .* scalar_proj)
    obstacle_to_ball = pball .- closest_pt
    distance = dot(obstacle_to_ball, obstacle_to_ball)

    if distance <= ball.radius^2
        # collision if the ball is not moving away
        velocity = (ball.xDot, ball.yDot)
        ball_to_obstacle = closest_pt .- pball

        angle = compute_angle(ball_to_obstacle, velocity)
        if angle > π
            angle = 2π - angle
        end

        if angle > (π / 1.99)
            return false
        end
        return true
    else
        return false
    end
end

function select_edge(intersect1::Tuple{Tuple{T,T},Tuple{T,T}}, intersect2::Tuple{Tuple{T,T},Tuple{T,T}}, ball::BallState{T}) where {T}
    velocity = (ball.xDot, ball.yDot)
    obstacle_vector1 = intersect1[2] .- intersect1[1]
    obstacle_vector2 = intersect2[2] .- intersect2[1]
    angle1 = compute_angle(velocity, obstacle_vector1)
    if angle1 > π
        angle1 -= π
    end

    angle2 = compute_angle(velocity, obstacle_vector2)
    if angle2 > π
        angle2 -= π
    end

    if abs(angle1 - π / 2.) < abs(angle2 - π / 2.)
        return intersect1
    else
        return intersect2
    end
end

# function circleShape(x, y, r)
# 	θ = LinRange(0, 2*π, 500)
# 	x .+ r*sin.(θ), y .+ r*cos.(θ)
# end

@recipe function f(env::PinBall)
	ballx, bally, bradius = env.ball.x, env.ball.y, env.ball.radius
	obstacles = env.obstacles
	tx, ty = env.config.target_pos
	tr = env.config.target_radius

	legend := false
	xlims := (0., 1.)
	ylims := (0., 1.)
	grid := false
	ticks := nothing
	foreground_color := :white
	aspect_ratio := 1.

	# obstacles
	for ob in obstacles
		@series begin
			seriestype := :shape
			seriescolor := :blue
			xpts = [p[1] for p in ob.points]
			ypts = [p[2] for p in ob.points]

			xpts, ypts
		end
	end

	# goal
	@series begin
		seriestype := :shape
		linecolor := nothing
		color := :green
		aspect_ratio := 1.
		fillalpha := 0.4

		circleShape(tx, ty, tr)
	end

	# ball
	@series begin
		seriestype := :shape
		linecolor := nothing
		color := :red
		aspect_ratio := 1.

		circleShape(ballx, bally, bradius)
	end

end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
"""
Bridges `CP.AllDifferent` to a series of `CP.DifferentFrom`.
"""
struct AllDifferent2DifferentFromBridge{T} <: MOIBC.AbstractBridge
    # An upper-triangular matrix (i.e. nothing if i < j, constraint if i >= j).
    # Standard sparse matrices cannot store anything that is not a Number
    # (more specifically, anything that does not implement `zero(T)`).
    # https://github.com/JuliaLang/julia/issues/30573
    cons::Dict{
        Tuple{Int, Int},
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T}},
    }
end

function MOIBC.bridge_constraint(
    ::Type{AllDifferent2DifferentFromBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.AllDifferent,
) where {T}
    return MOIBC.bridge_constraint(
        AllDifferent2DifferentFromBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{AllDifferent2DifferentFromBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.AllDifferent,
) where {T}
    f_scalars = MOIU.scalarize(f)
    dim = MOI.output_dimension(f)

    # Upper-triangular matrix of constraints: i >= j, i.e. d(d-1)/2 elements:
    #     \sum_{i=2}^{d} (n - i + 1) = d (d - 1) / 2
    cons = Dict{
        Tuple{Int, Int},
        MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T}},
    }()
    sizehint!(cons, dim * (dim - 1) / 2)

    for i in 1:dim
        for j in (i+1):dim
            cons[i, j] = MOI.add_constraint(
                model,
                f_scalars[i] - f_scalars[j],
                CP.DifferentFrom(zero(T))
            )
        end
    end

    return AllDifferent2DifferentFromBridge(cons)
end

function MOI.supports_constraint(
    ::Type{AllDifferent2DifferentFromBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.AllDifferent},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{AllDifferent2DifferentFromBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{AllDifferent2DifferentFromBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T}),
    ]
end

function MOI.get(
    b::AllDifferent2DifferentFromBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T},
    },
) where {T}
    return length(b.cons)
end

function MOI.get(
    b::AllDifferent2DifferentFromBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, CP.DifferentFrom{T},
    },
) where {T}
    return collect(values(b.cons))
end


text:
using LazyInverses
using LinearAlgebra
using BenchmarkTools

n = 128
x = randn(n)
y = randn(n)
A = randn(n, n)
A = A'A
C = cholesky(A)
invC = inverse(C)

println("ternary dot-product multiplication")
@btime dot($x, $invC, $x)
@btime dot($x, $invC, $y)
@btime dot($x, $C \ $x)
@btime dot($x, $C \ $y)

k = n
X = randn(k, n)
Y = randn(n, k)
println("ternary multiplication")
@btime *($X, $invC, $X')
@btime *($X, $invC, $Y)
@btime *($X, $C \ $X')
@btime *($X, $C \ $Y)


text:
# This is a markdown comment. We can use Markdown syntax, like calling this language `julia`, making stuff *italic* and **bold**, etc.
## This is a Julia comment.
#md # This will only show up in markdown output.
#nb # This will only show up in Jupyter output.

# This is some stuff from another file in this repo.

type TestType
    x::Any
    y::Real
    z::Integer
end

print("This Julia code should work no matter what.")

@assert 1 == 1

# This is a test of pre-processing. The date is: DATETIME.


text:
@testset "Constant" begin
    # Instantiation
    N1 = 10
    A1 = rand(Float64, N1) * 100
    A2 = rand(Float64, N1, N1) * 100
    A3 = rand(Float64, N1, N1, N1) * 100

    @testset "constructor" begin
        @test Constant() == Constant{Nearest}()
        for T in (Nearest, Previous, Next)
            it = Constant{T}()
            @test it isa Constant{T}
        end
    end

    @testset "Constant{Nearest}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i] == itp1(i + .3) == itp1(i - .3)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i,j] == itp2(i + .4, j - .3) == itp2(i - .4, j + .3)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i,j,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end

    @testset "Constant{Previous}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Previous}())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Previous}())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Previous}())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i] == itp1(i + .3) == itp1(i + .6)
                @test A1[i - 1] == itp1(i - .3) == itp1(i - .6)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i,j - 1] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i,j - 1,k] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end

    @testset "Constant{Next}" begin
        for (constructor, copier) in ((interpolate, x -> x), (interpolate!, copy))
            isinplace = constructor == interpolate!
            itp1 = @inferred(constructor(copier(A1), BSpline(Constant{Next}())))
            itp2 = @inferred(constructor(copier(A2), BSpline(Constant{Next}())))
            itp3 = @inferred(constructor(copier(A3), BSpline(Constant{Next}())))

            @test parent(itp1) === itp1.coefs
            @test Interpolations.lbounds(itp1) == (1,)
            @test Interpolations.ubounds(itp1) == (N1,)

            # Evaluation on provided data points
            for (itp, A) in ((itp1, A1), (itp2, A2), (itp3, A3))
                check_axes(itp, A, isinplace)
                check_inbounds_values(itp, A)
                check_oob(itp)
                can_eval_near_boundaries(itp)
            end

            # Evaluation between data points (tests constancy)
            for i in 2:N1 - 1
                @test A1[i + 1] == itp1(i + .3) == itp1(i + .6)
                @test A1[i] == itp1(i - .3) == itp1(i - .6)
            end
            # 2D
            for i in 2:N1 - 1, j in 2:N1 - 1
                @test A2[i + 1,j] == itp2(i + .4, j - .3) == itp2(i + .4, j - .6)
            end
            # 3D
            for i in 2:N1 - 1, j in 2:N1 - 1, k in 2:N1 - 1
                @test A3[i + 1,j,k + 1] == itp3(i + .4, j - .3, k + .1) == itp3(i + .4, j - .3, k + .2)
            end
        end
    end
end


text:
using Test

include("rotational-cipher.jl")

@testset "rotate function" begin
    @testset "rotate by n" begin
        @testset "no wrap" begin
            @test rotate(1, "a") == "b"
            @test rotate(1, 'a') == 'b'
            @test rotate(13, "m") == "z"
            @test rotate(13, 'm') == 'z'
        end
        @testset "wrap around" begin
            @test rotate(13, "n") == "a"
            @test rotate(13, 'n') == 'a'
        end
    end

    @testset "full rotation" begin
        @test rotate(26, "a") == "a"
        @test rotate(26, 'a') == 'a'
        @test rotate(0, "a") == "a"
        @test rotate(0, 'a') == 'a'
    end

    @testset "full strings" begin
        @test rotate(5, "OMG") == "TRL"
        @test rotate(5, "O M G") == "T R L"
        @test rotate(4, "Testing 1 2 3 testing") == "Xiwxmrk 1 2 3 xiwxmrk"
        @test rotate(21, "Let's eat, Grandma!") == "Gzo'n zvo, Bmviyhv!"
        @test rotate(13, "The quick brown fox jumps over the lazy dog.") == "Gur dhvpx oebja sbk whzcf bire gur ynml qbt."
    end
end

# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
# Additional exercises                                                        #
# Change @test_skip to @test for the optional bonus exercises from HINTS.md   #
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

# Bonus A
@testset "Bonus A: string literal R13" begin
    @test_skip R13"The quick brown fox jumps over the lazy dog." == "Gur dhvpx oebja sbk whzcf bire gur ynml qbt."
end

# Bonus B
@testset "Bonus B: string literals" begin
    @test_skip R5"OMG" == "TRL"
    @test_skip R4"Testing 1 2 3 testing" == "Xiwxmrk 1 2 3 xiwxmrk"
    @test_skip R21"Let's eat, Grandma!" == "Gzo'n zvo, Bmviyhv!"
    @test_skip R13"The quick brown fox jumps over the lazy dog." == "Gur dhvpx oebja sbk whzcf bire gur ynml qbt."
end


text:
@testset "Integer" begin
  @testset "range" begin

    r04 = fmpz(0):fmpz(4)
    r026 = fmpz(0):fmpz(2):fmpz(6)
    zlarge = fmpz(13)^100
    zjump = fmpz(11)^100
    rlarge = -zlarge:zlarge
    rlargejumps = -zlarge:zjump:zlarge

    @test fmpz(0):4 == 0:fmpz(4) == r04
    @test fmpz(0):2:6 == 0:2:fmpz(6) == 0:fmpz(2):6 == r026
    @test typeof(fmpz(0):4) == typeof(0:fmpz(4)) == typeof(r04)
    @test typeof(fmpz(0):fmpz(2):6) == typeof(0:fmpz(2):fmpz(6)) == typeof(0:fmpz(2):6) == typeof(r026)

    @test collect(r04) == fmpz[0, 1, 2, 3, 4]
    @test collect(r026) == fmpz[0, 2, 4, 6]
    @test collect(fmpz(-2):fmpz(5):fmpz(21)) == fmpz[-2, 3, 8, 13, 18]
    @test collect(fmpz(3):fmpz(-2):fmpz(-4)) == fmpz[3, 1, -1, -3]

    @test length(r04) == 5
    @test length(r026) == 4
    @test length(fmpz(-2):fmpz(5):fmpz(21)) == 5
    @test length(fmpz(3):fmpz(-2):fmpz(-4)) == 4

    # Julia assumes `length` to return an Integer, so if we want several base
    # functions to work, that should better be the case.
    @test length(r04) isa Integer
    @test length(r026) isa Integer
    @test length(rlarge) isa Integer
    @test length(rlargejumps) isa Integer

    @test r04[4] == r04[end-1] == 3
    @test r026[3] == r026[end-1] == 4
    @test rlarge[2zlarge] == rlarge[end-1] == rlarge[end-fmpz(1)] == zlarge - 1
    @test rlargejumps[2] == -zlarge + zjump
    @test rlargejumps[end] == zlarge - 2zlarge % zjump

    @test range(fmpz(0); step=2, stop=6) == fmpz(0):fmpz(2):fmpz(6)
    @test range(0; step=fmpz(2), stop=6) == fmpz(0):fmpz(2):fmpz(6)
    @test range(fmpz(0); step=2, length=4) == fmpz(0):fmpz(2):fmpz(6)
    VERSION >= v"1.7" && @test range(; stop=fmpz(0), length=3) == fmpz(-2):fmpz(0)

    @test 1 .+ r04 == fmpz(1):fmpz(5)
    @test fmpz(1) .+ r04 == fmpz(1):fmpz(5)
    @test fmpz(1) .+ (0:4) == fmpz(1):fmpz(5)
    @test 2 .* r04 == fmpz(0):fmpz(2):fmpz(8)
    @test fmpz(2) .* r04 == fmpz(0):fmpz(2):fmpz(8)
    @test 2 .* (0:4) == fmpz(0):fmpz(2):fmpz(8)

    @test mod(fmpz(7), r04) == fmpz(2)
    @test mod(7, r04) == 2
    @test mod(7, fmpz(10):fmpz(14)) == 12

    @test all(x -> x in r04, rand(r04, 5))
    @test all(x -> x in fmpz(0):fmpz(10):fmpz(100), rand(fmpz(0):fmpz(10):fmpz(100), 5))
    @test all(x -> x in rlarge, rand(rlarge, 20))
    @test all(x -> x in rlargejumps, rand(rlargejumps, 20))
  end
end


text:
module Tester

using Genie, App, Util, Genie.Configuration, Logger
SEARCHLIGHT_ON && eval(:(using Migration))

"""
    bootstrap_tests(cmd_args::String, config::Settings) :: Void

Sets up testing environment, includes test files, etc.
"""
function bootstrap_tests(cmd_args::String = "", config::Settings = App.config, resource::String = "") :: Void
  current_env = config.app_env

  set_test_env()

  include(abspath(joinpath(config.test_folder, "test_config.jl")))

  for (path, _, files) in (walkdir(abspath(joinpath(config.test_folder))) |> collect)
    for file_name in files
      isempty(resource) && endswith(file_name, "_test.jl") && include(joinpath(path, file_name))
      ! isempty(resource) && startswith(file_name, resource) && endswith(file_name, "_test.jl") && include(joinpath(path, file_name))
    end
  end

  App.config.app_env = current_env
  Logger.log("Switched app to >> $(uppercase(App.config.app_env)) << env", :debug)

  nothing
end

"""
    reset_db() :: Void

Prepares the test env DB running all migrations up.
"""
function reset_db() :: Void
  Migration.all_down()
  Migration.all_up()

  nothing
end

"""
    run_all_tests(cmd_args::String, config::Settings) :: Void

Runs all existing tests.
"""
function run_all_tests(cmd_args::String = "", config::Settings = App.config) :: Void
  bootstrap_tests(cmd_args, config)

  nothing
end

function run_tests() :: Void
  run_all_tests()
end
function run_tests(resource_name::Symbol) :: Void
  bootstrap_tests("", App.config, string(resource_name) |> lowercase)
end

"""
    set_test_env() :: Void

Switches Genie to the test env for the duration of the current execution.
"""
function set_test_env() :: Void
  # if ! is_test()
  #   Logger.log("You're attempting to run your test suite outside the TEST environment. This can lead to losing your production or development data, depending on your current/default environment.", :err, showst = false)
  # end
  if App.config.tests_force_test_env
    # Logger.log("Automatically switching to TEST environment to avoid data corruption. If you want to force running your test in a different environment, switch the `tests_force_test_env` variable to `false` in your env's config file.", :debug)
    App.config.app_env = TEST
    Logger.log("Switched app to >> $(uppercase(App.config.app_env)) << env", :debug)

    ! is_test() && error("Could not switch env")
  end

  nothing
end

end


text:
using Test
using PanelShift
using Dates

equalormi(x, y) = all([(ismissing(i) && ismissing(j)) || i == j for (i,j) in zip(x,y)])

@test equalormi(tlag([1;2;3], [4;5;6], 1), [missing;4;5])
@test equalormi(tlead([1;2;3], [4;5;6], 1), [5; 6; missing])

@test equalormi(tlag([1;2;3], [4;5;6], 2), [missing;missing;4])
@test equalormi(tlead([1;2;3], [4;5;6], 2), [6; missing; missing])

# unit-length vector
@test equalormi(tlag([1], [1]), [missing;])
@test equalormi(tlead([1], [1]), [missing;])

@test equalormi(tlag([1;3;5;6;7], [1;2;3;4;5], 2), [missing; 1; 2; missing; 3])
@test equalormi(tlag(float.([1;3;5;6;7]), [1;2;3;4;5], 2), [missing; 1; 2; missing; 3])

# non-numeric x and unequal gaps
@test equalormi(tlag([1;2;4;7;11], [:apple; :orange; :banana; :pineapple; :strawberry], 1), [missing; :apple; missing; missing; missing])
@test equalormi(tlag([1;2;4;7;11], [:apple; :orange; :banana; :pineapple; :strawberry], 2), [missing; missing; :orange; missing; missing])
@test equalormi(tlag([1;2;4;7;11], [:apple; :orange; :banana; :pineapple; :strawberry], 3), [missing; missing; :apple; :banana; missing])
@test equalormi(tlag([1;2;4;7;11], [:apple; :orange; :banana; :pineapple; :strawberry], 4), [missing; missing; missing; missing; :pineapple])
@test equalormi(tlead([1;2;4;7;11], [:apple; :orange; :banana; :pineapple; :strawberry], 4), [missing; missing; missing; :strawberry; missing])

# indexed by dates
@test equalormi(tlag([Date(2000,1,1), Date(2000, 1,2), Date(2000,1, 4)], [1,2,3], Day(1)), [missing; 1; missing])
@test equalormi(tlag([Date(2000,1,1), Date(2000, 1,2), Date(2000,1, 4)], [1,2,3], Day(2)), [missing; missing; 2])

# test shift
@test equalormi(tshift([1;2;3], [1;2;3], -1), tlead([1;2;3], [1;2;3], 1))
@test equalormi(tshift([1;2;3], [1;2;3], 1), tlag([1;2;3], [1;2;3], 1))

# safeguards
@test_throws ArgumentError tlag([1;2;2], [1,2,3])
@test_throws ArgumentError tlag([1;2;], [1,2,3])
@test_throws ArgumentError tlag([1;2;3], [1,2,3], 0)

tlead([1;NaN;3], [4;5;6], 1)
equalormi(tlead([1;2;3], [4;5;6], 1), [5; 6; missing])

using DataFrames

df = DataFrame(
    t = [1;2;3;4; 1;3;4; 1;4; 1],
    id = [1;1;1;1; 2;2;2; 3;3; 4],
    x = [1;2;3;4; 5;6;7; 8;9; 10]
)

panellag!(df, :id, :t, :x, :Lx)
@test equalormi(df.Lx, [missing; 1; 2; 3; missing; missing; 6; missing; missing; missing])


text:
using CowSay
using Documenter

DocMeta.setdocmeta!(CowSay, :DocTestSetup, :(using CowSay); recursive=true)

makedocs(;
    modules=[CowSay],
    authors="uncomfyhalomacro",
    repo="https://github.com/uncomfyhalomacro/CowSay.jl/blob/{commit}{path}#{line}",
    sitename="CowSay.jl",
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/uncomfyhalomacro/CowSay.jl",
    devbranch="main",
)


text:
using Coverage

cd(joinpath(@__DIR__, "..", "..")) do
    processed = process_folder()
    Codecov.submit_local(processed)
    covered_lines, total_lines = get_summary(processed)
    percentage = covered_lines / total_lines * 100
    println("($(percentage)%) covered")
end


text:
module TestSetProperty

include("preamble.jl")

struct ImmutableStruct
    a
    b
end

mutable struct MutableStruct
    a
    b
end

@testset begin
    @test setproperty!!((a=1, b=nothing), :b, 2) === (a=1, b=2)
    @test setproperty!!(ImmutableStruct(1, nothing), :b, 2) ==
        ImmutableStruct(1, 2)
    ms = MutableStruct(1, nothing)
    @test setproperty!!(ms, :b, 2) === ms
    @test ms.b == 2
end

end  # module


text:
using Brim
using Docile
using Lexicon

save("./doc/api.md", Brim)


text:
export BallOnPlate

"""
    BallOnPlate

The single-axis ball on a plate system example with a single actuator.

This system is a 2-state 1-input discrete-time system with bounds on both
the states and the inputs.

This example is from Section 9.6 of [1].

[1] S. Richter, ‘Computational complexity certification of gradient methods for real-time model predictive control’, Doctoral Thesis, ETH Zurich, 2012.
"""
struct BallOnPlate <: ControlBenchmarkProblem end

@addbenchmark BallOnPlate "lin01" "Ball on Plate (Single Axis)" "Single axis ball on plate" Set([:linear, :stable, :discrete])

function controlbenchmark( ::BallOnPlate )

    # System matrices
    A = [ 1.0 0.01;
          0.0 1.0 ]

    B = [-0.0004
         -0.0701]

    C = Matrix{Float64}( I, 2, 2 )

    D = fill( 0.0, (2, 1) )

    # State constraints
    stateBounds = Bounds( [ Bound( -0.2, 0.01 );
                            Bound( -0.1, 0.1 ) ] )

    # Input constraints
    inputBounds = Bounds( [ Bound( -0.0524, 0.0524 ) ] )

    return (; :sys => StateSpace( A, B, C, D, 0.01 ),
              :stateBounds => stateBounds,
              :inputBounds => inputBounds )

end


text:
"""
Trait container for dimension and array ordering

The default is `Forward()`, `Forward()`
"""
struct Order{D,A}
    dim::D
    array::A
end
Order() = Order(Forward(), Forward())

dimorder(order::Order) = order.dim
arrayorder(order::Order) = order.array

"""
Trait indicating that the array or dimension is in the normal forward order.
"""
struct Forward end

"""
Trait indicating that the array or dimension is in the reverse order.
Selector lookup or plotting will be reversed.
"""
struct Reverse end

Base.reverse(::Reverse) = Forward()
Base.reverse(::Forward) = Reverse()
Base.reverse(o::Order) = Order(revese(dimorder(o)), revese(arrayorder(o)))


text:
module EponymKeywordSyntax
export @eponym

using MacroTools: @capture, splitarg, MacroTools

function resolve_kwarg(kwarg)
    key, typ, splat, val = splitarg(kwarg)
    if splat
        kwarg
    elseif val == nothing
        :($key = $key)
    else
        kwarg
    end
end

"""

    @eponym(ex)

The `@eponym` macro infers the name of a keyword argument
from the name of its value. This is best explained by example:
```jldoctest
julia> using EponymKeywordSyntax

julia> f(args...; kw...) = (args=args, kw=kw)
f (generic function with 1 method)

julia> oh_dear_this_is_long_I_dont_want_to_write_it_twice=3
3
julia> f(oh_dear_this_is_long_I_dont_want_to_write_it_twice=oh_dear_this_is_long_I_dont_want_to_write_it_twice) == @eponym f(;oh_dear_this_is_long_I_dont_want_to_write_it_twice)
true

julia> a=1; b=2;

julia> f(a=a, b=b) == @eponym f(;a,b)
true

julia> (a=a, b=b) == @eponym (;a, b)
true

julia> (a=a, b=10) == @eponym (;a, b=10)
true

julia> f(1, a=a, b=b) == @eponym f(1;a, b)
true
```
"""
macro eponym(ex)
    esc(eponymmacro(ex))
end

function eponymmacro(ex)
    ret = if @capture ex (f_(args__; kw__))
        # function call
        kw = map(resolve_kwarg, kw)
        :($(f)($(args...); $(kw...)))
    elseif @capture ex (;kw__)
        kw = map(resolve_kwarg, kw)
        # NamedTuple
        :((;$(kw...)))
    elseif Meta.isexpr(ex, :call) || Meta.isexpr(ex, :tuple)
        ex
    else
        msg = """
        Macro `@eponym` can only be applied to function calls, or `NamedTuple` calls
        that begin with `;`.
        Got: $(ex)
        """
        throw(ArgumentError(msg))
    end
end

end # module


text:
## Julia 1.1.0

println("Hello World, Julia")

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
@info "Setting up for job..."
#JOB FILEPATHS
sib_wms_path = "/bench/PhD/NGS_binaries/BMI/sib_nuc_position_sequences.fa_wms.tr"
rys_wms_path = "/bench/PhD/NGS_binaries/BMI/rys_nuc_position_sequences.fa_wms.tr"

sib_df_binary = "/bench/PhD/NGS_binaries/BMI/sib_diff_positions"
rys_df_binary = "/bench/PhD/NGS_binaries/BMI/rys_diff_positions"
combined_df_binary = "/bench/PhD/NGS_binaries/BMI/combined_diff_positions"

sib_code_binary = "/bench/PhD/NGS_binaries/BMI/sib_diff_codes"
rys_code_binary = "/bench/PhD/NGS_binaries/BMI/rys_diff_codes"
combined_code_binary = "/bench/PhD/NGS_binaries/BMI/combined_diff_codes"

sib_diff_bg = "/bench/PhD/NGS_binaries/BMI/sib_diff_bg"
rys_diff_bg = "/bench/PhD/NGS_binaries/BMI/rys_diff_bg"
combined_diff_bg = "/bench/PhD/NGS_binaries/BMI/combined_diff_bg"

sib_e_pth = "/bench/PhD/NGS_binaries/BMI/sib_e"
rys_e_pth = "/bench/PhD/NGS_binaries/BMI/rys_e"
combined_e_pth = "/bench/PhD/NGS_binaries/BMI/combined_e"

#JOB CONSTANTS
const ensemble_size = 500
const no_sources = 8
const source_min_bases = 3
const source_max_bases = 10
@assert source_min_bases < source_max_bases
const source_length_range= source_min_bases:source_max_bases
const mixing_prior = .07
@assert mixing_prior >= 0 && mixing_prior <= 1

@info "Loading master libraries..."
using Distributed, Distributions, Serialization

@info "Adding workers..."
no_local_procs=2
worker_pool=addprocs(no_local_procs, topology=:master_worker)

@info "Loading libraries everywhere..."
@everywhere using BioMotifInference, Random
Random.seed!(myid()*10000)

@info "Assembling uninformative source priors..."
sib_ui_sp = BioMotifInference.assemble_source_priors(no_sources, Vector{Matrix{Float64}}())
rys_ui_sp = BioMotifInference.assemble_source_priors(no_sources, Vector{Matrix{Float64}}())
combined_ui_sp = BioMotifInference.assemble_source_priors(no_sources, Vector{Matrix{Float64}}())

@info "Loading BGHMM likelihood matrix binaries..."
sib_matrix=deserialize(sib_diff_bg)
rys_matrix=deserialize(rys_diff_bg)
combined_matrix=deserialize(combined_diff_bg)

@info "Loading coded observation sets..."
sib_obs = deserialize(sib_code_binary)
rys_obs = deserialize(rys_code_binary)
combined_obs = deserialize(combined_code_binary)

@info "Assembling sib IPM ensemble on uninformative priors..."
isfile(string(sib_e_pth,'/',"ens")) ? (sib_e = deserialize(string(sib_e_pth,'/',"ens"))) :
    (sib_e = BioMotifInference.IPM_Ensemble(worker_pool, sib_e_pth, ensemble_size, sib_ui_sp, (falses(0,0), mixing_prior), sib_matrix, sib_obs, source_length_range, posterior_switch=true); serialize(string(sib_e_pth,'/',"ens"),sib_e))
sib_e=[]; Base.GC.gc();

@info "Assembling rys IPM ensemble on uninformative priors..."
isfile(string(rys_e_pth,'/',"ens")) ? (rys_e = deserialize(string(rys_e_pth,'/',"ens"))) :
    (rys_e = BioMotifInference.IPM_Ensemble(worker_pool, rys_e_pth, ensemble_size, rys_ui_sp, (falses(0,0), mixing_prior), rys_matrix, rys_obs, source_length_range, posterior_switch=true); serialize(string(rys_e_pth,'/',"ens"),rys_e))
rys_e=[]; Base.GC.gc();

@info "Assembling combined IPM ensemble on uninformative priors..."
isfile(string(combined_e_pth,'/',"ens")) ? (combined_e = deserialize(string(combined_e_pth,'/',"ens"))) :
    (combined_e = BioMotifInference.IPM_Ensemble(worker_pool, combined_e_pth, ensemble_size, combined_ui_sp, (falses(0,0), mixing_prior), combined_matrix, combined_obs, source_length_range, posterior_switch=true); serialize(string(combined_e_pth,'/',"ens"),combined_e))
combined_e=[]; Base.GC.gc();

rmprocs(worker_pool)

@info "Job done!"


text:
"""
    CompositeBasisFunction

Union of several basis functions in the imaginary-time/real-frequency domain
domains
"""
struct CompositeBasisFunction
    polys::Vector{Any}
end

"""
    (::CompositeBasisFunction)(x::Real)

Evaluate basis function at position `x`
"""
function (obj::CompositeBasisFunction)(x::Real)
    return vcat(p(x) for p in obj.polys)
end

function (obj::CompositeBasisFunction)(x::Vector{T}) where {T<:Real}
    return vcat((p(x) for p in obj.polys)...)
end

"""
    CompositeBasisFunctionFT

Union of several basis functions in the imaginary-frequency domain
"""
struct CompositeBasisFunctionFT
    polys::Vector{Any}
end

"""
Evaluate basis function at frequency n
"""
function (obj::CompositeBasisFunctionFT)(n::Union{Int,Vector{Int}})
    return hcat(p(n) for p in obj.polys)
end

struct CompositeBasis <: AbstractBasis
    beta::Float64
    bases::Vector{AbstractBasis}
    u::Union{CompositeBasisFunction,Nothing}
    v::Union{CompositeBasisFunction,Nothing}
    uhat::Union{CompositeBasisFunctionFT,Nothing}
end

iswellconditioned(basis::CompositeBasis) = false

function _collect_polys(::Type{T}, polys) where {T}
    if any(isnothing, polys)
        return nothing
    else
        return T(polys)
    end
end

function CompositeBasis(bases::Vector{AbstractBasis})
    u = CompositeBasisFunction([b.u for b in bases])
    v = _collect_polys(CompositeBasisFunction, [b.v for b in bases])
    uhat = _collect_polys(CompositeBasisFunctionFT, [b.uhat for b in bases])
    return CompositeBasis(getbeta(first(bases)), bases, u, v, uhat)
end

function default_tau_sampling_points(basis::CompositeBasis)
    return sort(unique(vcat((default_tau_sampling_points(b) for b in basis.bases)...)))
end

function default_matsubara_sampling_points(basis::CompositeBasis; mitigate=true)
    return sort(
        unique(
            vcat((default_matsubara_sampling_points(b; mitigate) for b in basis.bases)...)
        ),
    )
end


text:
@testset "FBA with crowding constraints" begin
    model = load_model(model_paths["e_coli_core.json"])
    rid_weight = Dict(
        rid => 0.004 for rid in reactions(model) if
        !looks_like_biomass_reaction(rid) && !looks_like_exchange_reaction(rid)
    )

    sol = flux_balance_analysis_dict(
        model,
        Tulip.Optimizer;
        modifications = [
            change_optimizer_attribute("IPM_IterationsLimit", 1000),
            add_crowding_constraints(rid_weight),
            change_constraint("EX_glc__D_e"; lb = -6),
        ],
    )

    @test isapprox(
        sol["BIOMASS_Ecoli_core_w_GAM"],
        0.491026987015203,
        atol = TEST_TOLERANCE,
    )

    @test isapprox(sol["EX_ac_e"], 0.7084745257320869, atol = TEST_TOLERANCE)
end


text:
# This file is part of the Julia package ModularForms.jl
#
# Copyright (c) 2018-2019: Lynn Engelberts and Alexandru Ghitza.

"""
   ModularForms

A Julia package for computing with modular forms.

"""
module ModularForms

using Nemo

export delta_poly
export delta_qexp
export delta_k_qexp
export eisenstein_series_poly
export eisenstein_series_qexp
export eta_quotient
export hecke_operator_on_basis
export hecke_operator_on_qexp
export poly_to_power_series
export prime_range
export dim_Sk
export dim_Mk
export victor_miller_basis

include("delta.jl")
include("eis_series.jl")
include("eta_quotient.jl")
include("hecke.jl")
include("poly_to_power_series.jl")
include("prime_range.jl")
include("vm_basis.jl")

export func

"""
    func(x)

Return func
"""
func(x) = 2x + 1

end # module


text:
function assemble_system!(problem::AllenCahnProblem1D{TLBC, TRBC, TA}) where {TLBC <: NeumannBC, TRBC <: NeumannBC, TA <: SparseMatrixCSC}
    # Construct first row corresponding to uₓ(a) = vₗ
    problem.A[1,1] = -2
    problem.A[1,2] = 2

    # Construct rows corresponding to interior mesh points
    for i ∈ 2:problem.nₓ + 1
        problem.A[i, i - 1] = 1
        problem.A[i,i] = -2
        problem.A[i, i + 1] = 1
    end

    # Construct last row corresponding to uₓ(b) = vᵣ
    problem.A[end, end] = -2
    problem.A[end, end - 1] = 2

    @. problem.A *= 1/(problem.Δx)^2

    return problem
end

function assemble_system!(problem::AllenCahnProblem1D{TLBC, TRBC, TA}) where {TLBC <: PeriodicBC, TRBC <: PeriodicBC, TA <: SparseMatrixCSC}
    # Construct first row corresponding to u₀
    problem.A[1,1] = -2
    problem.A[1,2] = 1
    problem.A[1, end] = 1

    # Construct rows corresponding to interior mesh points
    for i ∈ 2:problem.nₓ
        problem.A[i, i - 1] = 1
        problem.A[i,i] = -2
        problem.A[i, i + 1] = 1
    end

    # Construct last row corresponding to uₙ
    problem.A[end, 1] = 1
    problem.A[end, end - 1] = 1
    problem.A[end, end] = -2

    @. problem.A *= 1/(problem.Δx)^2

    return problem
end


text:
# gradient estimation methods
## SGD

function update_gradient_sgd!(∇, recvbufs, epoch::Integer, repochs::Vector{<:Integer}; state=nothing, nreplicas, ncolumns::Integer, kwargs...)
    length(recvbufs) == length(repochs) || throw(DimensionMismatch("recvbufs has dimension $(length(recvbufs)), but repochs has dimension $(length(repochs))"))
    0 < nreplicas || throw(DomainError(nreplicas, "nreplicas must be positive"))
    epoch <= 1 || !isnothing(state) || error("expected state to be initiated for epoch > 1")
    nworkers = length(recvbufs)
    mod(nworkers, nreplicas) == 0 || throw(ArgumentError("nworkers must be divisible by nreplicas"))

    # record the epoch at which each partition was last updated
    if isnothing(state)
        uepochs = Dict{Int,Int}()
    else
        uepochs::Dict{Int,Int} = state
    end

    # add at most 1 replica of each partition to the overall gradient
    # the partitions are arranged sequentially, so if there are 2 partitions and 3 replicas, then
    # Vs is of length 6, and its elements correspond to partitions [1, 1, 1, 2, 2, 2]
    ∇ .= 0
    processed_nsamples = 0
    for worker_index in 1:nworkers

        # skip workers that we've never received anything from
        if repochs[worker_index] == 0
            continue
        end

        metadata = reinterpret(UInt16, metadata_view(recvbufs[worker_index]))
        if length(metadata) != 4
            @error "received incorrectly formatted metadata from the $(worker_index)-th worker in epoch $epoch: $metadata"
            continue
        end
        canary, worker_rank, nsubpartitions, subpartition_index = metadata
        if  canary != CANARY_VALUE
            @error "recieved incorrect canary value from the $(worker_index)-th worker in epoch $epoch: $canary"
            continue
        end
        if worker_rank != worker_index
            @error "unexpected rank for the $(worker_index)-th worker in epoch $epoch: $worker_rank"
            continue
        end
        if !(0 < subpartition_index <= nsubpartitions)
            @error "received incorrect sub-partition index from the $(worker_index)-th worker in epoch $epoch: $subpartition_index "
            continue
        end
        replica_index = ceil(Int, worker_index/nreplicas)
        partition_index = (replica_index-1)*nsubpartitions + subpartition_index

        # don't do anything if we didn't receive from this worker this epoch,
        # or if we've already updated that partition this epoch
        if repochs[worker_index] < epoch || (haskey(uepochs, partition_index) && uepochs[partition_index] == epoch)
            continue
        end

        # compute the number of samples that make up this partition, and add it to the total
        worker_nsamples = length(partition(ncolumns, div(nworkers, nreplicas), replica_index))
        subpartition_nsamples = length(partition(worker_nsamples, nsubpartitions, subpartition_index))
        processed_nsamples += subpartition_nsamples

        # add the sub-gradient computed by this worker
        uepochs[partition_index] = epoch
        Vw = reshape(data_view(recvbufs[worker_index]), size(∇)...)
        ∇ .+= Vw
    end

    # scale the (stochastic) gradient to make it unbiased estimate of the true gradient
    processed_fraction = processed_nsamples / ncolumns
    if !(processed_fraction ≈ 1)
        ∇ ./= processed_fraction
    end
    uepochs
end

## variance-reduced (table-based)

function update_gradient_vr!(∇, recvbufs, epoch::Integer, repochs::Vector{<:Integer}; state=nothing, nostale::Bool, nreplicas::Integer, nsubpartitions::Integer, kwargs...)
    length(recvbufs) == length(repochs) || throw(DimensionMismatch("recvbufs has dimension $(length(recvbufs)), but repochs has dimension $(length(repochs))"))
    0 < nreplicas || throw(DomainError(nreplicas, "nreplicas must be positive"))
    epoch <= 1 || !isnothing(state) || error("expected state to be initiated for epoch > 1")
    nworkers = length(recvbufs)
    mod(nworkers, nreplicas) == 0 || throw(ArgumentError("nworkers must be divisible by nreplicas"))
    npartitions = div(nworkers, nreplicas) * nsubpartitions

    # record the epoch at which each partition was last updated
    # store the previously computed partial gradients
    if isnothing(state)
        uepochs = zeros(Int, npartitions)
        worker_from_partition = zeros(Int, npartitions)
        sg = StochasticGradient(eltype(∇), npartitions, size(∇))
    else
        uepochs, worker_from_partition, sg = state
        worker_from_partition .= 0
    end

    # iterate over the received partial gradients to record which of them are newer than what we currently have
    for worker_index in 1:nworkers

        # skip workers that we've never received anything from
        if repochs[worker_index] == 0
            continue
        end

        metadata = reinterpret(UInt16, metadata_view(recvbufs[worker_index]))
        if length(metadata) != 4
            @error "received incorrectly formatted metadata from the $(worker_index)-th worker in epoch $epoch: $metadata"
            continue
        end
        canary, worker_rank, worker_nsubpartitions, subpartition_index = metadata
        if  canary != CANARY_VALUE
            @error "recieved incorrect canary value from the $(worker_index)-th worker in epoch $epoch: $canary"
            continue
        end
        worker_nsubpartitions == nsubpartitions || throw(ArgumentError("update_gradient_vr doesn't support dynamically changing the number of partitions"))
        if worker_rank != worker_index
            @error "unexpected rank for the $(worker_index)-th worker in epoch $epoch: $worker_rank"
            continue
        end
        if !(0 < subpartition_index <= nsubpartitions)
            @error "received incorrect sub-partition index from the $(worker_index)-th worker in epoch $epoch: $subpartition_index "
            continue
        end
        replica_index = ceil(Int, worker_index/nreplicas)
        partition_index = (replica_index-1)*nsubpartitions + subpartition_index

        # discard stale gradients if the nostale option is set
        if nostale && repochs[worker_index] != epoch
            continue
        end

        # skip updates that contain no new information
        if repochs[worker_index] <= uepochs[partition_index]
            continue
        end

        # record which worker has updates for which partition and how new it is
        worker_from_partition[partition_index] = worker_index
        uepochs[partition_index] = repochs[worker_index]
    end

    # update the gradient with the received sub-gradients
    partition_indices = [i for i in 1:npartitions if !iszero(worker_from_partition[i])]
    worker_indices = [worker_from_partition[i] for i in partition_indices]
    ∇s = [reshape(data_view(recvbufs[worker_index]), size(∇)...) for worker_index in worker_indices]
    update!(sg, zip(partition_indices, ∇s))

    # scale the gradient by the number of initialized sub-gradients
    f = initialized_fraction(sg)
    if isone(f)
        ∇ .= sg
    else
        ∇ .= sg ./ f
    end

    uepochs, worker_from_partition, sg
end

## variance-reduced (tree-based)

function setup_gradient_vrt(V, ;ncolumns::Integer, kwargs...)
    ∇i = zero(V)
    tg = TreeGradient(zero(V), ncolumns)
    ∇i, tg
end

function recvf_vrt(worker_index::Integer, epoch::Integer, repoch::Integer, recvbuf; state, nostale::Bool, nreplicas::Integer, ncolumns::Integer, nworkers::Integer, kwargs...)
    !isnothing(state) || throw(ArgumentError("state was not initialized"))
    gradient_state, _ = state
    ∇i, tg = gradient_state

    metadata = reinterpret(UInt16, metadata_view(recvbuf))
    canary, worker_rank, nsubpartitions, subpartition_index = metadata
    if  canary != CANARY_VALUE
        error("recieved incorrect canary value from the $(worker_index)-th worker in epoch $epoch: $canary")
    end
    if worker_rank != worker_index
        error("unexpected rank for the $(worker_index)-th worker in epoch $epoch: $worker_rank")
    end
    if !(0 < subpartition_index <= nsubpartitions)
        error("received incorrect sub-partition index from the $(worker_index)-th worker in epoch $epoch: $subpartition_index")
    end

    # discard stale gradients if the nostale option is set
    if nostale && repoch != epoch
        return
    end

    # compute which samples make up this partition
    replica_index = ceil(Int, worker_index/nreplicas)
    worker_samples = partition(ncolumns, div(nworkers, nreplicas), replica_index)
    worker_nsamples = length(worker_samples)
    subpartition_samples = first(worker_samples) .+ partition(worker_nsamples, nsubpartitions, subpartition_index) .- 1

    # insert the new gradient
    ∇i .= reshape(data_view(recvbuf), size(∇i)...)
    insert!(tg, first(subpartition_samples), last(subpartition_samples), ∇i)
    return
end

function update_gradient_vrt!(∇, recvbufs, epoch::Integer, repochs::Vector{<:Integer}; state, kwargs...)
    !isnothing(state) || throw(ArgumentError("state was not initialized"))
    ∇i, tg = state

    # scale the gradient by the number of initialized sub-gradients
    fraction_processed = tg.ninit / tg.n
    if isapprox(fraction_processed, 1)
        ∇ .= tg.∇
    else
        ∇ .= tg.∇ ./ fraction_processed
    end
    state
end


text:
module TestRand

using Phylo
using DataFrames

using Compat.Test

@testset "Nonultrametric()" begin
    # Create a 10 tip tree
    nu = Nonultrametric(10)
    @test eltype(nu) == NamedTree
    @test validate(rand(nu))
    @test Set(getleafnames(rand(nu))) == Set(getleafnames(rand(nu)))
    # Create a tree with named tips
    species = ["Dog", "Cat", "Human"]
    t = rand(Nonultrametric(species))
    @test ntrees(t) == 1
    @test validate(t)
    @test Set(getleafnames(t)) == Set(species)

    t2 = rand(Nonultrametric{BinaryTree{DataFrame, Vector{Float64}}}(species))
    @test length(getnoderecord(t2, species[1])) == 0
    t3 = rand(Nonultrametric{BinaryTree{DataFrame, Vector{String}}}(species))
    map(nodenameiter(t3)) do name
        setnoderecord!(t3, name, nodehistory(t3, name))
    end
    for name in nodenameiter(t3)
        @test all(getnoderecord(t3, name) .== nodehistory(t3, name))
    end
end

@testset "Ultrametric()" begin
    # Create a 10 tip tree
    u = Ultrametric(10)
    @test validate(rand(u))
    @test Set(getleafnames(rand(u))) == Set(getleafnames(rand(u)))
    tree = rand(u)
    @test ntrees(tree) == 1
    heights = map(x -> getheight(tree, x), getleafnames(tree))
    @test all(h -> h ≈ heights[1], heights)
    species = ["Dog", "Cat", "Human"]
    t = rand(Ultrametric(species))
    @test validate(t)
    @test Set(getleafnames(t)) == Set(species)

    numnodes = 50
    ul = Ultrametric{BinaryTree{DataFrame, Vector{Float64}}}(numnodes)
    u2 = rand(ul)
    @test eltype(ul) == BinaryTree{DataFrame, Vector{Float64}}
    @test length(nodefilter(isinternal, u2)) == numnodes - 2
end

@testset "TestSets" begin
    @test ntrees(rand(Ultrametric(10), 10)) == 10
    @test length(rand(Ultrametric(10), 10)) == 10
    names = ["One", "Two", "Three"]
    ts = rand(Nonultrametric(names), 20)
    @test length(ts) == 20
    @test ntrees(ts) == 20
    @test length(getleafnames(ts)) == length(names)
    @test length(getnodenames(ts)) == length(names) * 2 - 1
end
end


text:
@inline function discrete_derivatives(x::AbstractArray{T} where T<:Real;
        d::Int=1)
    @assert d ≥ 0 "discrete derivatives parameter must be greater than 0 (now it is $d)"

    res = zeros(eltype(x), length(x), d+1)

    @inbounds @simd for i in 1:length(x)
        res[i,1] = x[i]
    end

    @inbounds @simd for j in 2:(d+1)
        @inbounds @simd for i in 1:(length(x)-j+1)
            res[i,j] = res[i+1,j-1] - res[i,j-1]
        end
    end

    return res
end

@inline function discrete_derivative(x::AbstractArray{T} where T<:Real; d::Int=1)
    derivatives = discrete_derivatives(x, d=d)
    return derivatives[1:length(x)-d,d+1]
end


text:
function memalloc(bytesize::Int)
    ref = Ref{Cptr}()
    @apicall :cuMemAlloc (Ptr{Ptr{Cvoid}},Csize_t) ref bytesize
    ref[]
end

function memfree(ptr::Ptr)
    @apicall :cuMemFree (Ptr{Cvoid},) ptr
end

function meminfo()
    ref_free = Ref{Csize_t}()
    ref_total = Ref{Csize_t}()
    @apicall :cuMemGetInfo (Ptr{Csize_t},Ptr{Csize_t}) ref_free ref_total
    Int(ref_free[]), Int(ref_total[])
end


text:
"""
    install_with_pip()

    Uses the `Conda.pip` method to install
    napari using pyqt5
"""
function install_with_pip()
    # https://napari.org/docs/0.4.0/
    Conda.pip_interop(true)
    Conda.pip("install","napari[pyqt5]")
end

"""
    install_with_conda()

    Uses `Conda.add` to install napari from
    conda-forge
"""
function install_with_conda()
    # https://anaconda.org/conda-forge/napari
    Conda.add(["pyqt", "napari"]; channel = "conda-forge")
end


text:
export launch!

"""
    launch!(rocket::Rocket, Δt::Real) -> Float64
Simulate a rocket launch with with time steps `Δt` return max altitude
and time it was achieved
"""
function launch!(rocket::Rocket, Δt::Real)
    g = 9.80665  # acceleration caused by gravity
    v = 0        # velocity
    h = 0        # altitude

    for t in 0:Δt:1000
        m = mass(rocket)
        F = thrust(rocket) - m*g

        remaining = update!(rocket, t, Δt)

        # Any propellant and thrust left?
        if remaining == 0 || F <= 0
            return (t, h)
        end

        h += v*Δt
        a = F/m
        v += a*Δt
    end
end

"""
    launch!(callback::Function, rocket::Rocket, Δt::Real) -> Float64
Calls `callback(rocket, t, v, h)` on every iteration,
where `v` and `h` is the velocity and altitude at time `t`.
"""
function launch!(callback::Function, rocket::Rocket, Δt::Real)
    g = 9.80665  # acceleration caused by gravity
    v = 0        # velocity
    h = 0        # altitude

    for t in 0:Δt:1000
        callback(rocket, t, v, h)
        m = mass(rocket)
        F = thrust(rocket) - m*g

        remaining = update!(rocket, t, Δt)

        # Any propellant and thrust left?
        if remaining == 0 || F <= 0
            return (t, h)
        end

        h += v*Δt
        a = F/m
        v += a*Δt
    end
    callback(rocket, t, v, h)
    nothing
end


text:
using RNGPool

include("demo.jl")
include("lgModel.jl")

setRNGs(12345)

println("Using a linear Gaussian Model with n = 100, obsVariance = 1.0")
model, lM, _ = setupLGModel(100, 1.0)
visualizeCCSMC(model, lM, 128, :BS, true, false)
visualizeCCSMC(model, 128, true, false)
visualizeCCSMC(model, lM, 128, :AS, true, false)

visualizeCCSMC(model, lM, 128, :BS, true, true)
visualizeCCSMC(model, 128, true, true)
visualizeCCSMC(model, lM, 128, :AS, true, true)

println("\nUsing a linear Gaussian Model with n = 1000, obsVariance = 100.0")
model, lM, _ = setupLGModel(1000, 100.0)
visualizeCCSMC(model, lM, 512, :BS, true, false)
visualizeCCSMC(model, 512, true, false, 100)
visualizeCCSMC(model, lM, 512, :AS, true, false, 100)

visualizeCCSMC(model, lM, 512, :BS, true, true)
visualizeCCSMC(model, 512, true, true, 100)
visualizeCCSMC(model, lM, 512, :AS, true, true, 100)


text:
# Code Generator
# ==============

"""
Variable names used in generated code.

The following variable names may be used in the code.

- `p::Int`: current position of data
- `p_end::Int`: end position of data
- `p_eof::Int`: end position of file stream
- `ts::Int`: start position of token (tokenizer only)
- `te::Int`: end position of token (tokenizer only)
- `cs::Int`: current state
- `data::Any`: input data
- `mem::SizedMemory`: input data memory
- `byte::UInt8`: current data byte
"""
struct Variables
    p::Symbol
    p_end::Symbol
    p_eof::Symbol
    ts::Symbol
    te::Symbol
    cs::Symbol
    data::Symbol
    mem::Symbol
    byte::Symbol
end

struct CodeGenContext
    vars::Variables
    generator::Function
    checkbounds::Bool
    loopunroll::Int
    getbyte::Function
    clean::Bool
end

"""
    CodeGenContext(;
        vars=Variables(:p, :p_end, :p_eof, :ts, :te, :cs, :data, gensym(), gensym()),
        generator=:table,
        checkbounds=true,
        loopunroll=0,
        getbyte=Base.getindex,
        clean=false
    )

Create a code generation context.

Arguments
---------

- `vars`: variable names used in generated code
- `generator`: code generator (`:table`, `:inline` or `:goto`)
- `checkbounds`: flag of bounds check
- `loopunroll`: loop unroll factor (≥ 0)
- `getbyte`: function of byte access (i.e. `getbyte(data, p)`)
- `clean`: flag of code cleansing
"""
function CodeGenContext(;
        vars::Variables=Variables(:p, :p_end, :p_eof, :ts, :te, :cs, :data, gensym(), gensym()),
        generator::Symbol=:table,
        checkbounds::Bool=true,
        loopunroll::Integer=0,
        getbyte::Function=Base.getindex,
        clean::Bool=false)
    if loopunroll < 0
        throw(ArgumentError("loop unroll factor must be a non-negative integer"))
    elseif loopunroll > 0 && generator != :goto
        throw(ArgumentError("loop unrolling is not supported for $(generator)"))
    end
    # check generator
    if generator == :table
        generator = generate_table_code
    elseif generator == :inline
        generator = generate_inline_code
    elseif generator == :goto
        generator = generate_goto_code
    else
        throw(ArgumentError("invalid code generator: $(generator)"))
    end
    return CodeGenContext(vars, generator, checkbounds, loopunroll, getbyte, clean)
end

"""
    generate_init_code(context::CodeGenContext, machine::Machine)::Expr

Generate variable initialization code.
"""
function generate_init_code(ctx::CodeGenContext, machine::Machine)
    return quote
        $(ctx.vars.p)::Int = 1
        $(ctx.vars.p_end)::Int = 0
        $(ctx.vars.p_eof)::Int = -1
        $(ctx.vars.cs)::Int = $(machine.start_state)
    end
end

"""
    generate_exec_code(ctx::CodeGenContext, machine::Machine, actions=nothing)::Expr

Generate machine execution code with actions.
"""
function generate_exec_code(ctx::CodeGenContext, machine::Machine, actions=nothing)
    # make actions
    if actions == nothing
        actions = Dict{Symbol,Expr}()
    elseif actions == :debug
        actions = debug_actions(machine)
    elseif isa(actions, AbstractDict{Symbol,Expr})
        actions = Dict{Symbol,Expr}(collect(actions))
    else
        throw(ArgumentError("invalid actions argument"))
    end
    # generate code
    code = ctx.generator(ctx, machine, actions)
    if ctx.clean
        code = cleanup(code)
    end
    return code
end

function generate_table_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    action_dispatch_code, set_act_code = generate_action_dispatch_code(ctx, machine, actions)
    trans_table = generate_transition_table(machine)
    getbyte_code = generate_geybyte_code(ctx)
    set_cs_code = :(@inbounds $(ctx.vars.cs) = $(trans_table)[($(ctx.vars.cs) - 1) << 8 + $(ctx.vars.byte) + 1])
    eof_action_code = generate_eof_action_code(ctx, machine, actions)
    return quote
        $(ctx.vars.mem) = $(SizedMemory)($(ctx.vars.data))
        while $(ctx.vars.p) ≤ $(ctx.vars.p_end) && $(ctx.vars.cs) > 0
            $(getbyte_code)
            $(set_act_code)
            $(set_cs_code)
            $(action_dispatch_code)
            $(ctx.vars.p) += 1
        end
        if $(ctx.vars.p) > $(ctx.vars.p_eof) ≥ 0 && $(ctx.vars.cs) ∈ $(machine.final_states)
            $(eof_action_code)
            $(ctx.vars.cs) = 0
        elseif $(ctx.vars.cs) < 0
            $(ctx.vars.p) -= 1
        end
    end
end

function generate_transition_table(machine::Machine)
    trans_table = Matrix{Int}(undef, 256, length(machine.states))
    for j in 1:size(trans_table, 2)
        trans_table[:,j] .= -j
    end
    for s in traverse(machine.start), (e, t) in s.edges
        if !isempty(e.precond)
            error("precondition is not supported in the table-based code generator; try code=:inline or :goto")
        end
        for l in e.labels
            trans_table[l+1,s.state] = t.state
        end
    end
    return trans_table
end

function generate_action_dispatch_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    action_table = fill(0, (256, length(machine.states)))
    action_ids = Dict{Vector{Symbol},Int}()
    for s in traverse(machine.start)
        for (e, t) in s.edges
            if isempty(e.actions)
                continue
            end
            id = get!(action_ids, action_names(e.actions), length(action_ids) + 1)
            for l in e.labels
                action_table[l+1,s.state] = id
            end
        end
    end
    act = gensym()
    default = :()
    action_dispatch_code = foldr(default, action_ids) do names_id, els
        names, id = names_id
        action_code = rewrite_special_macros(ctx, generate_action_code(names, actions), false)
        return Expr(:if, :($(act) == $(id)), action_code, els)
    end
    action_code = :(@inbounds $(act) = $(action_table)[($(ctx.vars.cs) - 1) << 8 + $(ctx.vars.byte) + 1])
    return action_dispatch_code, action_code
end

function generate_inline_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    trans_code = generate_transition_code(ctx, machine, actions)
    eof_action_code = generate_eof_action_code(ctx, machine, actions)
    getbyte_code = generate_geybyte_code(ctx)
    return quote
        $(ctx.vars.mem) = $(SizedMemory)($(ctx.vars.data))
        while $(ctx.vars.p) ≤ $(ctx.vars.p_end) && $(ctx.vars.cs) > 0
            $(getbyte_code)
            $(trans_code)
            $(ctx.vars.p) += 1
        end
        if $(ctx.vars.p) > $(ctx.vars.p_eof) ≥ 0 && $(ctx.vars.cs) ∈ $(machine.final_states)
            $(eof_action_code)
            $(ctx.vars.cs) = 0
        elseif $(ctx.vars.cs) < 0
            $(ctx.vars.p) -= 1
        end
    end
end

function generate_transition_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    default = :($(ctx.vars.cs) = -$(ctx.vars.cs))
    return foldr(default, traverse(machine.start)) do s, els
        then = foldr(default, optimize_edge_order(s.edges)) do edge, els′
            e, t = edge
            action_code = rewrite_special_macros(ctx, generate_action_code(e.actions, actions), false)
            then′ = :($(ctx.vars.cs) = $(t.state); $(action_code))
            return Expr(:if, generate_condition_code(ctx, e, actions), then′, els′)
        end
        return Expr(:if, state_condition(ctx, s.state), then, els)
    end
end

function generate_goto_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    actions_in = Dict{Node,Set{Vector{Symbol}}}()
    for s in traverse(machine.start), (e, t) in s.edges
        push!(get!(actions_in, t, Set{Vector{Symbol}}()), action_names(e.actions))
    end
    action_label = Dict{Node,Dict{Vector{Symbol},Symbol}}()
    for s in traverse(machine.start)
        action_label[s] = Dict()
        if haskey(actions_in, s)
            for (i, names) in enumerate(actions_in[s])
                action_label[s][names] = Symbol("state_", s.state, "_action_", i)
            end
        end
    end

    blocks = Expr[]
    for s in traverse(machine.start)
        block = Expr(:block)
        for (names, label) in action_label[s]
            if isempty(names)
                continue
            end
            append_code!(block, quote
                @label $(label)
                $(rewrite_special_macros(ctx, generate_action_code(names, actions), false, s.state))
                @goto $(Symbol("state_", s.state))
            end)
        end
        append_code!(block, quote
            @label $(Symbol("state_", s.state))
            $(ctx.vars.p) += 1
            if $(ctx.vars.p) > $(ctx.vars.p_end)
                $(ctx.vars.cs) = $(s.state)
                @goto exit
            end
        end)
        default = :($(ctx.vars.cs) = $(-s.state); @goto exit)
        dispatch_code = foldr(default, optimize_edge_order(s.edges)) do edge, els
            e, t = edge
            if isempty(e.actions)
                if ctx.loopunroll > 0 && s.state == t.state && length(e.labels) ≥ 4
                    then = generate_unrolled_loop(ctx, e, t)
                else
                    then = :(@goto $(Symbol("state_", t.state)))
                end
            else
                then = :(@goto $(action_label[t][action_names(e.actions)]))
            end
            return Expr(:if, generate_condition_code(ctx, e, actions), then, els)
        end
        append_code!(block, quote
            @label $(Symbol("state_case_", s.state))
            $(generate_geybyte_code(ctx))
            $(dispatch_code)
        end)
        push!(blocks, block)
    end

    enter_code = foldr(:(@goto exit), machine.states) do s, els
        return Expr(:if, :($(ctx.vars.cs) == $(s)), :(@goto $(Symbol("state_case_", s))), els)
    end

    eof_action_code = rewrite_special_macros(ctx, generate_eof_action_code(ctx, machine, actions), true)

    return quote
        if $(ctx.vars.p) > $(ctx.vars.p_end)
            @goto exit
        end
        $(ctx.vars.mem) = $(SizedMemory)($(ctx.vars.data))
        $(enter_code)
        $(Expr(:block, blocks...))
        @label exit
        if $(ctx.vars.p) > $(ctx.vars.p_eof) ≥ 0 && $(ctx.vars.cs) ∈ $(machine.final_states)
            $(eof_action_code)
            $(ctx.vars.cs) = 0
        end
    end
end

function append_code!(block::Expr, code::Expr)
    @assert block.head == :block
    @assert code.head == :block
    append!(block.args, code.args)
    return block
end

function generate_unrolled_loop(ctx::CodeGenContext, edge::Edge, t::Node)
    # Generated code looks like this (when unroll=2):
    #     while p + 2 ≤ p_end
    #         l1 = $(getbyte)(data, p + 1)
    #         !$(generate_membership_code(:l1, e.labels)) && break
    #         l2 = $(getbyte)(data, p + 2)
    #         !$(generate_membership_code(:l2, e.labels)) && break
    #         p += 2
    #     end
    #     @goto ...
    @assert ctx.loopunroll > 0
    body = :(begin end)
    for k in 1:ctx.loopunroll
        l = Symbol(ctx.vars.byte, k)
        push!(
            body.args,
            quote
                $(generate_geybyte_code(ctx, l, k))
                $(generate_membership_code(l, edge.labels)) || begin
                    $(ctx.vars.p) += $(k-1)
                    break
                end
            end)
    end
    push!(body.args, :($(ctx.vars.p) += $(ctx.loopunroll)))
    quote
        while $(ctx.vars.p) + $(ctx.loopunroll) ≤ $(ctx.vars.p_end)
            $(body)
        end
        @goto $(Symbol("state_", t.state))
    end
end

function generate_eof_action_code(ctx::CodeGenContext, machine::Machine, actions::Dict{Symbol,Expr})
    return foldr(:(), machine.eof_actions) do s_as, els
        s, as = s_as
        action_code = rewrite_special_macros(ctx, generate_action_code(action_names(as), actions), true)
        Expr(:if, state_condition(ctx, s), action_code, els)
    end
end

function generate_action_code(list::ActionList, actions::Dict{Symbol,Expr})
    return generate_action_code(action_names(list), actions)
end

function generate_action_code(names::Vector{Symbol}, actions::Dict{Symbol,Expr})
    return Expr(:block, (actions[n] for n in names)...)
end

function generate_geybyte_code(ctx::CodeGenContext)
    return generate_geybyte_code(ctx, ctx.vars.byte, 0)
end

function generate_geybyte_code(ctx::CodeGenContext, varbyte::Symbol, offset::Int)
    code = :($(varbyte) = $(ctx.getbyte)($(ctx.vars.mem), $(ctx.vars.p) + $(offset)))
    if !ctx.checkbounds
        code = :(@inbounds $(code))
    end
    return code
end

function state_condition(ctx::CodeGenContext, s::Int)
    return :($(ctx.vars.cs) == $(s))
end

function generate_condition_code(ctx::CodeGenContext, edge::Edge, actions::Dict{Symbol,Expr})
    labelcode = generate_membership_code(ctx.vars.byte, edge.labels)
    precondcode = foldr(:(true), edge.precond) do p, ex
        name, value = p
        if value == BOTH
            ex1 = :(true)
        elseif value == TRUE
            ex1 = :( $(actions[name]))
        elseif value == FALSE
            ex1 = :(!$(actions[name]))
        else
            ex1 = :(false)
        end
        return Expr(:&&, ex1, ex)
    end
    return :($(labelcode) && $(precondcode))
end

function generate_membership_code(var::Symbol, set::ByteSet)
    min, max = minimum(set), maximum(set)
    @assert min isa UInt8 && max isa UInt8
    if max - min + 1 == length(set)
        # contiguous
        if min == max
            return :($(var) == $(min))
        else
            return :($(var) in $(min:max))
        end
    elseif max - min + 1 ≤ 64 && all(b - min ≥ max for b in 0x00:0xff if b < min)
        # storable in a 64-bit bitmap
        bitmap = UInt64(0)
        for x in set
            bitmap |= UInt64(1) << (x - min)
        end
        return :(($(UInt64(1)) << ($(var) - $(min))) & $(bitmap) != 0)
    else
        # fallback
        return foldr((range, cond) -> Expr(:||, :($(var) in $(range)), cond),
                     :(false),
                     sort(range_encode(set), by=length, rev=true))
    end
end

# Used by the :table and :inline code generators.
function rewrite_special_macros(ctx::CodeGenContext, ex::Expr, eof::Bool)
    args = []
    for arg in ex.args
        if isescape(arg)
            if !eof
                push!(args, quote
                    $(ctx.vars.p) += 1
                    break
                end)
            end
        elseif isa(arg, Expr)
            push!(args, rewrite_special_macros(ctx, arg, eof))
        else
            push!(args, arg)
        end
    end
    return Expr(ex.head, args...)
end

# Used by the :goto code generator.
function rewrite_special_macros(ctx::CodeGenContext, ex::Expr, eof::Bool, cs::Int)
    args = []
    for arg in ex.args
        if isescape(arg)
            if !eof
                push!(args, quote
                    $(ctx.vars.cs) = $(cs)
                    $(ctx.vars.p) += 1
                    @goto exit
                end)
            end
        elseif isa(arg, Expr)
            push!(args, rewrite_special_macros(ctx, arg, eof, cs))
        else
            push!(args, arg)
        end
    end
    return Expr(ex.head, args...)
end

function isescape(arg)
    return arg isa Expr && arg.head == :macrocall && arg.args[1] == Symbol("@escape")
end

function cleanup(ex::Expr)
    args = []
    for arg in ex.args
        if isa(arg, Expr)
            if arg.head == :line
                # pass
            elseif ex.head == :block && arg.head == :block
                append!(args, cleanup(arg).args)
            else
                push!(args, cleanup(arg))
            end
        else
            push!(args, arg)
        end
    end
    return Expr(ex.head, args...)
end

function debug_actions(machine::Machine)
    actions = Set{Symbol}()
    for s in traverse(machine.start)
        for (e, t) in s.edges
            union!(actions, a.name for a in e.actions)
        end
    end
    for as in values(machine.eof_actions)
        union!(actions, a.name for a in as)
    end
    function log_expr(name)
        return :(push!(logger, $(QuoteNode(name))))
    end
    return Dict{Symbol,Expr}(name => log_expr(name) for name in actions)
end

# Sort edges by its size in descending order.
function optimize_edge_order(edges)
    return sort!(copy(edges), by=e->length(e[1].labels), rev=true)
end

# Generic foldr.
function foldr(op::Function, x0, xs)
    function rec(xs, s)
        if s == nothing
            return x0
        else
            return op(s[1], rec(xs, iterate(xs, s[2])))
        end
    end
    return rec(xs, iterate(xs))
end


text:
"""
    IndicatorActiveOnFalseBridge{T}

The `IndicatorActiveOnFalseBridge` replaces an indicator constraint activated
on 0 with a variable ``z_0`` with the constraint activated on 1, with a variable ``z_1``.
It stores the added `variable_index` and added constraints:
- ``z_1 \\in \\mathbb{B}`` in `zero_one_cons`
- ``z_0 + z_1 == 1`` in `` in `disjunction_cons`
- The added `ACTIVATE_ON_ONE` indicator constraint in `indicator_cons_index`.
"""
struct IndicatorActiveOnFalseBridge{
    T,
    F<:MOI.AbstractVectorFunction,
    S<:MOI.AbstractScalarSet,
} <: AbstractBridge
    variable_index::MOI.VariableIndex
    zero_one_cons::MOI.ConstraintIndex{MOI.SingleVariable,MOI.ZeroOne}
    disjunction_cons::MOI.ConstraintIndex{
        MOI.ScalarAffineFunction{T},
        MOI.EqualTo{T},
    }
    indicator_cons_index::MOI.ConstraintIndex{
        F,
        MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE,S},
    }
end

function bridge_constraint(
    ::Type{IndicatorActiveOnFalseBridge{T,F,S}},
    model::MOI.ModelLike,
    f::MOI.VectorAffineFunction{T},
    s::IS,
) where {
    S<:MOI.AbstractScalarSet,
    T<:Real,
    F,
    IS<:MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO,S},
}
    f_scalars = MOIU.eachscalar(f)
    z2, zo_cons = MOI.add_constrained_variable(model, MOI.ZeroOne())
    # z1 + z2 == 1
    z1_z2 = MOIU.operate(+, T, f_scalars[1], MOI.SingleVariable(z2))
    dcons = MOIU.normalize_and_add_constraint(model, z1_z2, MOI.EqualTo(one(T)))
    f2 = MOIU.operate(vcat, T, MOI.SingleVariable(z2), f_scalars[2])
    ci = MOI.add_constraint(
        model,
        f2,
        MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(s.set),
    )
    return IndicatorActiveOnFalseBridge{T,F,S}(z2, zo_cons, dcons, ci)
end

function MOI.supports_constraint(
    ::Type{<:IndicatorActiveOnFalseBridge{T}},
    ::Type{<:MOI.VectorAffineFunction},
    ::Type{<:MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(
    ::Type{<:IndicatorActiveOnFalseBridge},
)
    return [(MOI.ZeroOne,)]
end
function MOIB.added_constraint_types(
    ::Type{IndicatorActiveOnFalseBridge{T,F,S}},
) where {T,F,S}
    return [
        (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
        (F, MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE,S}),
    ]
end

function concrete_bridge_type(
    ::Type{<:IndicatorActiveOnFalseBridge{T}},
    ::Type{F},
    ::Type{MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO,S}},
) where {T,F<:MOI.VectorAffineFunction,S<:MOI.AbstractScalarSet}
    return IndicatorActiveOnFalseBridge{T,F,S}
end

function concrete_bridge_type(
    ::Type{<:IndicatorActiveOnFalseBridge},
    ::Type{F},
    ::Type{MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO,S}},
) where {F<:MOI.VectorAffineFunction,S<:MOI.AbstractScalarSet}
    return IndicatorActiveOnFalseBridge{Float64,F,S}
end