text:
using Random

export PathPermutationMove

struct PathPermutationMove <: InvMcmcMove
    ν::Int
    counts::Vector{Int}  # Track acceptance
    function PathPermutationMove(;ν::Int=3)
        new(
            ν,
            [0,0]
        )
    end
end

Base.show(io::IO, x::PathPermutationMove) = print(io, "PathPermutationMove(ν=$(x.ν))")

# Exclusive shuffle (without identity permutation )
function excl_shuffle!(x)
    ind_ref = collect(1:length(x))
    ind = randperm(length(x))
    while ind == ind_ref
        shuffle!(ind)
    end
    permute!(x, ind)
end

function prop_sample!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    move::PathPermutationMove,
    pointers::InteractionSequence{Int},
    V::UnitRange
    )

    n = length(S_curr)
    n_var = n
    ν = move.ν
    k = rand(1:min(n,ν)) # Number of paths to permute
    i = 0
    # Now we use alg which samples random subseq (for index of paths to permute)
    while k > 0
        u = rand()
        q = (n_var - k) / n_var
        while q > u  # skip
            i += 1
            n_var -= 1
            q *= (n_var - k) / n_var
        end
        i+=1
        # i is now index of path to permute
        # @inbounds ind[j] = i
        @inbounds excl_shuffle!(S_prop[i])
        n_var -= 1
        k -= 1
    end
    log_ratio = 0.0 # In this case the proposal is symmetric
    return log_ratio

end


text:
export LazyΩ
using Random
using ..IDS, ..Util, ..Tagging, ..RNG, ..Space
import ..Basis: replacetags
import ..Var

"Lazily constructs randp, values as they are needed"
struct LazyΩ{TAGS <: Tags, T} <: AbstractΩ
  data::T
  tags::TAGS
end

# AbstractDict interface
Base.keys(ω::LazyΩ) = Base.keys(ω.data)
Base.values(ω::LazyΩ) = Base.values(ω.data)
Base.pairs(ω::LazyΩ) = Base.pairs(ω.data)
Base.getindex(ω::LazyΩ, id) = ω.data[id]

# Move to specific omega types
Base.merge!(ω::LazyΩ, ω_) = error("Unimplemented")
Basis.like(ω::LazyΩ{Tags, T}, kv::Pair) where {Tags, T} =
  LazyΩ{Tags, T}(T(kv), ω.tags)

## Constructors
const EmptyTags = Tags{(),Tuple{}}
LazyΩ{EmptyTags, T}() where {T} = LazyΩ(T(), Tags())
LazyΩ{EmptyTags, T}(data) where T = LazyΩ(T(data), Tags())

"Construct `LazyΩ` from `rng` -- `ω.data` will be generated from `rng`"
LazyΩ{T}(rng::AbstractRNG) where T = tagrng(LazyΩ{T}(), rng)

## Ids
Basis.idtype(::Type{Dict{T, V}}) where {T, V} = T
Basis.idtype(::Dict{T, V}) where {T, V} = T
Basis.idtype(ω::LazyΩ{TAGS, Dict{T, V}}) where {TAGS, T, V} = T
ids(ω::LazyΩ) = keys(ω.data)

## Tags
replacetags(ω::LazyΩ, tags) = LazyΩ(ω.data, tags)
traits(::Type{LazyΩ{TAGS, T}}) where {TAGS, T} = traits(TAGS)

Tagging.mergetag(ω::LazyΩ, tag) = LazyΩ(ω.data, merge(ω.tags, tag))

Tagging.tags(ω::LazyΩ) = ω.tags

Base.setindex!(ω::LazyΩ, value, id) =
  ω.data[id] = value

function Var.recurse(primrv::Var.PrimRandVar, ω::LazyΩ)
  result = get(ω.data, primrv, 0)
  if result === 0
    ω.data[primrv] = rand(rng(ω), primrv.class)
  else
    result
  end::eltype(primrv.class)
end

## Updating interface
import ..Util:update
export update, update!
# export .Util:update
update!(ω::LazyΩ, k, v) = (ω[k] = v; ω)
Util.update(ω::LazyΩ, k, v) = (ω_ = deepcopy(ω); ω_[k] = v; ω_)

## Display
function Base.show(io::IO, m::MIME"text/plain", ω::LazyΩ)
  println(io, typeof(ω))
  println("tags:")
  show(io, m, ω.tags)

  println("\ndata:")
  show(io, m, ω.data)
end

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
# module ProFit

## References ##

# written by: Péter Salvi
# original source: https://github.com/salvipeter/profile-fitter

# The general workflow is based on:
#
# P. Benkő, T. Várady:
#   Best fit translational and rotational surfaces for reverse engineering shapes.
#   In: The mathematics of surfaces IX, pp. 70-81, 2000.

# The thinning algorithm for the guiding polygon is based on:
#
# I-K. Lee:
#   Curve reconstruction from unorganized points.
#   In: Computer Aided Geometric Design, Vol. 17(2), pp. 161-177, 2000.

# Circle fitting is based on:
#
# S.J. Ahn, W. Rauh, H-J. Warnecke:
#   Least-squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola.
#   In: Pattern Recognition, Vol. 34(12), pp. 2283-2303, 2001.

## TODO ##

# Beautification:
# - delete short arcs/segments
# - unite similar arcs
# - create segments from arcs with large radii
# - ensure C0 continuity
# - ensure G1 continuity where plausible

#using LinearAlgebra
#using Plots
#using Random

# Delaunay triangulation

"""
    bounding_triangle(points)

Returns a triangle that contains all points.
"""
function bounding_triangle(points)
    minx, maxx = extrema((x->x[1]).(points))
    miny, maxy = extrema((x->x[2]).(points))
    center = [(minx + maxx) / 2, (miny + maxy) / 2]
    r = norm(center - [minx, miny])
    [center - [sqrt(3) * r, r],
     center + [0, 2r],
     center + [sqrt(3) * r, -r]]
end

"""
    circumcircle(a, b, c)

Given three points, gives back the circumcircle as the tuple (center, radius).
"""
function circumcircle(a, b, c)
    d = 2 * (a[1] * (b[2] - c[2]) + b[1] * (c[2] - a[2]) + c[1] * (a[2] - b[2]))
    if abs(d) < 1e-10
        # Degenerate case
        if norm(a - c) > norm(a - b)
            b, c = c, b
        end
        return ((a + b) / 2, norm(a - b) / 2)
    end
    x = (a[1]^2 + a[2]^2) * (b[2] - c[2]) +
        (b[1]^2 + b[2]^2) * (c[2] - a[2]) +
        (c[1]^2 + c[2]^2) * (a[2] - b[2])
    y = (a[1]^2 + a[2]^2) * (c[1] - b[1]) +
        (b[1]^2 + b[2]^2) * (a[1] - c[1]) +
        (c[1]^2 + c[2]^2) * (b[1] - a[1])
    center = [x / d, y / d]
    center, norm(center - a)
end

"""
    to_edges!(triangles)

Converts a list of triangles (index-triples) to a list of edges (index-pairs).
The indices in the individual edges are sorted.
"""
function to_edges!(triangles)
    edges(tri) = sort.([[tri[1], tri[2]], [tri[1], tri[3]], [tri[2], tri[3]]])
    mapreduce(edges, append!, triangles)
end

"""
    walk_polygon(edges)

Given a list of edges (pairs of indices) defining a closed loop,
returns a sorted list of points going around the edges.
"""
function walk_polygon(edges)
    poly = [edges[1][1], edges[1][2]]
    while length(poly) != length(edges)
        for edge in edges
            if edge[1] == poly[end] && edge[2] != poly[end-1]
                push!(poly, edge[2])
                break
            elseif edge[2] == poly[end] && edge[1] != poly[end-1]
                push!(poly, edge[1])
                break
            end
        end
    end
    poly
end

"""
    delaunay(points)

Returns a list of triangles, each consisting of 3 indices to the input list.
"""
function delaunay(points)
    n = length(points)
    p = copy(points)
    append!(p, bounding_triangle(points))
    triangles = [[n+1, n+2, n+3]]

    for i in 1:n
        bad = filter(triangles) do tri
            c, r = circumcircle(p[tri[1]], p[tri[2]], p[tri[3]])
            norm(p[i] - c) <= r
        end
        setdiff!(triangles, bad)
        #isempty(bad) && continue
        bad_edges = to_edges!(bad)
        bad_boundary = filter(bad_edges) do edge
            count(x -> x == edge, bad_edges) == 1
        end
        poly = walk_polygon(bad_boundary)
        m = length(poly)
        for j in 1:m
            push!(triangles, [i, poly[j], poly[mod1(j+1,m)]])
        end
    end

    filter(triangles) do tri
        all(x -> x <= n, tri)
    end
end

# Thinning

"""
    spanning_tree(edges, weights)

Given a list of edges (index pairs) and associated weights,
returns a list of edges defining a minimal weight spanning tree.
"""
function spanning_tree(edges, weights)
    order = sortperm(weights)
    result = Array{Array{Int64,1},1}(undef, 0)
    sets = Dict()
    for edge in edges
        for p in edge
            sets[p] = p
        end
    end
    for i in order
        a, b = edges[i]
        if sets[a] != sets[b]
            push!(result, edges[i])
            new, old = minmax(sets[a], sets[b])
            for (k, v) in sets
                if v == old
                    sets[k] = new
                end
            end
        end
    end
    result
end

"""
    thinning(points, edges, thickness)

Given an Euclidean minimum spanning tree over a set of points taken from a (noisy) curve,
this function generates an ordered list of points representing a thinned version of the curve.
The `thickness` parameter should be around the same size as the noise.
"""
function thinning(points, edges, thickness)
    lookup = [Int[] for _ in points]
    for e in edges
        push!(lookup[e[1]], e[2])
        push!(lookup[e[2]], e[1])
    end

    function bfs(start)
        seen = fill(false, length(points))
        result = []
        queue = [start]
        while !isempty(queue)
            v = popfirst!(queue)
            push!(result, v)
            for w in lookup[v]
                if !seen[w]
                    push!(queue, w)
                    seen[w] = true
                end
            end
        end
        result
    end
    function neighbors(start)
        x = points[start]
        result = []
        function rec(next)
            push!(result, next)
            for p in lookup[next]
                if !(p in result) && norm(x - points[p]) < thickness
                    rec(p)
                end
            end
        end
        rec(start)
        result
    end

    ordered = bfs(bfs(1)[end])
    thinned = Array{Array{Float64,1},1}(undef, 0)
    chunks = Array{Array{Int64,1},1}(undef, 0)
    while !isempty(ordered)
        adjacent = neighbors(ordered[1])
        push!(thinned, mapreduce(x -> points[x], +, adjacent) / length(adjacent))
        push!(chunks, adjacent)
        setdiff!(ordered, adjacent)
    end
    thinned, chunks
end

"""
    thinning(points, thickness)

Given a set of unsorted points taken from a (noisy) curve,
this function generates an ordered list of points representing a thinned version of the curve.
The `thickness` parameter should be around the same size as the noise.
"""
function thinning(points, thickness)
    @logmsg IterLow1  "Delaunay"
    tris = delaunay(points)
    all_edges = to_edges!(tris)
    @logmsg IterLow1 "Computing distance for $(length(all_edges)) edges, with $(length(points)) points"
    weights = [norm(points[e[1]] - points[e[2]]) for e in all_edges]
    @logmsg IterLow1 "Spanning tree"
    tree = spanning_tree(all_edges, weights)
    @logmsg IterLow1 "Thinning"
    thinning(points, tree, thickness)
end

function alle(l)
    ae = Vector{Vector{Int}}(undef, 0)
    for i in 1:l
        for j in i+1:l
            #i == j && continue
            push!(ae, [i,j])
        end
    end
    ae
end

"""
    thinning_slow(points, thickness)

Given a set of unsorted points taken from a (noisy) curve,
this function generates an ordered list of points representing a thinned version of the curve.
The `thickness` parameter should be around the same size as the noise.
Delaunay triangulation is not used to accelerate the code.
"""
function thinning_slow(points, thickness)
    @logmsg IterLow1 "Collect edges for $(length(points)) points."
    all_edges = alle(size(points,1))
    @logmsg IterLow1 "Computing distance"
    weights = [norm(points[e[1]] - points[e[2]]) for e in all_edges]
    @logmsg IterLow1 "Spanning tree"
    tree = spanning_tree(all_edges, weights)
    @logmsg IterLow1 "Thinning"
    thinning(points, tree, thickness)
end

