text:
# Tracing-based Automatic Differentiation
# =======================================
#
# Machine learning primarily needs [reverse-mode AD](./backandforth.ipynb), and
# tracing / operator overloading approaches are by far the most popular way to
# it; this is the technique used by ML frameworks from Theano to PyTorch. This
# notebook will cover the techniques used by those frameworks, as well as
# clarifying the distinction between the "static declaration"
# (Theano/TensorFlow) and "eager execution" (Chainer/PyTorch/Flux) approaches to
# AD.

include("utils.jl")

# Partial Evaluation
# ------------------
#
# Say we have a simple implementation of $x^n$ which we want to differentiate.

function pow(x, n)
  r = 1
  for i = 1:n
    r *= x
  end
  return r
end

pow(2, 3)

# We already know how to [differentiate Wengert lists](./intro.ipynb), but this
# doesn't look much like one of those. In fact, we can't write this program as a
# Wengert list at all, given that it contains control flow; and more generally
# our programs might have things like data structures or function calls that we
# don't know how to differentiate either.
#
# Though it's possible to generalise the Wengert list to handle these things,
# there's actually a simple and surprisingly effective alternative, called
# "partial evaluation". This means running some part of a program without
# running all of it. For example, given an expression like $x + 5 * n$ where we
# know $n = 3$, we can simplify to $x + 15$ even though we don't know what $x$
# is. This is a common trick in compilers, and Julia will often do it for you:

f(x, n) = x + 5 * n
g(x) = f(x, 3)

code_typed(g, Tuple{Int})[1]

# This suggests a solution to our dilemma above. If we know what $n$ is (say,
# $3$), we can write `pow(x, 3)` as $((1*x)*x)*x$, which _is_ a Wengert
# expression that we can differentiate. In effect, this is a kind of compilation
# from a complex language (Julia, Python) to a much simpler one.

# Static Declaration
# ------------------
#
# We want to trace all of the basic mathematical operations in the program,
# stripping away everything else. We'll do this using Julia's operator
# overloading; the idea is to create a new type which, rather than actually executing
# operations like $a + b$, records them into a Wengert list.

import Base: +, -

struct Staged
  w::Wengert
  var
end

a::Staged + b::Staged = Staged(w, push!(a.w, :($(a.var) + $(b.var))))

a::Staged - b::Staged = Staged(w, push!(a.w, :($(a.var) - $(b.var))))

# Actually, all of our staged definitions follow the same pattern, so we can
# just do them in a loop. We also add an extra method so that we can multiply
# staged values by constants.

for f in [:+, :*, :-, :^, :/]
  @eval Base.$f(a::Staged, b::Staged) = Staged(a.w, push!(a.w, Expr(:call, $(Expr(:quote, f)), a.var, b.var)))
  @eval Base.$f(a, b::Staged) = Staged(b.w, push!(b.w, Expr(:call, $(Expr(:quote, f)), a, b.var)))
end

# The idea here is to begin by creating a Wengert list (the "graph" in ML
# framework parlance), and create some symbolic variables which do not yet
# have numerical values.

w = Wengert()
x = Staged(w, :x)
y = Staged(w, :y)

# When we manipulate these variables, we'll get Wengert lists.

z = 2x + y
z.w |> Expr

# Crucially, this works with our original `pow` function!

w = Wengert()
x = Staged(w, :x)

y = pow(x, 3)
y.w |> Expr

# The rest is almost too easy! We already know how to derive this.

dy = derive_r(y.w, :x)
Expr(dy)

# If we dump the derived code into a function, we get code for the derivative
# of $x^3$ at any point (i.e. $3x^2$).

@eval dcube(x) = $(Expr(dy))

dcube(5)

# Congratulations, you just implemented TensorFlow.

# Eager Execution
# ---------------

# This approach has a crucial problem; because it works by stripping out control
# flow and parameters like $n$, it effectively freezes all of these things. We
# can get a specific derivative for $x^3$, $x^4$ and so on, but we can't get the
# general derivative of $x^n$ with a single Wengert list. This puts a severe
# limitation on the kinds of models we can express.$^1$
#
# The solution? Well, just re-build the Wengert list from scratch every time!

function D(f, x)
  x_ = Staged(w, :x)
  dy = derive(f(x_).w, :x)
  eval(:(let x = $x; $(Expr(dy)) end))
end

D(x -> pow(x, 3), 5)
#-
D(x -> pow(x, 5), 5)

# This gets us our gradients, but it's not going to be fast – there's a lot of overhead
# to building and evaluating the list/graph every time. There are two things we can
# do to alleviate this:
#
# 1. Interpret, rather compile, the Wengert list.
# 2. Fuse interpretation of the list (the numeric phase) with the building
#    and manipulation of the Wengert list (the symbolic phase).
#
# Implementing this looks a lot like the `Staged` object above. The key difference
# is that alongside the Wengert list, we store the numerical values of each variable
# and instruction as we go along. Also, rather than explicitly naming variables
# `x`, `y` etc, we generate names using `gensym()`.

gensym()
#-
struct Tape
  instructions::Wengert
  values
end

Tape() = Tape(Wengert(), Dict())

struct Tracked
  w::Tape
  var
end

function track(t::Tape, x)
  var = gensym()
  t.values[var] = x
  Tracked(t, var)
end

Base.getindex(x::Tracked) = x.w.values[x.var]

for f in [:+, :*, :-, :^, :/]
  @eval function Base.$f(a::Tracked, b::Tracked)
    var = push!(a.w.instructions, Expr(:call, $(Expr(:quote, f)), a.var, b.var))
    a.w.values[var] = $f(a[], b[])
    Tracked(a.w, var)
  end
  @eval function Base.$f(a, b::Tracked)
    var = push!(b.w.instructions, Expr(:call, $(Expr(:quote, f)), a, b.var))
    b.w.values[var] = $f(a, b[])
    Tracked(b.w, var)
  end
  @eval function Base.$f(a::Tracked, b)
    var = push!(a.w.instructions, Expr(:call, $(Expr(:quote, f)), a.var, b))
    a.w.values[var] = $f(a[], b)
    Tracked(a.w, var)
  end
end

# Now, when we call `pow` it looks a lot more like we are dealing with normal
# numeric values; but there is still a Wengert list inside.

t = Tape()
x = track(t, 5)

y = pow(x, 3)
y[]

y.w.instructions |> Expr

# Finally, we need to alter how we derive this list. The key insight is that
# since we already have values available, we don't need to explicitly build
# and evaluate the derivative code; instead, we can just evaluate each instruction
# numerically as we go along. We more-or-less just need to replace our symbolic
# functions like (`addm`) with the regular ones (`+`).
#
# This is, of course, not a particularly optimised implementation, and faster
# versions have many more tricks up their sleaves. But this gets at all the key
# ideas.

function derive(w::Tape, xs...)
  ds = Dict()
  val(x) = get(w.values, x, x)
  d(x) = get(ds, x, 0)
  d(x, Δ) = ds[x] = d(x) + Δ
  d(lastindex(w.instructions), 1)
  for v in reverse(collect(keys(w.instructions)))
    ex = w.instructions[v]
    Δ = d(v)
    if @capture(ex, a_ + b_)
      d(a, Δ)
      d(b, Δ)
    elseif @capture(ex, a_ * b_)
      d(a, Δ * val(b))
      d(b, Δ * val(a))
    elseif @capture(ex, a_^n_Number)
      d(a, Δ * n * val(a) ^ (n-1))
    elseif @capture(ex, a_ / b_)
      d(a, Δ * val(b))
      d(b, -Δ*val(a)/val(b)^2)
    else
      error("$ex is not differentiable")
    end
  end
  return map(x -> d(x.var), xs)
end

derive(y.w, x)

# With this we can implement a more general gradient function.

function gradient(f, xs...)
  t = Tape()
  xs = map(x -> track(t, x), xs)
  f(xs...)
  derive(t, xs...)
end

# Even with the limited set of gradients that we have, we're well on our way to
# differentiating more complex programs, like a custom `sin` function.

gradient((a, b) -> a*b, 2, 3)
#-
mysin(x) = sum((-1)^k/factorial(1.0+2k) * x^(1+2k) for k = 0:5)
#-
gradient(mysin, 0.5)
#-
cos(0.5)

# We can even take nested derivatives!

gradient(x -> gradient(mysin, x)[1], 0.5)
#-
-sin(0.5)

# Though the tracing approach has significant limitations, its power is in how
# easy it is to implement: one can build a fairly full-featured implementation,
# in almost any language, in a weekend. Almost all languages have the
# operator-overloading features required, and no matter how complex the host
# language, one ends up with a simple Wengert list.

# Note that we have not removed the need to apply our basic symbolic
# differentiation algorithm here. We are still looking up gradient definitions,
# reversing data flow and applying the chain rule – it's just interleaved with
# our numerical operations, and we avoid building the output into an explicit
# Wengert list.
#
# It's somewhat unusual to emphasise the symbolic side of AD, but I think it
# gives us an incisive way to understand the tradeoffs that different systems
# make. For example: TensorFlow-style AD has its numeric phase separate from
# Python's runtime, which makes it awkward to use. Conversely, PyTorch does run
# its numerical phase at runtime, but also its symbolic phase, making it
# impossible to optimise the backwards pass.
#
# We [observed](./forward.ipynb) that OO-based forward mode is particularly
# successful because it carries out its symbolic and numeric operations at
# Julia's compile and run time, respectively. In the [source to source reverse
# mode](./reverse.ipynb) notebook, we'll explore doing this for reverse mode as
# well.

# ### Footnotes

# $^1$ Systems like TensorFlow can also just provide ways to inject control flow
# into the graph. This brings us closer to a [source-to-source
# approach](./reverse.ipynb) where Python is used to build an expression in
# TensorFlows internal graph language.

# Fun fact: PyTorch and Flux's tapes are actually closer to the `Expr` format
# that we originally used, in which "tracked" tensors just have pointers to
# their parents (implicitly forming a graph/Wengert list/expression tree). A
# naive algorithm for backpropagation suffers from exponential runtime for the
# *exact* same reason that naive symbolic diff does; "flattening" this graph
# into a tree causes it to blow up in size.


text:
# olsvar.m
# Luca Brugnolini 2014
# Based on Kilian code

function olsvar(y,p)

# Input : data y, lag order p
# Output: coefficient estimate B, estimate of Cov(U), residual U

(T,K) = size(y);
Y = (y[p+1:T,:])';                     # LHS variable
X = ones(1,T-p);                       # constant term

for j=1:p                              # RHS variable
    Z = [y[p+1-j:T-j,:]]';
    X = [X;Z];

end

B = (Y*X')/(X*X');                     # OLS estimator
U = Y-B*X;                             # estimated residuals
SIGMA = U*U'/(T-p-p*K-1);              # estimated variance of disturbances

return B, SIGMA, U

end

ext: jl


text:
immutable MaxPooling{T,N} <: Transformation
    input::SumNode{T,N}
    output::OutputNode{T,N}
    stride::NTuple{N,Int}
end
MaxPooling(args...) = MaxPooling(Float64, args...)

function MaxPooling{T,N}(::Type{T}, sizein::NTuple{N,Int}, stride::NTuple{N,Int} = ntuple(i->2,N))
    sizeout = map(s->div(s,2), sizein)
    @show sizein, sizeout
    input = InputNode(T,sizein)
    output = OutputNode(T,sizeout)
    MaxPooling(input, output, stride)
end

function transform!{T}(pool::MaxPooling{T,2})
    nr,nc = pool.sizeout
    rstride,cstride = pool.stride
    for r=1:nr, c=1:nc
        # one before the start
        rstart = (r-1)*rstride
        cstart = (c-1)*cstride

        tile = view(pool.input.val, rstart+1:rstart+rstride, cstart+1:cstart+cstride)
        pool.output.val[r,c] = max(tile)
    end
    value(output)
