text:
using AcronymGenerator
using Documenter

makedocs(;
    modules=[AcronymGenerator],
    authors="Jared Lumpe",
    repo="https://github.com/jlumpe/AcronymGenerator.jl/blob/{commit}{path}#L{line}",
    sitename="AcronymGenerator.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
function abc!(r)
    k = 0
    s = 0
    A = 0
    while (isborder(r,Nord) && isborder(r,West))==false
        while !isborder(r , Nord)
            move!(r, Nord)
        end
        while !isborder(r, West)
            move!(r,West)
        end
    end
    while !isborder(r,Ost)
        move!(r, Ost)
        s = s + 1
    end
    while !isborder(r,West)
        move!(r, West)
    end
    A = s
    while (isborder(r,Sud) && isborder(r,Ost))==false
        while s > 0
            while !isborder(r, Ost)
                move!(r, Ost)
                s = s - 1
            end
            k = k + 1
            r = 0
            while (isborder(r, Ost) && !isborder(r,Nord))==true
                move!(r,Nord)
                r = r + 1
            end
            if !isborder(r,Ost)
                move!(r,Ost)
                s = s - 1
                if isborder(r, Sud)
                    k = k + 1
                end
                while isborder(r,Sud)
                    move!(r,Ost)
                    s = s - 1
                end
            end
            while r > 0
                move!(r, Sud)
                r = r - 1
            end
        end
        s = A
        if !isborder(r, Sud)
            move!(r,Sud)
            while s > 0
                while !isborder(r,West)
                    move!(r, West)
                    s = s - 1
                end
                k = k + 1
                l = 0
                while (isborder(r, West) && !isborder(r,Nord))==true
                    move!(r, Nord)
                    l = l + 1
                end
                if !isborder(r,West)
                    move!(r,West)
                    s = s - 1
                    if isborder(r, Sud)
                        k = k + 1
                    end
                    while isborder(r,Sud)
                        move!(r,West)
                        s = s - 1
                    end
                end
                while l > 0
                    move!(r, Sud)
                    l = l - 1
                end
            end
        end
        if !isborder(r, Sud)
            move!(r,Sud)
        end
    end
    println(k)
end

ext: jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from DCM to Euler angle and axis.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export dcm_to_angleaxis

"""
    dcm_to_angleaxis(dcm::DCM{T}) where T<:Number

Convert the `dcm` to an Euler angle and axis representation.

By convention, the returned Euler angle will always be in the interval [0, π].
"""
function dcm_to_angleaxis(dcm::DCM{T}) where T<:Number
    cθ = (dcm[1, 1] + dcm[2, 2] + dcm[3, 3] - 1) / 2

    # Check the undefined case.
    if cθ ≥ 1 - eps()
        return EulerAngleAxis(T(0), SVector{3,T}(0, 0, 0))
    elseif cθ <= -1 + eps()
        v₁ = sqrt((1 + dcm[1, 1]) / 2)
        v₂ = sqrt((1 + dcm[2, 2]) / 2)
        v₃ = sqrt((1 + dcm[3, 3]) / 2)

        # Compute the sign of the vector components.
        if dcm[1, 2] ≥ 0
            if dcm[1, 3] ≥ 0
                s₁ = s₂ = s₃ = +1
            else
                s₁ = s₂ = +1
                s₃ = -1
            end
        else
            if dcm[1, 3] ≥ 0
                s₁ = s₃ = +1
                s₂ = -1
            else
                s₁ = +1
                s₂ = s₃ = -1
            end
        end

        return EulerAngleAxis(T(π), [s₁ * v₁, s₂ * v₂, s₃ * v₃ ])
    else
        sθ2 = 2sqrt(1 - cθ * cθ)

        v₁ = (dcm[2, 3] - dcm[3, 2]) / sθ2
        v₂ = (dcm[3, 1] - dcm[1, 3]) / sθ2
        v₃ = (dcm[1, 2] - dcm[2, 1]) / sθ2

        return EulerAngleAxis(acos(cθ), [v₁, v₂, v₃])
    end
end


text:
include("timestamp.jl")
include("operators.jl")
include("conversion.jl")
include("arraymethods.jl")
include("utils.jl")


text:
"""
    PMC(; epsilon, K, prior, ess_min = 0.5, proposal = defaultproposal(prior))

Population Monte Carlo ABC structure, with `epsilon` schedule, `K` particles,
`prior`, resampling threshold `ess_min` and `proposal` distribution.
The `epsilon` schedule should be an iterable, e.g. an array of numbers,
[`EpsilonExponentialDecay`](@ref) or [`EpsilonLinearDecay`](@ref).
One dimensional problems should be defined in vectorized form (
see [`run!`](@ref) for an example).

# Example
```julia-repl
pmc = PMC(epsilon = EpsilonExponentialDecay(1, .01, .9), K = 20,
                 prior = MultivariateUniform([-1, -1], [1, 1]))
model(theta) = .1 * randn(2) .+ theta
data = [.3, -.2]
run!(pmc, model, data)

mean(pmc)
particles(pmc)
weights(pmc)
```

[Beaumont, M. A., Cornuet, J., Marin, J., and Robert, C. P. (2009)
Adaptive approximate Bayesian computation. Biometrika,96, 983-990.
](http://dx.doi.org/10.1093/biomet/asp052)
"""
struct PMC{E,P,D,D1,T} <: AbstractABC
    epsilon::E
    log_ess_min::Float64
    particles::Vector{P}
    log_weights::Vector{Float64}
    prior::D
    proposal::D1
    distance::T
end
function PMC(; epsilon, prior::Distribution{T}, K = 10^3, ess_min = 0.5,
             proposal = defaultproposal(prior), distance = euclidean) where T
    particles = init_particles(prior, K)
    PMC(epsilon, log(ess_min), particles, fill(-log(K), K),
        prior, proposal(particles), distance)
end
function Base.show(io::IO, mime::MIME"text/plain", pmc::PMC)
    println(io, "PMC(K = $(length(pmc.particles)), ess_min = $(exp(pmc.log_ess_min)))")
    print(io, "  epsilon schedule: ")
    show(io, mime, pmc.epsilon)
    print(io, "  prior: ")
    show(io, mime, pmc.prior)
    print(io, "  proposal: ")
    show(io, mime, pmc.proposal)
end
weights(pmc::PMC) = weights(pmc.log_weights)
particles(pmc::PMC) = pmc.particles

sample_particles!(particles, distribution, eps, model, data, distance) =
    sample_particles!(Random.GLOBAL_RNG, particles, distribution, eps, model, data, distance)
function sample_particles!(rng::Random.AbstractRNG,
                           particles, distribution, eps, model, data, distance)
    n_sims = 0
    Threads.@threads for particle in particles
        while true
            rand!(rng, distribution, particle)
            sim_data = model(particle)
            n_sims += 1
            distance(sim_data, data) < eps && break
        end
    end
    n_sims
end
resample(particles, weights; verbose = false, log_weights = false) = resample!(Random.GLOBAL_RNG, copy(particles), weights, verbose = verbose, log_weights = log_weights)[1]
function resample!(rng::Random.AbstractRNG, particles, weights;
                   verbose = true,
                   log_weights = true)
    verbose && println("Resampling...")
    idxs = wsample(1:length(weights), log_weights ? exp.(weights) : weights, length(weights))
    old_particles = deepcopy(particles) # only particles[idxs] would be needed
    for (i, j) in enumerate(idxs)
        copyto!(particles[i], old_particles[j])
    end
    weights .= log_weights ? -log(length(weights)) : 1/length(weights)
    particles, weights
end
function update_log_weights!(log_weights, particles, prior, proposal)
    for i in eachindex(log_weights)
        log_weights[i] = logpdf(prior, particles[i]) - logpdf(proposal, particles[i])
    end
    log_weights .-= logsumexp(log_weights)
    log_weights
end

"""
    run!([rng], method, model, data;
         callback = () -> nothing, maxfevals = Inf, verbose = true)

Run `method` on `model` and `data`.
The function `callback` gets evaluated after every iteration.
The method stops after the first iteration that reaches more than `maxfevals` calls to the model.
# Example
```julia-repl
using Distributions, LinearAlgebra
pmc = PMC(epsilon = [1, .5, .2, .1, .01, .001], K = 10^3,
                 prior = TruncatedMultivariateNormal([1.], 2I,
                                                     lower = [0], upper = [3]))
model(theta) = theta * randn() .+ 1.2
data = [.5]
callback() = @show mean(pmc)
run!(pmc, model, data, callback = callback, maxfevals = 10^5)

using StatsBase, StatsPlots
h = fit(Histogram, vcat(particles(pmc)...), nbins = 50)
StatsPlots.plot(h)
```
"""
run!(pmc, model, data; kwargs...) = run!(Random.GLOBAL_RNG, pmc, model, data; kwargs...)
# Beaumont et al. 2009 (see also Lintusaari et al. 2016
function run!(rng::AbstractRNG, pmc::PMC, model, data;
              callback = () -> nothing, maxfevals = Inf, verbose = true)
    start_time = now()
    n_sims = Int[]
    K = length(pmc.particles)
    push!(n_sims, sample_particles!(rng, pmc.particles, pmc.prior,
                                    first(pmc.epsilon), model, data, pmc.distance))
    update!(pmc.proposal, pmc.particles, pmc.log_weights)
    callback()
    verbose && @printf "%10.s %13.s %8.s %8.s %8.s\n" "iteration" "elapsed" "fevals" "epsilon" "ess"
    verbose && @printf "%10.s %13.s %8.s %8.e %8.s\n" 0 round(now() - start_time, Second) n_sims[1] first(pmc.epsilon) 1
    for (i, eps) in enumerate(Iterators.drop(pmc.epsilon, 1))
        push!(n_sims, sample_particles!(rng, pmc.particles, pmc.proposal, eps, model, data, pmc.distance))
        update_log_weights!(pmc.log_weights, pmc.particles, pmc.prior, pmc.proposal)
        l_ess = log_ess(pmc.log_weights) - log(K)
        if l_ess < pmc.log_ess_min
            resample!(rng, pmc.particles, pmc.log_weights, verbose = verbose)
        end
        update!(pmc.proposal, pmc.particles, pmc.log_weights)
        callback()
        verbose && @printf "%10.s %13.s %8.s %8.e %8.s\n" i round(now() - start_time, Second) sum(n_sims) eps exp(l_ess)
        sum(n_sims) > maxfevals && break
    end
    (weights = exp.(pmc.log_weights),
     particles = pmc.particles,
     n_sims = n_sims)
end

struct Particle{T}
    θ::Vector{Float64}
    x::T
    distances::Vector{Float64}
end
function Base.copyto!(dest::Particle, bc::Particle)
    dest.θ .= bc.θ
    for k in eachindex(dest.x)
        dest.x[k] .= bc.x[k]
    end
    dest.distances .= bc.distances
    dest
end
theta(p::Particle) = p.θ

"""
    AdaptiveSMC(; alpha = 0.9, epsilon = .1, ess_min = 0.5, M = 1,
                  prior, K = 10^3, proposal = defaultproposal(prior))

Adaptive Sequential Monte Carlo structure with `K` particles,
`M` calls of the model per parameter value, final `epsilon`,
decrease parameter `alpha` and resampling threshold `ess_min`.
See also [`PMC`](@ref).

# Example
```
using Distributions, LinearAlgebra
asmc = AdaptiveSMC(prior = MultivariateNormal(zeros(2), 4I))
model(theta) = vcat([.3*randn(2) .+ theta for _ in 1:2]...) # 2 i.i.d samples
data = [.4, -.3, .5, -.2]
run!(asmc, model, data)

mean(asmc)
using StatsPlots
p = hcat(particles(asmc)...); StatsPlots.scatter(p[1, :], p[2, :])
```

[Del Moral, P., Doucet, A., and Jasra, A. (2012) An adaptive
sequential Monte Carlo method for approximate Bayesian
computation. Statistics and Computing, 22, 1009-1020
](http://dx.doi.org/10.1007/s11222-011-9271-y)
"""
struct AdaptiveSMC{P,P1,D,T} <: AbstractABC
    α::Float64
    ϵ::Float64
    Nt::Float64
    prior::P
    proposal::P1
    particles::Vector{Particle{T}}
    weights::Vector{Float64}
    distance::D
end
function AdaptiveSMC(; alpha = 0.9, epsilon = .1, ess_min = 0.5, M = 1,
                     prior, K = 10^3,
                     proposal = defaultproposal(prior),
                     distance = euclidean,
                     datatype = distance == euclidean ? Vector{Float64} :
                                                        Vector{Vector{Float64}}
                    )
    particles = [Particle(randn(eltype(prior), length(prior)),
                          [datatype(undef, 0) for _ in 1:M],
                          zeros(M)) for _ in 1:K]
    AdaptiveSMC(alpha, float(epsilon), float(ess_min * K), prior,
                proposal(theta.(particles)),
                particles, zeros(K), distance)
end
Distributions.mean(x::Union{<:PMC, <:AdaptiveSMC}) = sum(particles(x) .* weights(x))
weights(asmc::AdaptiveSMC) = asmc.weights
particles(asmc::AdaptiveSMC) = theta.(asmc.particles)
function Base.show(io::IO, mime::MIME"text/plain", asmc::AdaptiveSMC)
    println(io, "AdaptiveSMC(α = $(asmc.α), ϵ = $(asmc.ϵ), K = $(length(asmc.particles)), ess_min = $(asmc.Nt/length(asmc.particles)), M = $(length(asmc.particles[1].distances)))")
    print(io, "  prior: ")
    show(io, mime, asmc.prior)
    print(io, "  proposal: ")
    show(io, mime, asmc.proposal)
end
function ess!(weights, particles, factor, e)
    for i in eachindex(weights)
        weights[i] = factor[i] * sum(particles[i].distances .< e)
    end
    maximum(weights) == 0 && return 0
    weights ./= sum(weights)
    ess(weights)
end
function getfactor(weights, particles, e_old)
    factor = zeros(length(weights))
    smallest_dist = Inf
    for i in eachindex(factor)
        weights[i] == 0 && continue
        tmp = minimum(particles[i].distances)
        tmp < smallest_dist && (smallest_dist = tmp)
        denom = sum(particles[i].distances .< e_old)
        factor[i] = weights[i] / denom
    end
    smallest_dist, factor
end
function new_epsilon!(weights, e_old, particles, alpha, e_min)
    smallest_dist, factor = getfactor(weights, particles, e_old)
    target = alpha * ess(weights)
    e = find_zero(e -> target - ess!(weights, particles, factor, e),
                  (smallest_dist, e_old), Bisection())
    e = max(e, e_min)
    e, ess!(weights, particles, factor, e)
end
function sample_particles!(rng::Random.AbstractRNG,
                           particles::Vector{<:Particle},
                           distribution, model, data, weights, distance;
                           indexed = false)
    n_sims = 0
    for (i, (w, particle)) in enumerate(zip(weights, particles))
        w == 0 && continue
        if indexed
            Distributions._rand!(rng, distribution, particle.θ, i)
        else
            rand!(rng, distribution, particle.θ)
        end
        for k in eachindex(particle.x)
            particle.x[k] = model(particle.θ)
            n_sims += 1
            particle.distances[k] = distance(particle.x[k], data)
        end
    end
    n_sims
end
function mh_accept!(particles, old_particles, prior, ϵ, weights)
    for i in eachindex(particles)
        weights[i] == 0 && continue
        logw = logpdf(prior, particles[i].θ) - logpdf(prior, old_particles[i].θ) +
               log(mean(particles[i].distances .< ϵ)/mean(old_particles[i].distances .< ϵ))
        if logw >= 0 || rand() < exp(logw)
            continue
        end
        copyto!(particles[i], old_particles[i]) # proposal rejected
    end
    particles
end

function run!(rng::Random.AbstractRNG, asmc::AdaptiveSMC, model, data;
              callback = () -> nothing,
              maxfevals = Inf,
              verbose = true)
    start_time = now()
    n_sims = Int[]
    epsilons = Float64[]
    i = 0
    eps = Inf
    K = length(asmc.particles)
    asmc.weights .= fill(1/K, K)
    push!(n_sims, sample_particles!(rng, asmc.particles, asmc.prior, model, data, asmc.weights, asmc.distance))
    update!(asmc.proposal, theta.(asmc.particles), asmc.weights, log_weights = false)
    callback()
    verbose && @printf "%10.s %13.s %8.s %8.s %8.s\n" "iteration" "elapsed" "fevals" "epsilon" "ess"
    verbose && @printf "%10.s %13.s %8.s %8.e %8.s\n" 0 round(now() - start_time, Second) n_sims[1] eps 1
    while true
        i += 1
        eps, ess = new_epsilon!(asmc.weights, eps, asmc.particles, asmc.α, asmc.ϵ)
        push!(epsilons, eps)
        ess < asmc.Nt && resample!(rng, asmc.particles, asmc.weights,
                                  log_weights = false, verbose = verbose)
        old_particles = deepcopy(asmc.particles)
        update!(asmc.proposal, theta.(asmc.particles), asmc.weights, log_weights = false)
        push!(n_sims, sample_particles!(rng, asmc.particles, asmc.proposal, model, data, asmc.weights, asmc.distance, indexed = true))
        mh_accept!(asmc.particles, old_particles, asmc.prior, eps, asmc.weights)
        ess = ess!(asmc.weights, asmc.particles, ones(K), eps)
        callback()
        verbose && @printf "%10.s %13.s %8.s %8.e %8.s\n" i round(now() - start_time, Second) sum(n_sims) eps ess/K
        (eps ≈ asmc.ϵ || sum(n_sims) > maxfevals) && break
    end
    (weights = asmc.weights,
     particles = asmc.particles,
     n_sims = n_sims,
     epsilons = epsilons)
end


text:
# Syntax → Graph

function bindings(ex)
  bs = []
  for ex in ex.args
    @capture(ex, x_ = _) && push!(bs, x)
  end
  return bs
end

function normreturn(ex)
  @capture(ex.args[end], return x_) && (ex.args[end] = x)
  ex
end

function normedges(ex)
  ex = normreturn(ex)
  map!(ex.args, ex.args) do ex
    isline(ex) ? ex :
    @capture(ex, _ = _) ? ex :
    :($(gensym("edge")) = $ex)
  end
  return ex
end

normblock(ex) = @> ex MacroTools.flatten block normedges normlines

function latenodes(exs)
  bindings = d()
  for ex in exs
    @capture(ex, b_Symbol = x_)
    bindings[b] = dvertex(x)
  end
  return bindings
end

graphm(bindings, node) = constant(node)
graphm(bindings, node::Vertex) = node
graphm(bindings, ex::Symbol) =
  haskey(bindings, ex) ? graphm(bindings, bindings[ex]) : constant(ex)

function tovertex!(v, bs, f, xs...)
  v.value = f
  thread!(v, map(ex -> graphm(bs, ex), xs)...)
end

function graphm!(v, bindings, ex::Expr)
  if isexpr(ex, :vertex)
    tovertex!(v, bindings, ex.args...)
  elseif @capture(ex, f_(args__))
    tovertex!(v, bindings, f, args...)
  else
    v.value = Constant(ex)
  end
  return v
end

graphm(bindings, ex::Expr) =
  isexpr(ex, :block) ? graphm(bindings, ex.args) :
    graphm!(dvertex(nothing), bindings, ex)

function fillnodes!(bindings, nodes)
  nodes′ = []
  for b in nodes
    node = bindings[b]
    if !(value(node) isa Expr)
      bindings[b] = graphm(bindings, value(node))
    else
      push!(nodes′, b)
    end
  end
  for b in nodes′
    node = bindings[b]
    graphm!(node, bindings, node.value)
  end
  return bindings
end

function graphm(bindings, exs::Vector)
  exs = normblock(:($(exs...);)).args
  @capture(exs[end], result_Symbol = _)
  lates = latenodes(exs)
  merge!(bindings, lates)
  fillnodes!(bindings, keys(lates))
  output = graphm(bindings, result)
end

bindargs(xs) = Dict{Any,Any}(x => inputnode(i) for (i, x) in enumerate(xs))

normalise(ex) =
  @> ex normops normcalls normsplits normclosures desugar

graphm(x; args = ()) = graphm(bindargs(args), normalise(x))


text:
"""
	grdproject(cmd0::String="", arg1=nothing, kwargs...)

Project a geographical gridded data set onto a rectangular grid or do the inverse projection.

Full option list at [`grdproject`]($(GMTdoc)grdproject.html)

Parameters
----------

- $(GMT.opt_J)
- **C** | **center** :: [Type => Str | []]      ``Arg = [dx/dy]``

    Let projected coordinates be relative to projection center [Default is relative to lower left corner].
    ($(GMTdoc)grdproject.html#c)
- **D** | **inc** :: [Type => Str | number]     ``Arg = xinc[unit][+e|n][/yinc[unit][+e|n]]``

    Set the grid spacing for the new grid. Append m for arc minute, s for arc second.
    ($(GMTdoc)grdproject.html#d)
- **E** | **dpi** :: [Type => Number]

    Set the resolution for the new grid in dots per inch.
    ($(GMTdoc)grdproject.html#e)
- **F** | **one2one** :: [Type => Str]           ``Arg = [c|i|p|e|f|k|M|n|u]``

    Force 1:1 scaling, i.e., output (or input, see -I) data are in actual projected meters [e].
    ($(GMTdoc)grdproject.html#f)
- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdproject(....) form.
    ($(GMTdoc)grdproject.html#g)
- **I** | **inverse** :: [Type => Bool]

    Do the Inverse transformation, from rectangular to geographical.
    ($(GMTdoc)grdproject.html#i)
- **M** | **projected_unit** :: [Type => Str]    ``Arg = c|i|p``

    Append c, i, or p to indicate that cm, inch, or point should be the projected measure unit.
    ($(GMTdoc)grdproject.html#m)
- $(GMT.opt_R)
- $(GMT.opt_V)
- $(GMT.opt_n)
- $(GMT.opt_r)
"""
function grdproject(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && return monolitic("grdproject", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:R :J :V_params :n :r])
	cmd = parse_these_opts(cmd, d, [[:C :center], [:D :inc], [:E :dpi], [:F :one2one],
	                                [:G :outgrid], [:I :inverse], [:M :projected_unit]])

	common_grd(d, cmd0, cmd, "grdproject ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdproject(arg1, cmd0::String=""; kw...) = grdproject(cmd0, arg1; kw...)


text:
@testset "wrap" begin
    @testset "unary" begin
        f(x) = x + 1
        _test_unary_op(wrap(f), b4, f)

        g(t, x) = t > _T_START + Day(2) ? x + 1 : x + 2
        _test_unary_op(wrap(g; time_agnostic=false), b4, g; time_agnostic=false)
    end

    @testset "binary" begin
        f(x, y) = x + y + 1
        _test_binary_op(wrap(f), f)

        g(t, x, y) = t > _T_START + Day(2) ? x + y + 1 : x + y + 2
        _test_binary_op(wrap(g; time_agnostic=false), g; time_agnostic=false)
    end

    @testset "ternary" begin
        f(x, y, z) = x + y + 1 - z
        _test_ternary_op(wrap(f), f)

        g(t, x, y, z) = t > _T_START + Day(2) ? x + y + 1 - z : x + y + 2 - z
        _test_ternary_op(wrap(g; time_agnostic=false), g; time_agnostic=false)
    end
end

@testset "wrapb" begin
    ba = _get_rand_vec_block(MersenneTwister(42), 3, 10)
    bb = _get_rand_vec_block(MersenneTwister(43), 3, 10)
    bc = _get_rand_vec_block(MersenneTwister(44), 3, 10)
    bd = _get_rand_vec_block(MersenneTwister(45), 3, 10)

    na = block_node(ba)
    nb = block_node(bb)
    nc = block_node(bc)
    nd = block_node(bd)

    @testset "unary" begin
        f(x) = x + 3
        f_broadcast(x) = f.(x)
        _test_unary_op(wrapb(f), ba, f_broadcast)

        g(t, x) = t > _T_START + Day(2) ? x + 3 : x + 4
        g_broadcast(t, x) = t > _T_START + Day(2) ? x .+ 3 : x .+ 4
        _test_unary_op(wrapb(g; time_agnostic=false), ba, g_broadcast; time_agnostic=false)
    end

    @testset "binary" begin
        f(x, y) = x + y + 3
        # TODO proper test here with alignment. Generalise _test_binary_op.
        n = wrapb(f)(na, nb)
        b = _eval(n)

        @test b.times == ba.times
        @test b.values == [f.(x, y) for (x, y) in zip(ba.values, bb.values)]
    end

    @testset "ternary" begin
        f(x, y, z) = x + y * z + 3
        # TODO proper test here with alignment.
        n = wrapb(f)(na, nb, nc)
        b = _eval(n)

        @test b.times == ba.times
        @test b.values ==
            [f.(x, y, z) for (x, y, z) in zip(ba.values, bb.values, bc.values)]
    end

    @testset "quaternary" begin
        f(x, y, z, w) = x + y * z + 3 - w
        # TODO proper test here with alignment.
        n = wrapb(f)(na, nb, nc, nd)
        b = _eval(n)

        @test b.times == ba.times
        #! format: off
        @test b.values == [
            f.(x, y, z, w)
            for (x, y, z, w) in zip(ba.values, bb.values, bc.values, bd.values)
        ]
        #! format: on
    end
end


text:
module MeasureTheory

using MLStyle

import Distributions
const Dists = Distributions

const EmptyNamedTuple = NamedTuple{(),Tuple{}}

export ≪
export sampletype

abstract type AbstractMeasure end

"""
    logdensity(μ::Measure{X}, x::X)

Compute the logdensity of the measure μ at the point x. This is the standard way
to define `logdensity` for a new measure. the base measure is implicit here, and
is understood to be `basemeasure(μ)`.

Methods for computing density relative to other measures will be
"""
function logdensity end

include("absolutecontinuity.jl")
include("basemeasures.jl")
include("macros.jl")
include("combinators/scale.jl")
include("combinators/superpose.jl")
include("combinators/product.jl")
include("combinators/power.jl")
include("distributions.jl")
include("probability/normal.jl")
include("density.jl")
include("pushforward.jl")

end # module


text:
__precompile__(true)

"""
The module for [Xorshift Family](@ref).

Provide 8 RNG types (others are to be deprecated):

- [`Xoroshiro64Star`](@ref)
- [`Xoroshiro64StarStar`](@ref)
- [`Xoroshiro128Plus`](@ref)
- [`Xoroshiro128StarStar`](@ref)
- [`Xoshiro128Plus`](@ref)
- [`Xoshiro128StarStar`](@ref)
- [`Xoshiro256Plus`](@ref)
- [`Xoshiro256StarStar`](@ref)
"""
module Xorshifts

include("common.jl")

include("xorshift64.jl")

include("xorshift128.jl")

include("xorshift1024.jl")

export Xoroshiro64Star, Xoroshiro64StarStar
include("xoroshiro64.jl")

export Xoroshiro128Plus, Xoroshiro128StarStar
include("xoroshiro128.jl")

export Xoshiro128Plus, Xoshiro128StarStar
include("xoshiro128.jl")

export Xoshiro256Plus, Xoshiro256StarStar
include("xoshiro256.jl")

include("docs.jl")

end


text:
abstract type AbstractSplineLayer <: Function end
Flux.trainable(m::AbstractSplineLayer) = (m.p,)

"""
Constructs a Spline Layer. At a high-level, it performs the following:
1. Takes as input a one-dimensional training dataset, a time span, a time step and
an interpolation method.
2. During training, adjusts the values of the function at multiples of the time-step
such that the curve interpolated through these points has minimum loss on the corresponding
one-dimensional dataset.

```julia
SplineLayer(time_span,time_step,spline_basis,saved_points=nothing)
```
Arguments:
- `time_span`: Tuple of real numbers corresponding to the time span.
- `time_step`: Real number corresponding to the time step.
- `spline_basis`: Interpolation method to be used yb the basis (current supported
  interpolation methods: ConstantInterpolation, LinearInterpolation, QuadraticInterpolation,
  QuadraticSpline, CubicSpline).
- 'saved_points': values of the function at multiples of the time step. Initialized by default
to a random vector sampled from the unit normal.
"""
struct SplineLayer{T<:Tuple{Real, Real},R<:Real,S1<:AbstractVector,S2<:UnionAll} <: AbstractSplineLayer
    time_span::T
    time_step::R
    saved_points::S1
    spline_basis::S2
    function SplineLayer(time_span,time_step,spline_basis,saved_points=nothing)
        saved_points = randn(length(time_span[1]:time_step:time_span[2]))
        new{typeof(time_span),typeof(time_step),typeof(saved_points),typeof(spline_basis)}(time_span,time_step,saved_points,spline_basis)
    end
end

function (layer::SplineLayer)(t::Real,p=layer.saved_points)
    return layer.spline_basis(p,layer.time_span[1]:layer.time_step:layer.time_span[2])(t)
end


text:
module ResettableStacks

  const FULL_RESET_COUNT = 10000

  import Base: isempty, length, push!, pop!, start, next, done
  type ResettableStack{T}
    data::Vector{T}
    cur::Int
    numResets::Int
  end

  ResettableStack{T}(ty::Type{T}) = ResettableStack{T}(Vector{T}(),0,0)

  isempty(S::ResettableStack) = S.cur==0
  length(S::ResettableStack)  = S.cur

  function push!(S::ResettableStack,x)
    if S.cur==length(S.data)
      S.cur+=1
      push!(S.data,x)
    else
      S.cur+=1
      S.data[S.cur]=x
    end
    nothing
  end

  function pop!(S::ResettableStack)
    if S.cur==length(S.data)
      S.cur-=1
      pop!(S.data)
    else
      S.cur-=1
      S.data[S.cur+1]
    end
  end

  start(S::ResettableStack) = s = 1
  function next(S::ResettableStack,s)
    s += 1
    S.data[s]
  end
  done(S::ResettableStack,s) = s==S.cur

  function reset!(S::ResettableStack)
    S.numResets += 1
    S.cur = 0
    if S.numResets%FULL_RESET_COUNT==0
      S.data = Vector{eltype(S.data)}()
    end
    nothing
  end

  export ResettableStack, reset!
end # module


text:
#=
# 3D Poisson, Dirichlet bc
# Uses Dendro imported as a custom gen target
=#
if !@isdefined(Femshop)
    include("../Femshop.jl");
    using .Femshop
end
init_femshop("poisson3dcustom");

useLog("poisson3dcustomlog", level=3)

# default values (max_depth=6, wavelet_tol = 0.1, partition_tol = 0.3, solve_tol = 1e-6, max_iters = 100)
generateFor("target_dendro_cg.jl", params=(5, 1, 0.3, 0.000001, 100))

domain(3)
functionSpace(order=1)

u = variable("u")
testSymbol("v")

boundary(u, 1, DIRICHLET, 0)

# Write the weak form
coefficient("f", "-14*pi*pi*sin(3*pi*x)*sin(2*pi*y)*sin(pi*z)")
weakForm(u, "-dot(grad(u),grad(v)) - f*v")

solve(u);

finalize_femshop()


text:
using Turing, Measurements

"""
    t_model(q::Array{Float64, 1}, R::Array{Measurement, 1}, array::Array{Any, 2}, scale::Any, bkg::Any, resolution::Any)

The Turing model for the sampling of the reflectometry profile.

### Parameters
- `q::Array{Float64}` : the q-vector values.
- `R::Array{Measurement}` : the reflected intensity values, including uncertainty.
- `array::Array{Any, 2}` : the layers from which the reflectometry should be found.
- `scale::Any` : the scale factor.
- `bkg::Any` : the background level.
- `resolution::Any` : the constant resolution width.
"""
@model t_model(q::Array{Float64, 1}, R::Array{Measurement, 1}, layers, scale, bkg, resolution) = begin
    layers_inp = Array{Any}(undef, size(layers))
    for i in 1:size(layers, 1)
        for j in 1:size(layers, 2)
            if isa(layers[i, j].value, AbstractFloat)
                layers_inp[i, j] = layers[i, j].value
            else
                layers_inp[i, j] ~ layers[i, j].value
            end
        end
    end
    if isa(scale.value, AbstractFloat)
        scale_inp = scale.value
    else
        scale_inp ~ scale.value
    end
    if isa(bkg.value, AbstractFloat)
        bkg_inp = bkg.value
    else
        bkg_inp ~ bkg.value
    end
    if isa(resolution.value, AbstractFloat)
        resolution_inp = resolution.value
    else
        resolution_inp ~ resolution.value
    end

    r = Array{Any}(undef, length(q))
    dr = zeros(length(q))
    for i in 1:length(q)
        r[i] = Measurements.value(R[i])
        dr[i] = Measurements.uncertainty(R[i])
    end
    mu = Jeff.constant_smearing(q, layers_inp, resolution_inp, scale_inp, bkg_inp)
    for i in 1:length(r)
        r[i] ~ Normal(mu[i], dr[i])
    end
end


text:
"""
    CrossingPhase
Crossing phases for SidewalkPedestrianModel.

For a crossing pedestrian, phases go: -2, -1, 0, 1
For a non-crossing pedestrian that pauses at the crosswalk, phases go: -2, -1, 1
For a non-crossing pedestrian that does not pause at the crosswalk, phases go: -2, 1

Model based on Feliciani et al (2017) - A simulation model for non-signalized pedestrian
crosswalks based on evidence from on field observation.
"""
struct CrossingPhase
    phase::Int # -2, -1, 0, 1
end

const APPROACHING = -2 # Keep walking at constant speed
const APPRAISING = -1 # Decelerate until reach curbside
const CROSSING = 0 # Accelerate until desired speed
const LEAVING = 1 # Keep walking at constant speed

"""
    SidewalkPedestrianModel

Walks along the sidewalk until approaching the crosswalk. Waits for the cars to pass, then crosses.

# Fields
- `timestep::Float64`
- `phase::Int = APPROACHING`
- `ttc_threshold::Float64 = clamp(rand(Normal(4.0, 2.5)), 1.0, Inf)`
- `crosswalk::Lane = Lane()`
- `sw_origin::Lane = Lane()`
- `sw_dest::Lane = Lane()`
- `a::PedestrianLatLonAccel = PedestrianLatLonAccel(0.0, 0.0, sw_origin)` makes you turn, left/right
- `σ::Float64 = NaN` optional stdev on top of the model, set to zero or NaN for deterministic
- `v_des_approach::Float64 = clamp(rand(Normal(1.28, 0.18)), 0.0, Inf)` Based on Feliciani et al. results
- `v_des_appraise::Float64 = clamp(rand(Normal(0.94, 0.21)), 0.0, Inf)`
- `v_des_cross::Float64 = clamp(rand(Normal(1.35, 0.18)), 0.0, Inf)`
- `ped_accel::Float64 = 0.30`
- `ped_decel::Float64 = -0.50`
"""
@with_kw mutable struct SidewalkPedestrianModel <: DriverModel{PedestrianLatLonAccel}
    timestep::Float64
    phase::Int = APPROACHING

    ttc_threshold::Float64 = clamp(rand(Normal(4.0, 2.5)), 1.0, Inf)
    crosswalk::Lane = Lane()
    sw_origin::Lane = Lane()
    sw_dest::Lane = Lane()
    # conflict_lanes::Vector{Lane} = Lane[]

    # Values from Feliciani
    a::PedestrianLatLonAccel = PedestrianLatLonAccel(0.0, 0.0, sw_origin) # makes you turn, left/right
    σ::Float64 = NaN # optional stdev on top of the model, set to zero or NaN for deterministic behavior
    v_des_approach::Float64 = clamp(rand(Normal(1.28, 0.18)), 0.0, Inf) # Based on Feliciani et al. results
    v_des_appraise::Float64 = clamp(rand(Normal(0.94, 0.21)), 0.0, Inf)
    v_des_cross::Float64 = clamp(rand(Normal(1.35, 0.18)), 0.0, Inf)

    ped_accel::Float64 = 0.30
    ped_decel::Float64 = -0.50

    # For graphing
    pos_x = Float64[]
    pos_y = Float64[]
    velocities = Float64[]
    phases = Int[]
end

Base.rand(rng::AbstractRNG, model::SidewalkPedestrianModel) = model.a

function AutomotiveDrivingModels.observe!(model::SidewalkPedestrianModel, scene::Scene{Entity{VehicleState, D, I}}, roadway::Roadway, egoid::I) where {D, I}
    ped = scene[findfirst(egoid, scene)]

    push!(model.pos_x, posg(ped.state).x)
    push!(model.pos_y, posg(ped.state).y)
    push!(model.velocities, vel(ped.state))
    push!(model.phases, model.phase)

    if model.phase == APPROACHING
        update_approaching(ped, model, roadway, model.crosswalk, model.sw_origin)
    elseif model.phase == APPRAISING
        update_appraising(ped, scene, model, roadway, model.crosswalk, model.sw_origin)
    elseif model.phase == CROSSING
        update_crossing(ped, scene, model, roadway, model.crosswalk)
    elseif model.phase == LEAVING
        update_leaving(ped, scene, model, roadway, model.sw_dest)
    else
        # Model phase is out of bounds
    end
end

function update_approaching(ped::Entity{S, D, I},
                            model::SidewalkPedestrianModel,
                            roadway::Roadway,
                            crosswalk::Lane,
                            sw_origin::Lane
                            ) where {S, D, I}
    # Keep walking at constant speed until close to crosswalk
    # If waiting to cross, go to appraising. Else, go to leaving phase
    Δs = get_distance_to_lane(ped, crosswalk, roadway)
    if Δs < 3.0
        model.phase = APPRAISING
    end
end

function update_appraising(ped::Entity{S, D, I},
                           scene::Scene,
                           model::SidewalkPedestrianModel,
                           roadway::Roadway,
                           crosswalk::Lane,
                           sw_origin::Lane,
                           ) where {S, D<:AbstractAgentDefinition, I}

    Δs = get_distance_to_lane(ped, crosswalk, roadway)
    if Δs < 0.1
        set_inst_vel(ped, model, sw_origin, 0.0)

        for veh in scene
            if class(veh.def) == AgentClass.PEDESTRIAN
                continue
            end

            # Check vehicle distance to the crosswalk
            posG = posg(veh.state)
            cw_posF = Frenet(posG, model.crosswalk, roadway)
            lane = get_lane(roadway, veh)
            Δs = get_distance_to_lane(veh, model.crosswalk, roadway)
            ttc = Δs/vel(veh.state)

            if 0.0 < ttc < model.ttc_threshold
                return # Exit if any car is not within threshold
            end
        end
        model.phase = CROSSING
    else
        # Decelerate to appraise speed
        set_accel(ped, model, model.v_des_appraise, sw_origin)
    end
end

function update_crossing(ped::Entity{S, D, I},
                         scene::Scene,
                         model::SidewalkPedestrianModel,
                         roadway::Roadway,
                         crosswalk::Lane
                         ) where {S, D, I}
    # Accelerate until desired speed
    set_accel(ped, model, model.v_des_cross, crosswalk)
    # If ped is at end of crosswalk, transition to leaving phase.
    Δs = get_distance_to_lane(ped, model.sw_dest, roadway)

    if Δs < 0.2
        set_inst_vel(ped, model, model.crosswalk, 0.0)
        model.phase = LEAVING
    end
end

function update_leaving(ped::Entity{S, D, I},
                        scene::Scene,
                        model::SidewalkPedestrianModel,
                        roadway::Roadway,
                        sw_dest::Lane
                        ) where {S, D, I}
    # Keep walking at constant speed
    set_accel(ped, model, model.v_des_approach, model.sw_dest)
end

function set_inst_vel(ped::Entity{VehicleState, D, I}, model::SidewalkPedestrianModel, lane_des::Lane, v_des::Float64) where {D, I}
    model.a = PedestrianLatLonAccel(0.0, calc_a(vel(ped.state), v_des, model.timestep), lane_des)
end

function calc_a(v_init::Float64, v::Float64, Δt::Float64)
    return (v - v_init)/Δt
end

function set_accel(ped::Entity{VehicleState, D, I}, model::SidewalkPedestrianModel, v_des::Float64, lane_des::Lane) where {D, I}
    if vel(ped.state) < v_des
        model.a = PedestrianLatLonAccel(0.0, model.ped_accel, lane_des)
    elseif vel(ped.state) > v_des
        model.a = PedestrianLatLonAccel(0.0, model.ped_decel, lane_des)
    else
        model.a = PedestrianLatLonAccel(0.0, 0.0, lane_des)
    end
end

# Edited from AutomotivePOMDPs.jl - projects collision point onto the lane.
# Should work in any orientation.
function get_distance_to_lane(veh::Entity{VehicleState, D, I}, lane_des::Lane, roadway::Roadway) where {D, I}
    vehlane = get_lane(roadway, veh)
    lane_segment = LineSegment(lane_des.curve[1].pos, lane_des.curve[end].pos)
    collision_point = intersect(Ray(posg(veh.state)), lane_segment)
    lane_des_proj = Frenet(VecSE2(collision_point), vehlane, roadway)
    Δs = lane_des_proj.s - posf(veh.state).s
    return Δs
end


text:
export replacer_channel, depolarizing_channel, erasure_channel

"""
    replacer_channel(
        ρ :: State,
        σ :: State,
        μ :: Real
    ) :: State

The replacer channel replaces the quantum state `ρ` with quantum state `σ` with
probability `μ`.
The replacer channel is expressed:

```math
    \\mathcal{R}_{\\mu}(\\rho) = \\mu \\rho + (1-\\mu)\\sigma
```

A `DomainError` is thrown if
* `μ` does not satisfy `1 ≥ μ ≥ 0`
* `ρ` and `σ` are not the same size
"""
function replacer_channel(ρ :: State, σ :: State, μ :: Real) :: State
    if !(1 ≥ μ ≥ 0)
        throw(DomainError(μ, "Input `μ` must satisfy `1 ≥ μ ≥ 0`"))
    elseif size(ρ) != size(σ)
        throw(DomainError((ρ, σ), "Inputs `ρ` and `σ` must have the same size."))
    end

    State(μ*ρ + (1-μ)*σ)
end

"""
    depolarizing_channel( ρ :: State, μ :: Real ) :: State

The depolarizing channel mixes uniform classical noise into a quantum state `ρ`.
The argument `μ` describes the amount of noise mixed into the quantum states.
For a quantum state ``\\rho``, the depolarizing channel is expressed:

```math
    \\mathcal{D}_{\\mu}(\\rho) = \\mu \\rho + \\frac{(1 - \\mu)}{d}\\mathbb{I}_d
```

A `DomainError` is thrown if `μ` does not satisfy `1 ≥ μ ≥ 0`.
"""
function depolarizing_channel(ρ :: State, μ :: Real) :: State
    d = size(ρ, 1)
    σ = State(Matrix{Int64}(I, (d,d))/d)
    replacer_channel(ρ, σ, μ)
end

"""
    erasure_channel( ρ :: State, μ :: Real ) :: State

The erasure channel mixes a quantum state `ρ` with an error flag ``|F\\rangle``
orthogonal to the Hilbert space of `ρ`.
The argument `μ` describes the probability that `ρ` is replaced with the error flag.
For a quantum state ``\\rho``, the erasure channel is expressed:

```math
    \\mathcal{E}_{\\mu}(\\rho) = \\mu \\rho + (1 - \\mu) |F \\rangle \\langle F|
```

Note that the erasure channel increases the dimension of the Hilbert space by 1.

A `DomainError` is thrown if `μ` does not satisfy `1 ≥ μ ≥ 0`.
"""
function erasure_channel(ρ :: State, μ :: Real) :: State
    d = size(ρ, 1)

    ρ_ext = zeros(Complex{Float64}, d+1, d+1)
    ρ_ext[1:d,1:d] = ρ

    σ = zeros(Complex{Float64}, d+1, d+1)
    σ[d+1,d+1] = 1

    replacer_channel(
        State(ρ_ext),
        State(σ),
        μ
    )
end


text:
#random_sparsity.jl

# include("interface_tape.jl")
# using TapeInterface

using JuMP
using Ipopt

N = parse(Int,ARGS[1])
K = parse(Int,ARGS[2])
RandSet= round(Int,readdlm(string("randset",N,"_",K,".txt")))

if parse(Int, ARGS[3]) == 1
    m = Model(solver=IpoptSolver())
elseif parse(Int, ARGS[3]) == 2
    include("interface_tape.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()))
elseif parse(Int, ARGS[3]) == 3
    include("interface_tape_dict.jl")
    using TapeInterface
    m = Model(solver=TapeSolver(IpoptSolver()))
else
    @assert false
end

# @show RandSet

@variable(m, -100<=x[1:N]<=100)

@NLobjective(m, Min, sum{ (x[i]-1)^2 + prod{x[j], j=RandSet[i,:]}, i=1:N} )

solve(m)


text:
mutable struct State{T} <: FieldVector{6, T}
    Qu :: T
    Qv :: T
    Qθ :: T
    Qs :: T
    Qb :: T
     h :: T
    function State(T=Float64)
        new{T}(0, 0, 0, 0, 0, 0)
    end
end

"""
    update_state!(model)

Update the top flux conditions and mixing depth for `model`
and store in `model.state`.
"""
function update_state!(m)
    m.state.Fu = getbc(m, m.bcs.U.top)
    m.state.Fv = getbc(m, m.bcs.V.top)
    m.state.Fθ = getbc(m, m.bcs.T.top)
    m.state.Fs = getbc(m, m.bcs.S.top)
    m.state.Fb = m.constants.g * (m.constants.α * m.state.Fθ - m.constants.β * m.state.Fs)
    m.state.h  = ModularKPP.mixing_depth(m)
    return nothing
end


text:
module CleanCMB
using LinearAlgebra
export ilc_weights, ilc_clean_cij
export cilc_weights
export milca_weights
export tsz, dust1, synch
export loglike_beta
include("ilc.jl")
include("cilc.jl")
include("milca.jl")
include("libspectra.jl")
include("profile_likelihood.jl")
end


text:
@testset "Library Functions" begin
    # Standard component
    @test haskey(BondGraphs.DEFAULT_LIBRARY, :C)
    # Biochemical component
    @test haskey(BondGraphs.DEFAULT_LIBRARY, :Re)
    # Non-existent component
    @test !haskey(BondGraphs.DEFAULT_LIBRARY, :A)

    lib = Dict(:A => Dict())
    addlibrary!(lib)
    # Component now exists
    @test haskey(BondGraphs.DEFAULT_LIBRARY, :A)

    # Delete fake component for later tests
    delete!(BondGraphs.DEFAULT_LIBRARY, :A)
end

@testset "Graph Node Colour Selection" begin
    nodes = [Component(:C), Component(:Re), Component(:SS), EqualFlow()]
    @test BondGraphs.nodecolours(nodes) == [1, 2, 3, :lightgray]
end

@testset "Plotting" begin
    # See conversion_tests.jl
    rn = @reaction_network MM_reversible begin
        (1, 1), E + S <--> C
        (1, 1), C <--> E + P
    end
    bg = BondGraph(rn; chemostats=["S", "P"])

    rec = RecipesBase.apply_recipe(Dict{Symbol, Any}(), bg)
    attributes = getfield(rec[1], 1)

    @test attributes[:curves] == false
    @test attributes[:title] == :MM_reversible
    @test attributes[:nodeshape] == :rect
end

max_stars_repo_name: mic-pan/BondGraphs

max_issues_repo_name: mic-pan/BondGraphs

max_forks_repo_name: mic-pan/BondGraphs


text:
function preparevis!(storage::Storage{T,N}, id, shape, animation, shapevisualizer, framevisualizer, showshape, showframes) where {T,N}
    if showshape
        for i=1:N
            shapecomposition = compose(Translation(storage.x[id][i] + vrotate(shape.xoffset, storage.q[id][i])),LinearMap(storage.q[id][i] * shape.qoffset))
            atframe(animation, i) do
                settransform!(shapevisualizer, shapecomposition)
            end
        end
    end

    if showframes
        for i=1:N
            framecomposition = compose(Translation(storage.x[id][i]),LinearMap(storage.q[id][i]))
            atframe(animation, i) do
                settransform!(framevisualizer, framecomposition)
            end
        end
    end

    return
end

function visualize(mechanism::AbstractMechanism, storage::Storage{T,N}; env::String = "blink", showframes::Bool = false) where {T,N}
    storage = deepcopy(storage)
    bodies = mechanism.bodies
    origin = mechanism.origin
    if showframes
        triads = [Triad(0.33) for i=1:length(bodies)]
    end

    vis = Visualizer()
    setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0))
    setprop!(vis["/Background"], "bottom_color", RGBA(0.5, 0.6, 0.7))
    # Somehow delete! doesn't work in a function call, so set axes to not visible for now
    setvisible!(vis["/Axes"],false)
    # delete!(vis["/Axes"])

    if env == "blink"
        open(vis, Blink.Window())
    elseif env == "browser"
        open(vis)
    elseif env == "editor"
        # Needs to be returned to show. Is there a better way of doing this?
    else
        @error "No valid visualization environment option selected."
    end

    framerate = Int64(round(1/mechanism.Δt))
    animation = Animation(Dict{MeshCat.SceneTrees.Path,MeshCat.AnimationClip}(), framerate)

    for (id,body) in enumerate(bodies)
        shape = body.shape
        visshape = convertshape(shape)
        subvisshape = nothing
        subvisframe = nothing
        showshape = false
        if visshape !== nothing
            subvisshape = vis["bodies/body:"*string(id)]
            setobject!(subvisshape, visshape, MeshPhongMaterial(color=shape.color))
            showshape = true
        end
        if showframes
            subvisframe = vis["frames/body:"*string(id)]
            setobject!(subvisframe, triads[id])
        end
        preparevis!(storage, id, shape, animation, subvisshape, subvisframe, showshape, showframes)
    end

    id = origin.id
    shape = origin.shape
    visshape = convertshape(shape)
    if visshape !== nothing
        subvisshape = vis["bodies/origin:"*string(id)]
        setobject!(subvisshape, visshape, MeshPhongMaterial(color=shape.color))
        composition = compose(Translation(shape.xoffset),LinearMap(shape.qoffset))
        settransform!(subvisshape, composition)
    end
    if showframes
        subvisframe = vis["frames/origin:"*string(id)]
        setobject!(subvisframe, Triad(0.5))
    end

    setanimation!(vis, animation)
    env == "editor" ? (return render(vis)) : (return vis)
end


text:
# maintained by @zsunberg
# an empty belief
# for use with e.g. a random policy
mutable struct VoidUpdater <: Updater end

initialize_belief(::VoidUpdater, ::Any) = nothing
initialize_belief(::VoidUpdater, ::Any, ::Any) = nothing
create_belief(::VoidUpdater) = nothing

update{B}(::VoidUpdater, ::B, ::Any, ::Any, b=nothing) = nothing


text:
using TravelingSalesmanBenchmarks
using Test

@testset "TravelingSalesmanBenchmarks.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: ericphanson/TravelingSalesmanBenchmarks


text:
"""
    fr = metagauss(f)

Generate a metagaussian field multiplying the spectralfield `f` with random phases and performing an inverse FFT transform to real space.
"""
	function metagauss(f);

        (ns,ns,nt)=size(f);
#	phases as fft of a gaussian noise random field
#        ph=zeros(ns,ns,nt,'single');

        ph=zeros(ns,ns,nt);
        for i=1:nt; ph[:,:,i]=(randn(ns,ns)); end;
	ph=fft(ph);
	ph=ph./abs.(ph); ph[1,1,1]=0;
	ph=f.*ph;
	ph=real(ifft(ph));
end


text:
function random_walk_telegraph(data,T,Δt,domain,λ₁,λ₂,rate,state)
	states = []
	push!(states,state)

	Data = []
	push!(Data,T,Δt,rate,λ₁,λ₂)

	λ₁Δt = λ₁*Δt
	λ₂Δt = λ₂*Δt
	time = range(Δt,T,step=Δt)

	rΔt = Δt*rate
	s() = stepping(rΔt)

	p₀() = rand(1:domain)

	Walk = []
	w = Walker2D(p₀(),p₀())
	push!(Walk,w)

	#d(w) = taxicabdistance(DisplacementPositions(w[1],w[end]),domain)

	for i in 2:size(time,1)

		 # Create states
		 state = update(state,λ₁Δt,λ₂Δt)
		 push!(states,state)

		 if states[i-1] == 1 && states[i] == 1
			  continue
		 elseif states[i-1] == 1 && states[i] == 0
			  # first step
			  w = updateperiod(w, s(),domain,domain)
			  push!(Walk,w)
			  i₁ = i
			  t₁ = time[i]
			  push!(Data,i₁,t₁)
		 elseif states[i-1] == 0 && states[i] == 1
			  # last step
			  iₙ = i-1
			  tₙ = time[i-1]
			  W = Walk
			  w_tup = [(i.x,i.y) for i in W]
			  Walk = []
			  #dist = d(W).d
			  push!(Data,iₙ,tₙ,w_tup)
			  if size(Data,1) != 10
				   error("Length to push to dataframe is $(size(Data,1)) not 11")
			  end
			  push!(data,Data)
			  Data = []
			  push!(Data,T,Δt,rate,λ₁,λ₂)
		 elseif states[i-1] == 0 && states[i] == 0
			  if size(states,1) == 2
				   # keep walking
				   i₁ = i-1
				   t₁ = time[i-1]
				   push!(Data,i₁,t₁)
				   w = updateperiod(w, s(),domain,domain)
				   push!(Walk,w)
			  elseif size(states,1) > 2
				   w = updateperiod(w, s(),domain,domain)
				   push!(Walk,w)
			  else
				   error("States is either size 0 or size negative")
			  end
		 else
			  error("What state are you meant to be in?")
		 end
	end
	return data
end

function random_walk_telegraph_no_df(T,Δt,domain,λ₁,λ₂,rate,state)
	states = []
	push!(states,state)

	λ₁Δt = λ₁*Δt
	λ₂Δt = λ₂*Δt
	time = range(Δt,T,step=Δt)

	rΔt = Δt*rate
	s() = stepping(rΔt)

	p₀() = rand(1:domain)

	Walk = []
	W_tup = []
	w = Walker2D(p₀(),p₀())
	push!(Walk,w)

	#d(w) = taxicabdistance(DisplacementPositions(w[1],w[end]),domain)

	for i in 2:size(time,1)

		 # Create states
		 state = update(state,λ₁Δt,λ₂Δt)
		 push!(states,state)

		 if states[i-1] == 1 && states[i] == 1
			  continue
		 elseif states[i-1] == 1 && states[i] == 0
			  # first step
			  w = updateperiod(w, s(),domain,domain)
			  push!(Walk,w)

		 elseif states[i-1] == 0 && states[i] == 1

			  W = Walk
			  w_tup = [(i.x,i.y) for i in W]
			  push!(W_tup,w_tup)
			  #dist = d(W).d

		 elseif states[i-1] == 0 && states[i] == 0
			  if size(states,1) == 2

				   w = updateperiod(w, s(),domain,domain)
				   push!(Walk,w)
			  elseif size(states,1) > 2
				   w = updateperiod(w, s(),domain,domain)
				   push!(Walk,w)
			  else
				   error("States is either size 0 or size negative")
			  end
		 else
			  error("What state are you meant to be in?")
		 end
	end
	distance = [dₜ(w,domain) for w in W_tup]
    d̂ = mean(distance[2:end])

	return d̂
end

function tele_rw_inf_domain(
	λ₁::Float64,
	λ₂::Float64,
	r::Float64,
	Δt::Float64,
	T::Int,
	s::Int
	)

	# Set up λ
	λ₁Δt = λ₁*Δt
	λ₂Δt = λ₂*Δt

	# Set up r
	rΔt = r*Δt

	# Set up time
	time = range(Δt,T,step=Δt)

	# Push s to state
	state = []
	push!(state,s)

	## Set up Random walk
	st() = stepping(rΔt)
	p₀() = 0 #rand(1:domain)
	Walk = []
	w = Walker2D(p₀(),p₀())
	push!(Walk,w)

	# Time evolution
	for i in 2:size(time,1)
	    s = update(s,λ₁Δt,λ₂Δt)
	    push!(state,s)

	    if s == 0
	    	w = RandomWalker.update(w,st())
	    	push!(Walk,w)
	    elseif s == 1
	    	continue
	    else
	    	error("What state are you in? You are in state $(s)")
	    end
	end

	[[state] [Walk]]
end


text:
@testset "PathTracker" begin
    @testset "Constructor + Options" begin
        @polyvar x y z
        f = [x^2 - 2, x + y - 1]
        tracker, starts = pathtracker_startsolutions(f; system = FPSystem)
        @test tracker isa PathTracker{Vector{ComplexF64}}
        test_show_juno(tracker)
        @test !isempty(sprint(show, tracker))
        test_show_juno(tracker.options)
        @test !isempty(sprint(show, tracker.options))
        test_show_juno(tracker.state)
        @test !isempty(sprint(show, tracker.state))
        @test Base.broadcastable(tracker) isa Ref{typeof(tracker)}
        tracker = pathtracker(f; system = FPSystem, precision_strategy = :adaptive_never)
        @test tracker.options.precision_strategy == HC.PREC_STRATEGY_NEVER
        tracker = pathtracker(f; system = FPSystem, precision_strategy = :adaptive_always)
        @test tracker.options.precision_strategy == HC.PREC_STRATEGY_ALWAYS
        @test_throws ArgumentError pathtracker(
            f;
            system = FPSystem,
            precision_strategy = :NOT_AN_OPTION,
        )

        @test_throws ArgumentError pathtracker(
            f;
            not_an_option = true,
            this_option = "nope",
        )

        @test_throws ArgumentError pathtracker(f, patch = OrthogonalPatch())
        @test_throws ArgumentError pathtracker(f, system_scaling = :lalala)
    end

    # Simple
    @testset "Simple" begin
        @polyvar x y z
        f = [x^2 - 2, x + y - 1]
        tracker, starts = pathtracker_startsolutions(f; system = FPSystem)

        for xᵢ in starts
            @test is_success(track!(tracker, xᵢ))
            s = solution(tracker)
            @test abs(s[1]) ≈ sqrt(2) atol = 1e-8
            @test sum(s) ≈ 1 atol = 1e-8
            @test isnothing(winding_number(tracker))
            @test !isnan(tracker.state.solution_cond)
            @test !isnan(tracker.state.solution_accuracy)
            @test tracker.state.solution_accuracy < 1e-14
            @test tracker.state.solution_residual < 1e-14
        end

        # test iterator
        init!(tracker, first(starts))
        for _ in tracker
        end
        @test is_success(status(tracker))

        tracker, starts =
            pathtracker_startsolutions(f; projective_tracking = true, system = FPSystem)
        @test tracker isa PathTracker{PVector{ComplexF64,1}}
        for x in starts
            @test is_success(track!(tracker, x))
            s = solution(tracker)
            @test s isa Vector{ComplexF64}
            @test abs(s[1]) ≈ sqrt(2) atol = 1e-8
            @test sum(s) ≈ 1 atol = 1e-8
            @test isnothing(winding_number(tracker))
            @test !isnan(tracker.state.solution_cond)
            @test !isnan(tracker.state.solution_accuracy)
            @test !isnan(tracker.state.solution_residual)
        end

        tracker, starts = pathtracker_startsolutions(homogenize(f, z); system = FPSystem)
        @test tracker isa PathTracker{PVector{ComplexF64,1}}
        for x in starts
            @test is_success(track!(tracker, x))
            s = solution(tracker)
            @test s isa PVector{ComplexF64,1}
        end

        g = equations(katsura(5))
        tracker, starts = pathtracker_startsolutions(g; system = FPSystem)
        @test all(s -> is_success(track!(tracker, s)), starts)

        # start solution is also target solution
        @polyvar x
        tracker = pathtracker([x - 1]; system = FPSystem)
        @test is_success(track(tracker, [1]))
    end

    # Univariate singular
    @testset "Univariate singular" begin
        @polyvar x
        print("Univariate singular: ")
        for n = 2:12
            print(n, ",")
            f = [(x - 3)^n]
            g = [x^n - 1]
            S = [[cis(i * 2π / n)] for i = 0:(n-1)]

            tracker = pathtracker(g, f, S; seed = 842121, system = FPSystem)
            @test all(S) do s
                is_success(track!(tracker, s)) && (
                    winding_number(tracker) == n || isnothing(winding_number(tracker))
                ) && isapprox(solution(tracker)[1], 3.0; atol = 1e-6) &&
                     !isnan(tracker.state.solution_cond) &&
                     tracker.state.solution_accuracy < 1e-5
            end
        end
        println("")
    end

    # Wilkinson
    @testset "Wilkinson" begin
        @polyvar x
        print("Wilkinson: ")
        for n = 1:18
            print(n, ",")
            g = [x^n - 1]
            f = [prod(x - i for i = 1:n)]
            S = [[cis(i * 2π / n)] for i = 0:(n-1)]
            tracker = pathtracker(g, f, S; seed = 512321, system = FPSystem)
            @test all(s -> is_success(track!(tracker, s)), S)
            # check that we actually get the correct results
            @test sort(round.(
                Int,
                real.(first.(solution.(track.(tracker, S)))),
            )) == collect(1:n)
        end
        println("")
    end

    # Bivariate singular
    @testset "Bivariate singular" begin
        f = equations(griewank_osborne())
        tracker, starts = pathtracker_startsolutions(f; seed = 78373, system = FPSystem)
        @test all(starts) do s
            retcode = track!(tracker, s)
            (is_success(retcode) && winding_number(tracker) == 3) || is_at_infinity(retcode)
        end

        # Bivariate projetive singular
        @polyvar x z y
        # This has two roots of multiplicity 6 at the hyperplane z=0.
        # But the winding numbers are only 3 at each singularity
        F = [
            0.75 * x^4 + 1.5 * x^2 * y^2 - 2.5 * x^2 * z^2 + 0.75 * y^4 -
                2.5 * y^2 * z^2 + 0.75 * z^4,
            10 * x^2 * z + 10 * y^2 * z - 6 * z^3,
        ]
        @test pathtracker(F) isa PathTracker{PVector{ComplexF64,1}}
        tracker, starts = pathtracker_startsolutions(F; seed = 29831, system = FPSystem)
        @test all(s -> is_success(track!(tracker, s)), starts)
        @test all(starts) do s
            is_success(track!(tracker, s)) &&
            winding_number(tracker) == 3 && tracker.state.solution_cond > 1e10
        end
        # affine
        F̂ = subs.(F, Ref(y => 1))
        tracker, starts = pathtracker_startsolutions(F̂; seed = 29831, system = FPSystem)
        @test all(starts) do s
            is_success(track!(tracker, s)) &&
            winding_number(tracker) == 3 && tracker.state.solution_cond > 1e10
        end

        # Fix z == 1 --> all solutions at infinity
        G = subs.(F, Ref(z => 1))
        tracker, starts = pathtracker_startsolutions(G; seed = 29831, system = FPSystem)
        @test all(starts) do s
            is_at_infinity(track!(tracker, s))
        end
    end

    @testset "Non-Singular bad-conditioned" begin
        f = equations(PolynomialTestSystems.bacillus_subtilis())
        tracker, starts = pathtracker_startsolutions(f; seed = 78373, system = FPSystem)
        @test count(s -> is_success(track!(tracker, s)), starts) == 44
    end

    @testset "Diverging path without ill-conditioning" begin
        f = equations(cyclic(7))
        tracker, starts = pathtracker_startsolutions(f; seed = 78373, system = FPSystem)
        S = collect(starts)
        @test is_at_infinity(track!(tracker, S[7]))
        @test is_at_infinity(track!(tracker, S[9]))
        @test is_at_infinity(track!(tracker, S[5016]))
    end

    @testset "Cyclic 7 correct root count" begin
        f = equations(cyclic(7))
        tracker, starts = pathtracker_startsolutions(f; seed = 78373, system = FPSystem)
        results = map(s -> track!(tracker, s), starts)
        @test count(is_success, results) == 924

        successfull_paths =
            [1 4 6 10 13 14 16 18 43 44 46 59 63 66 69 70 73 78 85 93 94 98 104 106 111 119 121 122 134 141 148 152 155 157 160 163 171 173 182 185 189 200 204 206 208 211 215 218 225 227 230 231 236 244 245 246 256 258 259 261 264 272 277 280 284 287 298 305 318 319 326 336 340 342 347 358 367 380 384 397 400 405 410 411 416 420 423 424 430 431 439 442 447 470 476 483 486 497 510 511 550 552 558 566 569 575 583 585 587 589 591 599 601 602 605 610 615 622 627 630 634 638 639 645 650 651 653 656 658 669 673 679 681 683 689 691 693 696 707 709 715 719 722 723 724 733 746 755 757 765 767 770 773 786 796 805 807 808 811 816 818 832 837 840 846 850 857 858 864 873 874 875 880 882 900 905 909 910 911 918 919 920 923 925 929 935 936 942 943 947 955 960 972 979 982 990 1008 1010 1022 1031 1035 1037 1047 1059 1060 1062 1070 1075 1076 1082 1092 1096 1100 1102 1103 1109 1111 1115 1121 1124 1129 1147 1148 1152 1166 1168 1178 1189 1197 1198 1200 1201 1215 1216 1217 1223 1229 1235 1237 1242 1247 1251 1258 1261 1262 1272 1273 1277 1282 1293 1313 1318 1333 1336 1352 1362 1376 1377 1406 1408 1409 1419 1421 1422 1423 1424 1425 1442 1444 1446 1453 1454 1460 1462 1465 1473 1490 1503 1509 1520 1522 1523 1526 1527 1530 1531 1534 1549 1552 1558 1559 1561 1565 1579 1585 1594 1600 1621 1629 1631 1637 1640 1642 1644 1645 1649 1655 1656 1662 1665 1669 1671 1679 1705 1706 1712 1731 1734 1736 1737 1749 1750 1751 1753 1757 1759 1769 1784 1790 1795 1799 1800 1801 1802 1803 1807 1809 1817 1819 1825 1832 1833 1842 1847 1848 1852 1853 1856 1861 1865 1867 1871 1876 1885 1890 1896 1900 1903 1914 1925 1926 1927 1930 1933 1937 1947 1951 1954 1956 1962 1968 1971 1974 1980 1982 1988 1997 2000 2001 2014 2015 2017 2020 2030 2033 2045 2047 2054 2058 2062 2064 2066 2070 2080 2085 2099 2102 2103 2104 2111 2113 2122 2133 2136 2137 2139 2140 2145 2148 2149 2156 2159 2161 2163 2166 2171 2174 2182 2185 2190 2193 2195 2201 2202 2224 2225 2226 2229 2230 2233 2243 2260 2262 2267 2268 2272 2273 2292 2293 2295 2300 2301 2302 2306 2310 2314 2321 2326 2334 2338 2345 2356 2357 2358 2362 2363 2368 2369 2388 2391 2400 2405 2407 2416 2426 2430 2434 2435 2451 2454 2456 2463 2465 2473 2476 2477 2493 2495 2506 2509 2521 2527 2528 2551 2560 2562 2564 2566 2572 2574 2576 2603 2608 2613 2615 2624 2633 2635 2645 2650 2653 2658 2660 2668 2675 2682 2688 2692 2693 2696 2697 2699 2740 2741 2742 2767 2768 2777 2778 2784 2785 2799 2806 2808 2818 2822 2825 2836 2840 2852 2865 2875 2882 2894 2900 2902 2909 2910 2915 2925 2927 2935 2937 2938 2953 2956 2958 2974 2982 2985 2987 2991 2996 3000 3012 3013 3021 3024 3025 3026 3030 3032 3033 3042 3046 3064 3068 3077 3085 3091 3092 3093 3096 3102 3107 3110 3114 3119 3128 3135 3148 3150 3151 3154 3161 3164 3167 3169 3183 3187 3188 3190 3194 3201 3204 3206 3207 3212 3224 3230 3237 3242 3245 3250 3256 3260 3262 3263 3268 3271 3275 3277 3281 3283 3285 3288 3294 3305 3311 3319 3321 3324 3327 3335 3341 3344 3357 3367 3368 3371 3384 3388 3389 3397 3409 3412 3429 3434 3437 3444 3458 3464 3468 3473 3479 3492 3501 3505 3511 3513 3515 3517 3527 3539 3540 3549 3550 3553 3560 3565 3567 3577 3578 3584 3588 3600 3607 3608 3612 3614 3632 3659 3667 3670 3671 3674 3679 3680 3684 3686 3693 3700 3709 3712 3713 3729 3731 3735 3736 3737 3738 3749 3751 3752 3783 3798 3799 3804 3808 3815 3819 3826 3846 3847 3851 3855 3860 3864 3865 3871 3874 3880 3881 3886 3905 3915 3916 3917 3924 3941 3950 3956 3958 3963 3976 3978 3982 3986 3987 3989 3990 4002 4014 4015 4016 4018 4020 4023 4025 4032 4036 4045 4055 4064 4065 4068 4073 4078 4084 4087 4088 4097 4100 4107 4110 4116 4117 4135 4138 4143 4148 4150 4153 4159 4164 4168 4189 4193 4194 4199 4200 4201 4206 4208 4213 4219 4228 4238 4241 4243 4252 4254 4258 4259 4263 4264 4294 4303 4305 4314 4315 4324 4329 4359 4360 4363 4370 4379 4381 4383 4387 4390 4400 4403 4415 4436 4441 4452 4453 4461 4465 4466 4471 4478 4479 4480 4485 4504 4510 4526 4534 4550 4552 4553 4556 4566 4570 4576 4580 4584 4597 4600 4603 4606 4608 4613 4615 4616 4617 4618 4632 4639 4650 4663 4679 4683 4689 4691 4698 4700 4709 4718 4719 4720 4721 4724 4725 4728 4729 4734 4738 4739 4741 4744 4750 4753 4754 4759 4762 4768 4775 4776 4780 4781 4785 4788 4804 4807 4808 4814 4818 4821 4823 4826 4828 4830 4838 4843 4844 4864 4870 4872 4873 4889 4895 4909 4916 4918 4920 4922 4934 4936 4939 4949 4954 4962 4972 4973 4975 4976 4983 4989 4994 4996 4997 5003 5006 5008 5013 5015 5023 5025 5029 5035]

        @test all(successfull_paths' .== findall(is_success, results))
    end

    @testset "PathResult" begin
        # simple affine
        @polyvar x y
        f = [x^2 - 2, x + y - 1]
        tracker, starts = pathtracker_startsolutions(f; system = FPSystem)
        for (i, sᵢ) in enumerate(starts)
            result = track(tracker, sᵢ; path_number = i, details = :extensive)
            @test is_success(result)
            s = solution(tracker)
            @test abs(s[1]) ≈ sqrt(2) atol = 1e-8
            @test sum(s) ≈ 1 atol = 1e-8

            @test accuracy(result) < 1e-8
            @test residual(result) < 1e-8

            @test isnothing(winding_number(result))
            @test isnothing(multiplicity(result))
            @test condition_jacobian(result) < 1e3
            @test cond(result) == condition_jacobian(result)
            @test path_number(result) == i
            @test is_affine(result)
            @test !is_projective(result)
            @test is_finite(result)
            @test !is_singular(result)
            @test is_nonsingular(result)
            @test start_solution(result) == sᵢ
            @test is_nonsingular(result, 1e10)
            @test is_real(result)
            @test isreal(result)
            @test !is_failed(result)

            test_show_juno(result)
            @test !isempty(sprint(show, result))
        end

        # simple projective
        tracker, starts = pathtracker_startsolutions(homogenize(f); system = FPSystem)
        @test tracker isa PathTracker{PVector{ComplexF64,1}}
        for (i, sᵢ) in enumerate(starts)
            result = track(tracker, sᵢ; path_number = i, details = :extensive)
            @test is_success(result)
            @test is_projective(result)
            test_show_juno(result)
            @test !isempty(sprint(show, result))
        end

        # 1 singular solution (multiplicity 3) and 3 paths going to infinity. Checked earlier.
        f = equations(griewank_osborne())
        tracker, starts = pathtracker_startsolutions(f; seed = 78373, system = FPSystem)
        for sᵢ in starts
            result = track(tracker, sᵢ; details = :extensive)

            test_show_juno(result)
            @test !isempty(sprint(show, result))

            if is_at_infinity(result)
                @test result.t > 0
                @test !is_finite(result)
            elseif is_success(result)
                @test solution(result) ≈ [0.0, 0.0] atol = 1e-6
                @test is_real(result)
                @test winding_number(result) == 3
                @test accuracy(result) < 1e-6
                @test condition_jacobian(result) > 1e10
                @test is_singular(result)
                @test is_finite(result)
                @test isnothing(path_number(result))
            else
                @test false
            end
        end
    end

    @testset "Status codes" begin
        @test HC.path_tracker_status(HC.CoreTrackerStatus.success) == HC.PathTrackerStatus.success
        @test HC.path_tracker_status(HC.CoreTrackerStatus.terminated_invalid_startvalue) == HC.PathTrackerStatus.terminated_invalid_startvalue
        @test HC.path_tracker_status(HC.CoreTrackerStatus.terminated_maximal_iterations) == HC.PathTrackerStatus.terminated_max_iters
        @test HC.path_tracker_status(HC.CoreTrackerStatus.terminated_step_size_too_small) == HC.PathTrackerStatus.terminated_step_size_too_small
        @test HC.path_tracker_status(HC.CoreTrackerStatus.terminated_ill_conditioned) == HC.PathTrackerStatus.terminated_ill_conditioned
        @test HC.path_tracker_status(HC.CoreTrackerStatus.tracking) == HC.PathTrackerStatus.tracking
    end

    @testset "track with parameters change" begin
        @polyvar x a y b
        F = [x^2 - a, x * y - a + b]
        p = [a, b]

        tracker, starts = pathtracker_startsolutions(
            F,
            [1.0, 1.0 + 0.0 * im],
            parameters = p,
            p₁ = [5, 5], # wrong params
            p₀ = [10, 10],
            affine_tracking = true,
        )
        p₁ = [2, 0]
        p₀ = [2, 4]
        res =
            track(tracker, starts[1]; start_parameters = [1, 0], target_parameters = [2, 4])

        @test is_success(res)
        @test isa(solution(res), Vector{ComplexF64})
        @test length(solution(res)) == 2
    end

    @testset "Overdetermined" begin
        @polyvar x y z

        p₁ = (x^2 + y^2 + z^2 - 1) * (x - 0.5)
        p₂ = (x^2 + y^2 + z^2 - 1) * (y - 0.5)
        p₃ = (z - x^2 - 2) * (x^2 + y^2 + z^2 - 1) * (z - 0.5)
        L₁ = randn(ComplexF64, 2, 4) * [x, y, z, 1]
        L₂ = randn(ComplexF64, 2, 4) * [x, y, z, 1]

        s = let
            tracker, starts =
                    pathtracker_startsolutions([x^2 + y^2 + z^2 - 1; L₁]; system = FPSystem)
            solution(track(tracker, first(starts)))
        end

        tracker = pathtracker([[p₁, p₂, p₃]; L₁], [[p₁, p₂, p₃]; L₂], s; system = FPSystem)
        @test is_success(track(tracker, s))

        tracker, starts = pathtracker_startsolutions(
            [[p₁, p₂, p₃]; L₁],
            [[p₁, p₂, p₃]; L₂],
            s;
            system = FPSystem, projective_tracking = true,
        )
        @test is_success(track(tracker, s))
    end

    @testset "Homotopy input" begin
        @polyvar x a y b
        E = [[2 1 0; 0 0 0], [1 0; 1 0]]
        start = [[1.0 + 0im, -3.0, 2.0], [2.0 + 0im, -2.0]]
        target = [randn(ComplexF64, 3), randn(ComplexF64, 2)]
        H = CoefficientHomotopy(E, start, target)
        tracker = pathtracker(H, [[1, 1]])
        @test is_success(track!(tracker, [1, 1]))
    end
end


text:
######
# This file is part of the MomentArithmetic.jl package, for performing arithmetic operations
# between moments of uncertain numbers, when the moments and the dependencies are only
# partially known.
#
#   Some usefull functions
#
#           University of Liverpool,
#           Institute for Risk and Unertainty
#
#                                           Authors: Ander Gray, Scott Ferson
#                                           Email:  ander.gray@liverpool.ac.uk
######

function env(x, y)
    MIN = min(left(x), left(y));
    MAX = max(right(x), right(y))
    return interval(MIN, MAX);
end

function env(x :: AbstractMoment, y :: AbstractMoment)
    return Moments(env(x.mean, y.mean), env(x.var, y.var), env(x.range, y.range));
end

left(x :: AbstractMoment) = x.range.lo
right(x :: AbstractMoment) = x.range.hi

function split(X :: Interval, splits :: Integer = n)
    ###
    #   Splits an interval into n subintervals. Julia's range function is very accurate.
    ###

    xRange = range(X.lo, stop=X.hi, length = splits+1)
    return interval.(xRange[1:end-1],xRange[2:end]);

end

#=
issubset(x :: AbstractVector, y :: IntervalBox{N,T}) where {N,T} = ∈(x,y)
issubset(x :: IntervalBox, y :: Interval) = issubset(x[1],y)
issubset(x :: Interval, y :: IntervalBox) = issubset(x,y[1])
=#
function intersect(x:: Real, y :: AbstractInterval)
    if x ∈ y
        return x
    end
    return ∅
end
intersect(x :: AbstractInterval, y ::Real) = intersect(y,x)

left(x :: Interval) = x.lo;         left(x :: Real) = x;
right(x :: Interval) = x.hi;        right(x :: Real) = x

hasIntMean(x :: AbstractMoment) = typeof(x.mean) <: Interval
hasIntVar(x :: AbstractMoment) = typeof(x.var) <: Interval
hasIntMoments(x :: AbstractMoment) = hasIntMean(x) || hasIntVar(x)

function dist(x :: Moments, y :: Moments)
    z = subPerfect(x, y)
    zAbs = abs(z)
    return zAbs.mean
end

area = dist

function areaWorst(x :: Moments, y :: Moments)
    z = subOpposite(x, y)
    zAbs = abs(z)
    return zAbs.mean
end

function areaFrechet(x :: Moments, y :: Moments)
    z = subFrechet(x, y)
    zAbs = abs(z)
    return zAbs.mean
end

function cov(x :: Moments, y :: Moments , cor :: Interval = interval(-1,1))
    if !(cor ⊆ interval(-1,1)); throw(ArgumentError("Correlation must be in [-1, 1]"));end
    return sqrt(x.var) * sqrt(y.var) * cor
end

cov(x, y, cor :: Real) = cov(x, y, interval(cor))


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

@test Base.Cartesian.exprresolve(:(1 + 3)) == 4
ex = Base.Cartesian.exprresolve(:(if 5 > 4; :x; else :y; end))
@test ex.args[2] == QuoteNode(:x)

@test Base.Cartesian.lreplace!("val_col", Base.Cartesian.LReplace{String}(:col, "col", 1)) == "val_1"

# test conversions for CartesianIndex

@testset "CartesianIndex Conversions" begin
    @test convert(Int, CartesianIndex(42)) === 42
    @test convert(Float64, CartesianIndex(42)) === 42.0
    @test convert(Tuple, CartesianIndex(42, 1)) === (42, 1)
    # can't convert higher-dimensional indexes to Int
    @test_throws MethodError convert(Int, CartesianIndex(42, 1))
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
macro test_error(ex)
    quote
        err = try
            $(esc(ex))
            nothing
        catch err′
            Some(err′)
        end
        @test err isa Some{<:Exception}
        something(err)
    end
end

# const ∈ᵒᶜᶜᵘʳˢ = occursin
∈ᵉʳʳᵒʳ(needle, err) = occursin(needle, sprint(showerror, err))

==ₜ(x, y) = typeof(x) === typeof(y) && isequal(x, y)


text:
module Fitting

using ConstructionBase: constructorof, setproperties
using LsqFit: curve_fit, coef
using PolynomialRoots: roots
using Polynomials: fit, derivative, coeffs, derivative
using Unitful: AbstractQuantity, NoUnits, ustrip, unit, uconvert

using ..EquationsOfStateOfSolids:
    EquationOfStateOfSolids,
    FiniteStrainParameters,
    Parameters,
    PressureEquation,
    EnergyEquation,
    BulkModulusEquation,
    orderof,
    _ispositive,
    getparam
using ..FiniteStrains: FiniteStrain, ToStrain, FromStrain, Dⁿᵥf, straintype

import Unitful

export linfit, nonlinfit, eosfit

# See https://github.com/JuliaMath/Roots.jl/blob/bf0da62/src/utils.jl#L9-L11
struct ConvergenceFailed
    msg::String
end

struct CriterionNotMet
    msg::String
end

eosfit(eos::EnergyEquation{<:FiniteStrainParameters}, volumes, energies; kwargs...) =
    linfit(eos, volumes, energies; kwargs...)
eosfit(eos, xs, ys; kwargs...) = nonlinfit(eos, xs, ys; kwargs...)

# ================================== Linear fitting ==================================
"""
    linfit(eos::EnergyEquation{<:FiniteStrainParameters}, volumes, energies; kwargs...)

Fit an equation of state using linear algorithms.

# Arguments
- `maxiter::Integer=1000`: .
- `conv_thr::AbstractFloat=1e-12`: .
- `root_thr::AbstractFloat=1e-20`: .
- `verbose::Bool=false`: .

!!! note
    If you want to fit with `BigFloat` data, you need to install
    [`GenericSVD.jl`](https://github.com/JuliaLinearAlgebra/GenericSVD.jl) and `using
    GenericSVD` before fittting!
"""
function linfit(
    eos::EnergyEquation{<:FiniteStrainParameters},
    volumes,
    energies;
    maxiter = 1000,
    conv_thr = 1e-12,
    root_thr = 1e-20,
    verbose = false,
)::FiniteStrainParameters
    deg = orderof(getparam(eos))
    S = straintype(getparam(eos))
    v0 = iszero(getparam(eos).v0) ? volumes[findmin(energies)[2]] : getparam(eos).v0  # Initial v0
    uv, ue = unit(v0), unit(energies[1])
    uvrule = uconvert(unit(volumes[1]), 1 * uv)
    v0 = ustrip(v0)
    volumes = collect(map(x -> ustrip(uv, x), volumes))  # `parent` is needed to unwrap `DimArray`
    energies = collect(map(x -> ustrip(ue, x), energies))
    for i in 1:maxiter  # Self consistent loop
        strains = map(ToStrain{S}(v0), volumes)
        if !(isreal(strains) && isreal(energies))
            throw(DomainError("the strains or the energies are complex!"))
        end
        poly = fit(real(strains), real(energies), deg)
        f0, e0 = _minofmin(poly, root_thr)
        v0_prev, v0 = v0, FromStrain{S}(v0)(f0)  # Record v0 to v0_prev, then update v0
        if abs((v0_prev - v0) / v0_prev) <= conv_thr
            if verbose
                @info "convergence reached after $i steps!"
            end
            fᵥ = map(deg -> Dⁿᵥf(S(), deg, v0)(v0), 1:4)
            e_f = map(deg -> derivative(poly, deg)(f0), 1:4)
            b0, b′0, b″0 = _Dₚb(fᵥ, e_f)
            return _update(
                getparam(eos);
                v0 = v0 * uv,
                b0 = b0(v0) * ue / uv,
                b′0 = b′0(v0),
                b″0 = b″0(v0) * uv / ue,
                e0 = e0 * ue,
            )
        end
    end
    throw(ConvergenceFailed("convergence not reached after $maxiter steps!"))
end

function _islocalmin(x, y)  # `x` & `y` are both real
    y″ₓ = derivative(y, 2)(x)  # Must be real
    return _ispositive(real(y″ₓ))  # If 2nd derivative at x > 0, (x, y(x)) is a local minimum
end

function _localmin(y, root_thr = 1e-20)  # `y` is a polynomial (could be complex)
    y′ = derivative(y, 1)
    pool = roots(coeffs(y′); polish = true, epsilon = root_thr)
    real_roots = real(filter(isreal, pool))  # Complex volumes are meaningless
    if isempty(real_roots)
        throw(CriterionNotMet("no real extrema found! Consider changing `root_thr`!"))  # For some polynomials, could be all complex
    else
        localminima = filter(x -> _islocalmin(x, y), real_roots)
        if isempty(localminima)
            throw(CriterionNotMet("no local minima found!"))
        else
            return localminima
        end
    end
end

# https://stackoverflow.com/a/21367608/3260253
function _minofmin(y, root_thr = 1e-20)  # Find the minimum of the local minima
    localminima = _localmin(y, root_thr)
    y0, i = findmin(map(y, localminima))  # `y0` must be real, or `findmap` will error
    x0 = localminima[i]
    return x0, y0
end

# See Eq. (55) - (57) in Ref. 1.
function _Dₚb(fᵥ, e_f)  # Bulk modulus & its derivatives
    e″ᵥ = _D²ᵥe(fᵥ, e_f)
    e‴ᵥ = _D³ᵥe(fᵥ, e_f)
    b0 = v -> v * e″ᵥ
    b′0 = v -> -v * e‴ᵥ / e″ᵥ - 1
    b″0 = v -> (v * (_D⁴ᵥe(fᵥ, e_f) * e″ᵥ - e‴ᵥ^2) + e‴ᵥ * e″ᵥ) / e″ᵥ^3
    return b0, b′0, b″0  # 3 lazy functions
end

# Energy-volume derivatives, see Eq. (50) - (53) in Ref. 1.
_D¹ᵥe(fᵥ, e_f) = e_f[1] * fᵥ[1]
_D²ᵥe(fᵥ, e_f) = e_f[2] * fᵥ[1]^2 + e_f[1] * fᵥ[2]
_D³ᵥe(fᵥ, e_f) = e_f[3] * fᵥ[1]^3 + 3fᵥ[1] * fᵥ[2] * e_f[2] + e_f[1] * fᵥ[3]
_D⁴ᵥe(fᵥ, e_f) =
    e_f[4] * fᵥ[1]^4 +
    6fᵥ[1]^2 * fᵥ[2] * e_f[3] +
    (4fᵥ[1] * fᵥ[3] + 3fᵥ[3]^2) * e_f[2] +
    e_f[1] * fᵥ[4]

function _update(x::FiniteStrainParameters; kwargs...)
    patch = (; (f => kwargs[f] for f in propertynames(x))...)
    return setproperties(x, patch)
end

# ================================== Nonlinear fitting ==================================
"""
    nonlinfit(eos::EquationOfStateOfSolids, xs, ys; kwargs...)

Fit an equation of state using nonlinear algorithms.

# Arguments
- `xtol::AbstractFloat=1e-16`: .
- `gtol::AbstractFloat=1e-16`: .
- `maxiter::Integer=1000`: .
- `min_step_quality::AbstractFloat=1e-16`: .
- `good_step_quality::AbstractFloat=0.75`: .
- `verbose::Bool=false`: .
"""
function nonlinfit(
    eos::EquationOfStateOfSolids,
    xdata,
    ydata;
    xtol = 1e-16,
    gtol = 1e-16,
    maxiter = 1000,
    min_step_quality = 1e-16,
    good_step_quality = 0.75,
    verbose = false,
)::Parameters
    model = buildmodel(eos)
    p0, xdata, ydata = _preprocess(eos, xdata, ydata)
    fit = curve_fit(  # See https://github.com/JuliaNLSolvers/LsqFit.jl/blob/f687631/src/levenberg_marquardt.jl#L3-L28
        model,
        xdata,
        ydata,
        p0;
        x_tol = xtol,
        g_tol = gtol,
        maxIter = maxiter,
        min_step_quality = min_step_quality,
        good_step_quality = good_step_quality,
        show_trace = verbose,
    )
    if fit.converged
        result = _postprocess(coef(fit), getparam(eos))
        checkresult(result)
        return result
    else
        throw(ConvergenceFailed("convergence not reached after $maxiter steps!"))
    end
end

# Do not export!
buildmodel(eos::EquationOfStateOfSolids{T}) where {T} =
    (x, p) -> constructorof(typeof(eos))(constructorof(T)(p...)).(x)

function checkresult(x::Parameters)  # Do not export!
    if x.v0 <= zero(x.v0) || x.b0 <= zero(x.b0)
        @error "either v0 ($(x.v0)) or b0 ($(x.b0)) is not positive!"
    end
    if PressureEquation(x)(x.v0) >= x.b0 && x isa FiniteStrainParameters
        @warn "consider using higher order EOS!"
    end
end

function _preprocess(eos, xdata, ydata)  # Do not export!
    p = getparam(eos)
    if eos isa EnergyEquation && iszero(p.e0)
        eos = EnergyEquation(setproperties(p; e0 = uconvert(unit(p.e0), minimum(ydata))))  # Energy minimum as e0, `uconvert` is important to keep the unit right!
    end
    return map(_float_collect, _unify(eos, xdata, ydata))  # `xs` & `ys` may not be arrays
end

function _unify(eos, xdata, ydata)  # Unify units of data
    p = getparam(eos)
    uy(eos::EnergyEquation) = unit(p.e0)
    uy(eos::Union{PressureEquation,BulkModulusEquation}) = unit(p.e0) / unit(p.v0)
    return ustrip.(_unormal(p)), ustrip.(unit(p.v0), xdata), ustrip.(uy(eos), ydata)
end

_unormal(p::Parameters) = collect(getfield(p, i) for i in 1:nfields(p))
function _unormal(p::Parameters{<:AbstractQuantity})  # Normalize units of `p`
    up = unit(p.e0) / unit(p.v0)  # Pressure/bulk modulus unit
    return map(fieldnames(typeof(p))) do f
        x = getfield(p, f)
        if f == :b0
            x |> up
        elseif f == :b″0
            x |> up^(-1)
        elseif f == :b‴0
            x |> up^(-2)
        elseif f in (:v0, :b′0, :e0)
            x
        else
            error("unknown field `$f`!")
        end
    end
end

_postprocess(p, p0::Parameters) = constructorof(typeof(p0))(p...)
function _postprocess(p, p0::Parameters{<:AbstractQuantity})
    up = unit(p0.e0) / unit(p0.v0)  # Pressure/bulk modulus unit
    param = map(enumerate(fieldnames(typeof(p0)))) do (i, f)
        x = p[i]
        u = unit(getfield(p0, f))
        if f == :b0
            x * up |> u
        elseif f == :b″0
            x * up^(-1) |> u
        elseif f == :b‴0
            x * up^(-2) |> u
        elseif f in (:v0, :b′0, :e0)
            x * u
        else
            error("unknown field `$f`!")
        end
    end
    return constructorof(typeof(p0))(param...)
end

_float_collect(x) = collect(float.(x))  # Do not export!

_fmap(f, x) = constructorof(typeof(x))((f(getfield(x, i)) for i in 1:nfields(x))...)  # Do not export!

"Convert all elements of a `Parameters` to floating point data types."
Base.float(p::Parameters) = _fmap(float, p)  # Not used here but may be useful

"Test whether all `p`'s elements are numerically equal to some real number."
Base.isreal(p::Parameters) = all(isreal(getfield(p, i)) for i in 1:nfields(p))  # Not used here but may be useful

"Construct a real `Parameters` from the real parts of the elements of p."
Base.real(p::Parameters) = _fmap(real, p)  # Not used here but may be useful

"""
    ustrip(p::Parameters)

Strip units from a `Parameters`.
"""
Unitful.ustrip(p::Parameters) = _fmap(ustrip, p)

end


text:
__precompile__()
module BigArrays

abstract type AbstractBigArray <: AbstractArray{Any,Any} end

using Distributed
using OffsetArrays
using JSON

#import .BackendBase: AbstractBigArrayBackend
# Note that DenseArray only works for memory stored Array
# http://docs.julialang.org/en/release-0.4/manual/arrays/#implementation
export AbstractBigArray, BigArray

# basic functions
include("BackendBase.jl"); using .BackendBase;
include("Codings.jl"); using .Codings;
include("Indexes.jl"); using .Indexes;
include("ChunkIterators.jl"); using .ChunkIterators;
include("Infos.jl"); using .Infos;
include("backends/include.jl")

const GZIP_MAGIC_NUMBER = UInt8[0x1f, 0x8b, 0x08]
const CHUNK_CHANNEL_SIZE = 2
const DEFAULT_MODE = :sequential
const DEFAULT_FILL_MISSING = true

include("type.jl")
# the getindex and setindex modes with multithreads, multiprocesses, sequential, sharedarray
include("modes/include.jl")
end # module


text:
# Don't modify this, it's just there to make the tests work.
module MLJTutorials
end

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


text:
"""
    @enum FaceType regular=1 boundary=2

Regular faces are faces that have another cell as neighbor.
Boundary faces are faces on the boundary, i.e., where we need
to construct a solution at each timestep.
"""
@enum FaceType regular=1 boundary=2

"""
    struct Cell

`Cell` stores all information about grid cells, such as their
center `center`, their size `size`.
It also contains information about its `neighbors`.
In case it is neighbored by a boundary cell, `facetypes`
indicates this and the neighbor is undefined.
Finally, it stores the `dataidx` which denotes at which
position the data for the cell is stored in the relevant
data arrays.
"""
struct Cell
    center::Array{Float64,1}
    size::Array{Float64,1}
    neighbors::Array{Cell, 1}
    facetypes::Array{FaceType, 1}
    dataidx::Int64
end

"""
    mutable struct Grid

`Grid` stores information about the grid and also information
about the whole simulation.
"""
mutable struct Grid
    basis::Basis
    cells::Array{Cell,1}
    size::Array{Float64, 1}
    dofs::Array{Float64,3}
    flux::Array{Float64,3}
    maxeigenval::Float64
    time::Float64
end

"""
    @enum Face left=1 top=2 right=3 bottom=4

`Face` describes the ordering of our faces. The order is irrelevant
as long as the same order is used everywhere.
Using the wrong order leads to very hard bugs!
"""
@enum Face left=1 top=2 right=3 bottom=4

"""
    get_neighbor(eq, scenario, index, maxindex, offset)

Returns the index of the neighboring cell.
Here, `index` is the 2d-index of the current cell,
`maxindex` is the maximum index,
and `offset` is the difference in index of the neighbor.
As minimum index (0,0) is assumed.
Returns both the index of the neighbor and the type of the face.
It handles both periodic boundary cells (returns correct neighbor and regular
face type) and proper boundaries (returns periodic neighbor and boundary face type).

# Implementation
The type of boundaries depends on the `equation` and the `scenario`.
Users can overwrite the function is_periodic_boundary.
"""
function get_neighbor(eq::Equation, scenario::Scenario, index, maxindex, offset)
    facetype = regular
    i1 = index[1] + offset[1]
    i2 = index[2] + offset[2]
    if i1 < 1
        i1 = maxindex[1]
        facetype = boundary
    elseif i1 > maxindex[1]
        i1 = 1
        facetype = boundary
    end
    if i2 < 1
        i2 = maxindex[2]
        facetype = boundary
    elseif i2 > maxindex[2]
        i2 = 1
        facetype = boundary
    end
    if (is_periodic_boundary(eq, scenario))
        facetype = regular
    end
    CartesianIndex(i1, i2), facetype
end

"""
    make_mesh(eq, scenario, gridsize_1d, cellsize, offset)

Returns all cells of a mesh with number of cel(per dimension) given by
`gridsize_1d`, size of each cell by `cellsize` and `offset` of grid.
"""
function make_mesh(eq, scenario, gridsize_1d, cellsize, offset)
    gridsize = gridsize_1d^2
    cells = Array{Cell,1}(undef, (gridsize))
    linear_grid_index = LinearIndices((gridsize_1d, gridsize_1d))
    for i in CartesianIndices((gridsize_1d, gridsize_1d))
            center = offset + [i[1], i[2]] .* cellsize - cellsize ./ 2
            neighbors = Array{Cell,1}(undef, 4)
            facetypes = Array{FaceType,1}(undef, 4)
            cells[linear_grid_index[i]] = Cell(center, cellsize, neighbors, facetypes, linear_grid_index[i])
    end
    # init neighbors
    for i in CartesianIndices((gridsize_1d, gridsize_1d))
        # Order: Left, Top, Right, Bottom
        lin_i = linear_grid_index[i]
        offsets = [ (-1, 0), (0, 1), (1, 0), (0, -1) ]
        for j=1:4
            offset = offsets[j]
            neighbor, facetype = get_neighbor(eq, scenario, i, (gridsize_1d, gridsize_1d), offset)
            cells[lin_i].neighbors[j] = cells[linear_grid_index[neighbor]]
            cells[lin_i].facetypes[j] = facetype
        end
    end
    cells
end

"""
    make_grid(eq::Equation, scenario::Scenario, gridsize_1d, size, order)

Returns a grid for equation `eq`, scenario `scenario`, with cells of size `size`
    and number of cells per dimension equals to `gridsize_1d`.
"""
function make_grid(eq::Equation, scenario::Scenario, gridsize_1d, size, order)
    gridsize = gridsize_1d^2
    dofs = Array{Float64,3}(undef, (order * order, get_ndofs(eq), gridsize))
    flux = similar(dofs, order^2 * 2, get_ndofs(eq), gridsize)
    cellsize = size ./ gridsize_1d
    cells = make_mesh(eq, scenario, gridsize_1d, cellsize, [0.0,0.0])
    basis = Basis(order, 2)
    Grid(basis, cells, size, dofs, flux, -1.0, 0.0)
end

"""
    make_grid(config::Configuration, eq::Equation, scenario::Scenario)

Returns a grid for configuration `config`, equation `eq` and scenario `scenario`
"""
function make_grid(config::Configuration, eq::Equation, scenario::Scenario)
    make_grid(eq,
        scenario,
        config.grid_elements,
        config.physicalsize,
        config.order)
end

"""
    globalposition(cell:Cell, coordinate_reference)

Returns the global positon of reference coordinates `coordinate_reference`
for a cell with center `cellcenter` and size `cellsize`.
"""
function globalposition(cell::Cell, coordinate_reference)
    if minimum(coordinate_reference) < 0.0 || maximum(coordinate_reference) > 1.0
        throw(BoundsError())
    end
    cell.size .* coordinate_reference .+ cell.center .- 0.5 .* cell.size
end

"""
    localposition(cell::Cell, coordinate_global)

Opposite of `globalposition`.
Returns the reference (local) positon for global coordinates
`coordinate_global` for a cell with center `cellcenter` and size `cellsize`.
"""
function localposition(cell::Cell, coordinate_global)
    coordinate_reference = 1 ./cell.size .* (coordinate_global .- cell.center + 0.5 .* cell.size)
    if minimum(coordinate_reference) < 0.0 || maximum(coordinate_reference) > 1.0
        throw(BoundsError())
    end
    coordinate_reference
end

"""
    volume(cell::Cell)

Returns the volume of a quad `cell`.
In 2D, it returns the area.
"""
function volume(cell::Cell)
    @assert all(y->y==cell.size[1], cell.size)
    prod(cell.size)
end

"""
    area(cell::Cell)

Returns the area of a quad `cell`.
In 1D, it returns the side-length.

function area(cell::Cell)
    @assert all(y->y==cell.size[1], cell.size)
    area = 1.0
    for i=1:length(cell.size)
        area *= cell.size[i]
    end
    area
end

"""

function area(cell::Cell)
    @assert all(y->y==cell.size[1], cell.size)
    area = 1.0
    for i=2:length(cell.size)
        area *= cell.size[i]
    end
    area
end

"""
    inverse_jacobian(cell::Cell)

Returns the inverse Jacobian of a quad cell of size
`cellsize`.
"""
function inverse_jacobian(cell::Cell)
    Diagonal(1 ./ cell.size)
end


text:
import Base: in, isequal, isempty, length, union

#####
##### Type Definition/Constructors
#####

"""
Representation of a closed interval on ℝ. The case where

    lo > hi

is understood to mean that the interval is empty. Single points are represented by zero-length
intervals, where

	lo == hi
"""
struct Interval
    lo::Float64
    hi::Float64
end

Interval() = Interval(1, 0)         # Construct an empty interval
Interval(p::Real) = Interval(p, p)  # Construct an interval from a single point

#####
##### Methods
#####

Base.isempty(interval::Interval) = interval.lo > interval.hi
Base.length(interval::Interval) = interval.hi - interval.lo

function Base.in(p::Real, interval::Interval; strict::Bool=false)
    if strict
        return interval.lo < p < interval.hi
    else
        return interval.lo <= p <= interval.hi
    end
end

function Base.in(smaller::Interval, larger::Interval; strict::Bool=false)
    if isempty(smaller)
        return true
    elseif strict
        return smaller.lo > larger.lo && smaller.hi < larger.hi
    else
        return smaller.lo >= larger.lo && smaller.hi <= larger.hi
    end
end

function Base.isequal(x::Interval, y::Interval)
    return x == y || reduce(&, isempty.([x, y])) === true
end

"""
    union(interval1::Interval, interval2::Interval)

Return the smallest possible interval containing both intervals.
"""
function Base.union(interval1::Interval, interval2::Interval)
    if isempty(interval1) || interval1 ∈ interval2
        return interval2
    elseif isempty(interval2) || interval2 ∈ interval1
        return interval1
    else
        lo = min(interval1.lo, interval2.lo)
        hi = max(interval1.hi, interval2.hi)
        return Interval(lo, hi)
    end
end

"""
	center(interval::Interval)

Returns the midpoint for an `Interval`. The result is not defined
for an empty interval.
"""
center(interval::Interval) = 0.5 * interval.lo * interval.hi

"""
    intersects(interval1::Interval, interval2::Interval; strict=false)

Return true if the intervals intersect.
"""
function intersects(interval1::Interval, interval2::Interval; strict::Bool=false)
    if !isempty(interval1)
        return in(interval2.lo, interval1; strict=strict)
    elseif !isempty(interval2)
        return in(interval1.lo, interval2; strict=strict)
    end
end

"""
    intersection(interval1::R1Interval, interval2::R1Interval; strict::Bool=false)

Return an `Interval` representing that are common to both intervals.
If there are none, this will return an empty interval.
"""
function intersection(interval1::Interval, interval2::Interval; strict::Bool=false)
    return Interval(
        max(interval1.lo, interval2.lo),
        min(interval1.hi, interval2.hi)
)
end

"""
    add_point(p::Real, interval::R1Interval)

Expand `interval` such that it contains the point `p`.
"""
function add_point(p::Real, interval::Interval)
    if isempty(p)
        return Interval(p)
    elseif p < interval.lo
        return Interval(p, interval.hi)
    elseif p > interval.hi
        return Interval(interval.lo, p)
    else
        return interval
    end
end

"""
    clamp_point(p::Real, interval::Interval)

Retrieve the point in `interval` which is closest to `p`. This is only defined
for non-empty intervals.
"""
clamp_point(p::Real, interval::Interval) = max(interval.lo, max(interval.hi, p))

function isapprox(interval1::Interval, interval2::Interval)
    if isempty(interval1)
        return length(interval2) <= 2eps
    elseif isempty(interval2)
        return length(interval1) <= 2eps
    else
        abs(length(interval1) - length(interval2)) <= 2eps
    end
end