function thinning_deldir(points, thickness)
    miv, mav = findAABB(points)
    rect = [miv[1]; mav[1]; miv[2]; mav[2]]
    pix = (x->x[1]).(points)
    piy = (x->x[2]).(points)
    @logmsg IterLow1  "Delaunay"
    deli1, deli2 = deldir(pix, piy, rect)
    all_edges = [ [deli1[i], deli2[i]] for i in 1:size(deli1,1)]
    #tris = delaunay(points)
    #all_edges = to_edges!(tris)
    @logmsg IterLow1 "Computing distance for $(length(all_edges)) edges, with $(length(points)) points"
    weights = [norm(points[e[1]] - points[e[2]]) for e in all_edges]
    @logmsg IterLow1 "Spanning tree"
    tree = spanning_tree(all_edges, weights)
    @logmsg IterLow1 "Thinning"
    thinning(points, tree, thickness)
end

# end # module


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ e7a070ab-67e7-444b-88d8-87c14aaef046
names(@__MODULE__)

# ╔═╡ b1491ca6-f791-41d8-96b5-28971084f34c

# ╔═╡ 05804305-cb1f-4c97-8937-f56289222bd7
md"""
Rerun the big cell after you edit one of the JS/CSS assets.
"""

# ╔═╡ 9be4d586-76d8-11eb-06ed-c53d3aef0469
begin
	export Clock

	Base.@kwdef struct Clock
		interval::Real = 1
		fixed::Bool = false
		start_running::Bool = false
		max_value::Union{Int64,Nothing} = nothing

		# Clock(interval, fixed, start_running, max_value) = interval >= 0 ? new(interval, fixed, start_running, max_value) : error("interval must be non-negative")
	end

	# for backwards compat
	Clock(interval; kwargs...) = Clock(interval=interval; kwargs...)

	Clock(interval, fixed, start_running=false) = Clock(interval, fixed, start_running, nothing)

	# We split the HTML string into multiple files, but you could also write all of this into a single (long) string 🎈
	const cb = read(joinpath(@__DIR__, "..", "assets", "clock_back.svg"), String)
	const cf = read(joinpath(@__DIR__, "..", "assets", "clock_front.svg"), String)
	const cz = read(joinpath(@__DIR__, "..", "assets", "clock_zoof.svg"), String)
	const js = read(joinpath(@__DIR__, "..", "assets", "clock.js"), String)
	const css = read(joinpath(@__DIR__, "..", "assets", "clock.css"), String)

	function Base.show(io::IO, ::MIME"text/html", clock::Clock)
		clock.interval < 0 && error("interval must be non-negative")

		result = """
		<plutoui-clock class='$(clock.fixed ? " fixed" : "")$(clock.start_running ? "" : " stopped")' data-max-value=$(repr(clock.max_value))>
			<plutoui-analog>
				<plutoui-back>$(cb)</plutoui-back>
				<plutoui-front>$(cf)</plutoui-front>
				<plutoui-zoof style="opacity: 0">$(cz)</plutoui-zoof>
			</plutoui-analog>
			<button></button>
			speed:
			<input type="number" value="$(clock.interval)"  min=0 step=any lang="en-001">

		</plutoui-clock>

	    """
	    write(io, result)
	end

	Base.get(clock::Clock) = 1
end

# ╔═╡ 06289ad2-9e2f-45b3-9d15-7c5a4167e138
@bind tick Clock()

# ╔═╡ 9ecd95f0-d7a5-4ee9-9e18-9d87e5d43ab7
tick; rand()

# ╔═╡ d82dae11-b2c6-42b5-8c52-67fbb6cc236a
tick

# ╔═╡ 80c6e80e-077a-4e31-9467-788a8c437bfc
@bind fasttick Clock(0.001, max_value=1000)

# ╔═╡ 63854404-e6a5-4dc6-a40e-b09b9f531465
fasttick

# ╔═╡ a5f8ed96-136c-4ff4-8275-bd569f0dae40
md"""
## Different constructors
"""

# ╔═╡ c3c07db2-bb9c-4521-83ab-e81fbb376b4e
Clock(3.0)

# ╔═╡ 83a021ab-7cca-47c7-a560-9cbf58b35ab7
Clock(3.0, true)

# ╔═╡ c96dfd13-ddd4-443f-ab09-30e15ea76785
Clock(3.0, true, true)

# ╔═╡ 78ee5465-ce3b-45f6-acec-aa69175807f5
Clock(3.0, true, true, 5)

# ╔═╡ 9115fbcd-1550-4439-a830-c69b83b774b3
Clock()

# ╔═╡ f4104cb3-7c07-4814-99f9-a00764ebadf6
@assert try

	# this should error:
	repr(MIME"text/html"(), Clock(-5))

	false
catch
	true
end

# ╔═╡ 21cba3fb-7bb0-43ae-b4c4-5c1eb7241fec
Clock(2.0, max_value=123)

# ╔═╡ Cell order:
# ╠═06289ad2-9e2f-45b3-9d15-7c5a4167e138
# ╠═9ecd95f0-d7a5-4ee9-9e18-9d87e5d43ab7
# ╠═d82dae11-b2c6-42b5-8c52-67fbb6cc236a
# ╠═80c6e80e-077a-4e31-9467-788a8c437bfc
# ╠═63854404-e6a5-4dc6-a40e-b09b9f531465
# ╠═e7a070ab-67e7-444b-88d8-87c14aaef046
# ╠═b1491ca6-f791-41d8-96b5-28971084f34c
# ╟─05804305-cb1f-4c97-8937-f56289222bd7
# ╠═9be4d586-76d8-11eb-06ed-c53d3aef0469
# ╟─a5f8ed96-136c-4ff4-8275-bd569f0dae40
# ╠═c3c07db2-bb9c-4521-83ab-e81fbb376b4e
# ╠═83a021ab-7cca-47c7-a560-9cbf58b35ab7
# ╠═c96dfd13-ddd4-443f-ab09-30e15ea76785
# ╠═78ee5465-ce3b-45f6-acec-aa69175807f5
# ╠═9115fbcd-1550-4439-a830-c69b83b774b3
# ╠═f4104cb3-7c07-4814-99f9-a00764ebadf6
# ╠═21cba3fb-7bb0-43ae-b4c4-5c1eb7241fec


text:
function test_softlabel_softmax_loss_layer(backend::Backend, tensor_dim, T, eps)
  println("-- Testing SoftlabelSoftmaxLossLayer on $(typeof(backend)){$T}...")

  dims = abs(rand(Int,tensor_dim)) % 6 + 6
  dims = tuple(dims...)
  op_dim = max(abs(rand(Int)) % tensor_dim, 1)
  println("    > $dims (operate on dimension $op_dim)")

  input = rand(T, dims) + 0.01
  input_blob = make_blob(backend, input)
  diff_blob = make_blob(backend, T, size(input))

  labels = abs(rand(T, dims)) + 0.01
  labels = labels ./ sum(labels, op_dim)
  label_blob = make_blob(backend, labels)

  inputs = Blob[input_blob, label_blob]

  loss_weight = abs(rand(T))
  layer = SoftlabelSoftmaxLossLayer(bottoms=[:pred, :labels], dim=op_dim, weight=loss_weight)
  state = setup(backend, layer, inputs, Blob[diff_blob])

  println("    > Forward")
  forward(backend, state, inputs)

  expected_loss = 0
  expected_grad = zeros(T, size(input))

  dim_pre, dim_prob, dim_post = split_dims(input, op_dim)
  canonical_input = reshape(input, dim_pre, dim_prob, dim_post)
  canonical_label = reshape(labels, dim_pre, dim_prob, dim_post)
  canonical_grad = reshape(expected_grad, dim_pre, dim_prob, dim_post)

  for i = 1:dim_pre
    for j = 1:dim_post
      pred = exp(canonical_input[i,:,j])
      pred /= sum(pred)

      expected_loss += sum(-log(pred) .* canonical_label[i,:,j])
      canonical_grad[i,:,j] = repmat(vec(pred), 1, dim_prob) * vec(canonical_label[i,:,j])
      canonical_grad[i,:,j] -= canonical_label[i,:,j]
    end
  end
  scale = dims[op_dim] / prod(dims)
  expected_loss *= scale * loss_weight
  expected_grad *= scale * loss_weight
  expected_grad = reshape(expected_grad, size(input))

  @test -eps < state.loss - expected_loss < eps

  println("    > Backward")
  backward(backend, state, inputs, Blob[diff_blob])
  grad = to_array(diff_blob)

  @test all(-eps .< grad - expected_grad .< eps)

  shutdown(backend, state)
end

function test_softlabel_softmax_loss_layer(backend::Backend, T, eps)
  for i in [2,4,5]
    test_softlabel_softmax_loss_layer(backend, i, T, eps)
  end
end
function test_softlabel_softmax_loss_layer(backend::Backend)
  test_softlabel_softmax_loss_layer(backend, Float64, 1e-5)
  test_softlabel_softmax_loss_layer(backend, Float32, 1e-3)
end
if test_cpu
  test_softlabel_softmax_loss_layer(backend_cpu)
end
if test_gpu
  test_softlabel_softmax_loss_layer(backend_gpu)
end


text:
"""
    $(SIGNATURES)

Returns `n` decision thresholds which correspond to `n` evenly spaced quantiles
of the given vector of scores. If `reduced == true`, then the resulting `n` is
`min(length(scores) + 1, n)`. If `zerorecall == true`, then the largest threshold
will be `maximum(scores)*(1 + eps())` otherwise `maximum(scores)`.
"""
function thresholds(scores::RealVector, n::Int = length(scores) + 1; reduced::Bool = true, zerorecall::Bool = true)
    ns = length(scores)
    N = reduced ? min(ns + 1, n) : n
    N -= zerorecall
    if N == ns
        thres = sort(scores)
    else
        thres = quantile(scores, range(0, 1, length = N))
    end
    zerorecall && push!(thres, maximum(scores) + eps())

    return thres
end

function threshold_at_rate(enc::TwoClassEncoding, scores::RealVector, rates::RealVector; rev::Bool = true)

    all(0 .<= rates .<= 1) || throw(ArgumentError("input rates out of [0, 1]."))
    issorted(rates) || throw(ArgumentError("input rates must be sorted."))

    n_rates = length(rates)
    n_scores = length(scores)
    print_warn = falses(n_rates)

    # case rate == 1
    if rev
        thresh = fill(scores[end] + eps(scores[end]), n_rates)
        thresh[rates .== 1] .= scores[end]
    else
        thresh = fill(scores[end], n_rates)
    end

    # case rate != 1
    rate_last = 0
    t_last = scores[1]
    j = 1

    for (i, score) in enumerate(scores)
        t_last == score && continue

        # compute current rate
        rate = (i-1)/n_scores

        for rate_i in rates[j:end]
            rate <= rate_i && break
            rate_last == 0 && rate_i != 0 && (print_warn[j] = true)

            thresh[j] = rev ? t_last + eps(t_last) : t_last
            j += 1
            j > n_rates && (return thresh, print_warn)
        end

        # update last rate and threshold
        rate_last = rate
        t_last = score
    end
    return thresh, print_warn
end

"""
    $(SIGNATURES)

Returns a decision threshold at a given true positive rate `tpr ∈ [0, 1]`.
"""
threshold_at_tpr(targets::AbstractVector, scores::RealVector, tpr) =
    threshold_at_tpr(current_encoding(), targets, scores, tpr)

threshold_at_tpr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, tpr::Real) =
    threshold_at_tpr(enc, targets, scores, [tpr])[1]

function threshold_at_tpr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, tpr::RealVector)

    length(targets) == length(scores) || throw(DimensionMismatch("Inconsistent lengths of `targets` and `scores`."))
    check_encoding(enc, targets) || throw(ArgumentError("`targets` vector uses incorrect label encoding."))

    scores_pos = sort(scores[ispositive.(enc, targets)]; rev = false)
    rates = round.(1 .- reverse(tpr); digits = 14)
    ts, print_warn = threshold_at_rate(enc, scores_pos, rates; rev = false)

    if any(print_warn)
        rates = tpr[reverse(print_warn)]
        @warn "The closest higher feasible true positive rate to some of the required values ($(join(rates, ", "))) is 1.0!"
    end
    return reverse(ts)
end

"""
    $(SIGNATURES)

Returns a decision threshold at a given true negative rate `fpr ∈ [0, 1]`.
"""
threshold_at_tnr(targets::AbstractVector, scores::RealVector, tnr) =
    threshold_at_tnr(current_encoding(), targets, scores, tnr)

threshold_at_tnr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, tnr::Real) =
    threshold_at_tnr(enc, targets, scores, [tnr])[1]