end

function grad!{T}(pool::MaxPooling{T,2})
    # TODO: pass the full gradient to the input with the maximum value in the tile
    #   note: we should really store this index during the forward pass... maybe have a matrix of indices??
end


text:
export @device, @ctrl, @measure

"""
    @device [options] <generic circuit definition>

Entry for defining a generic quantum program. A generic quantum program is a function takes
a set of classical arguments as input and return a quantum program that can be furthur compiled
into pulses or other quantum instructions. The quantum program can return classical values from
device if `return` statement is declared explicitly, or it always return nothing, and mutates the
quantum register.

# Supported Semantics

- [`@ctrl`](@ref): Keyword for controlled gates in quantum circuit.
- [`@measure`](@ref): Keyword for measurement in quantum circuit.

The function marked by `@device` can be multiple dispatched like other Julia function. The only difference
is that it always returns a quantum circuit object that should be runable on quantum device by feeding it
the location of qubits and the pointer to quantum register.

# Options

- `target`, compilation target, default is `:julia`, see **Compilation Targets** for details.

## Compilation Targets

- `:julia`, default target, compiles the program to Julia program.
- `:qasm`, compiles the program to [openQASM](https://github.com/Qiskit/openqasm).

# Example

We can define a Quantum Fourier Transformation in the following recursive way

```julia
@device function qft(n::Int)
    1 => H
    for k in 2:n
        @ctrl k 1=>shift(2π/2^k)
    end

    if n > 1
        2:n => qft(n-1)
    end
end
```

This will give us a generic quantum circuit `qft` with 1 method.
"""
macro device end

macro device(ex)
    return esc(device_m(__module__, ex))
end

macro device(args...)
    options = args[1:end-1]
    ex = args[end]

    kwargs = Pair[]
    for each in options
        if (each isa Expr) && (each.head === :(=))
            each.args[2] isa QuoteNode || throw(ParseError("expect a Symbol, got $(each.args[2])"))
            push!(kwargs, each.args[1] => each.args[2].value)
        else
            throw(Meta.ParseError("Invalid Syntax, expect a compile option, got $each"))
        end
    end

    return esc(device_m(__module__, ex; kwargs...))
end

"""
    @ctrl k <gate location>

Keyword for controlled gates in quantum circuit. It must be used inside `@device`. See also [`@device`](@ref).
"""
macro ctrl end

"""
    @measure <location> [operator] [configuration]

Keyword for measurement in quantum circuit. It must be used inside `@device`. See also [`@device`](@ref).

# Arguments

- `<location>`: a valid `Locations` argument to specifiy where to measure the register
- `[operator]`: Optional, specifiy which operator to measure
- `[configuration]`: Optional, it can be either:
    - `remove=true` will remove the measured qubits
    - `reset_to=<bitstring>` will reset the measured qubits to given bitstring
"""
macro measure end


text:
using ValidatedNumerics
using .DynamicDefinition, .PwDynamicDefinition
using .Contractors

using .DynamicDefinition: derivative, plottable

using TaylorSeries: Taylor1

struct Iterate <: Dynamic
    D::PwMap
    n::Int
end

Base.show(io::IO, D::Iterate) = print(io, "$(D.n)-times iterate of: ", D.D)

function (D::Iterate)(x::Taylor1)
    y = x
    for i = 1:D.n
        y = (D.D)(y)
    end
    return y
end

DynamicDefinition.nbranches(D::Iterate) = nbranches(D.D)^D.n
DynamicDefinition.is_full_branch(D::Iterate) = is_full_branch(D.D)
DynamicDefinition.domain(D::Iterate) = domain(D.D)

"""
compose_endpoints(D, endpointsf)

Find discontinuity points of f ∘ D, where f has discontinuity points endpointsf
and D is a given PwMap.
"endpoints" always include the extrema of the given domain.
"""
function compose_endpoints(D, endpointsf)
    v = [D.endpoints[1]]
    for k = 1:nbranches(D)
        preims = [preim(D, k, x) for x in endpointsf[2:end-1]]
        if orientation(D, k) > 0
            # TODO: still trouble if the branches are not full
            append!(v, preims)
        else
            append!(v, preims[end:-1:1])
        end
        append!(v, [D.endpoints[k+1]])
    end
    return v
end

function DynamicDefinition.endpoints(D::Iterate)
    @assert D.n >= 1
    endpoints = D.D.endpoints
    for k = 2:D.n
        endpoints = compose_endpoints(D.D, endpoints)
    end
    return endpoints
end

"""
unpack(k, b, D)

Convert an integer k∈[1,b^n] into a tuple ∈[1,b]^n bijectively

This is used to index preimages: the k'th of the b^n preimages of an Iterate
corresponds to choosing the v[i]'th branch when choosing the i'th preimage, for i = 1..k,
where v = unpack(k, b, n)
"""
function unpack(k, b, n)
    @assert 1 ≤ k ≤ b^n
    v = fill(0, n)
    k = k-1
    for i = 1:n
        (k, v[n+1-i]) = divrem(k, b)
    end
    return v .+ 1
end

function evaluate_branch(D::Iterate, k, x)
    @assert 1 ≤ k ≤ nbranches(D)
    n = D.n
    v = unpack(k, nbranches(D.D), n)
    for i = 1:n
        x = branch(D.D, v[i])(x)
    end
    return x
end

DynamicDefinition.branch(D::Iterate, k) = x -> evaluate_branch(D, k, x)

using LinearAlgebra: Bidiagonal

function Jac(fprime, v::Vector{T}) where {T}
    dv = fprime.(v)
    ev = -ones(T, length(v)-1)
    return Bidiagonal{T}(dv, ev, :U)
end

"""
Compute the preimage of an Iterate D in its k'th branch
"""
function DynamicDefinition.preim(D::Iterate, k, y, ϵ=1e-15; max_iter = 100)
    @assert 1 <= k <= nbranches(D)

    n = D.n

    v = unpack(k, nbranches(D.D), n)

    fs = D.D.Ts[v]
    S = [hull(D.D.endpoints[v[i]], D.D.endpoints[v[i]+1]) for i in 1:n]
    return nthpreimage!(y, fs, S)[1]
end

function DynamicDefinition.orientation(D::Iterate, k)
    v = unpack(k, nbranches(D.D), D.n)
    return prod([orientation(D.D, k) for k in v])
end

function DynamicDefinition.plottable(D::Iterate, x)
	@assert 0 <= x <= 1
    for k = 1:D.n
        x = DynamicDefinition.plottable(D.D, x)
        x = clamp(x, 0, 1)
    end
    return x
end

using RecipesBase
@recipe f(::Type{Iterate}, D::Iterate) = x -> plottable(D, x)


text:
n=10
g10 = CompleteGraph(n)
z = copy(g10)
for k=2:5
    z = blkdiag(z, g10)
    add_edge!(z, (k-1)*n, k*n)
    c, ch = label_propagation(z)
    a = collect(n:n:k*n)
    a = Int[div(i-1,n)+1 for i=1:k*n]
    # check the number of community
    @test length(unique(a)) == length(unique(c))
    # check the partition
    @test a == c
end


text:
using aoc2019.Fuel: total, correct_total
using DelimitedFiles

input = joinpath(@__DIR__, "input")
weights = readdlm(input, Int)

p1 = total(weights)
@assert p1 == 3420719

p2 = correct_total(weights)
@assert p2 == 5128195

print("-----------------------------------------------------------------------\n")
print("the tyranny of the rocket equation -- part one\n    fuel required : $p1\n")
print("the tyranny of the rocket equation -- part two\n    fuel required : $p2\n")
print("-----------------------------------------------------------------------\n")


text:
# ---------------------------------------------------------------------------- #
#
#   Elas.jl
#
#   Struct for elasticity problems
#
#   λυτέος
#   Fall 2017
#
#   Max Opgenoord
#
# ---------------------------------------------------------------------------- #

include("material.jl")

"""
    Elas

Elas type:
Used as placeholder for information about structural elasticity problems,
e.g name, source function, and boundary conditions.
"""
struct Elas{S<:Function,B<:Function} <: Problem

  name::String     # Name of problem

  mat::Material    # Material used

  source::S               # Rhs of problem

  bctype::Vector{Int64}   # Relates each boundary to a boundary function type
  bcnorm::Bool            # Whether or not boundary functions are defined normal to the surface
  bcfunc::Vector{B} # Boundary function for each boundary face

end

"""
    Elas(src::Function, bcf::Array{Function}, bctype::Array{Int64},
         mat::Material; name = "Answer to Life Universe and Everything",
         bcnorm = false)

Constructor for `Problem` type. `src` is a function that returns the source
function for the problem. `bcf` is an Array of Functions, which specifies the
boundary condition on each boundary surface.
"""
function Elas(src::Function, bcf::Vector{Function}, bctype::Vector{Int64},
              mat::Material; name = "Answer to Life Universe and Everything",
              bcnorm = false)

  if length(bctype) != length(bcf)
    error("Problem:: bctype needs have the same size as bcf")
  end

  Elas(name, src, bctype, bcnorm, bcf)

end


text:
using Blink
w = Blink.AtomShell.Window()
loadurl(w, string("file://" ,pwd(),"/index.html"))

max_stars_repo_name: wolfiex/RopaJL

max_issues_repo_name: wolfiex/RopaJL

max_forks_repo_name: wolfiex/RopaJL


text:
using Cairo
using Random

# shape functions,

function randpos(n,w::Real,h::Real)
    @static if VERSION >= v"0.7"
        Random.seed!(141413)
    else
        srand(141413)
    end
    px = rand(n)*w
    py = rand(n)*h
    return (px,py)
end

function clear_bg(c::CairoContext,w::Real,h::Real)
    save(c)
    set_source_rgb(c,1.0,1.0,1.0)
    rectangle(c, 0,0,w,h)
    paint(c)
    restore(c)
end