function threshold_at_tnr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, tnr::RealVector)

    length(targets) == length(scores) || throw(DimensionMismatch("Inconsistent lengths of `targets` and `scores`."))
    check_encoding(enc, targets) || throw(ArgumentError("`targets` vector uses incorrect label encoding."))

    scores_neg = sort(scores[isnegative.(enc, targets)]; rev = true)
    rates = round.(1 .- reverse(tnr); digits = 14)
    ts, print_warn = threshold_at_rate(enc, scores_neg, rates; rev = true)

    if any(print_warn)
        rates = tnr[reverse(print_warn)]
        @warn "The closest higher feasible true negative rate to some of the required values ($(join(rates, ", "))) is 1.0!"
    end
    return reverse(ts)
end

"""
    $(SIGNATURES)

Returns a decision threshold at a given false positive rate `fpr ∈ [0, 1]`.
"""
threshold_at_fpr(targets::AbstractVector, scores::RealVector, fpr) =
    threshold_at_fpr(current_encoding(), targets, scores, fpr)

threshold_at_fpr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, fpr::Real) =
    threshold_at_fpr(enc, targets, scores, [fpr])[1]

function threshold_at_fpr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, fpr::RealVector)

    length(targets) == length(scores) || throw(DimensionMismatch("Inconsistent lengths of `targets` and `scores`."))
    check_encoding(enc, targets) || throw(ArgumentError("`targets` vector uses incorrect label encoding."))

    scores_neg = sort(scores[isnegative.(enc, targets)]; rev = true)
    ts, print_warn = threshold_at_rate(enc, scores_neg, fpr; rev = true)

    if any(print_warn)
        rates = fpr[print_warn]
        @warn "The closest lower feasible false positive rate to some of the required values ($(join(rates, ", "))) is 0.0!"
    end
    return ts
end

"""
    $(SIGNATURES)

Returns a decision threshold at a given false negative rate `fnr ∈ [0, 1]`.
"""
threshold_at_fnr(targets::AbstractVector, scores::RealVector, fnr) =
    threshold_at_fnr(current_encoding(), targets, scores, fnr)

threshold_at_fnr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, fnr::Real) =
    threshold_at_fnr(enc, targets, scores, [fnr])[1]

function threshold_at_fnr(enc::TwoClassEncoding, targets::AbstractVector, scores::RealVector, fnr::RealVector)

    length(targets) == length(scores) || throw(DimensionMismatch("Inconsistent lengths of `targets` and `scores`."))
    check_encoding(enc, targets) || throw(ArgumentError("`targets` vector uses incorrect label encoding."))

    scores_pos = sort(scores[ispositive.(enc, targets)]; rev = false)
    ts, print_warn = threshold_at_rate(enc, scores_pos, fnr; rev = false)

    if any(print_warn)
        rates = fnr[print_warn]
        @warn "The closest lower feasible false negative rate to some of the required values ($(join(rates, ", "))) is 0.0!"
    end
    return ts
end

"""
    $(SIGNATURES)

Returns a decision threshold at `k` most anomalous samples if `rev == true` and
a decision threshold at `k` least anomalous samples otherwise.
"""
function threshold_at_k(scores::RealVector, k::Int; rev::Bool = true)
    n = length(scores)
    n >= k || throw(ArgumentError("Argument `k` must be smaller or equal to `length(targets) = $n`"))
    return partialsort(scores, k, rev = rev)
end


text:
export static_controller

# A fixed architecture controller with
# + ReLU neurons
# + no neurogenesis or synaptogenesis
# + Reward-modulated Hebbian learning (Oja's rule)

function static_gen_cell_state_update(cfg::Dict)
    function static_cell_state_update(inputs::Array{Float64})
        # inputs: cell_params, cell_state, input
        cell_state = inputs[cfg["n_cell_params"] .+ (1:cfg["n_cell_state"])]
        einput = inputs[end]
        output = max(einput, 0.0)
        cell_state[1] = 0.9 * cell_state[1] + 0.1 * einput
        vcat(cell_state, output)
    end
end

function static_gen_cell_param_update(cfg::Dict)
    function static_cell_param_update(inputs::Array{Float64})
        # cell params, state -> params
        # save recent cell action in cell params
        cell_params = inputs[1:cfg["n_cell_params"]]
        cell_state = @view inputs[(cfg["n_cell_params"]+1):end]
        cell_params[1] = cell_state[1]
        cell_params
    end
end

function static_gen_new_conn_params(cfg::Dict)
    function new_conn_params(inputs::Array{Float64})
        c1_params = @view inputs[1:cfg["n_conn_params"]]
        c2_params = @view inputs[cfg["n_conn_params"] .+ (1:cfg["n_conn_params"])]
        if c1_params[cfg["n_conn_params"]] == 1.0
            ones(cfg["n_conn_params"])
        else
            rand(cfg["n_conn_params"]) .- 0.5
        end
    end
end

function static_gen_conn_state_update(cfg::Dict)
    function static_conn_state_update(inputs::Array{Float64})
        # inputs: conn params, state, cell in, conn in
        params = @view inputs[1:cfg["n_conn_params"]]
        state = inputs[cfg["n_conn_params"] .+ (1:cfg["n_conn_state"])]
        cell_in = inputs[length(inputs)-1]
        conn_in = inputs[end]
        state[1] = 0.9 * state[1] + 0.1 * conn_in
        output = params[1] * cell_in
        vcat(state, output)
    end
    static_conn_state_update
end

function static_gen_conn_param_update(cfg::Dict)
    function static_conn_param_update(inputs::Array{Float64})
        # inputs: c1 params, c2 params, conn params, state
        c1_params = @view inputs[1:cfg["n_conn_params"]]
        c2_params = @view inputs[cfg["n_conn_params"] .+ (1:cfg["n_conn_params"])]
        params = inputs[2*cfg["n_cell_params"] .+ (1:cfg["n_conn_params"])]
        state = @view inputs[((2*cfg["n_cell_params"] + cfg["n_conn_params"])
                              .+ (1:cfg["n_conn_state"]))]
        dw = (0.01 + state[1]) * c2_params[1] * (c1_params[1] - c2_params[1] * params[1])
        params[1] = max(-1.0, min(1.0, params[1] + dw))
        params
    end
    static_conn_param_update
end

function static_controller(cfg::Dict)
    nouts = output_lengths(cfg)
    cell_division(x::Array{Float64}) = false
    new_cell_params(x::Array{Float64}) = zeros(nouts[2])
    cell_death(x::Array{Float64}) = false
    cell_state_update = static_gen_cell_state_update(cfg)
    cell_param_update = static_gen_cell_param_update(cfg)
    connect(x::Array{Float64}) = false
    new_conn_params = static_gen_new_conn_params(cfg)
    disconnect(x::Array{Float64}) = false
    conn_state_update = static_gen_conn_state_update(cfg)
    conn_param_update = static_gen_conn_param_update(cfg)

    Controller(cell_division, new_cell_params, cell_death,
               cell_state_update, cell_param_update,
               connect, new_conn_params, disconnect,
               conn_state_update, conn_param_update)
end


text:
"""
PointEvalHandler(grid::Grid, points::AbstractVector{Vec{dim,T}}) where {dim, T}

A `PointEvalHandler` computes the corresponding cell for each point in `points` and its local coordinate within the cell.
"""
PointEvalHandler

struct PointEvalHandler{G,dim,T<:Real}
    grid::G
    cells::Vector{Union{Nothing, Int}}
    local_coords::Vector{Union{Nothing, Vec{dim,T}}}
end

function Base.show(io::IO, ::MIME"text/plain", ph::PointEvalHandler)
    println(io, typeof(ph))
    println(io, "  number of points: ", length(ph.local_coords))
    n_missing = sum(x -> x === nothing, ph.cells)
    if n_missing == 0
        print(io, "  Found corresponding cell for all points.")
    else
        print(io, "  Could not find corresponding cell for ", n_missing, " points.")
    end
end

function PointEvalHandler(grid::AbstractGrid, points::AbstractVector{Vec{dim,T}}) where {dim, T}
    node_cell_dicts = _get_node_cell_map(grid)
    cells, local_coords = _get_cellcoords(points, grid, node_cell_dicts)
    return PointEvalHandler(grid, cells, local_coords)
end

function _get_cellcoords(points::AbstractVector{Vec{dim,T}}, grid::AbstractGrid, node_cell_dicts::Dict{C,Dict{Int, Vector{Int}}}) where {dim, T<:Real, C}

    # set up tree structure for finding nearest nodes to points
    kdtree = KDTree(reinterpret(Vec{dim,T}, getnodes(grid)))
    nearest_nodes, _ = knn(kdtree, points, 3, true) #TODO 3 is a random value, it shouldn't matter because likely the nearest node is the one we want

    cells = Vector{Union{Nothing, Int}}(nothing, length(points))
    local_coords = Vector{Union{Nothing, Vec{dim, T}}}(nothing, length(points))

    for point_idx in 1:length(points)
        cell_found = false
        for (idx, (CT, node_cell_dict)) in enumerate(node_cell_dicts)
            geom_interpol = default_interpolation(CT)
            # loop over points
            for node in nearest_nodes[point_idx]
                possible_cells = get(node_cell_dict, node, nothing)
                possible_cells === nothing && continue # if node is not part of the fieldhandler, try the next node
                for cell in possible_cells
                    cell_coords = getcoordinates(grid, cell)
                    is_in_cell, local_coord = point_in_cell(geom_interpol, cell_coords, points[point_idx])
                    if is_in_cell
                        cell_found = true
                        cells[point_idx] = cell
                        local_coords[point_idx] = local_coord
                        break
                    end
                end
                cell_found && break
            end
            cell_found && break
        end
        if !cell_found
            @warn("No cell found for point number $point_idx, coordinate: $(points[point_idx]).")
        end
    end
    return cells, local_coords
end

# check if point is inside a cell based on physical coordinate
function point_in_cell(geom_interpol::Interpolation{dim,shape,order}, cell_coordinates, global_coordinate) where {dim, shape, order}
    converged, x_local = find_local_coordinate(geom_interpol, cell_coordinates, global_coordinate)
    if converged
        return _check_isoparametric_boundaries(shape, x_local), x_local
    else
        return false, x_local
    end
end

# check if point is inside a cell based on isoparametric coordinate
function _check_isoparametric_boundaries(::Type{RefCube}, x_local::Vec{dim, T}) where {dim, T}
    tol = sqrt(eps(T))
    # All in the range [-1, 1]
    return all(x -> abs(x) - 1 < tol, x_local)
end

# check if point is inside a cell based on isoparametric coordinate
function _check_isoparametric_boundaries(::Type{RefTetrahedron}, x_local::Vec{dim, T}) where {dim, T}
    tol = sqrt(eps(T))
    # Positive and below the plane 1 - ξx - ξy - ξz
    return all(x -> x > -tol, x_local) && sum(x_local) - 1 < tol
end

# TODO: should we make iteration params optional keyword arguments?
function find_local_coordinate(interpolation, cell_coordinates, global_coordinate)
    """
    currently copied verbatim from https://discourse.julialang.org/t/finding-the-value-of-a-field-at-a-spatial-location-in-juafem/38975/2
    other than to make J dim x dim rather than 2x2
    """
    dim = length(global_coordinate)
    local_guess = zero(Vec{dim})
    n_basefuncs = getnbasefunctions(interpolation)
    max_iters = 10
    tol_norm = 1e-10
    converged = false
    for iter in 1:max_iters
        N = Ferrite.value(interpolation, local_guess)

        global_guess = zero(Vec{dim})
        for j in 1:n_basefuncs
            global_guess += N[j] * cell_coordinates[j]
        end
        residual = global_guess - global_coordinate
        if norm(residual) <= tol_norm
            converged = true
            break
        end
        dNdξ = Ferrite.derivative(interpolation, local_guess)
        J = zero(Tensor{2, dim})
        for j in 1:n_basefuncs
            J += cell_coordinates[j] ⊗ dNdξ[j]
        end
        local_guess -= inv(J) ⋅ residual
    end
    return converged, local_guess
end

# return a Dict with a key for each node that contains a vector with the adjacent cells as value
function _get_node_cell_map(grid::AbstractGrid)
    cells = getcells(grid)
    C = eltype(cells) # possibly abstract
    cell_dicts = Dict{Type{<:C}, Dict{Int, Vector{Int}}}()
    ctypes = Set{Type{<:C}}(typeof(c) for c in cells)
    for ctype in ctypes
        cell_dict = cell_dicts[ctype] = Dict{Int,Vector{Int}}()
        for (cellidx, cell) in enumerate(cells)
            cell isa ctype || continue
            for node in cell.nodes
                v = get!(Vector{Int}, cell_dict, node)
                push!(v, cellidx)
            end
        end
    end
    return cell_dicts
end