# ddotsx: fill random dots with different methods, disc only (no ring)
function ddots1(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    new_path(cr)
    for i=1:n
       move_to(cr,px[i],py[i])
       rel_move_to(cr,radius,0)
       arc(cr, px[i], py[i], radius, 0, 2*pi)
    end
    close_path(cr)

    set_source_rgb(cr, 0, 0, 1.0)
    fill(cr)
end

function ddots2(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    for i=1:n
        new_path(cr)
        move_to(cr,px[i],py[i])
        rel_move_to(cr,radius,0)
        arc(cr, px[i], py[i], radius, 0, 2*pi)
        close_path(cr)
        set_source_rgb(cr, 0, 0, 1.0)
        fill(cr)
    end
end

function ddots3(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    new_path(cr)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_cap(cr,Cairo.CAIRO_LINE_CAP_ROUND)
    set_line_width(cr,radius*2.0)

    for i=1:n
       move_to(cr,px[i],py[i])
       rel_line_to(cr,0,0)
       stroke(cr)
    end

    close_path(cr)

end

function ddots4(cr::CairoContext,
    rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    cc = radius + 1

    rectangle(cr,0,0,2*cc,2*cc)
    clip(cr)
    push_group(cr)
    arc(cr,radius,radius,radius,0,2*pi)
    set_source_rgb(cr,0,0,1.0)
    fill(cr)
    p = pop_group(cr)
    reset_clip(cr)

    for i=1:n
        save(cr)
        translate(cr,px[i]-cc,py[i]-cc)
        set_source(cr,p)
        paint(cr)
        restore(cr)
    end
end

function ddots5(cr::CairoContext,
    rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    cc = radius + 1

    s1 = Cairo.CairoARGBSurface(cc*2,cc*2)
    c1 = Cairo.CairoContext(s1)
    rectangle(c1,0,0,2*cc,2*cc)
    set_source_rgba(c1,0,0,0,0)
    paint(c1)
    new_path(c1)
    arc(c1,radius,radius,radius,0,2*pi)
    set_source_rgb(c1,0,0,1.0)
    fill(c1)

    p = Cairo.CairoPattern(s1)

    for i=1:n
        save(cr)
        translate(cr,px[i]-cc,py[i]-cc)
        set_source(cr,p)
        paint(cr)
        restore(cr)
    end
end

# rdotsx: fill random dots with different methods
function rdots1(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    new_path(cr)
    for i=1:n
       move_to(cr,px[i],py[i])
       rel_move_to(cr,radius,0)
       arc(cr, px[i], py[i], radius, 0, 2*pi)
    end
    close_path(cr)

    set_source_rgb(cr, 0, 0, 1.0)
    fill_preserve(cr)
    set_line_width(cr,1.0)
    set_source_rgb(cr, 0, 0, 0)
    stroke(cr)
end

function rdots2(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    for i=1:n
        new_path(cr)
        move_to(cr,px[i],py[i])
        rel_move_to(cr,radius,0)
        arc(cr, px[i], py[i], radius, 0, 2*pi)
        close_path(cr)
        set_source_rgb(cr, 0, 0, 1.0)
        fill_preserve(cr)
        set_line_width(cr,1.0)
        set_source_rgb(cr, 0, 0, 0)
        stroke(cr)
    end
end

function rdots3(cr::CairoContext, rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    new_path(cr)

    set_line_cap(cr,Cairo.CAIRO_LINE_CAP_ROUND)

    set_source_rgb(cr, 0, 0, 0)
    set_line_width(cr,radius*2 + 2.0)
    for i=1:n
       move_to(cr,px[i],py[i])
       rel_line_to(cr,0,0)
       stroke(cr)
    end

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,radius*2)
    for i=1:n
       move_to(cr,px[i],py[i])
       rel_line_to(cr,0,0)
       stroke(cr)
    end

end

function rdots4(cr::CairoContext,
    rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    cc = radius + 3.0

    rectangle(cr,0,0,2*cc,2*cc)
    clip(cr)
    push_group(cr)
    arc(cr,radius,radius,radius,0,2*pi)
    set_source_rgb(cr,0,0,1.0)
    fill_preserve(cr)
    set_source_rgb(cr,0,0,0)
    set_line_width(cr,1.0)
    stroke(cr)
    p = pop_group(cr)
    reset_clip(cr)

    for i=1:n
        save(cr)
        translate(cr,px[i]-cc,py[i]-cc)
        set_source(cr,p)
        paint(cr)
        restore(cr)
    end
end

function rdots5(cr::CairoContext,
    rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    cc = radius + 3.0

    s1 = Cairo.CairoARGBSurface(cc*2,cc*2)
    c1 = Cairo.CairoContext(s1)
    rectangle(c1,0,0,2*cc,2*cc)
    set_source_rgba(c1,0,0,0,0)
    paint(c1)
    new_path(c1)
    arc(c1,radius+1,radius+1,radius,0,2*pi)
    set_source_rgb(c1,0,0,1.0)
    fill_preserve(c1)
    set_source_rgb(c1,0,0,0)
    set_line_width(c1,1.0)
    stroke(c1)

    p = Cairo.CairoPattern(s1)

    for i=1:n
        save(cr)
        translate(cr,px[i]-cc,py[i]-cc)
        set_source(cr,p)
        paint(cr)
        restore(cr)
    end
end

# alpha transparency using paint_with_alpha
function rdots6(cr::CairoContext,
    rect_width::Real, rect_height::Real, radius::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    cc = radius + 3.0

    s1 = Cairo.CairoARGBSurface(cc*2,cc*2)
    c1 = Cairo.CairoContext(s1)
    rectangle(c1,0,0,2*cc,2*cc)
    set_source_rgba(c1,0,0,0,0)
    paint_with_alpha(c1, 0.5)
    new_path(c1)
    arc(c1,radius+1,radius+1,radius,0,2*pi)
    set_source_rgb(c1,0,0,1.0)
    fill_preserve(c1)
    set_source_rgb(c1,0,0,0)
    set_line_width(c1,1.0)
    stroke(c1)

    p = Cairo.CairoPattern(s1)

    for i=1:n
        save(cr)
        translate(cr,px[i]-cc,py[i]-cc)
        set_source(cr,p)
        paint_with_alpha(cr, 1 - i/n)
        restore(cr)
    end
end

# lines0, random x,y lines
function lines0(cr::CairoContext, rect_width::Real, rect_height::Real, width::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    new_path(cr)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,width)

    new_path(cr)
    move_to(cr,px[1],py[1])
    for i=2:n
       line_to(cr,px[i],py[i])
    end
    stroke(cr)
end

# lines1, x sorted, line a plot
function lines1(cr::CairoContext, rect_width::Real, rect_height::Real, width::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    px = sort(px)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,width)

    new_path(cr)
    move_to(cr,px[1],py[1])
    for i=2:n
        line_to(cr,px[i],py[i])
    end
    stroke(cr)
end

# lines2, x sorted, independent lines per coord
function lines2(cr::CairoContext, rect_width::Real, rect_height::Real, width::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    px = sort(px)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,width)

    new_path(cr)
    for i=1:n-1
        move_to(cr,px[i],py[i])
        line_to(cr,px[i+1],py[i+1])
        stroke(cr)
    end
end

# lines3, x sorted, in clusters of 100
function lines3(cr::CairoContext, rect_width::Real, rect_height::Real, width::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    px = sort(px)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,width)
    new_path(cr)
    for c = 1:100:n
        move_to(cr,px[c],py[c])
        for i=c:min(c+100,n)
            line_to(cr,px[i],py[i])
        end
        stroke(cr)
    end
end

# lines4, x sorted, in clusters of 1000
function lines4(cr::CairoContext, rect_width::Real, rect_height::Real, width::Real, n::Int)

    clear_bg(cr,rect_width,rect_height)
    px,py = randpos(n,rect_width,rect_height)

    px = sort(px)

    set_source_rgb(cr, 0, 0, 1.0)
    set_line_width(cr,width)
    new_path(cr)
    for c = 1:1000:n
        move_to(cr,px[c],py[c])
        for i=c:min(c+1000,n)
            line_to(cr,px[i],py[i])
        end
        stroke(cr)
    end
end


text:
using ITensors
using ITensors.NDTensors
using Test

@testset "NDTensors compatibility" begin
  i = Index([QN(0) => 1, QN(1) => 1])

  T = BlockSparseTensor(Float64, [Block(1, 1)], (i', dag(i)))
  @test nnzblocks(T) == 1
  @test nzblocks(T) == [Block(1, 1)]

  T = BlockSparseTensor(Float64, [Block(1, 1)], [i', dag(i)])
  @test nnzblocks(T) == 1
  @test nzblocks(T) == [Block(1, 1)]

  T = BlockSparseTensor(Float64, [Block(1, 1)], IndexSet(i', dag(i)))
  @test nnzblocks(T) == 1
  @test nzblocks(T) == [Block(1, 1)]

  @testset "blockdim" begin
    i = Index(2)
    @test_throws ErrorException blockdim(i, Block(1))
    @test_throws ErrorException blockdim(i, 1)
    @test_throws ErrorException blockdim(1, Block(1))
    @test_throws ErrorException blockdim(1, 1)
  end
end


text:
abstract type TraitParam end

abstract type TraitParamCons{A, B<:TraitParam} <: TraitParam
end

abstract type TraitParamNil <: TraitParam
end

abstract type Trait{P <: TraitParam} end

abstract type TraitInstance end

struct MethodInterface
    name     :: Symbol
    argsty   :: Any
    retty    :: Any
    fresh    :: Vector{Any}
    infer    :: Set{Symbol} # including type class parameters inferred in callsites.
end

struct DefaultMethod
    name :: Symbol
    impl :: Expr
end
ClassMethod = Union{LineNumberNode, MethodInterface, DefaultMethod}

function extract_default(impl::LineNumberNode) end
function extract_default(@nospecialize(impl))
    function extract_default_inner(call)
        @when (:($f($(_...))) || f) && if f isa Symbol end = call begin
            f
        @when :($f :: $_) = call
            extract_default_inner(f)
        @when :($f where {$(_...)}) = call
            extract_default_inner(f)
        @otherwise
            nothing
        end
    end
    name = @match impl begin
        Expr(:function, call, _)    => extract_default_inner(call)
        Expr(:(=), call, _)         => extract_default_inner(call)
        Expr(:macrocall, _..., arg) => extract_default_inner(arg)
        _ => nothing
    end
    name === nothing && return nothing
    (name, impl)
end

function collect_default!(impl::LineNumberNode, syms :: Set{Symbol}) end
function collect_default!(@nospecialize(impl), syms :: Set{Symbol})
    found_func = extract_default(impl)
    found_func === nothing && error("Not a method definition: $impl.")
    push!(syms, found_func[1])
    nothing
end

function extract_impl(impl::LineNumberNode) end
function extract_impl(@nospecialize(impl))
    function extract_impl_inner!(call)
        @when (:($f($(_...)))) && if f isa Symbol end = call begin
            # fix mutual references when implementing interfaces.
            call.args[1] = gensym(f)
            f
        @when f::Symbol = call
            f
        @when :($f :: $_) = call
            extract_impl_inner!(f)
        @when :($f where {$(_...)}) = call
            extract_impl_inner!(f)
        @otherwise
            nothing
        end
    end
    name = @match impl begin
        Expr(:function, call, _)    => extract_impl_inner!(call)
        Expr(:(=), call, _)         => extract_impl_inner!(call)
        Expr(:macrocall, _..., arg) => extract_impl_inner!(arg)
        _ => nothing
    end
    name === nothing && return nothing
    (name, impl)
end

extract_method_interface(f::Symbol, argtys::AbstractArray, retty::Any, fresh::AbstractArray) =
    begin
    argsty = :($Tuple{$(argtys...)})
    infer  = occurred_sym(:($argsty where {$(fresh...)}))
    MethodInterface(f, argsty, retty, fresh, infer)
    end

extract_method(stmt::Union{Expr, LineNumberNode})::ClassMethod =
    @match stmt begin
        :($f :: [$(args...)] where {$(fresh...)} => $ret) => extract_method_interface(f, args, ret, fresh)
        :($f :: [$(args...)] => $ret)                     => extract_method_interface(f, args, ret, Any[])
        :($f :: $arg where {$(fresh...)} => $ret)         => extract_method_interface(f, [arg], ret, fresh)
        :($f :: $arg => $ret)                             => extract_method_interface(f, [arg], ret, Any[])
        ::LineNumberNode              => stmt
        _                             => begin
            def = extract_default(stmt)
            def === nothing && error("Malformed method interface or default method $stmt.")
            DefaultMethod(def[1], def[2])
        end
    end

extract_tvars(var :: Union{Symbol, Expr})::Symbol =
    @match var begin
        :($a <: $_) => a
        :($a >: $_) => a
        :($_ >: $a >: $_) => a
        :($_ <: $a <: $_) => a
        a::Symbol         => a
    end

struct TypeParams
    as_arguments  :: Vector{Expr}
    as_where      :: Vector{Symbol}
end

TypeParams(syms :: Syms) where Syms <: AbstractArray{Symbol} =
    TypeParams(
       [:(::Type{$sym}) for sym in syms],
       syms
    )

function occurred_sym(sig)::Set{Symbol}
    @match sig begin
        Expr(hd)              => Set(Symbol[])
        :[$(args...)]         => union!(map(occurred_sym, args)...)
        :($a where {$(b...)}) => setdiff!(occurred_sym(a), map(extract_tvars, b))
        Expr(_, tl...)        => union!(map(occurred_sym, tl)...)
        a :: Symbol           => Set([a])
        _                     => Set(Symbol[])
    end
end

function get_functional_dependency(expr)
    error("Malformed functional dependency $expr, example expected form is '{T2 = f(T1), T3 = g(T1)}'")
end


text:
@testset "717.1-bit-and-2-bit-characters.jl" begin
    @test is_one_bit_character([1, 1, 1, 1, 0]) == true
    @test is_one_bit_character([1, 0]) == false
    @test is_one_bit_character([1, 1, 1, 0]) == false
    @test is_one_bit_character([1, 0, 1, 1, 0]) == true
end


text:
export Half

struct Half{M} <: AbstractMeasure
    parent::M
end

function Base.show(io::IO, μ::Half)
    print(io, "Half")
    show(io, μ.parent)
end

unhalf(μ::Half) = μ.parent

function basemeasure(μ::Half)
    inbounds(x) = x > 0
    constℓ = logtwo
    varℓ() = 0.0
    base = basemeasure(unhalf(μ))
    FactoredBase(inbounds, constℓ, varℓ, base)
end

function Base.rand(rng::AbstractRNG, T::Type, μ::Half)
    return abs(rand(rng, T, unhalf(μ)))
end

logdensity(μ::Half, x) = logdensity(unhalf(μ), x)


text:
module StateSpaceModelsEM

# Packages
using Printf
using LinearAlgebra
using Statistics
using Random

include("utils/utils.jl")

include("ACQR/ACQR.jl")
export
	ACQR_kfilter, ACQR_kfilter_s,
	ACQR_ksmoother, ACQR_ksmoother_s,
	ACQR_em, ACQR_em_s,
	ACQR_em_starting

include("ACQRS/ACQRS.jl")
export ACQRS_sub

include("ABCDQR/ABCDQR.jl")
export
	ABCDQR_kfilter, ABCDQR_kfilter_s,
	ABCDQR_em, ABCDQR_em_s,
	ABCDQR_em_starting

include("acqr/acqr.jl")
export
	acqr,acqrEM,acqrKal,
	acqr_kfilter, acqr_kfilter_s,
	acqr_ksmoother, acqr_ksmoother_s,
	acqr_em, acqr_em_s,
	acqr_kalman

end # module


text:
#
# Löwner-John ellipsoid example
#
# Compute the inner and outer ellipsoids of a polytope.
#
# Adapter from example in the MOSEK distro.
#
# The inner ellipsoidal approximation to a polytope
#
#    S = { x \in R^n | Ax < b }.
#
# maximizes the volume of the inscribed ellipsoid,
#
#    { x | x = C*u + d, || u ||_2 <= 1 }.
#
# The volume is proportional to det(C)^(1/n), so the
# problem can be solved as
#
#   maximize         t
#   subject to       t       <= det(C)^(1/n)
#               || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
#                 C is PSD
#
# which is equivalent to a mixed conic quadratic and semidefinite
# programming problem.
#
#
# The outer ellipsoidal approximation to a polytope given
# as the convex hull of a set of points
#
#     S = conv{ x1, x2, ... , xm }
#
# minimizes the volume of the enclosing ellipsoid,
#
#   { x | || P*(x-c) ||_2 <= 1 }
#
# The volume is proportional to det(P)^{-1/n}, so the problem can
# be solved as
#
#   minimize         t
#   subject to       t       >= det(P)^(-1/n)
#               || P*xi + c ||_2 <= 1,  i=1,...,m
#                 P is PSD.
#
# References:
# [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
#

using JuMP

# Models the convex set
#
#   S = { (x, t) \in R^n x R | x >= 0, t <= (x1 * x2 * ... * xn)^(1/n) }
#
function geometric_mean(m, x, t)
    n = length(x)
    if n == 1
        @constraint(m,t-x[1] <= 0)
    else
        t2 = @variable(m)
        @constraint(m,1.0*[t2; x[n]; t] in MOI.PowerCone(1-1.0/n))
        geometric_mean(m,x[1:n-1], t2)
    end
end

# Purpose: Models the hypograph of the n-th power of the
# determinant of a positive definite matrix. See [1,2] for more details.
#
#   The convex set (a hypograph)
#
#   C = { (X, t) \in S^n_+ x R |  t <= det(X)^{1/n} },
#
#   can be modeled as the intersection of a semidefinite cone
#
#   [ X, Z; Z^T Diag(Z) ] >= 0
#
#   and a number of rotated quadratic cones and affine hyperplanes,
#
#   t <= (Z11*Z22*...*Znn)^{1/n}  (see geometric_mean).
function det_root(m,X,t)
    n = floor(Int,sqrt(length(X)))

    @variable(m,Y[1:2*n,1:2*n],PSD)
    @constraint(m,reshape(Y[n+1:2*n,1:n] - Y[n+1:2*n,n+1:2*n],n*n) in MOI.Zeros(n*n))
    @constraint(m,reshape(X - Y[1:n,1:n], n*n) in MOI.Zeros(n*n))

    geometric_mean(m,[Y[i+n,i+n] for i in 1:n],t)
end

#  The inner ellipsoidal approximation to a polytope
#
#     S = { x \in R^n | Ax < b }.
#
#  maximizes the volume of the inscribed ellipsoid,
#
#     { x | x = C*u + d, || u ||_2 <= 1 }.
#
#  The volume is proportional to det(C)^(1/n), so the
#  problem can be solved as
#
#    maximize         t
#    subject to       t       <= det(C)^(1/n)
#                || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
#                C is PSD
#
#  which is equivalent to a mixed conic quadratic and semidefinite
#  programming problem.
function lownerjohn_inner(solver :: MOI.AbstractOptimizer,
                          A      :: Array{Float64,2},
                          b      :: Vector{Float64})
    m,n = size(A)

    model = direct_model(solver)

    @variable(model, t >= 0)
    @variable(model, C[1:n,1:n], PSD)
    @variable(model, d[1:n])

    # [ b-A*d ; A*C ] in QCone

    for i in 1:m
        @constraint(model, [ b[i] - (A[i,:]' * d) ; (A[i,:]'*C)'] in MOI.SecondOrderCone(n+1))
    end

    det_root(model,C,t)

    @objective(model, Max, t)

    optimize(model)

    reshape([JuMP.resultvalue(item) for item in C],(n,n)),[JuMP.resultvalue(item) for item in d]
end

# The outer ellipsoidal approximation to a polytope given
# as the convex hull of a set of points
#
#   S = conv{ x1, x2, ... , xm }
#
# minimizes the volume of the enclosing ellipsoid,
#
#   { x | || P*x-c ||_2 <= 1 }
#
# The volume is proportional to det(P)^{-1/n}, so the problem can
# be solved as
#
#   maximize         t
#   subject to       t       <= det(P)^(1/n)
#               || P*xi - c ||_2 <= 1,  i=1,...,m
#               P is PSD.
function lownerjohn_outer(solver :: MOI.AbstractOptimizer,
                          x      :: Array{Float64,2})
    m,n = size(x)
    model = direct_model(solver)

    @variable(model,t >= 0)
    @variable(model,P[1:n,1:n],PSD)
    @variable(model,c[1:n])

    # (1, Px-c) in cone
    for i in 1:m
        @constraint(model,[1 ; (x[i,:]' * P)' - c] in MOI.SecondOrderCone(n+1))
    end

    # t <= det(P)^{1/n}
    det_root(model,P,t)

    @objective(model,Max,t)

    optimize(model)

    [ JuMP.resultvalue(item) for item in P ], [ JuMP.resultvalue(item) for item in c ]
end

# Probably do something else here:
using MosekTools
solver = Mosek.Optimizer

# Test

# Vertices of a pentagon in 2D
p = [0.0  0.0;
     1.0  3.0;
     5.5  4.5;
     7.0  4.0;
     7.0  1.0;
     3.0 -2.0]

n = size(p,1)
# The hyperplane representation of the same polytope
A = [ (-p[:,2] + [p[n,2];p[1:n-1,2]])  (p[:,1] - [p[n,1];p[1:n-1,1]]) ]
b = A[:,1] .* p[:,1] + A[:,2] .* p[:,2]

Po, co = lownerjohn_outer(solver(), p)
Ci, di = lownerjohn_inner(solver(), A, b)
Poinv = Po ^ -1

println("Solution:")
println("  Outer ellipsoid: Po = $Po, co = $co")
println("  Inner ellipsoid: Ci = $Ci, di = $di")

if Pkg.installed("Plots") != nothing && Pkg.installed("GR") != nothing
    println("Plotting solution...")
    import Plots
    Plots.gr()
    Plots.plot([p[:,1] ; p[1,1]], [p[:,2] ; p[1,2]])

    Plots.plot!(t -> Ci[1,:]' * [ cos(t),sin(t) ] + di[1],
                t -> Ci[2,:]' * [ cos(t),sin(t) ] + di[2],
                0,2*pi)

    Plots.plot!(t -> Poinv[1,:]' * ( [cos(t),sin(t)] + c ),
                t -> Poinv[2,:]' * ( [cos(t),sin(t)] + c ),
                0,2*pi)

    println("Plot done. Press enter to continue.")
    readline()
else
    println("Install Plots package to plot solution")
end


text:
@testset "Recruitment tests" begin

piseed = floor( Int, pi * 1_000_000 )
eseed = floor( Int, exp( 1 ) * 1_000_000 )
s2seed = floor( Int, sqrt( 2 ) * 1_000_000 )

mpSim = ManpowerSimulation( "sim" )
setSimulationLength!( mpSim, 300 )

attribute = Attribute( "attribute A" )
addInitialAttributeValue!( attribute, "value A", 1 )
addSimulationAttribute!( mpSim, attribute )

node = BaseNode( "node A" )
setNodeTarget!( node, 50 )
addNodeRequirement!( node, "attribute A", "value A" )
addSimulationBaseNode!( mpSim, node )

@testset "Basic recruitment test" begin
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig=false )
    report = nodeFluxReport( mpSim, 12, :in, "active" )["active"]
    @test all( report[:, Symbol( "external => active" )] .== 10 ) &&
        all( report[:, Symbol( "recruitment A: external => node A" )] .== 10 )
end  # @testset "Basic recruitment test"

@testset "Recruitment schedule" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig=false )
    report = nodeFluxReport( mpSim, 12, :in, "active" )["active"]
    @test all( report[:, Symbol( "external => active" )] .==
        repeat( [10, 0], 13 ) )

    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig=false )
    report = nodeFluxReport( mpSim, 12, :in, "active" )["active"]
    @test all( report[:, Symbol( "external => active" )] .==
        repeat( [0, 10], 13 ) )

end  # @testset "Recruitment schedule"

@testset "Adaptive recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 12 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentAdaptiveRange!( recruitment, 5, 10 )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig=false )
    report = nodeFluxReport( mpSim, 12, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    targetRec = vcat( fill( 10, 5 ), fill( 5, 21 ) )
    @test all( report .== targetRec )
end  # @testset "Adaptive recruitment test"

@testset "Disc random recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 1 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentDist!( recruitment, :disc, Dict( 5 => 1.0, 10 => 2.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )[[1, end]]
    pval = pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) )
    @test pval > 0.05
    @test all( reportCounts .== [93, 208] )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )[[1, end]]
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) == pval
    @test all( reportCounts .== [93, 208] )

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )[[1, end]]

    if reportCounts == [93, 208]
        @warn "Statistical fluke, despite seed randomisation recruitment is the same."
    else
        @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) != pval
        @test reportCounts != [93, 208]
    end

    run( mpSim, saveConfig=false, seed=eseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )[[1, end]]
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) > 0.05
    @test all( reportCounts .== [105, 196] )

    run( mpSim, saveConfig=false, seed=s2seed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )[[1, end]]
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) > 0.05
    @test all( reportCounts .== [90, 211] )
end  # @testset "Disc random recruitment test"

@testset "PUnif random recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 1 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentDist!( recruitment, :pUnif,
        Dict( 5 => 1.0, 8 => 2.0, 11 => 1.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    p = [1.0, 1.0, 1.0, 2.0, 2.0, 2.0]
    p /= sum( p )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    pval = pvalue( ChisqTest( reportCounts, p ) )
    @test pval > 0.05
    @test all( reportCounts .== [28, 33, 38, 73, 62, 67] )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) == pval
    @test all( reportCounts .== [28, 33, 38, 73, 62, 67] )

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) != pval
    @test reportCounts != [28, 33, 38, 73, 62, 67]

    run( mpSim, saveConfig=false, seed=eseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) > 0.05
    @test all( reportCounts .== [38, 27, 36, 49, 77, 74] )

    run( mpSim, saveConfig=false, seed=s2seed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) > 0.05
    @test all( reportCounts .== [38, 27, 30, 68, 61, 77] )
end  # @testset "PUnif random recruitment test"

@testset "PLin random recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 1 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentDist!( recruitment, :pLin, Dict( 5 => 1.0, 10 => 11.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    p = [1.0, 3.0, 5.0, 7.0, 9.0, 11.0]
    p /= sum( p )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    pval = pvalue( ChisqTest( reportCounts, p ) )
    @test pval > 0.05
    @test all( reportCounts .== [5, 30, 34, 71, 72, 89] )

    run( mpSim, saveConfig=false, seed=piseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) == pval
    @test all( reportCounts .== [5, 30, 34, 71, 72, 89] )

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) != pval
    @test reportCounts != [5, 30, 34, 71, 72, 89]

    run( mpSim, saveConfig=false, seed=eseed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) > 0.05
    @test all( reportCounts .== [9, 28, 40, 76, 56, 92] )

    run( mpSim, saveConfig=false, seed=s2seed )
    report = nodeFluxReport( mpSim, 1, :in, "active" )["active"]
    report = report[:, Symbol( "external => active" )]
    reportCounts = counts( Vector{Int}( report ), 5:10 )
    @test pvalue( ChisqTest( reportCounts, p ) ) > 0.05
    @test all( reportCounts .== [10, 20, 35, 62, 80, 94] )
end  # @testset "PLin random recruitment test"

setSimulationLength!( mpSim, 12 )
subpop = Subpopulation( "All" )

@testset "Fixed age recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    setRecruitmentAgeFixed!( recruitment, 60 )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    run( mpSim, saveConfig=false )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    @test report[1, Symbol( 60.0 )] == 1000
end  # @testset "Fixed Age recruitment test"

@testset "Disc random age recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    setRecruitmentAgeDist!( recruitment, :disc,
        Dict( 60.0 => 1.0, 120.0 => 2.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )

    run( mpSim, saveConfig=false, seed=piseed )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    reportCounts = Vector{Int}( [report[1, Symbol( 60.0 )],
        report[1, Symbol( 120.0 )]] )
    pval = pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) )
    @test pval > 0.05
    @test all( reportCounts .== [344, 656] )

    for age in 72.0:12.0:108.0
        @test report[1, Symbol( age )] == 0
    end  # for age in 72.0:12.0:108.0

    run( mpSim, saveConfig=false, seed=piseed )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    reportCounts = Vector{Int}( [report[1, Symbol( 60.0 )],
        report[1, Symbol( 120.0 )]] )
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) == pval
    @test all( reportCounts .== [344, 656] )

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    reportCounts = Vector{Int}( [report[1, Symbol( 60.0 )],
        report[1, Symbol( 120.0 )]] )
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) != pval
    @test reportCounts != [344, 656]

    run( mpSim, saveConfig=false, seed=eseed )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    reportCounts = Vector{Int}( [report[1, Symbol( 60.0 )],
        report[1, Symbol( 120.0 )]] )
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) > 0.05
    @test all( reportCounts .== [346, 654] )

    run( mpSim, saveConfig=false, seed=s2seed )
    report = subpopulationAgeReport( mpSim, [0.0], 12, :age, subpop )["All"]
    reportCounts = Vector{Int}( [report[1, Symbol( 60.0 )],
        report[1, Symbol( 120.0 )]] )
    @test pvalue( ChisqTest( reportCounts, [1/3, 2/3] ) ) > 0.05
    @test all( reportCounts .== [323, 677] )
end  # @testset "Disc random age recruitment test"

@testset "PUnif random age recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    setRecruitmentAgeDist!( recruitment, :pUnif,
        Dict( 60.0 => 1.0, 96.0 => 3.0, 120.0 => 1.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    testDistribution1 = Uniform( 60, 96 )
    testDistribution2 = Uniform( 96, 120 )
    testDistribution = MixtureModel( [testDistribution1, testDistribution2],
        [0.25, 0.75] )

    run( mpSim, saveConfig=false, seed=piseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    pval = pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) )
    @test pval > 0.05

    run( mpSim, saveConfig=false, seed=piseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges,
        testDistribution ) ) == pval

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges,
        testDistribution ) ) != pval

    run( mpSim, saveConfig=false, seed=eseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05

    run( mpSim, saveConfig=false, seed=s2seed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05
end  # @testset "PUnif random age recruitment test"

@testset "PLin random age recruitment test" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    setRecruitmentAgeDist!( recruitment, :pLin,
        Dict( 60.0 => 1.0, 120.0 => 6.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    testDistribution = Truncated( SymTriangularDist( 120, 72 ), 60, 120 )

    run( mpSim, saveConfig=false, seed=piseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    pval = pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) )
    @test pval > 0.05

    run( mpSim, saveConfig=false, seed=piseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges,
        testDistribution ) ) == pval

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges,
        testDistribution ) ) != pval

    run( mpSim, saveConfig=false, seed=eseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05

    run( mpSim, saveConfig=false, seed=s2seed + 1 )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05
end  # @testset "PLin random age recruitment test"

@testset "PLin random age recruitment test 2" begin
    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    setRecruitmentAgeDist!( recruitment, :pLin,
        Dict( 60.0 => 1.0, 90.0 => 6.0, 120.0 => 1.0 ) )
    addSimulationRecruitment!( mpSim, recruitment )

    @test verifySimulation!( mpSim )
    testDistribution = Truncated( SymTriangularDist( 90, 36 ), 60, 120 )

    run( mpSim, saveConfig=false, seed=piseed - 1 )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    pval = pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) )
    @test pval > 0.05

    run( mpSim, saveConfig=false, seed=piseed - 1 )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) ==
        pval

    run( mpSim, saveConfig=false, seed=piseed - 1, sysEnt=true )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges,
        testDistribution ) ) != pval

    run( mpSim, saveConfig=false, seed=eseed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05

    run( mpSim, saveConfig=false, seed=s2seed )
    recruitmentAges = Vector{Float64}( DataFrame( DBInterface.execute(
        mpSim.simDB, "SELECT * FROM Personnel_Sim" ) )[:, :ageAtRecruitment] )
    @test pvalue( ExactOneSampleKSTest( recruitmentAges, testDistribution ) ) >
        0.05
end  # @testset "PLin random age recruitment test 2"

@testset "Attribute initialisation test" begin
    attribute = Attribute( "attribute B" )
    addInitialAttributeValue!( attribute, ("A", 1.0), ("B", 2.0), ("C", 1.0) )
    addSimulationAttribute!( mpSim, attribute )

    clearSimulationRecruitment!( mpSim )
    recruitment = Recruitment( "recruitment A" )
    setRecruitmentSchedule!( recruitment, 24 )
    setRecruitmentTarget!( recruitment, "node A" )
    setRecruitmentFixed!( recruitment, 1000 )
    addSimulationRecruitment!( mpSim, recruitment )

    subpop1 = Subpopulation( "B/A" )
    subpop2 = Subpopulation( "B/B" )
    subpop3 = Subpopulation( "B/C" )
    addSubpopulationCondition!.( [subpop1, subpop2, subpop3],
        MPcondition.( "attribute B", ==, ["A", "B", "C"] ) )
    setSubpopulationSourceNode!.( [subpop1, subpop2, subpop3], "node A" )

    @test verifySimulation!( mpSim )

    run( mpSim, saveConfig=false, seed=piseed )
    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3 )
    reportCounts = Int.( [report[1, Symbol( "B/A" )],
        report[1, Symbol( "B/B" )], report[1, Symbol( "B/C" )]] )
    @test pvalue( ChisqTest( reportCounts, [1/4, 1/2, 1/4] ) ) > 0.05
    @test all( reportCounts .== [250, 486, 264] )

    run( mpSim, saveConfig=false, seed=piseed )
    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3 )
    reportCounts = Int.( [report[1, Symbol( "B/A" )],
        report[1, Symbol( "B/B" )], report[1, Symbol( "B/C" )]] )
    @test all( reportCounts .== [250, 486, 264] )

    run( mpSim, saveConfig=false, seed=piseed, sysEnt=true )
    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3 )
    reportCounts = Int.( [report[1, Symbol( "B/A" )],
        report[1, Symbol( "B/B" )], report[1, Symbol( "B/C" )]] )
    @test reportCounts != [250, 486, 264]

    run( mpSim, saveConfig=false, seed=eseed )
    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3 )
    reportCounts = Int.( [report[1, Symbol( "B/A" )],
        report[1, Symbol( "B/B" )], report[1, Symbol( "B/C" )]] )
    @test pvalue( ChisqTest( reportCounts, [1/4, 1/2, 1/4] ) ) > 0.05
    @test all( reportCounts .==[243, 496, 261] )

    run( mpSim, saveConfig=false, seed=s2seed )
    report = subpopulationPopReport( mpSim, 12, subpop1, subpop2, subpop3 )
    reportCounts = Int.( [report[1, Symbol( "B/A" )],
        report[1, Symbol( "B/B" )], report[1, Symbol( "B/C" )]] )
    @test pvalue( ChisqTest( reportCounts, [1/4, 1/2, 1/4] ) ) > 0.05
    @test all( reportCounts .== [248, 500, 252] )