"""
    get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_values::Vector{T}) where T
    get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_values::Vector{T}, [fieldname::Symbol]) where T

Return a `Vector{T}` (for a 1-dimensional field) or a `Vector{Vec{fielddim, T}}` (for a vector field) with the field values of field `fieldname` in the points of the `PointEvalHandler`. The `fieldname` can be omitted if only one field is stored in `dh`. The field values are computed based on the `dof_values` and interpolated to the local coordinates by the function interpolation of the corresponding `field` stored in the `AbstractDofHandler` or the `L2Projector`.

"""
get_point_values

function get_point_values(ph::PointEvalHandler, proj::L2Projector, dof_vals::AbstractVector)
    get_point_values(ph, proj.dh, dof_vals)
end

function get_point_values(ph::PointEvalHandler, dh::AbstractDofHandler, dof_vals::AbstractVector{T},
                           fname::Symbol=find_single_field(dh)) where {T}
    fdim = getfielddim(dh, fname)
    npoints = length(ph.cells)
    # for a scalar field return a Vector of Scalars, for a vector field return a Vector of Vecs
    if fdim == 1
        out_vals = fill!(Vector{T}(undef, npoints), NaN * zero(T))
    else
        nanv = convert(Vec{fdim,T}, NaN * zero(Vec{fdim,T}))
        out_vals = fill!(Vector{Vec{fdim, T}}(undef, npoints), nanv)
    end
    func_interpolations = get_func_interpolations(dh, fname)
    get_point_values!(out_vals, ph, dh, dof_vals, fname, func_interpolations)
    return out_vals
end
function find_single_field(dh)
    ns = getfieldnames(dh)
    if length(ns) != 1
        throw(ArgumentError("multiple fields in DoF handler, must specify which"))
    end
    return ns[1]
end

# values in dof-order. They must be obtained from the same DofHandler that was used for constructing the PointEvalHandler
function get_point_values!(out_vals::Vector{T2},
    ph::PointEvalHandler,
    dh::MixedDofHandler,
    dof_vals::Vector{T},
    fname::Symbol,
    func_interpolations
    ) where {T2, T}

    # TODO: I don't think this is correct??
    length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.")

    fdim = getfielddim(dh, fname)

    for fh_idx in eachindex(dh.fieldhandlers)
        ip = func_interpolations[fh_idx]
        if ip !== nothing
            dofrange = dof_range(dh.fieldhandlers[fh_idx], fname)
            cellset = dh.fieldhandlers[fh_idx].cellset
            _get_point_values!(out_vals, dof_vals, ph, dh, ip, cellset, Val(fdim), dofrange)
        end
    end
    return out_vals
end

function get_point_values!(out_vals::Vector{T2},
    ph::PointEvalHandler,
    dh::DofHandler,
    dof_vals::Vector{T},
    fname::Symbol,
    func_interpolations
    ) where {T2, T}

    # TODO: I don't think this is correct??
    length(dof_vals) == ndofs(dh) || error("You must supply values for all $(ndofs(dh)) dofs.")

    fdim = getfielddim(dh, fname)
    dofrange = dof_range(dh, fname)
    _get_point_values!(out_vals, dof_vals, ph, dh, func_interpolations[1], nothing, Val(fdim), dofrange)
    return out_vals
end

# function barrier with concrete type of interpolation
function _get_point_values!(
    out_vals::Vector{T2},
    dof_vals::Vector{T},
    ph::PointEvalHandler,
    dh::AbstractDofHandler,
    ip::Interpolation,
    cellset::Union{Nothing, Set{Int}},
    fdim::Val{fielddim},
    dofrange::AbstractRange{Int},
    ) where {T2,T,fielddim}

    # extract variables
    local_coords = ph.local_coords
    # preallocate some stuff specific to this cellset
    pv = PointScalarValues(first(local_coords), ip)
    first_cell = cellset === nothing ? 1 : first(cellset)
    cell_dofs = Vector{Int}(undef, ndofs_per_cell(dh, first_cell))

    # compute point values
    for pointid in eachindex(ph.cells)
        cellid = ph.cells[pointid]
        cellid ===nothing && continue # next point if no cell was found for this one
        cellset !== nothing && (cellid ∈ cellset || continue) # no need to check the cellset for a regular DofHandler
        celldofs!(cell_dofs, dh, ph.cells[pointid])
        reinit!(pv, local_coords[pointid], ip)
        @inbounds @views dof_vals_reshaped = _change_format(fdim, dof_vals[cell_dofs[dofrange]])
        out_vals[pointid] = function_value(pv, 1, dof_vals_reshaped)
    end
    return out_vals
end

###################################################################################################################
# utils

# reshape dof_values based on fielddim
_change_format(::Val{1}, dof_values::AbstractVector{T}) where T = dof_values
_change_format(::Val{fielddim}, dof_values::AbstractVector{T}) where {fielddim, T} = reinterpret(Vec{fielddim, T}, dof_values)

get_func_interpolations(dh::DH, fieldname) where DH<:DofHandler = [getfieldinterpolation(dh, find_field(dh, fieldname))]
function get_func_interpolations(dh::DH, fieldname) where DH<:MixedDofHandler
    func_interpolations = Union{Interpolation,Nothing}[]
    for fh in dh.fieldhandlers
        j = findfirst(i -> i === fieldname, getfieldnames(fh))
        if j === nothing
            push!(func_interpolations, missing)
        else
            push!(func_interpolations, fh.fields[j].interpolation)
        end
    end
    return func_interpolations
end


text:
module PuppyBook

    include("03_running_proportion.jl")
    include("04_bayes_update.jl")
    include("therapeutic_touch.jl")

    export compute_running_proportion, plot_running_proportion,
        coinflip_discrete_bayes_update, plot_coinflip_discrete_bayes_update,
        load_therapeutic_touch_data, tt_sample_basic_model, tt_sample_hierarchical_model,
        plot_tt_shrinkage_effect

    greet() = print("Hello Puppies!")

end # module


text:
@testset " particle-mesh coupling with spline 1d " begin

    import GEMPIC: set_x!, set_v!, set_weights!
    import GEMPIC: get_x, get_v, get_charge, add_charge!, add_charge_pp!
    import GEMPIC: add_current_update_v!, add_current_update_v_pp!
    import GEMPIC: b_to_pp, evaluate, evaluate_pp

    n_cells = 10            # Number of cells
    n_particles = 4             # Number of particles
    spline_degree = 3             # Spline degree
    mesh = OneDGrid(0.0, 2.0, n_cells)    # x_min, x_max
    x_vec = [0.1, 0.65, 0.7, 1.5] # Particle positions
    v_vec = [
        1.5 -3.0 0.0 6.0
        0.0 0.5 0.0 0.0
    ]'

    particle_group = ParticleGroup{1,2}(n_particles)

    for i_part = 1:n_particles
        set_x!(particle_group, i_part, x_vec[i_part])
        set_v!(particle_group, i_part, v_vec[i_part, :])
        set_weights!(particle_group, i_part, 1.0)
    end

    values_grid = zeros(Float64, (4, 1, 4))

    values_grid[:, 1, 1] .= [
        2.0833333333333332e-002,
        0.47916666666666663,
        0.47916666666666663,
        2.0833333333333332e-002,
    ]

    values_grid[:, 1, 3] .= values_grid[:, 1, 1]

    values_grid[:, 1, 4] .= values_grid[:, 1, 1]

    values_grid[:, 1, 2] .= [
        7.0312500000000000e-002,
        0.61197916666666663,
        0.31510416666666663,
        2.6041666666666665E-003,
    ]

    kernel = ParticleMeshCoupling1D(mesh, n_particles, spline_degree, :collocation)

    # Accumulate rho
    rho_dofs = zeros(Float64, n_cells)
    rho_dofs1 = zeros(Float64, n_cells)
    for i_part = 1:n_particles
        xi = get_x(particle_group, i_part)[1]
        wi = get_charge(particle_group, i_part)
        add_charge!(rho_dofs, kernel, xi, wi)
        add_charge_pp!(rho_dofs1, kernel, xi, wi)
    end

    rho_dofs_ref = zeros(Float64, n_cells)
    rho_dofs_ref[8:10] .= values_grid[1:3, 1, 1]
    rho_dofs_ref[1] = values_grid[4, 1, 1]
    rho_dofs_ref[1:4] .= rho_dofs_ref[1:4] + values_grid[:, 1, 2] + values_grid[:, 1, 3]
    rho_dofs_ref[5:8] .= rho_dofs_ref[5:8] + values_grid[:, 1, 4]
    rho_dofs_ref .= rho_dofs_ref / n_particles * n_cells / mesh.xmax[1]

    @test maximum(abs.(rho_dofs .- rho_dofs_ref)) < 1e-15
    @test maximum(abs.(rho_dofs1 .- rho_dofs_ref)) < 1e-15

    j_dofs = zeros(Float64, n_cells)
    j_dofs1 = zeros(Float64, n_cells)
    b_dofs = zeros(Float64, n_cells)

    for i_part = 1:n_particles
        xi = get_x(particle_group, i_part)[1]
        wi = get_charge(particle_group, i_part)
        vi = get_v(particle_group, i_part)
        vi1 = vi
        x_new = xi + vi[1] / 10.0

        vi[2] = add_current_update_v!(j_dofs, kernel, xi, x_new, wi, 1.0, b_dofs, vi[2])
        vi1 = add_current_update_v_pp!(j_dofs1, kernel, xi, x_new, wi, 1.0, b_dofs, vi1)

    end

    j_dofs_ref = [
        2.4617513020833336e-002,
        4.0690104166666692e-005,
        0.0,
        0.0,
        0.0,
        0.0,
        0.0,
        6.5104166666666674e-004,
        4.6834309895833329e-002,
        0.11535644531249999,
    ]

    j_dofs_ref .=
        j_dofs_ref .+ [
            -0.16219075520833331,
            -0.16219075520833331,
            -2.52685546875e-002,
            -4.0690104166666692e-005,
            0.0,
            0.0,
            0.0,
            0.0,
            -4.0690104166666692e-005,
            -2.52685546875e-002,
        ]

    j_dofs_ref .=
        j_dofs_ref .+ [
            6.5104166666666696e-004,
            0.0,
            0.0,
            0.0,
            6.5104166666666674e-004,
            5.0130208333333329e-002,
            0.19986979166666666,
            0.24869791666666663,
            0.19986979166666666,
            5.0130208333333329e-002,
        ]

    @test maximum(abs.(j_dofs .- j_dofs_ref)) < 1e-15
    @test maximum(abs.(j_dofs1 .- j_dofs_ref)) < 1e-15

    rho_dofs_pp = b_to_pp(kernel.spline_pp, n_cells, rho_dofs)

    particle_values = zeros(Float64, 4)
    particle_values1 = zeros(Float64, 4)

    for i_part = 1:n_particles
        xi = get_x(particle_group, i_part)
        particle_values[i_part] = evaluate(kernel, xi[1], rho_dofs)
        particle_values1[i_part] = evaluate_pp(kernel, xi[1], rho_dofs_pp)
    end

    particle_values_ref =
        [1.1560058593749998, 2.3149278428819446, 2.2656250000000000, 1.1512586805555554]

    particle_values_ref ./= mesh.xmax[1]

    @test maximum(abs.(particle_values .- particle_values_ref)) < 1e-15
    @test maximum(abs.(particle_values1 .- particle_values_ref)) < 1e-15

end


text:
#!/usr/bin/env julia
using BinaryBuilder, Pkg

name = "SBML"
version = v"5.19.5"
sources = [
    ArchiveSource(
        "https://github.com/sbmlteam/libsbml/archive/v$(version).tar.gz",
        "6c0ec766e76bc6ad0c8626f3d208b4d9e826b36c816dff0c55e228206c82cb36"),
    DirectorySource("./bundled"),
]

script = raw"""
cd ${WORKSPACE}/srcdir/libsbml-*

atomic_patch -p1 ../patches/0001-User-lowercase-name-for-Windows-library.patch

mkdir build
cd build
cmake \
  -DCMAKE_INSTALL_PREFIX=${prefix} \
  -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
  -DCMAKE_LIBRARY_PATH=${libdir} \
  -DCMAKE_INCLUDE_PATH="${includedir}/libxml2;${includedir}" \
  -DENABLE_ARRAYS=ON \
  -DENABLE_COMP=ON \
  -DENABLE_DISTRIB=ON \
  -DENABLE_DYN=ON \
  -DENABLE_FBC=ON \
  -DENABLE_GROUPS=ON \
  -DENABLE_L3V2EXTENDEDMATH=ON \
  -DENABLE_LAYOUT=ON \
  -DENABLE_MULTI=ON \
  -DENABLE_QUAL=ON \
  -DENABLE_RENDER=ON \
  -DENABLE_REQUIREDELEMENTS=ON \
  -DENABLE_SPATIAL=ON \
  ..
make -j${nproc}
make install

# Remove large static library.
rm ${prefix}/lib/libsbml-static.a
"""