end  # @testset "Attribute initialisation test"

end  # @testset "Recruitment tests"


text:
include("$(@__FILE__)/../src/utils/containers.jl")

c1 = Container("hello", [1, 5])
c2 = Container("myboy", ["hjo", "hi"])

c3 = append!(c1, c2)
c4 = append!(c2, c1)


text:
# Default values shared across instruments
const default_line_width_mps = predict_intrinsic_stellar_line_width(5780,v_rot=1800)  # m/s
default_chunk_size_factor = 3       # For default_calc_chunk_width TODO: Figure out what value to use.  Ask Alex
default_calc_chunk_width = ChunkWidthFixedΔlnλ(default_chunk_size_factor*default_line_width_mps/speed_of_light_mps)  # Used by read_mask_espresso and read_vald for assigning lambda_lo and lambda_hi
default_min_chunk_Δv = 20000           # m/s  for ChunkWidthFixedΔlnλ

default_Δv_to_avoid_tellurics = 15000.0  # m/s

 # TODO: OPT: Make above const once settle on good values

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
module AbstractGPs

    using ChainRulesCore
    using Distributions
    using FillArrays
    using LinearAlgebra
    using Reexport
    @reexport using KernelFunctions
    using Random
    using Statistics

    using KernelFunctions: ColVecs, RowVecs

    export GP, mean, cov, std, cov_diag, mean_and_cov, marginals, rand,
        logpdf, elbo, dtc, posterior, approx_posterior, VFE, DTC, sampleplot,
        update_approx_posterior, LatentGP, ColVecs, RowVecs

    # Various bits of utility functionality.
    include(joinpath("util", "common_covmat_ops.jl"))

    # AbstractGP interface and FiniteGP interface.
    include(joinpath("abstract_gp", "abstract_gp.jl"))
    include(joinpath("abstract_gp", "finite_gp.jl"))

    # Basic GP object.
    include(joinpath("gp", "mean_function.jl"))
    include(joinpath("gp", "gp.jl"))

    # Efficient exact and approximate posterior GP implementations.
    include(joinpath("posterior_gp", "posterior_gp.jl"))
    include(joinpath("posterior_gp", "approx_posterior_gp.jl"))

    # LatentGP object to accomodate GPs with non-gaussian likelihoods.
    include(joinpath("latent_gp", "latent_gp.jl"))

    # Plotting utilities.
    include(joinpath("util", "plotting.jl"))
end # module


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/subsolver_config/config/ipopt.jl
# Configuration adjustment subroutines for Ipopt.
#############################################################################

function set_default_config!(::Val{:ipopt}, ext::ExtensionType, m::T) where {T <: MOI.AbstractOptimizer}

    MOI.set(m, MOI.RawParameter("max_iter"),3000)
    MOI.set(m, MOI.RawParameter("acceptable_tol"), 1E30)
    MOI.set(m, MOI.RawParameter("acceptable_iter"), 300)
    MOI.set(m, MOI.RawParameter("constr_viol_tol"), 0.000001)
    MOI.set(m, MOI.RawParameter("acceptable_compl_inf_tol"), 0.000001)
    MOI.set(m, MOI.RawParameter("acceptable_dual_inf_tol"), 1.0)
    MOI.set(m, MOI.RawParameter("acceptable_constr_viol_tol"), 0.000001)
    MOI.set(m, MOI.RawParameter("print_level"), 0)

    return nothing
end


text:
using Plots, Polyhedra
using LazySets: center, HalfSpace

include("simplify_set.jl")
include("analyze_sequence.jl")
include("sample.jl")

function create_plot(; layout=1)
    return plot(aspect_ratio=1, layout=layout)
end

function plot_results(x0::AbstractVector, ε::Number, P_ua::LazySet,
                      P_oa::LazySet)
    P_full = □(x0, ε)
    return plot_results(P_full, P_ua, P_oa)