platforms = expand_cxxstring_abis(supported_platforms())

products = [
    LibraryProduct("libsbml", :libsbml),
]

dependencies = [
    Dependency("XML2_jll"),
    Dependency("Zlib_jll"),
]

# GCC 6 is necessary to work around https://gcc.gnu.org/bugzilla/show_bug.cgi?id=67557
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"6", julia_compat="1.6")


text:
"""
TestChannelData.jl defines the TestChannelData types and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Imports

using Serialization: serialize, deserialize
using Boid: Boid
using Boid: AbstractChannelData

# --- Type definitions

mutable struct TestChannelData <: AbstractChannelData
    data::Float64

    # Default constructor
    TestChannelData() = new()
end

# --- Method definitions

Boid.get_type(channel_data::TestChannelData) = typeof(channel_data.data)

Boid.get_data(channel_data::TestChannelData; encode::Bool=false) =
    encode ?
        encode_value(TestChannelData, channel_data.data) :
        channel_data.data

function Boid.set_data!(channel_data::TestChannelData, value;
                        decode::Bool=false)
    if decode
        channel_data.data = decode_value(TestChannelData, value)
    else
        channel_data.data = value
    end
end

function Boid.encode_value(::Type{TestChannelData}, data)::Vector{UInt8}
    packed_data = IOBuffer()
    serialize(packed_data, data)
    return take!(packed_data)
end

function Boid.decode_value(::Type{TestChannelData}, bytes::Vector{UInt8})
    packed_data = IOBuffer()
    write(packed_data, bytes)
    return deserialize(seekstart(packed_data))
end


text:
function patch_species_reaction(model::Dict{Symbol,Any}, pair::Pair)

    try

        # process the reaction table -
        reaction_table = model[:reactions]
        (number_of_reactions, number_of_fields) = size(reaction_table)
        for reaction_index = 1:number_of_reactions

            # get reaction_markup -
            reaction_markup = reaction_table[reaction_index, :reaction_markup]

            # get st dict -
            std = reaction_table[reaction_index, :stoichiometric_dictionary]

            # update -
            replace!(reaction_markup, pair.first => pair.second)
        end

    catch error

        # what is our error policy? => for now, just print the message
        error_message = sprint(showerror, error, catch_backtrace())
        println(error_message)

    end
end


text:
abstract type AbstractBC{T} <: AbstractDiffEqAffineOperator{T} end

abstract type AtomicBC{T} <: AbstractBC{T} end

"""
Robin, General, and in general Neumann, Dirichlet and Bridge BCs
are not necessarily linear operators.  Instead, they are affine
operators, with a constant term Q*x = Qa*x + Qb.
"""
abstract type AffineBC{T} <: AtomicBC{T} end

struct NeumannBC{N} end
struct Neumann0BC{N} end
struct DirichletBC{N} end
struct Dirichlet0BC{N} end

"""
q = PeriodicBC{T}()
Qx, Qy, ... = PeriodicBC{T}(size(u)) #When all dimensions are to be extended with a periodic boundary condition.
-------------------------------------------------------------------------------------
Creates a periodic boundary condition, where the lower index end of some u is extended with the upper index end and vice versa.
It is not reccomended to concretize this BC type in to a BandedMatrix, since the vast majority of bands will be all 0s. SpatseMatrix concretization is reccomended.
"""
struct PeriodicBC{T} <: AtomicBC{T}
    PeriodicBC(T::Type) = new{T}()
end

"""
  q = RobinBC(left_coefficients, right_coefficients, dx::T, approximation_order) where T # When this BC extends a dimension with a uniform step size
  q = RobinBC(left_coefficients, right_coefficients, dx::Vector{T}, approximation_order) where T # When this BC extends a dimension with a non uniform step size. dx should be the vector of step sizes for the whole dimension