end

function plot_results(P_full::LazySet, P_ua::LazySet, P_oa::LazySet)
    p = create_plot()
    plot!(p, P_full, color=:red)
    plot!(p, P_oa, color=:yellow)
    plot!(p, simplify_overapprox(P_oa), opacity=0.2, color=:yellow)
    plot!(p, P_ua, color=:green)
    plot!(p, simplify_underapprox(P_ua), opacity=0.2, color=:green)
    return p
end

function plot_results(ε::Number, P_uas::AbstractVector, P_oas::AbstractVector,
                      xs::AbstractVector)
    p = create_plot()
    m = length(xs)
    p = plot_results(xs[2], ε, P_uas[1], P_oas[1])
    plot!(p, tickfontsize=20)
    plots = [p]
    p = plot(plots...)
    return p
end

function plot_results(ε::Number, P_uas::AbstractVector, P_oas::AbstractVector,
                      xs::AbstractVector, samples_sat_list::AbstractVector,
                      samples_unsat_list::AbstractVector)
    p = create_plot()
    plot!(p, xticks=[round(xs[1][1]; digits=2)],
             yticks=[round(xs[1][2]; digits=2)])
    plot!(p, □(xs[1], ε), color=:yellow)
    plot!(p, □(xs[1], ε), color=:green)
    plots = [p]
    m = length(xs)
    for j in 2:m
        p = plot_results(xs[j], ε, P_uas[j-1], P_oas[j-1])
        plot!(p, xticks=[round(xs[j][1]; digits=2)],
                 yticks=[round(xs[j][2]; digits=2)])
        push!(plots, p)
    end
    for (samples_sat, samples_unsat) in zip(samples_sat_list,
                                            samples_unsat_list)
        p = plot_results(xs[m], ε, P_uas[m-1], P_oas[m-1])
        plot!(p, xticks=[round(xs[end][1]; digits=2)],
                 yticks=[round(xs[end][2]; digits=2)])
        for (samples, color) in [(samples_sat, :green), (samples_unsat, :red)]
            x = [s[1] for s in samples]
            y = [s[2] for s in samples]
            plot!(p, x, y, seriestype=:scatter, color=color, label=nothing)
        end
        push!(plots, p)
    end
    p = plot(plots...)
    return p
end

function analyze_sequence_example(α, npieces; log::Bool=false)
    A = [0.0 1.0    ; -1.0 0.0]
    B = [0.0 1.0 - α; -1.0 0.0]

    As = [A for i in 1:npieces]
    Bs = [B for i in 1:npieces]
    δ = 4*π / npieces
    times = [i * δ for i in 1:npieces]

    x0 = [1.0, 1.0]
    ε = 0.1
    m = 10
    ζ = ε / 100

    answer, definite_answer, P_uas, P_oas, xs =
        analyze_sequence(As, Bs, times, x0, ε, m, ζ; log=log, return_sets=true)

    if npieces == 1
        P_uas = [P_uas[1]]
        P_oas = [P_oas[1]]
        samples_sat = samples_unsat = []
        p = plot_results(ε, P_uas, P_oas, xs)
    else
        nsamples = 1000
        samples_sat, samples_unsat = sample_points(As, Bs, xs, P_oas[end],
                                                   times, ε, nsamples; log=log)
        p = plot_results(ε, P_uas, P_oas, xs, [samples_sat], [samples_unsat])
    end

    return p
end

p = analyze_sequence_example(0.01, 1)
savefig("membership_different_1.pdf")

p = analyze_sequence_example(0.01, 23)
savefig("membership_different_23.pdf")

p = analyze_sequence_example(0.0, 1)
savefig("membership_same_1.pdf")

p = analyze_sequence_example(0.0, 23)
savefig("membership_same_23.pdf")

max_stars_repo_name: HySynth/HySynth

max_issues_repo_name: HySynth/HySynth

max_forks_repo_name: HySynth/HySynth


text:
# Rohan McLure 2019 (C), Australian National University, licensed under MIT.

"""
Construct permitting synchronisation to be occuring simultaneously with a number of different processing preparing to wrote to our local array.
"""
mutable struct Handshake
    pid_events::Vector{Union{Channel{Nothing}, Nothing}}

    function Handshake(participants::Vector{Int64})
        largest_proc = maximum(procs())
        events = map(1:largest_proc) do pid
            if pid in participants
                return Channel{Nothing}(0)
            else
                return nothing
            end
        end
        new(events)
    end
end

"""
    put!(hs::Handshake, put_as::Int64)

Synchronisation primatives for custom synchronisation type.
"""
function put!(hs::Handshake, put_as::Int64)
    put!(hs.pid_events[put_as], nothing)
end

"""
    take!(hs::Handshake, take_from::Int64)

Synchronisation primatives for custom synchronisation type.
"""
function take!(hs::Handshake, take_from::Int64)
    take!(hs.pid_events[take_from])
end

"""
    ArrayChannel(T::Type, participants::Vector{Int64}, dims::Int64...)

Channel construct that associates itself with an array at each of the specified 'participating' processes.
Allows for synchronous, in place communication between bufferes on different processes by means of [`put!`](@ref), [`take!`](@ref), [scatter!](@ref), [gather!](@ref) or [reduce!](@ref).

# Examples
```jldoctest
julia> A = ArrayChannel(Float64, workers(), 2, 2)

Nothing displayed as master process is not part of workers(), and as such
does not host any local data.
```

"""
mutable struct ArrayChannel
    lock::Union{Nothing, ReentrantLock}
    cond_take::Union{Nothing, Handshake}
    cond_put::Union{Nothing, Channel{Nothing}}
    buffer::Union{Nothing, InPlaceArray}
    scratch1::Union{Nothing, InPlaceArray}
    scratch2::Union{Nothing, InPlaceArray}
    participants::Union{Nothing, Vector{Int64}}
    rrid::RRID

    # Null constructor
    function ArrayChannel(id::RRID)
        ch = new(nothing, nothing, nothing, nothing, nothing, nothing, nothing, id)
        channels[id] = ch
        return ch
    end

    # Headless constructor
    function ArrayChannel(A::InPlaceArray{T,N}, ps::Vector{Int64}, id::RRID, s1=nothing, s2=nothing) where {T,N}
        ps = sort(ps)
        ch = new(ReentrantLock(), Handshake(ps), Channel{Nothing}(1), A, s1, s2, ps, id)
        channels[id] = ch
        return ch
    end

    # Main constructor
    function ArrayChannel(T::Type, participants::Vector{Int64}, dims...)
        id = RRID()
        participants = sort(participants)
        ch = new(ReentrantLock(), Handshake(participants), Channel{Nothing}(1), nothing, nothing, nothing, participants, id)
        @sync for proc in participants
            if proc == myid()
                @async begin
                    ipa = InPlaceArray(Array{T}(undef, dims...), id)
                    s1 = InPlaceArray(Array{T}(undef, dims...), id, 1)
                    s2 = InPlaceArray(Array{T}(undef, dims...), id, 2)
                    ch.buffer = ipa
                    ch.scratch1 = s1
                    ch.scratch2 = s2
                    channels[id] = ch
                end
            else
                @spawnat proc begin
                    ipa = InPlaceArray(Array{T}(undef, dims...), id)
                    s1 = InPlaceArray(Array{T}(undef, dims...), id, 1)
                    s2 = InPlaceArray(Array{T}(undef, dims...), id, 2)
                    ArrayChannel(ipa, participants, id, s1, s2)
                end
            end
        end
        return ch
    end
end

function get_arraychannel(id::RRID)
    return channels[id]
end

"""
    put!(ac::ArrayChannel, send_to::Int64, [tag::Union{Distributed.RRID, Nothing}])

put! initiates two blocking remotecalls for each worker in the workerpool. The first waits on the receiver to authorises the buffer to be overwritten, the second writes the data.
"""
function put!(ac::ArrayChannel, send_to::Int64, tag::Union{ArrayChannel, RRID, Nothing}=nothing)
    if tag isa ArrayChannel
        tag = tag.rrid
    end
    same_channel = tag === nothing
    lock(ac.lock) do
        id = same_channel ? ac.rrid : (tag::RRID)
        place = ac.buffer
        remotecall_wait(send_to, id, myid()) do rrid, my_pid
            # From the rrid, get the ArrayChannel reference, and wait on cond_take
            X = get_arraychannel(rrid)
            put!(X.cond_take, my_pid)
        end
        if same_channel
            remotecall_wait(send_to, id, place) do id, payload
                # Serialise the InPlaceArray, but do nothing with it.
                X = get_arraychannel(id)
                put!(X.cond_put, nothing)
            end
        else
            # Target another channel's buffer by anonymously altering the buffer metadata
            target(place, send_to, 0, tag)
        end
    end
end

"""
    take!(ac::ArrayChannel, recv_from::Int64)

take! signals to the other owners of the ArrayChannel the intention to overwrite the buffer with new array data, and then waits for it to be written.
"""
function take!(ac::ArrayChannel, recv_from::Int64)
    # Accept any code
    take!(ac.cond_take, recv_from)
    lock(ac.lock) do
        take!(ac.cond_put)
    end
    return ac
end

function serialize(S::AbstractSerializer, ac::ArrayChannel)
    writetag(S.io, OBJECT_TAG)
    serialize(S, typeof(ac))
    serialize(S, ac.rrid)
end

function deserialize(S::AbstractSerializer, t::Type{<:ArrayChannel}) where {T,N}
    id = deserialize(S) :: RRID
    ch = ac_get_from(id)
    if ch ≠ nothing
        return ch
    end
    return ArrayChannel(id)
end

"""
    reduce!(op, ac::ArrayChannel, root::Int64)
Participate in a reduction on all processes participating in this ArrayChannel, where only the 'root' while be affected by the result.

Blocks until this process' role in the reduction is complete.

# Example
```jldoctest
Reduce by taking the vector sum, with root at process one
julia> reduce!(+, A, 1)
```
"""
function reduce!(op, ac::ArrayChannel, root::Int64)
    lock(ac.lock) do
        peers = ac.participants
        n = length(peers)
        pow_2 = 2^Int(ceil(log(2,n))) # Smallest power of two ≤
        ac_id = ac.rrid
        idx = indexin(myid(), peers)[1]
        root_idx = indexin(root, peers)[1]
        # Reshape around root
        pos = (idx < root_idx ? pow_2 - n : pow_2) - (idx - root_idx)
        z = pos
        lowest_z = pow_2 - n + 1
        k = 0
        leaf = z % 2 != 0 || z == lowest_z
        while z % 2 == 0
            if pos - 2^k >= lowest_z
                recv_idx = mod(idx + 2^k, n)
                recv_idx = recv_idx == 0 ? n : recv_idx
                recv_from = peers[recv_idx]
                take!(ac.cond_take, recv_from)
                take!(ac.cond_put)

                # Perform computation
                acc_buffer = myid() == root ? ac.buffer : ac.scratch1
                if k == 0
                    broadcast!(op, acc_buffer, ac.buffer, ac.scratch2)
                else
                    broadcast!(op, acc_buffer, acc_buffer, ac.scratch2)
                end
            end
            z ÷= 2
            k += 1
        end

        if myid() != root
            sender_idx = mod(idx - 2^k, n)
            sender_idx = sender_idx == 0 ? n : sender_idx
            send_to = peers[sender_idx]

            remotecall_wait(send_to, ac_id, peers[idx]) do rrid, sender_id
                X = get_arraychannel(rrid)
                put!(X.cond_take, sender_id)
            end

            if leaf
                target(ac.buffer, send_to, 2)
            else
                target(ac.scratch1, send_to, 2)
            end
        end
    end
end

"""
Since InPlaceArrays are serialised in place (subject to their buffer reference), use this method to anonymously alter the buffer reference before sendint to a buffer of choice.

Buffers are counted as follows:
The reference to the main buffer for the ArrayChannel is given a code zero.
Currently only to scratch buffers are needed for `reduce!`, however a positive integer will reference these.
"""
function target(ipa::InPlaceArray{T,N}, proc_id::Int64, buf::Int64, channel::Union{RRID, Nothing}=nothing) where {T,N}
    new_id = channel === nothing ? ipa.rrid : (channel::RRID)
    backup_id = ipa.rrid
    backup_buf_no = ipa.buf_no
    ipa.buf_no = buf
    ipa.rrid = new_id
    remotecall_wait(proc_id, ipa) do place
        X = get_arraychannel(place.rrid)
        put!(X.cond_put, nothing)
    end
    ipa.buf_no = backup_buf_no
    ipa.rrid = backup_id
end

# Required to show InPlaceArray objects
function show(S::IO, ac::ArrayChannel)
    if ac.buffer == nothing
        return
    end
	invoke(show, Tuple{IO, InPlaceArray}, S, ac.buffer::InPlaceArray)
end

function copy!(AC::ArrayChannel, src::AbstractArray)
    copy!(AC.buffer.src, src)
end

function copy!(dest::AbstractArray, AC::ArrayChannel)
    copy!(dest, AC.buffer.src)
end

function axes(AC::ArrayChannel)
    buffer = AC.buffer.src
    buffer === nothing ? nothing : axes(buffer)
end

function get_local_data(AC::ArrayChannel)
    return AC.buffer.src
end

function getindex(ac::ArrayChannel, keys...)
    getindex(ac.buffer.src, keys...)
end

function setindex!(ac::ArrayChannel, v, keys...)
    setindex!(ac.buffer.src, v, keys...)
end

function fill!(ac::ArrayChannel, v)
    fill!(ac.buffer.src, v)
end

function length(ac::ArrayChannel)
    length(ac.buffer.src)
end

function iterate(ac::ArrayChannel)
    iterate(ac.buffer.src)
end

function iterate(ac::ArrayChannel, state)
    iterate(ac.buffer.src)
end

function broadcast(op, As::ArrayChannel...)
    broadcast(op, map(x->x.buffer.src, As))
end

function broadcast!(op, dest::ArrayChannel, As::AbstractArray...)
    broadcast!(op, dest.buffer.src, As)
end

function broadcast!(op, dest::AbstractArray, As::ArrayChannel...)
    broadcast!(op, dest, map(x->x.buffer.src, As))
end

function ac_get_from(id::RRID)
    if id in keys(channels)
        return channels[id]
    end
    nothing
end

global channels = Dict{RRID, ArrayChannel}()


text:
module Registers

  using NativeDebugger
  export RegisterValue, RegisterSet
  abstract RegisterSet
  import Base: convert, getindex, -, +

  abstract Architecture

  """
  Represents the value of a machine register.

  T is the architecturally defined maximum size of the register.
  We keep both the register's value as well as a bitmask of which
  bits in the value we know to be valid.
  """
  immutable RegisterValue{T}
      value::T
      mask::T
  end
  convert{T<:Integer}(::Type{T}, r::RegisterValue) = convert(T,r[])
  convert{T}(::Type{RegisterValue{T}},value::T) = RegisterValue{T}(value, -1%T)
  getindex(v::RegisterValue) = v.value & v.mask
  isvalid{T}(v::RegisterValue{T}) = v.mask == (-1%T)
  invalidated{T}(v::RegisterValue{T}) = RegisterValue{T}(v.value, 0)
  for op in (:+, :-)
      @eval ($op)(r::RegisterValue, x::Integer) = RegisterValue(($op)(r.value, x), r.mask)
      @eval ($op)(x::Integer, r::RegisterValue) = RegisterValue(($op)(x, r.value), r.mask)
  end

  function Base.show(io::IO, r::RegisterValue)
      value = hex(r[], 2sizeof(r[]))
      print_with_color(:green, io, "0x")
      for (i,c) in enumerate(value)
          shift = (8sizeof(r[]) - 4i)
          b = (r.mask & (UInt64(0xf) << shift)) >> shift
          print_with_color((b == 0xf) ? :green :
                           (b == 0x0) ? :red :
                           :yellow, io, string(c))
      end
  end

  import NativeDebugger: getarch

  # RegisterSets should ideally support the following operations
  export invalidate_regs!, set_sp!, set_ip!, set_dwarf!, get_dwarf,
    getarch, intptr

  # This will be provided automatically
  export get_syscallarg

  function ip end
  function invalidate_regs! end
  function set_sp! end
  function set_ip! end
  function set_dwarf! end
  function get_dwarf end
  function intptr end

  function get_syscallarg(regs, idx)
     get_syscallarg(getarch(regs), regs, idx)
  end

  # Convenience fallback
  function get_dwarf(RC, sym::Symbol)
     get_dwarf(getarch(RC), RC, sym)
  end

  function set_dwarf!(RC, sym::Symbol, val)
     set_dwarf!(getarch(RC), RC, sym, val)
  end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""ClusterRole is a cluster level, logical grouping of PolicyRules that can be referenced as a unit by a RoleBinding or ClusterRoleBinding.

    IoK8sApiRbacV1ClusterRole(;
        aggregationRule=nothing,
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        rules=nothing,
    )

    - aggregationRule::IoK8sApiRbacV1AggregationRule : AggregationRule is an optional field that describes how to build the Rules for this ClusterRole. If AggregationRule is set, then the Rules are controller managed and direct changes to Rules will be stomped by the controller.
    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata.
    - rules::Vector{IoK8sApiRbacV1PolicyRule} : Rules holds all the PolicyRules for this ClusterRole
"""
mutable struct IoK8sApiRbacV1ClusterRole <: SwaggerModel
    aggregationRule::Any # spec type: Union{ Nothing, IoK8sApiRbacV1AggregationRule } # spec name: aggregationRule
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiRbacV1PolicyRule} } # spec name: rules

    function IoK8sApiRbacV1ClusterRole(;aggregationRule=nothing, apiVersion=nothing, kind=nothing, metadata=nothing, rules=nothing)
        o = new()
        validate_property(IoK8sApiRbacV1ClusterRole, Symbol("aggregationRule"), aggregationRule)
        setfield!(o, Symbol("aggregationRule"), aggregationRule)
        validate_property(IoK8sApiRbacV1ClusterRole, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiRbacV1ClusterRole, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiRbacV1ClusterRole, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiRbacV1ClusterRole, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        o
    end
end # type IoK8sApiRbacV1ClusterRole

const _property_map_IoK8sApiRbacV1ClusterRole = Dict{Symbol,Symbol}(Symbol("aggregationRule")=>Symbol("aggregationRule"), Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("rules")=>Symbol("rules"))
const _property_types_IoK8sApiRbacV1ClusterRole = Dict{Symbol,String}(Symbol("aggregationRule")=>"IoK8sApiRbacV1AggregationRule", Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("rules")=>"Vector{IoK8sApiRbacV1PolicyRule}")
Base.propertynames(::Type{ IoK8sApiRbacV1ClusterRole }) = collect(keys(_property_map_IoK8sApiRbacV1ClusterRole))
Swagger.property_type(::Type{ IoK8sApiRbacV1ClusterRole }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiRbacV1ClusterRole[name]))}
Swagger.field_name(::Type{ IoK8sApiRbacV1ClusterRole }, property_name::Symbol) =  _property_map_IoK8sApiRbacV1ClusterRole[property_name]

function check_required(o::IoK8sApiRbacV1ClusterRole)
    true
end

function validate_property(::Type{ IoK8sApiRbacV1ClusterRole }, name::Symbol, val)
end


text:
##

using ACE, ACEbase
using Printf, Test, LinearAlgebra, ACE.Testing, Random
using ACE: evaluate, evaluate_d, SymmetricBasis, PIBasis,
           grad_config, grad_params, O3
using ACEbase.Testing: fdtest

##

@info("Basic test of LinearACEModel construction and evaluation")

# construct the 1p-basis
maxdeg = 6
ord = 3
Bsel = SimpleSparseBasis(ord, maxdeg)

B1p = ACE.Utils.RnYlm_1pbasis(; maxdeg=maxdeg)

# generate a configuration
nX = 10
Xs = rand(PositionState{Float64}, B1p.bases[1], nX)
cfg = ACEConfig(Xs)

φ = ACE.Invariant()
basis = SymmetricBasis(φ, B1p, O3(), Bsel)

BB = evaluate(basis, cfg)
c = rand(length(BB)) .- 0.5
naive = ACE.LinearACEModel(basis, c, evaluator = :naive)
standard = ACE.LinearACEModel(basis, c, evaluator = :standard)

## FIO

@info("Check FIO")
using ACEbase.Testing: test_fio
println(@test(all(test_fio(naive; warntype = false))))
println(@test(all(test_fio(standard; warntype = false))))

##

evaluate(naive, cfg) ≈  evaluate(standard, cfg)
grad_params(naive, cfg) ≈ grad_params(standard, cfg)
grad_config(naive, cfg) ≈ grad_config(standard, cfg)

evaluate_ref(basis, cfg, c) = sum(evaluate(basis, cfg) .* c)

grad_config_ref(basis, cfg, c) = permutedims(evaluate_d(basis, cfg)) * c

grad_params_ref(basis, cfg, c) = evaluate(basis, cfg)

grad_params_config_ref(basis, cfg, c) = evaluate_d(basis, cfg)

(fun, funref, str) = (ACE.grad_params_config, grad_params_config_ref, "grad_params_config")

for (fun, funref, str) in [
         (evaluate, evaluate_ref, "evaluate"),
         (ACE.grad_config, grad_config_ref, "grad_config"),
         (ACE.grad_params, grad_params_ref, "grad_params"),
         (ACE.grad_params_config, grad_params_config_ref, "grad_params_config"),
      ]
   @info("Testing `$str` for different model evaluators")
   for ntest = 1:30
      local c, cfg
      cfg = rand(PositionState{Float64}, B1p.bases[1], nX) |> ACEConfig
      c = rand(length(basis)) .- 0.5
      ACE.set_params!(naive, c)
      ACE.set_params!(standard, c)
      val = funref(basis, cfg, c)
      val_naive = fun(naive, cfg)
      val_standard = fun(standard, cfg)
      print_tf(@test( val ≈ val_naive ≈ val_standard ))
   end
   println()
end

##


text:
abstract type AbstractChain end

# Chain interface
Base.getindex(ch::AbstractChain, i::Integer) = throw(MethodError(getindex, (typeof(ch),Integer)))
Base.length(ch::AbstractChain) = throw(MethodError(length, (typeof(ch),)))
Base.eltype(ch::AbstractChain) = throw(MethodError(eltype, (typeof(ch),)))
function Base.iterate(ch::AbstractChain, state=1)
    state > length(ch) && return nothing
    return (ch[state], state+1)
end

"""Empty Chain"""
struct EmptyChain <: AbstractChain end
Base.length(ch::EmptyChain) = 0

""" k-Chain is a formal linear combination of k-cells.

    Given a set E, we can construct a free R-module M that has E ⊆ M as a basis M, the module M is the module of the formal linear combinations of elements of E, or free module over E, R^{(E)}, s.t. given a finite subset {X_1, ..., X_n} of E, a formal linear combination of X_1, ..., X_n is

    a_1 X_1 + ··· + a_n X_n,

    where the a_i ∈ R.