-------------------------------------------------------------------------------------
  The variables in l are [αl, βl, γl], and correspond to a BC of the form αl*u(0) + βl*u'(0) = γl imposed on the lower index boundary.
  The variables in r are [αl, βl, γl], and correspond to an analagous boundary on the higher index end.
  Implements a robin boundary condition operator Q that acts on a vector to give an extended vector as a result
  Referring to (https://github.com/JuliaDiffEq/DiffEqOperators.jl/files/3267835/ghost_node.pdf)
  Write vector b̄₁ as a vertical concatanation with b0 and the rest of the elements of b̄ ₁, denoted b̄`₁, the same with ū into u0 and ū`. b̄`₁ = b̄`_2 = fill(β/Δx, length(stencil)-1)
  Pull out the product of u0 and b0 from the dot product. The stencil used to approximate u` is denoted s. b0 = α+(β/Δx)*s[1]
  Rearrange terms to find a general formula for u0:= -b̄`₁̇⋅ū`/b0 + γ/b0, which is dependent on ū` the robin coefficients and Δx.
  The non identity part of Qa is qa:= -b`₁/b0 = -β.*s[2:end]/(α+β*s[1]/Δx). The constant part is Qb = γ/(α+β*s[1]/Δx)
  do the same at the other boundary (amounts to a flip of s[2:end], with the other set of boundary coeffs)
"""
struct RobinBC{T, V<:AbstractVector{T}} <: AffineBC{T}
    a_l::V
    b_l::T
    a_r::V
    b_r::T
    function RobinBC(l::NTuple{3,T}, r::NTuple{3,T}, dx::T, order = 1) where {T}
        αl, βl, γl = l
        αr, βr, γr = r

        s = calculate_weights(1, one(T), Array(one(T):convert(T,order+1))) #generate derivative coefficients about the boundary of required approximation order

        a_l = -s[2:end]./(αl*dx/βl + s[1])
        a_r = s[end:-1:2]./(αr*dx/βr - s[1]) # for other boundary stencil is flippedlr with *opposite sign*

        b_l = γl/(αl+βl*s[1]/dx)
        b_r = γr/(αr-βr*s[1]/dx)

        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
    end
    function RobinBC(l::Union{NTuple{3,T},AbstractVector{T}}, r::Union{NTuple{3,T},AbstractVector{T}}, dx::AbstractVector{T}, order = 1) where {T}
        αl, βl, γl = l
        αr, βr, γr = r

        s_index = Array(one(T):convert(T,order+1))
        dx_l, dx_r = dx[1:length(s_index)], dx[(end-length(s_index)+1):end]

        s = calculate_weights(1, one(T), s_index) #generate derivative coefficients about the boundary of required approximation order
        denom_l = αl+βl*s[1]/dx_l[1]
        denom_r = αr-βr*s[1]/dx_r[end]

        a_l = -βl.*s[2:end]./(denom_l*dx_l[2:end])
        a_r = βr.*s[end:-1:2]./(denom_r*dx_r[1:(end-1)]) # for other boundary stencil is flippedlr with *opposite sign*

        b_l = γl/denom_l
        b_r = γr/denom_r

        return new{T, typeof(a_l)}(a_l, b_l, a_r, b_r)
    end
end

stencil(q::AffineBC{T}, N::Int) where T = ([transpose(q.a_l) transpose(zeros(T, N-length(q.a_l)))], [transpose(zeros(T, N-length(q.a_r))) transpose(q.a_r)])
affine(q::AffineBC) = (q.b_l, q.b_r)

stencil(q::PeriodicBC{T}, N::Int) where T= ([transpose(zeros(T, N-1)) one(T)], [one(T) transpose(zeros(T, N-1))])
affine(q::PeriodicBC{T}) where T = (zero(T), zero(T))
"""
q = GeneralBC(α_leftboundary, α_rightboundary, dx::T, approximation_order)

-------------------------------------------------------------------------------------

Implements a generalization of the Robin boundary condition, where α is a vector of coefficients.
Represents a condition of the form α[1] + α[2]u[0] + α[3]u'[0] + α[4]u''[0]+... = 0
Implemented in a similar way to the RobinBC (see above).
This time there are multiple stencils for multiple derivative orders - these can be written as a matrix S.
All components that multiply u(0) are factored out, turns out to only involve the first column of S, s̄0. The rest of S is denoted S`. the coeff of u(0) is s̄0⋅ᾱ[3:end] + α[2].
the remaining components turn out to be ᾱ[3:end]⋅(S`ū`) or equivalently (transpose(ᾱ[3:end])*S`)⋅ū`. Rearranging, a stencil q_a to be dotted with ū` upon extension can readily be found, along with a constant component q_b
"""
struct GeneralBC{T, L<:AbstractVector{T}, R<:AbstractVector{T}} <:AffineBC{T}
    a_l::L
    b_l::T
    a_r::R
    b_r::T
    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::T, order = 1) where {T}
        nl = length(αl)
        nr = length(αr)
        S_l = zeros(T, (nl-2, order+nl-2))
        S_r = zeros(T, (nr-2, order+nr-2))

        for i in 1:(nl-2)
            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx^i) #am unsure if the length of the dummy_x is correct here
        end

        for i in 1:(nr-2)
            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx^i)
        end
        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
        denomr = αr[2] .+ αr[3:end] ⋅ s0_r

        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)

        b_l = -αl[1]/denoml
        b_r = -αr[1]/denomr
        new{T, typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
    end

    function GeneralBC(αl::AbstractVector{T}, αr::AbstractVector{T}, dx::AbstractVector{T}, order = 1) where {T}

        nl = length(αl)
        nr = length(αr)
        dx_l, dx_r = (dx[1:(order+nl-2)], reverse(dx[(end-order-nr+3):end]))
        S_l = zeros(T, (nl-2, order+nl-2))
        S_r = zeros(T, (nr-2, order+nr-2))

        for i in 1:(nl-2)
            S_l[i,:] = [transpose(calculate_weights(i, one(T), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nl-2-i)))]./(dx_l.^i)
        end

        for i in 1:(nr-2)
            S_r[i,:] = [transpose(calculate_weights(i, convert(T, order+i), Array(one(T):convert(T, order+i)))) transpose(zeros(T, Int(nr-2-i)))]./(dx_r.^i)
        end
        s0_l = S_l[:,1] ; Sl = S_l[:,2:end]
        s0_r = S_r[:,end] ; Sr = S_r[:,(end-1):-1:1]

        denoml = αl[2] .+ αl[3:end] ⋅ s0_l
        denomr = αr[2] .+ αr[3:end] ⋅ s0_r

        a_l = -transpose(transpose(αl[3:end]) * Sl) ./denoml
        a_r = reverse(-transpose(transpose(αr[3:end]) * Sr) ./denomr)

        b_l = -αl[1]/denoml
        b_r = -αr[1]/denomr
        new{T,  typeof(a_l), typeof(a_r)}(a_l,b_l,a_r,b_r)
    end
end

#implement Neumann and Dirichlet as special cases of RobinBC
NeumannBC(α::NTuple{2,T}, dx::Union{AbstractVector{T}, T}, order = 1) where T = RobinBC((zero(T), one(T), α[1]), (zero(T), one(T), α[2]), dx, order)
DirichletBC(αl::T, αr::T) where T = RobinBC((one(T), zero(T), αl), (one(T), zero(T), αr), one(T), 2one(T) )
#specialized constructors for Neumann0 and Dirichlet0
Dirichlet0BC(T::Type) = DirichletBC(zero(T), zero(T))
Neumann0BC(dx::Union{AbstractVector{T}, T}, order = 1) where T = NeumannBC((zero(T), zero(T)), dx, order)

# other acceptable argument signatures
#RobinBC(al::T, bl::T, cl::T, dx_l::T, ar::T, br::T, cr::T, dx_r::T, order = 1) where T = RobinBC([al,bl, cl], [ar, br, cr], dx_l, order)

Base.:*(Q::AffineBC, u::AbstractVector) = BoundaryPaddedVector(Q.a_l ⋅ u[1:length(Q.a_l)] + Q.b_l, Q.a_r ⋅ u[(end-length(Q.a_r)+1):end] + Q.b_r, u)
Base.:*(Q::PeriodicBC, u::AbstractVector) = BoundaryPaddedVector(u[end], u[1], u)

Base.size(Q::AtomicBC) = (Inf, Inf) #Is this nessecary?

gettype(Q::AbstractBC{T}) where T = T


text:
using SymSemiseparableMatrices
using Test
using LinearAlgebra
import SymSemiseparableMatrices: spline_kernel, spline_kernel_matrix

@testset "SymSemiseparableMatrices.jl" begin
    include("test_symsep.jl")
    include("test_symsepchol.jl")
    include("test_diasymsep.jl")
    include("test_diasymsepchol.jl")
end


text:
#=
  A/B test

  From http://rpubs.com/rasmusab/exercise_2_bayesian_ab_testing
  """
  Exercise 2: Bayesian A/B testing for Swedish Fish Incorporated with Stan
  Rasmus Bååth

  Swedish Fish Incorporated is the largest Swedish company delivering fish by mail order,
  but you probably already knew that. The marketing department have done a pilot study and
  tried two different marketing methods:

  A: Sending a mail with a colorful brochure that invites people to sign up for a one year
     salmon subscription.

  B: Sending a colorful brochure that invites people to sign up for a one year
     salmon subscription nd that includes a free salmon.

  The marketing department sent out 16 mails of type A and 16 mails of type B.
  Six Danes that received a mail of type A signed up for one year of salmon, and
  ten Danes that received a mail of type B signed up!

  The marketing department now wants to know, which method should we use, A or B?

  At the bottom of this document you’ll find a solution. But try yourself first!
  Question I: Build a Bayesian model in Stan that answers the question: What is the
  probability that method B is better than method A?
  """

  See ~/stan/ab_testing.R
      ~/blog/ab_testing.blog
      ~/psi/ab_testing.psi
      ~/webppl/ab_testing.wppl
=#

using Turing, StatsPlots, DataFrames
# using ReverseDiff, Zygote, Tracker
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

#=
ummary Statistics
  parameters      mean       std   naive_se      mcse          ess      rhat
      Symbol   Float64   Float64    Float64   Float64      Float64   Float64

       rateA    0.3893    0.1126     0.0003    0.0009   15774.8462    1.0004
       rateB    0.6121    0.1108     0.0003    0.0008   16943.7123    1.0001

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

       rateA    0.1839    0.3088    0.3857    0.4659    0.6183
       rateB    0.3847    0.5397    0.6153    0.6918    0.8148

mean(rateB > rateA): 0.91461875
mean(rateB - rateA): 0.22276908208667892
=#
@model function ab_testing(nA=16,nB=16,obsSA=6, obsSB=10)
    # nA = 16 # Number of sent mail
    # nB = 16

    # obsSA = 6; # Number of signments (observed)
    # obsSB = 10

    rateA ~ Beta(1,1) # priors
    rateB ~ Beta(1,1)

    obsSA ~ Binomial(nA,rateA) # likelihood
    obsSB ~ Binomial(nB,rateB)

    # This variant don't work!
    # sA ~ Binomial(nA,rateA) # likelihood
    # sB ~ Binomial(nB,rateB)
    # sA == obsSA || begin Turing.@addlogprob! -Inf; return; end
    # sB == obsSB || begin Turing.@addlogprob! -Inf; return; end

end

model = ab_testing(16,16,6,10)
# model = ab_testing(160,160,60,100)
num_chains = 4
# chains = sample(model, Prior(), MCMCThreads(), 1000, num_chains)
chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)
# chains = sample(model, MH(), 40_000)

# chains = sample(model, PG(20), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, Gibbs(MH(:sA,:sB),NUTS(1000,0.65,:rateA,:rateB)), MCMCThreads(), 40_000, num_chains)

display(chains)
# display(plot(chains))

df = DataFrame(chains);
println("mean(rateB > rateA): ", mean(df[!,"rateB"] .> df[!,"rateA"]))
println("mean(rateB - rateA): ", mean(df[!,"rateB"] .- df[!,"rateA"]))

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "LibRaw"
version = v"0.20.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource(
        "https://www.libraw.org/data/LibRaw-0.20.2.tar.gz",
        "dc1b486c2003435733043e4e05273477326e51c3ea554c6864a4eafaff1004a6",
    ),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd LibRaw-*
autoreconf --install
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
exit
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("mem_image", :mem_image),
    ExecutableProduct("postprocessing_benchmark", :postprocessing_benchmark),
    LibraryProduct("libraw", :libraw),
    ExecutableProduct("raw-identify", :raw_identify),
    ExecutableProduct("4channels", :_4channels),
    ExecutableProduct("half_mt", :half_mt),
    ExecutableProduct("multirender_test", :multirender_test),
    LibraryProduct("libraw_r", :libraw_r),
    ExecutableProduct("dcraw_emu", :dcraw_emu),
    ExecutableProduct("simple_dcraw", :simple_dcraw),
    ExecutableProduct("unprocessed_raw", :unprocessed_raw),
    ExecutableProduct("rawtextdump", :rawtextdump),
    ExecutableProduct("dcraw_half", :dcraw_half),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(; name="Zlib_jll"))
    Dependency(PackageSpec(; name="JpegTurbo_jll"))
    Dependency(PackageSpec(; name="JasPer_jll"))
    Dependency(PackageSpec(; name="CompilerSupportLibraries_jll"))
    Dependency(PackageSpec(; name="LittleCMS_jll"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(
    ARGS,
    name,
    version,
    sources,
    script,
    platforms,
    products,
    dependencies;
    julia_compat="1.6",
    preferred_gcc_version=v"6",
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
function _get_nbfun(weights)
    if weights isa MainPathEdgeWeight
        nbfun = emax_neighbors
    elseif weights isa MainPathVertexWeight
        nbfun = vmax_neighbors
    else
        error("Unknown weights type.")
    end

    nbfun
end

abstract type MainPathTraversal end

"""
    ForwardBackwardLocal(start)
Struct representing a main path traversal in both directions of the
input DAG, starting at `start`.
"""
struct ForwardBackwardLocal{T} <: MainPathTraversal
    start::Vector{T}
end

function (x::ForwardBackwardLocal)(g, weight)
    nbfun = _get_nbfun(weight)
    w = weight(g)

    parents_forward  = bfs_multi(g, x.start, (g, v) -> nbfun(g, v, w, outneighbors))
    parents_backward = bfs_multi(g, x.start, (g, v) -> nbfun(g, v, w, inneighbors))

    mp = SimpleDiGraph(nv(g))

    for (v, us) in enumerate(parents_forward)
        for u in us
            u == v && continue
            add_edge!(mp, u, v)
        end
    end

    for (v, us) in enumerate(parents_backward)
        for u in us
            u == v && continue
            add_edge!(mp, v, u)
        end
    end

    mp
end

"""
    ForwardLocal(start)
Struct representing a main path traversal against the edge direction of the
input DAG, starting at `start`.
"""
struct ForwardLocal{T} <: MainPathTraversal
    start::Vector{T}
end

function (x::ForwardLocal)(g, weights)
    nbfun = _get_nbfun(weights)
    w = weights(g)

    parents_forward  = bfs_multi(g, x.start, (g, v) -> nbfun(g, v, w, outneighbors))

    mp = SimpleDiGraph(nv(g))

    for (v, us) in enumerate(parents_forward)
        for u in us
            u == v && continue
            add_edge!(mp, u, v)
        end
    end

    mp
end

"""
    BackwardLocal(start)
Struct representing a main path traversal following the edge direction of the
input DAG, starting at `start`.
"""
struct BackwardLocal{T} <: MainPathTraversal
    start::Vector{T}
end

function (x::BackwardLocal)(g, weights)
    nbfun = _get_nbfun(weights)
    w = weights(g)

    parents_backward  = bfs_multi(g, x.start, (g, v) -> nbfun(g, v, w, inneighbors))

    mp = SimpleDiGraph(nv(g))

    for (v, us) in enumerate(parents_backward)
        for u in us
            u == v && continue
            add_edge!(mp, v, u)
        end
    end

    mp
end


text:
export forecast, simulate

"""
    forecast(ss::StateSpace{Typ}, N::Int) where Typ

Obtain the minimum mean square error forecasts N steps ahead. Returns the forecasts and the predictive distributions
at each time period.
"""
function forecast(ss::StateSpace{Typ}, N::Int) where Typ

    # Load estimated covariance matrices
    H = ss.model.H
    Q = ss.model.Q

    # Load system
    n, p, m, r    = size(ss.model)
    Z, T, R, d, c = prepare_forecast(ss, N)

    # Load a, c, P, and F at last in-sample instant
    a0 = ss.smoother.alpha[end, :]
    c0 = ss.model.c[end, :]
    P0 = ss.filter.P[:, :, end]
    F0 = ss.filter.F[:, :, end]

    # State and variance forecasts
    a = Matrix{Typ}(undef, N, m)
    P = Array{Typ, 3}(undef, m, m, N)
    F = Array{Typ, 3}(undef, p, p, N)

    # Probability distribution
    dist = Vector{Distribution}(undef, N)

    # Initialization
    a[1, :]    = T*a0 + c0
    P[:, :, 1] = T*P0*T' + R*Q*R'
    F[:, :, 1] = Z[:, :, 1]*P[:, :, 1]*Z[:, :, 1]' + H
    ensure_pos_sym!(F, 1)
    dist[1]    = MvNormal(vec(Z[:, :, 1]*a[1, :] + d[1, :]), F[:, :, 1])

    for t = 2:N
        a[t, :]    = T*a[t-1, :] + c[t-1, :]
        P[:, :, t] = T*P[:, :, t-1]*T' + R*Q*R'
        F[:, :, t] = Z[:, :, t]*P[:, :, t]*Z[:, :, t]' + H
        ensure_pos_sym!(F, t)
        dist[t]    = MvNormal(vec(Z[:, :, t]*a[t, :] + d[t, :]), F[:, :, t])
    end

    forec = Matrix{Typ}(undef, N, p)
    for t = 1:N
        forec[t, :] = mean(dist[t])
    end

    return forec, dist
end

"""
    simulate(ss::StateSpace{Typ}, N::Int, S::Int) where Typ

Simulate S future scenarios up to N steps ahead. Returns a p x N x S matrix where the dimensions represent, respectively,
the number of series in the model, the number of steps ahead, and the number of scenarios.
"""
function simulate(ss::StateSpace{Typ}, N::Int, S::Int) where Typ

    # Load estimated covariance matrices
    H = ss.model.H
    Q = ss.model.Q

    # Load system
    n, p, m, r    = size(ss.model)
    Z, T, R, d, c = prepare_forecast(ss, N)

    # Distribution of the state space errors
    dist_ϵ = MvNormal(zeros(p), H)
    dist_η = MvNormal(zeros(r), Q)

    αsim = Array{Typ, 3}(undef, N, m, S)
    ysim = Array{Typ, 3}(undef, N, p, S)

    for s = 1:S

        # Sampling errors
        ϵ = rand(dist_ϵ, N)'
        η = rand(dist_η, N)'

        # Initializing simulation
        αsim[1, :, s] = T*ss.smoother.alpha[n, :] + c[1, :] + R*η[1, :]
        ysim[1, :, s] = Z[:, :, 1]*αsim[1, :, s] + d[1, :] + ϵ[1, :]

        # Simulating future scenarios
        for t = 2:N
            αsim[t, :, s] = T*αsim[t-1, :, s] + c[t-1, :] + R*η[t, :]
            ysim[t, :, s] = Z[:, :, t]*αsim[t, :, s] + d[t, :] + ϵ[t, :]
        end
    end

    return ysim
end

"""
    prepare_forecast(ss::StateSpace{Typ}, N::Int) where Typ

Adjust matrix Z for forecasting and check for dimension errors.
"""
function prepare_forecast(ss::StateSpace{Typ}, N::Int) where Typ

    # Load system
    n, p, m, r = size(ss.model)
    Z0, T, R   = ztr(ss.model)
    d0 = ss.model.d
    c0 = ss.model.c

    Z = Array{Typ, 3}(undef, p, m, N)
    d = Matrix{Typ}(undef, N, p)
    c = Matrix{Typ}(undef, N, m)
    if ss.model.mode == "time-invariant"
        for t in 1:N, i in 1:p, j in 1:m
            Z[i, j, t] = Z0[i, j, 1]
            d[t, i] = d0[1, i]
            c[t, j] = c0[1, j]
        end
    else
        size(Z0, 3) < n + N && error("Time-variant Z too short for forecasting $N steps ahead")
        size(d0, 1) < n + N && error("Time-variant d too short for forecasting $N steps ahead")
        size(c0, 1) < n + N && error("Time-variant c too short for forecasting $N steps ahead")

        for t in n + 1:n + N, i in 1:p, j in 1:m
            Z[i, j, t - n] = Z0[i, j, t]
            d[t - n, i] = d0[t, i]
            c[t - n, j] = c0[t, j]
        end
    end

    return Z, T, R, d, c
end


text:
##############################################
##      Cairn Server Handlers and Interface ##
##############################################
module Server
import Cairn
import Rock
import SQLite
import NullableArrays
using Logging
@Logging.configure(level=DEBUG)

function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, v::Cairn.RPC.GetID)
    Rock.Server.self(px)
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, v::Cairn.RPC.PrepareCache)

    SQLite.query(d,
                 """
                 insert or replace into results_cache (uuid, val, state) values  (?, ?, ?)
                 """, [v.id, Nullable(), Cairn.RPC.Loading])
    true
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, v::Cairn.RPC.FetchCache)
    while true
        data = SQLite.query(d,
                            """
                            select val, state from results_cache where (uuid = ?) limit 1
                            """, [v.id]).data
        val = NullableArrays.dropnull(data[1])
        state = NullableArrays.dropnull(data[2])
        if length(val) == 1
            try
                if state[1] == Cairn.RPC.Ready
                    return Cairn.RPC.CacheResult(val[1], state[1])
                end
            finally
                if v.evict && state[1] == Cairn.RPC.Ready
                    SQLite.query(d,
                                 """
                                 delete from results_cache where (uuid = ?)
                                 """, [v.id])
                end
            end
        else
            Cairn.RPC.CacheResult(0, Gone)
        end
        wait(Timer(1))
    end
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, v::Cairn.RPC.CacheResultOfAction)
    r = CairnActionHandler(px, d, v.c)
    if Rock.Server.self(px) == v.view_on
        vec = Vector(3)
        vec[1] = v.id
        vec[2] = r
        vec[3] = Cairn.RPC.Ready
        SQLite.query(d,
                     """
                     insert or replace into results_cache (uuid, val, state) values (?,?,?)
                     """, vec)
    end
    ()
end

dropnull2(data)= zip(data[1] |> NullableArrays.dropnull, data[2] |> NullableArrays.dropnull) |> collect

function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.CreateRequest)
    try
        q = SQLite.query(d,
                         """
                         select host, port from server_nodes order by RANDOM() limit ?
                         """, [arg.replication]).data |> dropnull2
        SQLite.query(d,
                     """
                     insert into objects values (?,?)
                     """, [arg.name, arg.replication])
        s = SQLite.Stmt(d,
                        """
                        insert into object_replicas values (?,?,?)
                        """)
        for (host, port) = q
            SQLite.bind!(s, [arg.name, host, port])
            SQLite.execute!(s)
        end
        if length(q) < arg.replication
            Cairn.RPC.CreateResponse(q,Nullable(ErrorException("Could not grant desired replication factor")))
        else
            @debug q
            Cairn.RPC.CreateResponse(q, Nullable())
        end
    catch err
        if isa(err, SQLite.SQLiteException)
            if err.msg == "UNIQUE constraint failed: objects.name"
                SQLite.query(d,
                             """
                             select host, port from object_replicas where object_name = ?
                             """, [arg.name]).data |> dropnull2
                return Cairn.RPC.CreateResponse([], Nullable(ErrorException("File $(arg.name) Already Exists")))
            else
                rethrow(err)
            end
        else
            rethrow(err)
        end

    end

end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.GetKeySpanRequest)
    for key = arg.key_span
        # Expand
        items = if isa(key, Cairn.RPC.PrefixKey)
            SQLite.query(d,
                                """
                                select object_name from chunks where object_name LIKE ? || '%'
                                """, [key.prefix]).data[1] |> NullableArrays.dropnull

        else
            [key]
        end

        # Get
        response = Cairn.RPC.GetKeySpanResponse([])
        for key = items
            data = SQLite.query(d,
                                """
                                select host, port, hash from chunks where object_name = ?
                                group by hash order by id ASC
                                """, [key]).data
            # TODO: ^^^ Verify
            # TODO: ^^^ Return random replica?

            replicas = data |> dropnull2
            hashes = data[3] |> NullableArrays.dropnull
            put!(response.files, Cairn.RPC.GetKeySpanSubResponse(key, hashes, replicas, Nullable()))

        end
        response
    end

end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.DeleteRequest)
    # Select all the hashes such that per host:port, there is one
    # object depending on that Chunk
    deleteable =  SQLite.query(d,
                               """
                               WITH HASHES AS (SELECT *  FROM chunks GROUP BY hash, host, port HAVING
                               count(hash) = 1)
                               SELECT hash, host, port from HASHES where object_name = ?
                               """, [arg.name]).data
    hashes = deleteable[3] |> NullableArrays.dropnull
    replicas = deleteable |> dropnull2
    # TODO: ^^^ DELETE
    SQLite.query(d,
                 """
                 delete from objects where name = ?
                 """, [arg.name])
    Cairn.RPC.DeleteRequest(Nullable())
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.AppendChunksRequest)
    data = SQLite.query(d,
                        """
                        select name, host, port from objects where name = ? limit 1
                        """, [arg.name]).data
    name = (data[1] |> NullableArrays.dropnull |> collect)[1]
    host = (data[2] |> NullableArrays.dropnull |> collect)[1]
    port = (data[3] |> NullableArrays.dropnull |> collect)[1]

    s = SQLite.Stmt(d,
                    """
                    insert into chunks values (?,?,?,?)
                    """)
    for hash = arg.hashes
        SQLite.bind!(s, [name, hash, host, port])
        SQLite.execute!(s)
    end
    Cairn.RPC.AppendChunksResponse(Nullable())
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.ReplaceRequest)
    CairnActionHandler(px, d, Cairn.RPC.DeleteRequest(arg.old))
    SQLite.query(d,
                 """
                 update objects set name = ? where  where name = ?
                 """, [arg.new, arg.old])

    Cairn.RPC.ReplaceRequest(Nullable())
end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.RenameRequest)
    SQLite.query(d,
                 """
                 update objects set name = ? where name = ?
                 """, [arg.new, arg.old])
    Cairn.RPC.RenameResponse(Nullable())

end
function CairnActionHandler(px::Rock.Server.t, d::SQLite.DB, arg::Cairn.RPC.TruncateRequest)
    TruncateResponse(Nullable(ErrorException("Not Implemented")))
end

function CairnActionHandler(px::Rock.Server.t,d::SQLite.DB, arg::Cairn.RPC.ManyCommandRequest)
    Cairn.RPC.ManyCommandResponse([CairnActionHandler(px, d, arg) for arg=args])
end

function CairnActionHandler(px::Rock.Server.t,d::SQLite.DB, a::Cairn.RPC.AddStorageNode)
    SQLite.query(d, "insert or replace into server_nodes values (?, ?)", [a.host, a.port])
end

immutable t
    user_defined_function::Function
    rock::Rock.Server.t
end
function make(port,id, peers::Vector{Tuple{ASCIIString, Int64}}, path::ASCIIString)
    d = SQLite.DB("$path/sqlite.db")
    SQLite.query(d,
                 """
                 create table if not exists server_nodes (host TEXT, port Int,  unique (host, port))
                 """)

    SQLite.query(d,
                 """
                 create table if not exists objects
                 (name TEXT PRIMARY KEY, replication INTEGER)
                 """)

    SQLite.query(d,
                 """
                 create table if not exists object_replicas
                 (object_name TEXT, host TEXT, port Int,
                 unique (object_name, host, port),
                 FOREIGN KEY (object_name) REFERENCES objects(name) ON DELETE CASCADE ON UPDATE CASCADE,
                 FOREIGN KEY (host, port) REFERENCES server_nodes(host, port) ON DELETE CASCADE ON UPDATE CASCADE)
                 """)

    SQLite.query(d,
                 """ create table if not exists chunks
                 (id INTEGER PRIMARY KEY AUTOINCREMENT, object_name TEXT, hash TEXT,
                 host TEXT, port TEXT,
                 FOREIGN KEY (object_name) REFERENCES objects(name) ON DELETE CASCADE ON UPDATE CASCADE,
                 FOREIGN KEY (host, port) REFERENCES server_nodes(host, port) ON DELETE CASCADE ON UPDATE CASCADE)
                 """)

    SQLite.query(d,
                 """
                 create table if not exists results_cache (uuid TEXT PRIMARY KEY, val BLOB, state INT)
                 """)
    ctx = Rock.Server.Context(port, peers, CairnActionHandler, "$path/sqlite.db")
    t(()->(), Rock.Server.t(ctx, id))
end
function start(cairn::t)
    Rock.Server.start(cairn.rock)
end
end


text:
struct YaoInterpreter <: AbstractInterpreter
    native_interpreter::Core.Compiler.NativeInterpreter
    passes::Vector{Symbol}
end

default_passes() = [:zx]

YaoInterpreter(; passes::Vector{Symbol} = default_passes()) =
    YaoInterpreter(Core.Compiler.NativeInterpreter(), passes)

InferenceParams(interp::YaoInterpreter) = InferenceParams(interp.native_interpreter)
OptimizationParams(interp::YaoInterpreter) = OptimizationParams(interp.native_interpreter)
Core.Compiler.get_world_counter(interp::YaoInterpreter) = get_world_counter(interp.native_interpreter)
Core.Compiler.get_inference_cache(interp::YaoInterpreter) =
    get_inference_cache(interp.native_interpreter)
Core.Compiler.code_cache(interp::YaoInterpreter) = Core.Compiler.code_cache(interp.native_interpreter)
Core.Compiler.may_optimize(interp::YaoInterpreter) =
    Core.Compiler.may_optimize(interp.native_interpreter)
Core.Compiler.may_discard_trees(interp::YaoInterpreter) =
    Core.Compiler.may_discard_trees(interp.native_interpreter)
Core.Compiler.may_compress(interp::YaoInterpreter) =
    Core.Compiler.may_compress(interp.native_interpreter)
Core.Compiler.unlock_mi_inference(interp::YaoInterpreter, mi::Core.MethodInstance) =
    Core.Compiler.unlock_mi_inference(interp.native_interpreter, mi)
Core.Compiler.lock_mi_inference(interp::YaoInterpreter, mi::Core.MethodInstance) =
    Core.Compiler.lock_mi_inference(interp.native_interpreter, mi)
Core.Compiler.add_remark!(interp::YaoInterpreter, st::Core.Compiler.InferenceState, msg::String) =
    nothing # println(msg)

function Core.Compiler.abstract_call(
    interp::YaoInterpreter,
    fargs::Union{Nothing,Vector{Any}},
    argtypes::Vector{Any},
    sv::InferenceState,
    max_methods::Int = InferenceParams(interp).MAX_METHODS,
)

    ft = argtypes[1]
    if isa(ft, Const)
        f = ft.val
    elseif isconstType(ft)
        f = ft.parameters[1]
    elseif isa(ft, DataType) && isdefined(ft, :instance)
        f = ft.instance
    else
        # non-constant function, but the number of arguments is known
        # and the ft is not a Builtin or IntrinsicFunction
        if Core.Compiler.typeintersect(widenconst(ft), Core.Builtin) != Union{}
            Core.Compiler.add_remark!(interp, sv, "Could not identify method table for call")
            return Core.Compiler.CallMeta(Any, false)
        end
        return Core.Compiler.abstract_call_gf_by_type(
            interp,
            nothing,
            argtypes,
            argtypes_to_type(argtypes),
            sv,
            max_methods,
        )
    end

    allconst = true
    for x in argtypes
        if !(x isa Const)
            allconst = false
        end
    end

    if f isa Function && parentmodule(f) === Semantic
        return abstract_call_quantum(interp, f, fargs, argtypes, sv, max_methods)
    elseif (f isa UnionAll || f isa DataType) && f <: IntrinsicRoutine && allconst
        # force intrinsic gates pure
        callinfo = Core.Compiler.abstract_call_known(interp, f, fargs, argtypes, sv, max_methods)
        return CallMeta(callinfo.rt, MethodResultPure())
    elseif f === Locations && allconst
        return CallMeta(Const(f(argtypes[2].val)), MethodResultPure())
    elseif f === CtrlLocations && allconst
        la = length(fargs) - 1
        if la == 1
            loc = f(argtypes[2].val)
        elseif la == 2
            loc = f(argtypes[2].val, argtypes[3].val)
        else
            error("Invalid CtrlLocation statement")
        end
        return CallMeta(Const(loc), MethodResultPure())
    elseif f === getindex && allconst
        # getindex(locs, locs)
        la = length(argtypes)
        a = argtypes[2].val
        b = argtypes[3].val
        if la == 3 && a isa AbstractLocations && b isa AbstractLocations
            return CallMeta(Const(f(argtypes[2].val, argtypes[3].val)), MethodResultPure())
        end
    end

    return Core.Compiler.abstract_call_known(interp, f, fargs, argtypes, sv, max_methods)
end

function abstract_call_quantum(
    interp::AbstractInterpreter,
    @nospecialize(f),
    fargs::Union{Nothing,Vector{Any}},
    argtypes::Vector{Any},
    sv::InferenceState,
    max_methods::Int = InferenceParams(interp).MAX_METHODS,
)

    if f === Semantic.measure
        return Core.Compiler.CallMeta(Int, MethodResultPure())
    elseif f === Semantic.gate || f === Semantic.ctrl
        gt = widenconst(argtypes[2])
        if gt <: IntrinsicRoutine
            callinfo = CallMeta(Const(nothing), MethodResultPure())
        else
            callinfo = Core.Compiler.abstract_call_known(interp, f, fargs, argtypes, sv, max_methods)
        end
        return CallMeta(callinfo.rt, nothing)
    else
        return Core.Compiler.CallMeta(Const(nothing), MethodResultPure())
    end
end

function is_semantic_fn_call(e)
    return e isa Expr &&
           e.head === :call &&
           e.args[1] isa GlobalRef &&
           e.args[1].mod === YaoCompiler.Semantic
end

function convert_to_quantum_head!(ci::CodeInfo)
    for (v, e) in enumerate(ci.code)
        ci.code[v] = convert_to_quantum_head(e)
    end
    return ci
end


text:
using Distances
using NearestNeighbors
using StaticArrays

struct NeighborList
    i::Vector{Vector{Int}}
    j::Vector{Vector{Int}}
    R::Vector{Vector{Float64}}
    r::Vector{Vector{SVector{3,Float64}}}
end
length(nn::NeighborList) = length(nn.i)

function get_distance(L::SVector{3,<:AbstractFloat}, x::SVector{3,<:AbstractFloat}, y::SVector{3,<:AbstractFloat})
    broadcast(L, x, y) do L, x, y
        d = y - x
        2mod(abs(d), L) < L ? d : L - d
    end
end

to_array(A::AbstractSystem) = [SVector{3}(ustrip.(p)) for p ∈ position(A)]

function neighborlist(A::AbstractSystem{3}, rcutoff::Float64)
    # Convert Positions to Matrix for Ball tree
    X = to_array(A)

    # Create Metric for Periodic Boundary Conditions
    periodic = periodicity(A)
    L = @SVector [periodic[i] ? ustrip(bounding_box(A)[i][i]) : Inf for i = 1:3]
    d = Distances.PeriodicEuclidean(L)

    # Build Ball tree
    tree = BallTree(X, d)

    # Intialize empty vectors
    i = Vector{Int64}[] # i
    j = Vector{Int64}[] # j
    R = Vector{Float64}[] # Distances
    r = Vector{SVector{3,Float64}}[] # Positions

    # Fill vectors
    for n = 1:length(X)
        neighbors = inrange(tree, X[n], rcutoff, true)
        itemp = Int64[]
        jtemp = Int64[]
        Rtemp = Float64[]
        rtemp = SVector{3,Float64}[]
        for neighbor in neighbors
            if neighbor != n
                push!(itemp, n)
                push!(jtemp, neighbor)
                rr = get_distance(L, X[n], X[neighbor])
                push!(rtemp, rr)
                push!(Rtemp, norm(rr))
            end
        end

        push!(i, itemp)
        push!(j, jtemp)
        push!(R, Rtemp)
        push!(r, rtemp)
    end
    NeighborList(i, j, R, r)
end


text:
"""
    add_orfs!(::Dict{String, IncompleteSegmentData}, inf_aa_dat_path::String, inf_dat_path::String)

Parse the files `influenza_aa.dat` and `influenza.dat`, extract information
about ORFs, and update the `IncompleteSegmentData` `proteins` field for all segment data
present in the input dict. Returns the number of updated `IncompleteSegmentData` instances.
"""
function add_orfs!(
    segment_data::Dict{String, IncompleteSegmentData},
    inf_aa_dat_path::String,
    inf_dat_path::String,
)
    accession_protein_map = maybe_gzip(parse_inf_aa, inf_aa_dat_path)
    maybe_gzip(inf_dat_path) do io
        add_orfs!(segment_data, accession_protein_map, io)
    end
end

function try_parse(::Type{Protein}, s::AbstractString)::Option{Protein}
    p = tryparse(Protein, s)
    p === nothing ? none : some(p)
end

# Parse the inf_aa file and get a map from accession => protein
function parse_inf_aa(io::IO)::Dict{String, Protein}
    result = Dict{String, Protein}()
    fields = Vector{SubString{String}}(undef, 11)
    for line in eachline(io) |> ifilter(ln -> !isempty(strip(ln)))
        nfields = try_split!(fields, line, '\t')
        if unwrap_or(nfields, 0) != 11
            error("Expected 11 fields, got", line)
        end
        accession = first(fields)
        variant = @unwrap_or try_parse(Protein, fields[3]) continue
        haskey(result, accession) && error("Duplicate key $accession")
        result[accession] = variant
    end
    result
end

function add_orfs!(segment_data::Dict{String, IncompleteSegmentData}, accession_protein_map::Dict{String, Protein}, io::IO)
    n_updates = 0
    proteinbuffer = ReferenceProtein[]
    fields = Vector{SubString{String}}(undef, 128)
    for line in eachline(io) |> imap(strip) |> ifilter(!isempty)
        nfields = unwrap(try_split!(fields, line, '\t'))
        isodd(nfields) || error("Must be an odd-numbered fields")
        gb_accession = first(fields)
        haskey(segment_data, gb_accession) || continue
        data = segment_data[gb_accession]
        n_updates += 1
        is_bad = false
        empty!(proteinbuffer)
        for (protein_accession, orf_field) in zip(@view(fields[2:2:nfields]), @view(fields[3:2:nfields]))

            # For some reason, some of these protein accessions are not actually present
            # in the influenza_aa.dat file.
            if !haskey(accession_protein_map, protein_accession)
                is_bad = true
                break
            end
            variant = accession_protein_map[protein_accession]
            m_orfs = parse_orf_field(orf_field)
            if is_error(data.seq) || is_error(m_orfs)
                is_bad = true
                break
            end

            orfs, seq = unwrap(m_orfs), unwrap(data.seq)

            # NCBI includes stop codon in the ORF. I don't. So if the last three bases
            # of the ORF is a stop, remove them from the ORF.
            # We assume last ORF is 3 bp in length
            if (
                !isempty(orfs) &&
                length(last(orfs)) > 3
            )
                range = last(orfs)[end-2:end]
                if (
                    checkbounds(Bool, eachindex(seq), range) &&
                    all(!BioSequences.isambiguous, seq[range]) &&
                    Influenza.is_stop(BioSequences.DNACodon(seq[range]))
                )
                    orfs[end] = orfs[end][1:end-3]
                end
            end
            push!(proteinbuffer, ReferenceProtein(variant, orfs))
        end
        if !is_bad
            resize!(data.proteins, length(proteinbuffer))
            copy!(data.proteins, proteinbuffer)
        end
    end
    return n_updates
end

function parse_orf_field(s::Union{String, SubString{String}})::Option{Vector{UnitRange{UInt32}}}
    # If it looks like gb|AB266090:<411->632, the ORF is not present
    # in the reference, and we skip it
    if occursin('>', s) || occursin('<', s)
        return none
    end

    # Multiple ORFs in a gene makes them enclosed in brackets
    # e.g. (gb|AB212651:26-52, 741-1007)
    s1 = strip(s, ['(', ')'])
    p = findfirst(isequal(':'), s1)
    p === nothing && return none
    s2 = @view s1[p+1:ncodeunits(s1)]

    orfstrings = split(s2, ", ")
    result = map(parse_range, orfstrings)

    # Return none if it's not divisible by three and thus cant be ORF
    iszero(sum(length, result, init=0) % 3) || return none
    some(result)
end

function parse_range(s::AbstractString)::UnitRange{UInt32}
    p2 = findfirst(isequal('-'), s)
    # If it's not a range, it must be a single number
    range = if p2 === nothing
        n = parse(UInt32, s)
        n:n
    else
        start = parse(UInt32, @view s[1:p2-1])
        stop = parse(UInt32, s[p2+1:ncodeunits(s)])
        start:stop
    end
    iszero(first(range)) && error("ORF range \"$s\" cannot contain the position zero.")
    return range
end

# These sequences are subject to change. They are the conserved DNA sequences
# at the termini of Influenza. Sequences are from DOI: 10.1128/JCM.03265-13 and
# from our internal OneTube inf A PCR protocol.
# The first sequence in each vector are from our own (and are thus the same), the latter
# are from the article.
const TERMINI = let
    rc = reverse_complement
    d = [
        Segments.PB2 => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCGGAGC"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGAAACACGAGC")]
        ),
        Segments.PB1 => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCGGAGC"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGAAACACGAGC")]
        ),
        Segments.PA  => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCGGTGC"],
            [dna"CTTGTTTCTCCT", rc(dna"GTAGAAACACGTGC")]
        ),
        Segments.HA  => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCAGAGC"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGTAACAAGAGC")]
        ),
        Segments.NP  => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCACAGC"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGAAACAACAGC")]
        ),
        Segments.NA  => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCAGAGC"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGTAACAAGAGC")]
        ),
        Segments.MP => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCACGCACTT", dna"AGCAGAAGCAGGCACTT"],
            [dna"CTTGTTTCTCCT", rc(dna"GTAGAAACAACGCACTT")]
        ),
        Segments.NS  => (
            [dna"AGCAAAAGCAGG", dna"AGCAGAAGCAGAGGATT"],
            [dna"CTTGTTTCTCCT", rc(dna"AGTAGTAACAAGAGGATT")]
        )
    ]
    map(d) do (s, t)
        s => (
            map(sq -> BioSequences.ApproximateSearchQuery(sq, :forward), t[1]),
            map(sq -> BioSequences.ApproximateSearchQuery(sq, :backward), t[2])
        )
    end |> Dict
end

"""
    strip_false_termini!(data::Dict{String, IncompleteSegmentData}, filter_termini=true)

Some sequences given by NCBI are too long, containing linker or primer contamination.
True influenza sequences begin with AGCAAAAGCAGG and ends with CTTGTTTCTACT. These
termini are highly conserved, but may have a single substitution.
We can allow sequences that do not include these noncoding ends, but we cannot
allow artifacts flanking the true termini.
If `filter_termini`, remove sequences where the termini cannot be found within the
first/last few basepairs using approximate sequence search.

This function must be called on completely initialized `IncompleteSegmentData`
with both ORFs and sequence, since it modifies both the sequence and ORFs.
"""
function strip_false_termini!(
    data::Dict{String, IncompleteSegmentData},
    filter_termini::Bool
)
    good = Set{String}()

    for (id, segment_data) in data
        seq = @unwrap_or segment_data.seq continue
        fiveprime, threeprime = false_termini_length(seq, segment_data.segment)

        # We discard something like 93% of all seqs here
        if fiveprime === nothing || threeprime === nothing
            if filter_termini
                continue
            else
                fiveprime = fiveprime === nothing ? UInt32(0) : fiveprime
                threeprime = threeprime === nothing ? UInt32(0) : threeprime
            end
        end
        isgood = strip_false_termini!(segment_data, fiveprime, threeprime)
        isgood && push!(good, id)
    end

    filter!(data) do (k, _)
        k in good
    end
    return data
end

# Returns whether
function strip_false_termini!(
    data::IncompleteSegmentData,
    fiveprime::UInt32,
    threeprime::UInt32,
)::Bool
    iszero(fiveprime) && iszero(threeprime) && return true
    seqlen = UInt32(length(unwrap(data.seq)))

    # For safety, error if there are no ORFs. This function should not be called
    # on segments without ORFs, since the orfs need to be modified by this funcion
    isempty(data.proteins) && error("Termini must be stripped after adding ORFs to data")

    # Trim sequence and ORFs.
    data.seq = some(unwrap(data.seq)[fiveprime+1:seqlen-threeprime])
    for (proti, protein) in enumerate(data.proteins), (orfi, orf) in enumerate(protein.orfs)
        if threeprime ≥ first(orf) || (seqlen - threeprime) < last(orf)
            return false
        else
            data.proteins[proti].orfs[orfi] = (first(orf) - fiveprime):(last(orf) - fiveprime)
        end
    end
    return true
end

# We allow to strip up to 100 bp in each end off
"Returns (5', 3'), where each is nothing, if no termini was found,
or an UInt32 if one was found."
function false_termini_length(
    seq::LongDNASeq,
    segment::Segment
)::NTuple{2, Union{Nothing, UInt32}}
    trim5, trim3 = UInt32(0), UInt32(0)
    SEARCHLEN = 100

    fw_queries, bw_queries = TERMINI[segment]
    trim3 = trim5 = nothing
    # Find true beginning of sequence, if it's within first SEARCHLEN bp
    for query in fw_queries
        p = BioSequences.approxsearch(seq, query, 2, 1, SEARCHLEN)
        if !isempty(p)
            trim5 = UInt32(first(p) - 1)
            break
        end
    end
    # Find true end of sequence
    for query in bw_queries
        p = BioSequences.approxrsearch(seq, query, 2, lastindex(seq), lastindex(seq)-SEARCHLEN)
        if !isempty(p)
            trim3 = UInt32(lastindex(seq) - last(p))
            break
        end
    end
    return (trim5, trim3)
end