"""
mutable struct Chain{R} <: AbstractChain
    dim::Int
    coefs::Vector{R}   # coefficients from ring R
    elems::Vector{Int} # indexes to elements of E
end
Chain(coefs::Vector{R}, elems::Vector{Int}) where {R} = Chain{R}(0, coefs, elems)
Chain(dim::Int, ::Type{R}) where {R} = Chain{R}(dim, R[], Int[])
Chain(::Type{R}) where {R} = Chain(0,R)
dim(ch::Chain) = ch.dim
setdim!(ch::Chain, dim::Int) = (ch.dim = dim)

# implement interface
Base.getindex(ch::Chain, i::Integer) = (ch.coefs[i], ch.elems[i])
Base.length(ch::Chain) = length(ch.coefs)
Base.eltype(ch::Chain{R}) where {R} = R

function Base.show(io::IO, ch::Chain)
    if length(ch.elems) == 0
        print(io, "0")
    else
        for (i,(c,el)) in enumerate(zip(ch.coefs, ch.elems))
            i != 1 && print(io, " + ")
            print(io, "$c[$el]")
        end
    end
end

function Base.push!(ch::Chain{R}, coef::R, idx::Int) where {R}
    push!(ch.coefs, coef)
    push!(ch.elems, idx)
    return ch
end
Base.push!(ch::Chain{R}, e::Tuple{R,Int}) where {R} = push!(ch, e[1], e[2])
Base.push!(ch::Chain{R}, e::Pair{R,Int}) where {R}  = push!(ch, e[1], e[2])
function Base.append!(a::Chain, b::Chain)
    append!(a.coefs, b.coefs)
    append!(a.elems, b.elems)
    return a
end
+(ch::Chain{R}, e::Tuple{R,Int}) where {R} = push!(ch, e[1], e[2])

+(a::Chain{R}, b::Chain{R}) where {R} = Chain{R}(dim(a), vcat(a.coefs, b.coefs), vcat(a.elems, b.elems))
-(a::Chain{R}, b::Chain{R}) where {R} = Chain{R}(dim(a), vcat(a.coefs, -b.coefs), vcat(a.elems, b.elems))

*(ch::Chain{R}, r::R) where {R} = Chain{R}(dim(ch), ch.coefs*r, ch.elems)
*(r::R, c::Chain{R}) where {R} = c*r

function simplify(ch::Chain)
    R = eltype(ch)
    tmpChain = Dict{Int,R}()
    for (c,el) in ch
        if haskey(tmpChain, el)
            tmpChain[el] += c
        else
            tmpChain[el] = c
        end
    end
    cc = Chain(dim(ch), R)
    for (el,c) in tmpChain
        c != zero(R) && push!(cc, (c,el))
    end
    return cc
end


text:
include("simAn.jl")

"""
    assemble!(ata_model::AbstractModel; solver="jump", starting_design=Matrix{Float64}(undef, 0, 0), results_folder="results", start_temp=0.1, geom_temp=0.1, n_item_sample=1, n_test_sample=1, opt_feas=0.0, n_fill=1, max_time=1000.00, max_conv=2, feas_nh=0, opt_nh=5, verbosity=2, optimizer_constructor="GLPK", optimizer_attributes=[("tm_lim", 1000)])

# Description

Assemble the tests.

# Arguments

- **`ata_model::AbstractModel`** : Required. The model built with ATA fuctions.
- **`solver`** : Optional. Default: `"jump"`. Values: `"jump"`, `"siman"`. The solving interface to be used (JuMP or internal solver based on Simulated Annealing).
- **`starting_design`** : Optional. Default: `Matrix{Float64}(undef, 0, 0)`. The starting design matrix. Must be a `Matrix{Float64}`.
- **`results_folder`** : Optional. Default: `"results"`. The folder in which the output is stored.

## siman arguments

  - **`start_temp`** : Optional. Default: `0.1`. Values: `[0, Inf]`. Starting temperature, set to minimum for short journeys (if 0 worse solutions will never be accepted).
  - **`geom_temp`** : Optional. Default: `0.1`. Values: `[0, Inf)`. Decreasing geometric factor.
  - **`n_item_sample`** : Optional. Default: `1`. Values: `[1, Inf]`. Number of items to alter. Set to minimum for a shallow analysis, set to maximum for a deep analysis of the neighbourhoods.
  - **`n_test_sample`** : Optional. Default: `1`. Values: `[1, Inf]`. Number of tests to alter. Set to minimum for a shallow analysis, set to maximum for a deep analysis of the neighbourhoods.
  - **`opt_feas`** : Optional. Default: `0.0`. Values: `[0, Inf)`. Optimality/feasibility balancer, if 0 only feasibility of solution is analysed. Viceversa, if 1, only optimality is considered (uncontrained model). All the other values are accepted but produce uninterpretable results.
  - **`n_fill`** : Optional. Default: `1`. Values: `[0, Inf)`. Number of fill-up phases, usually 1 is sufficient, if start_temp is high it can be higher. If a starting_design is supplied, it can be set to 0.
  - **`verbosity`** : Optional. Default: `2`. Values: `1` (minimal), `2` (detailed). Verbosity level. In the console '+' stands for improvement, '_' for accepting worse solution. The dots are the fill-up improvement steps.

    * Termination criteria

      - **`max_time`** : Optional. Default: `1000.0`. Values: `[0, Inf)`. Time limit in seconds.
      - **`max_conv`** : Optional. Default: `2`. Values: `[1, Inf)`. Maximum convergence, stop when, after max_conv rounds no improvements have been found. Set to minimum for shallow analysis, increase it for deep analysis of neighbourhoods.
      - **`feas_nh`** : Optional. Default: `0`. Values: `[1, Inf)`. Maximum number of Feasibility neighbourhoods to explore, set to the minimum if the model is small or not highly constrained.
      - **`opt_nh`** : Optional. Default: `5`. Values: `[1, Inf)`. Maximum number of Optimality neighbourhoods to explore, set to the minimum if the model is highly constrained.

## jump arguments

  - **`optimizer_constructor`** : Optional. Default: `"GLPK"`. Values: `"GLPK"`, `"Knitro"`, `"Gurobi"`, `"Cbc"`, `"CPLEX"`, `"Xpress"`, `"SCIP"`, `"Juniper"`. JuMP solver selection. Remember to load the required package before assemble!.
  - **`optimizer_attributes`** : Optional. Default: `[("tm_lim", 1000)]`. Values: An array of pairs `(attribute, value)`. Attributes and related values for the JuMP solver.

## other keyword arguments
  - **`kwargs...`** : Optional.
"""
function assemble!(
    ata_model::AbstractModel;
    solver = "jump",
    starting_design = Matrix{Float64}(undef, 0, 0),
    results_folder = "results",
    start_temp = 0.1,
    geom_temp = 0.1,
    n_item_sample = 1,
    n_test_sample = 1,
    opt_feas = 0.0,
    n_fill = 1,
    max_time = 1000.0,
    max_conv = 2,
    feas_nh = 0,
    opt_nh = 5,
    verbosity = 2,
    optimizer_constructor = "GLPK",
    optimizer_attributes = [("tm_lim", 1000)],
    kwargs...,
)
    _group_by_friends!(ata_model)
    if solver == "siman"
        siman!(
            ata_model;
            results_folder = results_folder,
            starting_design = starting_design,
            start_temp = start_temp,
            geom_temp = geom_temp,
            max_time = max_time,
            n_item_sample = n_item_sample,
            n_test_sample = n_test_sample,
            max_conv = max_conv,
            verbosity = verbosity,
            opt_feas = opt_feas,
            n_fill = n_fill,
            feas_nh = feas_nh,
            opt_nh = opt_nh,
            kwargs...,
        )
    elseif solver == "jump"
        jump!(
            ata_model;
            results_folder = results_folder,
            starting_design = starting_design,
            optimizer_constructor = optimizer_constructor,
            optimizer_attributes = optimizer_attributes,
            kwargs...,
        )
    else
        error!(ata_model, "only \"siman\" and \"jump\" are supported.")
    end
    return nothing
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

module KrigingEstimators

using Meshes
using GeoStatsBase
using Variography

using LinearAlgebra: Factorization, Symmetric
using LinearAlgebra: bunchkaufman, cholesky, issuccess, ⋅
using Combinatorics: multiexponents
using Unitful

import GeoStatsBase: fit, predict, status
import GeoStatsBase: solve

include("estimators.jl")

export
  # estimators
  KrigingEstimator,
  SimpleKriging,
  OrdinaryKriging,
  UniversalKriging,
  ExternalDriftKriging,
  fit, predict, status,
  weights, combine

end


text:
using BayesianDataFusion
using Base.Test

using DataFrames

A = randn(15,2);
B = randn(20,2);
C = randn(2,2);

X = Float64[ sum(A[i,:].*B[j,:].*C[k,:]) for i in 1:size(A,1), j in 1:size(B,1), k in 1:size(C,1)]

df = DataFrame(A=Int64[], B=Int64[], C=Int64[], v=Float64[])
for i=1:size(A,1), j=1:size(B,1), k=1:size(C,1)
  push!(df, Any[i, j, k, X[i,j,k]])
end

rd = RelationData(df)
@test rd.entities[1].name == "A"
@test rd.entities[2].name == "B"
@test rd.entities[3].name == "C"

assignToTest!(rd.relations[1], 10)
result = macau(rd, burnin=50, psamples=10, num_latent=2, verbose=false)

## checking prediction for whole matrix for tensors
Yhat = pred_all(rd.relations[1])
@test size(Yhat) == size(X)
yprod   = rd.entities[1].model.sample[:,4]
yprod .*= rd.entities[2].model.sample[:,2]
yprod .*= rd.entities[3].model.sample[:,1]
@test_approx_eq Yhat[4,2,1] sum(yprod)+rd.relations[1].model.mean_value