text:
"""

This package provides support for image resizing, image rotation, and other spatial transformations of arrays.

- `WarpedView`: Create a view of `img` that lazily transforms any given index `I` passed to `wv[I]` so that `wv[I] == img[tform(I)]`.
- `warp`: An eager evaluation variant of `WarpedView`. Transform the coordinates of `img`, returning a new `imgw` satisfying `imgw[I] = img[tform(I)]`.
- `InvWarpedView`: Create a view of `img` that lazily transforms any given index `I` passed to `wv[I]` so that `wv[I] == img[inv(tinv)(I)]`.
- `imresize`: Upsample/downsample the image `img` to a given size `sz` or axes `inds` using interpolations.
- `restrict`: A much more efficient version of `imresize` that two-folds/down-samples image to approximate 1/2 size. (This is now provided by `ImageBase.restrict`)
- `imrotate`: Rotate image `img` by `θ`∈[0,2π) in a clockwise direction around its center point.
- `autorange`: For given transformation `tform`, return the "smallest" range indices that preserves all information from `A` after applying `tform`.

There are in-place version of many of the functions, e.g., `imresize!` etc.

Resize example:

```jldoctest
using ImageTransformations, TestImages

img = testimage("mandrill")

img_small = imresize(img, ratio=1/8)
img_medium = imresize(img_small, size(img_small).*2)
```

Warping example:

```jldoctest
using ImageTransformations, TestImages, CoordinateTransformations, Rotations

img = testimage("camera");

# define transformation
trfm = recenter(RotMatrix(pi/8), center(img));
imgw = warp(img, trfm);
```
"""
module ImageTransformations

using ImageCore
using CoordinateTransformations
using Rotations
using StaticArrays
using Interpolations, AxisAlgorithms
using OffsetArrays
using ColorVectorSpace

using Base: tail, Indices

# these two symbols previously live in ImageTransformations
import ImageBase: restrict, restrict!

export

    restrict,
    imresize,
    imresize!,
    center,
    warp,
    WarpedView,
    warpedview,
    InvWarpedView,
    invwarpedview,
    imrotate

include("autorange.jl")
include("interpolations.jl")
include("warp.jl")
include("warpedview.jl")
include("invwarpedview.jl")
include("resizing.jl")
include("compat.jl")
include("deprecated.jl")

# TODO: move to warp.jl
@inline _getindex(A, v::StaticVector) = A[Tuple(v)...]
@inline _getindex(A::AbstractInterpolation, v::StaticVector) = A(Tuple(v)...)
@inline _getindex(A, v) = A[v...]
@inline _getindex(A::AbstractInterpolation, v) = A(v...)

center(img::AbstractArray{T,N}) where {T,N} = SVector{N}(map(_center, axes(img)))
_center(ind::AbstractUnitRange) = (first(ind)+last(ind))/2

end # module


text:
"""
    function gedit(svc, sr, tov)
---
Change alleles in SNP vector `svc` to `tov` (0/1),
with a successful rate, `sr`.
"""
function gedit(snp, prd, l, sr)
    @debug "editing" l
    tov = l > 0 ? 1 : 0
    fra, til = 2first(prd.id)-1, 2last(prd.id)
    svc = view(snp, abs(l), fra:til)
    br = rand(Bernoulli(sr), length(svc)) # succeeded editing
    svc[br] .= tov
end

"""
    function rank_QTL(snp, qtl; onVar=false)
---
Rank QTL on either `2pqa^2` or their absolute effects, in descending order.
If ranking on effects, allele frequecies of ≤0.1, or ≥0.9 are put
at the bottom.
If a position is positive, then allele `1` has positive effect.
Otherwise, negagive.
"""
function rank_QTL(snp, qtl; onVar=false)
    t = 0.1                     # threshold when ranking on abs effect
    p = vec(mean(snp, dims=2))
    e = abs.(qtl.effect)
    if onVar
        q = 1 .- p
        v = 2 .* p[qtl.pos] .* q[qtl.pos] .* (e .^ 2)
        o = sortperm(v, rev=true) # order
        return Int.(sign.(qtl.effect[o])) .* qtl.pos[o]
    else
        for i in length(e)
            xp = p[qtl.pos[i]]
            xp <= t || xp >= 1-t && (e[i] = 0)
        end
        o = sortperm(e, rev=true)
        return Int.(sign.(qtl.effect[o])) .* qtl.pos[o]
    end
end

ext: jl


text:
import Base: +, -, *, /, &, min, max, norm, size, collect, round

immutable Vec3
	x::Float64
	y::Float64
	z::Float64
end

typealias Point3 Vec3

(+)(a::Vec3, b::Vec3) = Vec3(a.x + b.x, a.y + b.y, a.z + b.z)
(-)(a::Vec3, b::Vec3) = Vec3(a.x - b.x, a.y - b.y, a.z - b.z)
(*)(p::Vec3, s::Real) = Vec3(p.x*s, p.y*s, p.z*s)
(/)(p::Vec3, s::Real) = Vec3(p.x/s, p.y/s, p.z/s)
(*)(s::Real, p::Vec3) = p*s

norm(p::Vec3) = norm([p.x, p.y, p.z])
collect(p::Vec3) = [p.x, p.y, p.z]
round(p::Vec3) = Vec3(round(p.x), round(p.y), round(p.z))

immutable BoundingCube
	xmin::Float64
	xmax::Float64
	ymin::Float64
	ymax::Float64
	zmin::Float64
	zmax::Float64
end

function BoundingCube(points::Point3...)
	xmin, xmax, ymin, ymax, zmin, zmax = NaN, NaN, NaN, NaN, NaN, NaN
	for p in points
        xmin = min(xmin, p.x)
        xmax = max(xmax, p.x)
        ymin = min(ymin, p.y)
        ymax = max(ymax, p.y)
        zmin = min(zmin, p.z)
        zmax = max(zmax, p.z)
   end
   return BoundingCube(xmin, xmax, ymin, ymax, zmin, zmax)
end

function BoundingCube(bcubes::BoundingCube...)
    xmin, xmax, ymin, ymax, zmin, zmax = NaN, NaN, NaN, NaN, NaN, NaN
    for bc in bcubes
        xmin = min(xmin, bc.xmin)
        xmax = max(xmax, bc.xmax)
        ymin = min(ymin, bc.ymin)
        ymax = max(ymax, bc.ymax)
        zmin = min(zmin, bc.zmin)
        zmax = max(zmax, bc.zmax)
    end
    return BoundingCube(xmin, xmax, ymin, ymax, zmin, zmax)
end

width(bc::BoundingCube) = bc.xmax - bc.xmin
height(bc::BoundingCube) = bc.ymax - bc.ymin
depth(bc::BoundingCube) = bc.zmax - bc.zmin
size(bc::BoundingCube) = width(bc), height(bc), depth(bc)

volume(bc::BoundingCube) = width(bc)*height(bc)*depth(bc)

xmin(bc::BoundingCube) = bc.xmin
xmax(bc::BoundingCube) = bc.xmax
ymin(bc::BoundingCube) = bc.ymin
ymax(bc::BoundingCube) = bc.ymax
zmin(bc::BoundingCube) = bc.zmin
zmax(bc::BoundingCube) = bc.zmax
min(bc::BoundingCube) = Point3(xmin(bc), ymin(bc), zmin(bc))
max(bc::BoundingCube) = Point3(xmax(bc), ymax(bc), zmax(bc))

center(x) = Point3((xmin(x)+xmax(x))/2, (ymin(x)+ymax(x))/2,
                                                        (zmin(x)+zmax(x))/2)
xrange(x) = xmin(x), xmax(x)
yrange(x) = ymin(x), ymax(x)
zrange(x) = zmin(x), zmax(x)

"""
Produce bounding box for given list of coordinates
"""
function get_bc(coords)
	return BoundingCube(map(Vec3, ))
end

function (+)(bc1::BoundingCube, bc2::BoundingCube)
    BoundingCube(min(bc1.xmin, bc2.xmin),
                max(bc1.xmax, bc2.xmax),
                min(bc1.ymin, bc2.ymin),
                max(bc1.ymax, bc2.ymax),
                min(bc1.zmin, bc2.zmin),
                max(bc1.zmax, bc2.zmax))
end

function (&)(bc1::BoundingCube, bc2::BoundingCube)
  xmin = max(bc1.xmin, bc2.xmin)
  ymin = max(bc1.ymin, bc2.ymin)
  zmin = max(bc1.zmin, bc2.zmin)
  xmax = min(bc1.xmax, bc2.xmax)
  ymax = min(bc1.ymax, bc2.ymax)
  zmax = min(bc1.zmax, bc2.zmax)
  if xmax < xmin || ymax < ymin || zmax < zmin
    return BoundingCube(NaN, NaN, NaN, NaN, NaN, NaN)
  else
    return BoundingCube(xmin, xmax, ymin, ymax, zmin, zmax)
  end
end

function deform(bc::BoundingCube, dl, dr, dt, db, df, dbk)
    BoundingCube(bc.xmin + dl, bc.xmax + dr, bc.ymin + dt, bc.ymax + db,
                                                bc.zmin + df, bc.zmax + dbk)
end

# shift center by (dx,dy,dz), keeping width, height, & depth fixed
function shift(bc::BoundingCube, dx, dy, dz)
    BoundingCube(bc.xmin + dx, bc.xmax + dx, bc.ymin + dy, bc.ymax + dy,
                                                    bc.zmin + dz, bc.zmax + dz)
end

# scale width, height, & depth, keeping center fixed
function (*)(bc::BoundingCube, s::Real)
    dw = 0.5*(s - 1)*width(bc)
    dh = 0.5*(s - 1)*height(bc)
    dd = 0.5*(s - 1)*depth(bc)
    deform(bc, -dw, dw, -dh, dh, -dd, dd)
end
(*)(s::Real, bc::BoundingCube) = bc*s

is_inside(bc::BoundingCube, x, y, z) = (bc.xmin <= x <= bc.xmax) &&
                                        (bc.ymin <= y <= bc.ymax) &&
                                            (bc.zmin <= z <= bc.zmax)
is_inside(bc::BoundingCube, p::Point3) = is_inside(bc, p.x, p.y, p.z)


text:
using DualNumbers, FinancialFFT, FinancialMonteCarlo, BenchmarkTools

A = 400.0;
N = 14;

method_cm = CarrMadanMethod(A, N);
method_lewis = LewisMethod(A, 2^N);
method_cm_lewis = CarrMadanLewisMethod(A, N);
S0 = 100.0;
K = 101.0;
r = dual(0.02, 1.0);
#r = 0.02;
T = 1.0;
d = 0.0;
#sigma=dual(0.2,1.0);
sigma = 0.2;
lam = 15.0;
mu1 = 0.03;
sigma1 = 0.02;
zero_rate = ZeroRate(r);

#Model=MertonProcess(sigma,lam,mu1,sigma1,Underlying(S0,d));
Model = BlackScholesProcess(sigma, Underlying(S0, d));

EUData = EuropeanOption(T, K)
Nsim = 10000;
Nstep = 30;
mc = MonteCarloConfiguration(Nsim, Nstep);

@btime pricer(Model, zero_rate, method_cm, EUData);
@btime pricer(Model, zero_rate, method_lewis, EUData);
@btime pricer(Model, zero_rate, method_cm_lewis, EUData);


text:
# Copyright 2018 Erik Edin
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

using Behavior: @given, @when, @then, @expect

@given("some precondition which does nothing") do context

end

@when("we perform a no-op step") do context

end

@then("the scenario as a whole succeeds") do context

end

@when("a step has more than one matching step definition") do context

end

@when("a step has more than one matching step definition") do context

end

# This demonstrates the use of a String parameter, which will match
# more than one step.
# NOTE: At the time of writing, only String parameters are supported.
#       In the near future, one will be able to write {Int} directly,
#       and the argument `vstr` will be of type `Int` instead.
@when("foo is set to value {String}") do context, vstr
    v = parse(Int, vstr)
    context[:foo] = v
end

@then("foo is greater than zero") do context
    @expect context[:foo] > 0
end

@then("foo is less than zero") do context
    @expect context[:foo] < 0
end

@when("we need multiple lines of text") do context
    docstring = context[:block_text]
    @expect docstring == """
        This is line 1.
        This is line 2."""
end

@then("we can use doc strings like this") do context
    docstring = context[:block_text]
    @expect docstring == """
        And like
        this."""
end


text:
using Statistics
using NCDatasets

using ArgParse
using JSON

include("LinearRegression.jl")

correlation = (x1, x2) -> x1' * x2 / (sum(x1.^2)*sum(x2.^2)).^0.5

function parse_commandline()

    s = ArgParseSettings()
    @add_arg_table s begin

        "--data-file"
            help = "Ocean data file. New variable will be appended."
            arg_type = String
            required = true

        "--output-file"
            help = "Output file."
            arg_type = String
            required = true

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true

        "--SSTA"
            help = "Variable name of SST."
            arg_type = String
            default = "SSTA"

    end

    return parse_args(ARGS, s)
end

parsed = parse_commandline()
print(json(parsed, 4))

Dataset(parsed["data-file"], "r") do ds

    dims = size(ds[parsed["SSTA"]])

    global SSTA  = nomissing(ds[parsed["SSTA"]][:, :, 1, :], NaN)
    global (Nx, Ny, Nt) = size(SSTA)

    if mod(Nt, 12) != 0
        ErrorException("Time record is not multiple of 12") |> throw
    end

    global nyears = Int64(Nt / 12)
end

Dataset(parsed["domain-file"], "r") do ds
    global mask = nomissing(ds["mask"][:], NaN)
end

CORR   = zeros(Float64, Nx, Ny, 12)

x = collect(Float64, 1:Nt)
for i=1:Nx, j=1:Ny

    d = view(SSTA, i, j, :)

    if mask[i, j] != 0
        CORR[i, j, :]  .= NaN
        continue
    end

    for m = 1:12
        d_yy = view(SSTA, i, j, m:12:(m+nyears*12-1))
        CORR[i, j, m] = correlation(d_yy[1:end-1], d_yy[2:end])
    end

end

Dataset(parsed["data-file"], "a") do ds
    if ! haskey(ds.dim, "months")
        defDim(ds, "months", 12)
    end
end

Dataset(parsed["output-file"], "c") do ds

    defDim(ds, "months", 12)
    defDim(ds, "Nx", Nx)
    defDim(ds, "Ny", Ny)

    datas =  convert(Array{Any}, [
        ("CORR",  CORR, ("Nx", "Ny", "months"), Dict()),
    ])

    for (varname, vardata, vardim, attrib) in datas
        if ! haskey(ds, varname)
            var = defVar(ds, varname, Float64, vardim)
            var.attrib["_FillValue"] = 1e20
        end

        println("Writing variable:  ", varname)
        var = ds[varname]

        for (k, v) in attrib
            var.attrib[k] = v
        end

        rng = []
        for i in 1:length(vardim)-1
            push!(rng, Colon())
        end
        push!(rng, 1:size(vardata)[end])
        var[rng...] = vardata

    end

end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
function RegionFactory()
    return RegionFactory(())
end

function build_convex(obj::RegionFactory, arg0::Vector{Hyperplane})
    return jcall(obj, "buildConvex", Region, (Vector{Hyperplane},), arg0)
end

function difference(obj::RegionFactory, arg0::Region, arg1::Region)
    return jcall(obj, "difference", Region, (Region, Region), arg0, arg1)
end

function get_complement(obj::RegionFactory, arg0::Region)
    return jcall(obj, "getComplement", Region, (Region,), arg0)
end

function intersection(obj::RegionFactory, arg0::Region, arg1::Region)
    return jcall(obj, "intersection", Region, (Region, Region), arg0, arg1)
end

function union(obj::RegionFactory, arg0::Region, arg1::Region)
    return jcall(obj, "union", Region, (Region, Region), arg0, arg1)
end

function xor(obj::RegionFactory, arg0::Region, arg1::Region)
    return jcall(obj, "xor", Region, (Region, Region), arg0, arg1)
end


text:
"""
    nearest_neighbor_count(mc, ϵ = 1e-6)

Determines the number of nearest neighbors by bond distance.
"""
function nearest_neighbor_count(mc, ϵ = 1e-6)
    l = lattice(mc)
    dirs = directions(l)
    @assert dirs[1] == [0.0, 0.0]

    sq_norm = dot(dirs[2], dirs[2])
    for i in 3:length(l)
        if !(sq_norm * (1.0 - ϵ)  <= dot(dirs[i], dirs[i]) <= sq_norm * (1.0 + ϵ))
            return i-2
        end
    end
    return length(l)-1
end

################################################################################
### Measurement constructors
################################################################################

function checkflavors(model, N=2)
    if nflavors(model) != N
        @warn(
            "$N flavors are required, but $(nflavors(model)) have been found"
        )
    end
    nothing
end

# This has lattice_iteratorator = Nothing, because it straight up copies G
function greens_measurement(
        mc::DQMC, model::Model, greens_iterator = Greens();
        capacity = _default_capacity(mc), eltype = geltype(mc),
        obs = let
            N = length(lattice(model)) * nflavors(model)
            LogBinner(zeros(eltype, (N, N)), capacity=capacity)
        end, kwargs...
    )
    Measurement(
        mc, model, greens_iterator, nothing, greens_kernel,
        obs = obs; kwargs...
    )
end

function occupation(
        mc::DQMC, model::Model; wrapper = nothing,
        lattice_iterator = EachSiteAndFlavor(), kwargs...
    )
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    Measurement(mc, model, Greens(), li, occupation_kernel; kwargs...)
end

function charge_density(
        mc::DQMC, model::Model, greens_iterator;
        wrapper = nothing, lattice_iterator = EachSitePairByDistance(), kwargs...
    )
    checkflavors(model)
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    Measurement(mc, model, greens_iterator, li, cdc_kernel; kwargs...)
end

charge_density_correlation(mc, m; kwargs...) = charge_density(mc, m, Greens(); kwargs...)
charge_density_susceptibility(mc, m; kwargs...) = charge_density(mc, m, CombinedGreensIterator(mc); kwargs...)

"""
    magnetization(mc, model, dir[; lattice_iterator, kwargs...])

Returns the x, y or z magnetization measurement given `dir = :x`, `:y` or `:z`
respectively.

NOTE:
We're skipping the multiplication by `-1im` during the measurement of the y
magnetization. To get the correct result, multiply the final result by `-1im`.
"""
function magnetization(
        mc::DQMC, model::Model, dir::Symbol;
        wrapper = nothing, lattice_iterator = EachSite(), kwargs...
    )
    checkflavors(model)
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    if dir == :x;
        return Measurement(mc, model, Greens(), li, mx_kernel; kwargs...)
    elseif dir == :y;
        return Measurement(mc, model, Greens(), li, my_kernel; kwargs...)
    elseif dir == :z;
        return Measurement(mc, model, Greens(), li, mz_kernel; kwargs...)
    else throw(ArgumentError("`dir` must be :x, :y or :z, but is $dir"))
    end

end

function spin_density(
        dqmc, model, dir::Symbol, greens_iterator;
        wrapper = nothing, lattice_iterator = EachSitePairByDistance(), kwargs...
    )
    checkflavors(model)
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    dir in (:x, :y, :z) || throw(ArgumentError("`dir` must be :x, :y or :z, but is $dir"))
    if     dir == :x
        return Measurement(dqmc, model, greens_iterator, li, sdc_x_kernel; kwargs...)
    elseif dir == :y
        return Measurement(dqmc, model, greens_iterator, li, sdc_y_kernel; kwargs...)
    elseif dir == :z
        return Measurement(dqmc, model, greens_iterator, li, sdc_z_kernel; kwargs...)
    else
        error("$greens_iterator not recongized")
    end
end
spin_density_correlation(args...; kwargs...) = spin_density(args..., Greens(); kwargs...)
spin_density_susceptibility(mc, args...; kwargs...) = spin_density(mc, args..., CombinedGreensIterator(mc); kwargs...)

function pairing(
        dqmc::DQMC, model::Model, greens_iterator;
        K = 1 + nearest_neighbor_count(dqmc), wrapper = nothing,
        lattice_iterator = EachLocalQuadByDistance(1:K),
        kernel = pc_combined_kernel, kwargs...
    )
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    Measurement(dqmc, model, greens_iterator, li, kernel; kwargs...)
end
pairing_correlation(mc, m; kwargs...) = pairing(mc, m, Greens(); kwargs...)
pairing_susceptibility(mc, m; kwargs...) = pairing(mc, m, CombinedGreensIterator(mc); kwargs...)

"""
    current_current_susceptibility(dqmc, model[; K = #NN + 1, kwargs...])

Creates a measurement recordering the current-current susceptibility
`OBS[dr, dr1, dr2] = Λ_{dr1, dr2}(dr) = ∫ 1/N \\sum_ ⟨j_{dr2}(i+dr, τ) j_{dr1}(i, 0)⟩ dτ`.

By default this measurement synchronizes `dr1` and `dr2`, simplifying the
formula to
`OBS[dr, dr'] = Λ_{dr'}(dr) = ∫ 1/N \\sum_ ⟨j_{dr'}(i+dr, τ) j_{dr'}(i, 0)⟩ dτ`.
The direction/index `dr'` runs from `1:K` through the `directions(mc)` which are
sorted. By default we include on-site (1) and nearest neighbors (2:K).

The current density `j_{dr}(i, τ)` with `j = i + dr` is given by
`j_{j-i}(i, τ) = i T_{ij} c_i^†(τ) c_j(0) + h.c.`. Note that the Hermitian
conjugate generates a term in `-dr` direction with prefactor `-i`. If `+dr` and
`-dr` are included in the calculation of the superfluid density (etc) you may
see larger values than expected.

"""
function current_current_susceptibility(
        dqmc::DQMC, model::Model;
        K = 1 + nearest_neighbor_count(dqmc),
        greens_iterator = CombinedGreensIterator(dqmc), wrapper = nothing,
        lattice_iterator = EachLocalQuadByDistance(2:K), kwargs...
    )
    @assert is_approximately_hermitian(hopping_matrix(dqmc, model)) "CCS assumes Hermitian matrix"
    li = wrapper === nothing ? lattice_iterator : wrapper(lattice_iterator)
    Measurement(dqmc, model, greens_iterator, li, cc_kernel; kwargs...)
end

function boson_energy_measurement(dqmc, model; kwargs...)
    Measurement(dqmc, model, Nothing, nothing, energy_boson; kwargs...)
end

function noninteracting_energy(dqmc, model; kwargs...)
    Measurement(dqmc, model, Greens(), nothing, nonintE_kernel; kwargs...)
end

# These require the model to implement intE_kernel
function interacting_energy(dqmc, model; kwargs...)
    Measurement(dqmc, model, Greens(), nothing, intE_kernel; kwargs...)
end

function total_energy(dqmc, model; kwargs...)
    Measurement(dqmc, model, Greens(), nothing, totalE_kernel; kwargs...)
end

################################################################################
### Wicks expanded Kernel functions
################################################################################

# Notes on math and stuff:
#=
- Use Wicks theorem to go from one expectation value with 2N operators to
  N expectation values with 2 operators each
- use G_{ij} = ⟨c_i c_j^†⟩
- use ⟨c_i^† c_j⟩ = δ_{ij} - ⟨c_j c_i^†⟩ = (I - G)_{ji} = swapop(::GreensMatrix)[i, j]
- we define G[i,j] as spin up, G[i+N,j+N] as spin down (N = number of sites)
- we used (I - G)[i, j+N] = G[i, j+N] etc when possible
- Gl0_{i,j} = ⟨c_i(l) c_j(0)^†⟩
- ⟨c_i(0) c_j(l)^†⟩ = ⟨(c_j(l) c_i(0)^†)^†⟩ = (⟨c_j(l) c_i(0)^†⟩)^†
                    = Gl0_{j,i}^† = Gl0_{j,i}^*
=#

"""
    greens_kernel(mc, model, G::GreensMatrix)

Returns the unprocessed Greens function `greens(mc) = {⟨cᵢcⱼ^†⟩}`.

* Lattice Iterators: `nothing` (zero index)
* Greens Iterators: `Greens` or `GreensAt`
"""
greens_kernel(mc, model, G::GreensMatrix) = G.val

"""
    occupation_kernel(mc, model, i::Integer, G::GreensMatrix)

Returns the per index occupation `⟨nᵢ⟩`.

* Lattice Iterators: `EachSiteAndFlavor`, `EachSite`
* Greens Iterators: `Greens` or `GreensAt`
"""
occupation_kernel(mc, model, i::Integer, G::GreensMatrix) = 1 - G[i, i]

"""
    cdc_kernel(mc, model, ij::NTuple{2, Integer}, G::GreensMatrix)
    cdc_kernel(mc, model, ij::NTuple{2, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair charge density `⟨nᵢ(τ) nⱼ(0)⟩ - ⟨nᵢ(τ)⟩⟨nⱼ(0)⟩` where
`τ = 0` for the first signature.

* Lattice Iterators: `OnSite`, `EachSitePair` or `EachSitePairByDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function cdc_kernel(mc, model, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    N = length(lattice(mc))
    # ⟨n↑n↑⟩
    swapop(G)[i, i] * swapop(G)[j, j] +
    swapop(G)[i, j] * G[i, j] +
    # ⟨n↑n↓⟩
    swapop(G)[i, i] * swapop(G)[j+N, j+N] +
    swapop(G)[i, j+N] * G[i, j + N] +
    # ⟨n↓n↑⟩
    swapop(G)[i+N, i+N] * swapop(G)[j, j] +
    swapop(G)[i+N, j] * G[i+N, j] +
    # ⟨n↓n↓⟩
    swapop(G)[i+N, i+N]  * swapop(G)[j+N, j+N] +
    swapop(G)[i+N, j+N] * G[i+N, j+N]
end

function cdc_kernel(mc, model, ij::NTuple{2}, packed_greens::NTuple{4})
    i, j = ij
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(mc))
    # ⟨n↑(l)n↑(0)⟩
    swapop(Gll)[i, i] * swapop(G00)[j, j] +
    swapop(G0l)[i, j] * Gl0[i, j] +
    # ⟨n↑(l)n↓(0)⟩
    swapop(Gll)[i, i] * swapop(G00)[j+N, j+N] +
    swapop(G0l)[i, j+N] * Gl0[i, j + N] +
    # ⟨n↓(l)n↑(0)⟩
    swapop(Gll)[i+N, i+N] * swapop(G00)[j, j] +
    swapop(G0l)[i+N, j] * Gl0[i+N, j] +
    # ⟨n↓(l)n↓(0)⟩
    swapop(Gll)[i+N, i+N] * swapop(G00)[j+N, j+N] +
    swapop(G0l)[i+N, j+N] * Gl0[i+N, j+N]
end

"""
    mx_kernel(mc, model, i::Integer, G::GreensMatrix)

Returns the per-site x-magnetization `⟨cᵢ↑^† cᵢ↓ + cᵢ↓^† cᵢ↑⟩`.

* Lattice Iterators: `EachSite`
* Greens Iterators: `Greens` or `GreensAt`
"""
function mx_kernel(mc, model, i, G::GreensMatrix)
    N = length(lattice(model))
    -G[i+N, i] - G[i, i+N]
end

"""
    my_kernel(mc, model, i::Integer, G::GreensMatrix)

Returns the per-site y-magnetization `-⟨cᵢ↑^† cᵢ↓ - cᵢ↓^† cᵢ↑⟩` without the
imaginary prefactor.

* Lattice Iterators: `EachSite`
* Greens Iterators: `Greens` or `GreensAt`
"""
function my_kernel(mc, model, i, G::GreensMatrix)
    N = length(lattice(model))
    G[i+N, i] - G[i, i+N]
end

"""
    mz_kernel(mc, model, i::Integer, G::GreensMatrix)

Returns the per-site z-magnetization `⟨nᵢ↑ - nᵢ↓⟩`.

* Lattice Iterators: `EachSite`
* Greens Iterators: `Greens` or `GreensAt`
"""
function mz_kernel(mc, model, i, G::GreensMatrix)
    N = length(lattice(model))
    G[i+N, i+N] - G[i, i]
end

"""
    sdc_x_kernel(mc, model, ij::NTuple{2, Integer}, G::GreensMatrix)
    sdc_x_kernel(mc, model, ij::NTuple{2, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair x-spin density `⟨mxᵢ(τ) mxⱼ(0)⟩ - ⟨mxᵢ(τ)⟩⟨mxⱼ(0)⟩`
where `τ = 0` for the first signature.

* Lattice Iterators: `OnSite`, `EachSitePair` or `EachSitePairByDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function sdc_x_kernel(mc, model, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    N = length(lattice(model))
    G[i+N, i] * G[j+N, j] + G[i+N, i] * G[j, j+N] +
    G[i, i+N] * G[j+N, j] + G[i, i+N] * G[j, j+N] +
    swapop(G)[i, j+N] * G[i+N, j] + swapop(G)[i, j] * G[i+N, j+N] +
    swapop(G)[i+N, j+N] * G[i, j] + swapop(G)[i+N, j] * G[i, j+N]
end
function sdc_x_kernel(mc, model, ij::NTuple{2}, packed_greens::NTuple{4})
    i, j = ij
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    Gll[i+N, i] * G00[j+N, j] + Gll[i+N, i] * G00[j, j+N] +
    Gll[i, i+N] * G00[j+N, j] + Gll[i, i+N] * G00[j, j+N] +
    swapop(G0l)[i, j+N] * Gl0[i+N, j] + swapop(G0l)[i, j] * Gl0[i+N, j+N] +
    swapop(G0l)[i+N, j+N] * Gl0[i, j] + swapop(G0l)[i+N, j] * Gl0[i, j+N]
end

"""
    sdc_y_kernel(mc, model, ij::NTuple{2, Integer}, G::GreensMatrix)
    sdc_y_kernel(mc, model, ij::NTuple{2, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair x-spin density `⟨myᵢ(τ) myⱼ(0)⟩ - ⟨myᵢ(τ)⟩⟨myⱼ(0)⟩`
where `τ = 0` for the first signature.

* Lattice Iterators: `OnSite`, `EachSitePair` or `EachSitePairByDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function sdc_y_kernel(mc, model, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    N = length(lattice(model))
    - G[i+N, i] * G[j+N, j] + G[i+N, i] * G[j, j+N] +
      G[i, i+N] * G[j+N, j] - G[i, i+N] * G[j, j+N] +
    - swapop(G)[i, j+N] * G[i+N, j] + swapop(G)[i, j] * G[i+N, j+N] +
      swapop(G)[i+N, j+N] * G[i, j] - swapop(G)[i+N, j] * G[i, j+N]
end
function sdc_y_kernel(mc, model, ij::NTuple{2}, packed_greens::NTuple{4})
    i, j = ij
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    - Gll[i+N, i] * G00[j+N, j] + Gll[i+N, i] * G00[j, j+N] +
      Gll[i, i+N] * G00[j+N, j] - Gll[i, i+N] * G00[j, j+N] +
    - swapop(G0l)[i, j+N] * Gl0[i+N, j] + swapop(G0l)[i, j] * Gl0[i+N, j+N] +
      swapop(G0l)[i+N, j+N] * Gl0[i, j] - swapop(G0l)[i+N, j] * Gl0[i, j+N]
end

"""
    sdc_z_kernel(mc, model, ij::NTuple{2, Integer}, G::GreensMatrix)
    sdc_z_kernel(mc, model, ij::NTuple{2, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair x-spin density `⟨mzᵢ(τ) mzⱼ(0)⟩ - ⟨mzᵢ(τ)⟩⟨mzⱼ(0)⟩`
where `τ = 0` for the first signature.

* Lattice Iterators: `OnSite`, `EachSitePair` or `EachSitePairByDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function sdc_z_kernel(mc, model, ij::NTuple{2}, G::GreensMatrix)
    i, j = ij
    N = length(lattice(model))
    swapop(G)[i, i]     * swapop(G)[j, j] -
    swapop(G)[i, i]     * swapop(G)[j+N, j+N] -
    swapop(G)[i+N, i+N] * swapop(G)[j, j] +
    swapop(G)[i+N, i+N] * swapop(G)[j+N, j+N] +
    swapop(G)[i, j] * G[i, j]     - swapop(G)[i, j+N] * G[i, j+N] -
    swapop(G)[i+N, j] * G[i+N, j] + swapop(G)[i+N, j+N] * G[i+N, j+N]
end
function sdc_z_kernel(mc, model, ij::NTuple{2}, packed_greens::NTuple{4})
    i, j = ij
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    swapop(Gll)[i, i] * swapop(G00)[j, j] -
    swapop(Gll)[i, i] * swapop(G00)[j+N, j+N] -
    swapop(Gll)[i+N, i+N] * swapop(G00)[j, j] +
    swapop(Gll)[i+N, i+N] * swapop(G00)[j+N, j+N] +
    swapop(G0l)[i, j] * Gl0[i, j]     - swapop(G0l)[i, j+N] * Gl0[i, j+N] -
    swapop(G0l)[i+N, j] * Gl0[i+N, j] + swapop(G0l)[i+N, j+N] * Gl0[i+N, j+N]
end

"""
    pc_kernel(mc, model, ij::NTuple{4, Integer}, G::GreensMatrix)
    pc_kernel(mc, model, ij::NTuple{4, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair pairing `⟨Δᵢⱼ(τ) Δₖₗ^†(0)⟩` where `τ = 0` for the
first signature and `Δᵢⱼ(τ) = 0.5 cᵢ↑(τ) cⱼ↓(τ)`. The Delta operator is a
deferred version of the  pair field operator `Δᵢ = 0.5 ∑ₐ f(a) cᵢ↑ cᵢ₊ₐ↓` which
leaves the execution of the sum for after the simulation.

* Lattice Iterators: `EachLocalQuadByDistance` or `EachLocalQuadBySyncedDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function pc_kernel(mc, model, sites::NTuple{4}, G::GreensMatrix)
    pc_kernel(mc, model, sites, (G, G, G, G))
end
function pc_kernel(mc, model, sites::NTuple{4}, packed_greens::NTuple{4})
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    # Δ_v(src1, trg1)(τ) Δ_v^†(src2, trg2)(0)
    # G_{i, j}^{↑, ↑}(τ, 0) G_{i+d, j+d'}^{↓, ↓}(τ, 0) -
    # G_{i, j+d'}^{↑, ↓}(τ, 0) G_{i+d, j}^{↓, ↑}(τ, 0)
    Gl0[src1, src2] * Gl0[trg1+N, trg2+N] - Gl0[src1, trg2+N] * Gl0[trg1+N, src2]
end

"""
    pc_alt_kernel(mc, model, ij::NTuple{4, Integer}, G::GreensMatrix)
    pc_alt_kernel(mc, model, ij::NTuple{4, Integer}, Gs::NTuple{4, GreensMatrix})

Returns the per-site-pair pairing `⟨Δᵢⱼ^†(τ) Δₖₗ(0)⟩` where `τ = 0` for the
first signature and `Δᵢⱼ(τ) = 0.5 cᵢ↑(τ) cⱼ↓(τ)`. The Delta operator is a
deferred version of the  pair field operator `Δᵢ = 0.5 ∑ₐ f(a) cᵢ↑ cᵢ₊ₐ↓` which
leaves the execution of the sum for after the simulation.

* Lattice Iterators: `EachLocalQuadByDistance` or `EachLocalQuadBySyncedDistance`
* Greens Iterators: `Greens`, `GreensAt` or `CombinedGreensIterator`
"""
function pc_alt_kernel(mc, model, sites::NTuple{4}, G::GreensMatrix)
    pc_alt_kernel(mc, model, sites, (G, G, G, G))
end
function pc_alt_kernel(mc, model, sites::NTuple{4}, packed_greens::NTuple{4})
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    # Δ_v^†(src1, trg1)(τ) Δ_v(src2, trg2)(0)
    # (I-G)_{j, i}^{↑, ↑}(0, τ) (I-G)_{j+d', i+d}^{↓, ↓}(0, τ) -
    # (I-G)_{j, i+d}^{↑, ↓}(0, τ) G_{j+d', i}^{↓, ↑}(0, τ)
    swapop(G0l)[trg1+N, trg2+N] * swapop(G0l)[src1, src2] -
    G0l[src2, trg1+N] * G0l[trg2+N, src1]
end

function pc_combined_kernel(mc, model, sites::NTuple{4}, G)
    # Δ^† Δ + Δ Δ^†
    # same as in https://arxiv.org/pdf/1912.08848.pdf
    pc_kernel(mc, model, sites, G) + pc_alt_kernel(mc, model, sites, G)
end

function pc_ref_kernel(mc, model, sites::NTuple{4}, G::GreensMatrix)
    # Δ^† Δ + Δ Δ^† but ↑ and ↓ are swapped
    pc_ref_kernel(mc, model, sites, (G, G, G, G))
end
function pc_ref_kernel(mc, model, sites::NTuple{4}, packed_greens::NTuple{4})
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    Gl0[src1+N, src2+N] * Gl0[trg1, trg2] -
    Gl0[src1+N, trg2] * Gl0[trg1, src2+N] +
    (I[trg2, trg1] - G0l[trg2, trg1]) * (I[src2, src1] - G0l[src2+N, src1+N]) -
    (I[src2, trg1] - G0l[src2+N, trg1]) * (I[trg2, src1] - G0l[trg2, src1+N])
end

function cc_kernel(mc, model, sites::NTuple{4}, packed_greens::NTuple{4})
    # Computes
    # ⟨j_{t2-s2}(s2, l) j_{t1-s1}(s1, 0)⟩
    # where t2-s2 (t1-s1) is usually a NN vector/jump, and
    # j_{t2-s2}(s2, l) = i \sum_σ [T_{ts} c_t^†(l) c_s(τ) - T_{st} c_s^†(τ) c_t(τ)]
    src1, trg1, src2, trg2 = sites
	G00, G0l, Gl0, Gll = packed_greens
    N = length(lattice(model))
    T = mc.stack.hopping_matrix
    output = zero(eltype(G00))

    # Iterate through (spin up, spin down)
    for σ1 in (0, N), σ2 in (0, N)
        s1 = src1 + σ1; t1 = trg1 + σ1
        s2 = src2 + σ2; t2 = trg2 + σ2
        # Note: if H is real and Hermitian, T can be pulled out and the I's cancel
        # Note: This matches crstnbr/dqmc if H real, Hermitian
        # Note: I for G0l and Gl0 do not always cancel
        # I have a tex document for this now
        output += (
            (T[s2, t2] * (I[t2, s2] - Gll[t2, s2]) - T[t2, s2] * (I[t2, s2] - Gll[s2, t2])) *
            (T[t1, s1] * (I[s1, t1] - G00[s1, t1]) - T[s1, t1] * (I[s1, t1] - G00[t1, s1])) +
            - T[t2, s2] * T[t1, s1] * swapop(G0l)[t2, s1] * Gl0[s2, t1] +
            + T[t2, s2] * T[s1, t1] * swapop(G0l)[t2, t1] * Gl0[s2, s1] +
            + T[s2, t2] * T[t1, s1] * swapop(G0l)[s2, s1] * Gl0[t2, t1] +
            - T[s2, t2] * T[s1, t1] * swapop(G0l)[s2, t1] * Gl0[t2, s1]
        )
    end

    output
end

@inline function nonintE_kernel(mc, model, G::GreensMatrix)
    # <T> = \sum Tji * (Iij - Gij) = - \sum Tji * (Gij - Iij)
    T = mc.stack.hopping_matrix
    nonintE(T, G.val)
end

# TODO should this be moved/work differently?
nonintE(T::AbstractArray, G::GreensMatrix) = nonintE(T, G.val)
function nonintE(T::AbstractArray, G::AbstractArray)
    output = zero(eltype(G))
    for i in axes(G, 1), j in axes(G, 2)
        output += T[j, i] * (I[i, j] - G[i, j])
    end
    # 2 because we're using spin up/down symmetry
    2.0 * output
end
function nonintE(T::BlockDiagonal{X, N}, G::BlockDiagonal{X, N}) where {X, N}
    output = zero(eltype(G))
    @inbounds n = size(T.blocks[1], 1)
    @inbounds for i in 1:N
        t = T.blocks[i]
        g = G.blocks[i]
        @turbo for k in 1:n, l in 1:n
            output += t[k,l] * (ifelse(k==l, 1.0, 0.0) - g[k,l])
        end
    end
    output
end

function totalE_kernel(mc, model, G::GreensMatrix)
    nonintE_kernel(mc, model, G) + intE_kernel(mc, model, G)
end


text:
module TestDDataFrame
	using Base.Test
	using DataArrays
	using DataFrames

	const datafile = joinpath(dirname(@__FILE__), "data", "distributed", "test.csv")
	const nloops = 10

	function load_pkgs()
	    println("loading packages...")
	    @everywhere using Blocks
	    @everywhere using DataFrames
	end

	if nprocs() < 4
	    addwrkrs = 4 - nprocs()
	    println("adding $addwrkrs more processors...")
	    addprocs(addwrkrs)
	end
	println("\tnprocs: $(nprocs())")
	load_pkgs()

	df = dreadtable(datafile, header=false)
	names(df)
	names!(df, ["c1","c2","c3","c4","c5","c6","c7","c8","c9","c10"])

	sum_result = df+df
	mul_result = 2*df
	eq_result = (sum_result .== mul_result)
	@assert all(eq_result)

	df1 = dreadtable(open(datafile), 1000)
	@assert nrow(df1) == nrow(df)
	@assert ncol(df1) == ncol(df)
	@assert isapprox(sum(matrix(colsums(df1))), sum(matrix(colsums(df))))
end


text:
export
ruleVBGaussianMixtureM,
ruleVBGaussianMixtureW,
ruleVBGaussianMixtureZBer,
ruleVBGaussianMixtureZCat,
ruleVBGaussianMixtureOut

function ruleVBGaussianMixtureM(dist_out::ProbabilityDistribution,
                                dist_switch::ProbabilityDistribution,
                                dist_factors::Vararg{Union{Nothing, ProbabilityDistribution{Univariate}}})
    # Univariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])
    k = findfirst(dist_means .== nothing) # Find factor
    z_bar = clamp.(unsafeMeanVector(dist_switch), tiny, 1.0 - tiny)

    return Message(Univariate, GaussianMeanPrecision, m=unsafeMean(dist_out), w=z_bar[k]*unsafeMean(dist_precs[k]))
end

function ruleVBGaussianMixtureM(dist_out::ProbabilityDistribution,
                                dist_switch::ProbabilityDistribution,
                                dist_factors::Vararg{Union{Nothing, ProbabilityDistribution{Multivariate}, ProbabilityDistribution{MatrixVariate}}})
    # Multivariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])
    k = findfirst(dist_means .== nothing) # Find factor
    z_bar = clamp.(unsafeMeanVector(dist_switch), tiny, 1.0 - tiny)

    return Message(Multivariate, GaussianMeanPrecision, m=unsafeMean(dist_out), w=z_bar[k]*unsafeMean(dist_precs[k]))
end

function ruleVBGaussianMixtureW(dist_out::ProbabilityDistribution,
                                dist_switch::ProbabilityDistribution,
                                dist_factors::Vararg{Union{Nothing, ProbabilityDistribution{Univariate}}})
    # Univariate update
    dist_means = collect(dist_factors[1:2:end]) # TODO: make more efficient use of conversions to compute m and v in one go
    dist_precs = collect(dist_factors[2:2:end])
    k = findfirst(dist_precs .== nothing) # Find factor
    (m_mean_k, v_mean_k) = unsafeMeanCov(dist_means[k])
    (m_out, v_out) = unsafeMeanCov(dist_out)
    z_bar = unsafeMeanVector(dist_switch)

    return Message(Univariate, Gamma,
        a = 1.0 + 0.5*z_bar[k],
        b = 0.5*z_bar[k]*(v_out + v_mean_k + (m_out - m_mean_k)^2))
end

function ruleVBGaussianMixtureW(dist_out::ProbabilityDistribution,
                                dist_switch::ProbabilityDistribution,
                                dist_factors::Vararg{Union{Nothing, ProbabilityDistribution{Multivariate}, ProbabilityDistribution{MatrixVariate}}})
    # Multivariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])
    k = findfirst(dist_precs .== nothing) # Find factor
    (m_mean_k, v_mean_k) = unsafeMeanCov(dist_means[k])
    (m_out, v_out) = unsafeMeanCov(dist_out)
    z_bar = unsafeMeanVector(dist_switch)
    d = dims(dist_means[1])

    return Message(MatrixVariate, Wishart,
        nu = 1.0 + z_bar[k] + d,
        v  = cholinv(z_bar[k]*( v_out + v_mean_k + (m_out - m_mean_k)*(m_out - m_mean_k)')))
end

function softmax(v::Vector{Float64})
    r = v .- maximum(v)
    clamp!(r, -100.0, 0.0)
    exp.(r)./sum(exp.(r))
end

function ruleVBGaussianMixtureZBer( dist_out::ProbabilityDistribution,
                                    dist_switch::Any,
                                    dist_m1::ProbabilityDistribution,
                                    dist_w1::ProbabilityDistribution,
                                    dist_m2::ProbabilityDistribution,
                                    dist_w2::ProbabilityDistribution)
    # Uni- and Multivariate update
    U = Vector{Float64}(undef, 2)
    U[1] = averageEnergy(GaussianMeanPrecision, dist_out, dist_m1, dist_w1)
    U[2] = averageEnergy(GaussianMeanPrecision, dist_out, dist_m2, dist_w2)

    return Message(Univariate, Bernoulli, p=softmax(-U)[1])
end

function ruleVBGaussianMixtureZCat( dist_out::ProbabilityDistribution,
                                    dist_switch::Any,
                                    dist_factors::Vararg{ProbabilityDistribution})
    # Uni- and Multivariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])

    n_factors = length(dist_means)
    U = Vector{Float64}(undef, n_factors)
    for k = 1:n_factors
        U[k] = averageEnergy(GaussianMeanPrecision, dist_out, dist_means[k], dist_precs[k])
    end

    return Message(Univariate, Categorical, p=softmax(-U))
end

function ruleVBGaussianMixtureOut(  dist_out::Any,
                                    dist_switch::ProbabilityDistribution,
                                    dist_factors::Vararg{ProbabilityDistribution{Univariate}})
    # Univariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])
    z_bar = unsafeMeanVector(dist_switch)

    w  = 0.0
    xi = 0.0
    for k = 1:length(z_bar)
        w  += z_bar[k]*unsafeMean(dist_precs[k])
        xi += unsafeMean(dist_precs[k])*unsafeMean(dist_means[k])*z_bar[k]
    end

    return Message(Univariate, GaussianWeightedMeanPrecision, xi=xi, w=w)
end

function ruleVBGaussianMixtureOut(  dist_out::Any,
                                    dist_switch::ProbabilityDistribution,
                                    dist_factors::Vararg{Union{ProbabilityDistribution{Multivariate}, ProbabilityDistribution{MatrixVariate}}})
    # Multivariate update
    dist_means = collect(dist_factors[1:2:end])
    dist_precs = collect(dist_factors[2:2:end])
    z_bar = unsafeMeanVector(dist_switch)
    d = dims(dist_means[1])

    w = Diagonal(zeros(d))
    xi = zeros(d)
    for k = 1:length(z_bar)
        w  += z_bar[k]*unsafeMean(dist_precs[k])
        xi += unsafeMean(dist_precs[k])*unsafeMean(dist_means[k])*z_bar[k]
    end

    return Message(Multivariate, GaussianWeightedMeanPrecision, xi=xi, w=w)
end


text:
# This file is derived from parts of Julia's math tests from https://github.com/JuliaLang/julia/blob/master/test/math.jl

# Test math functions. We compare to BigFloat instead of hard-coding
# values, assuming that BigFloat has an independent and independently
# tested implementation.
for T in (Float32, Float64)
    x = T(1//3)
    y = T(1//2)
    yi = 4
    # Test random values
    @test xacos(x) ≈ acos(big(x))
    @test xacos_u1(x) ≈ acos(big(x))
    @test xacosh(1+x) ≈ acosh(big(1+x))
    @test xasin(x) ≈ asin(big(x))
    @test xasin_u1(x) ≈ asin(big(x))
    @test xasinh(x) ≈ asinh(big(x))
    @test xatan(x) ≈ atan(big(x))
    @test xatan_u1(x) ≈ atan(big(x))
    @test xatan2(x,y) ≈ atan2(big(x),big(y))
    @test xatan2_u1(x,y) ≈ atan2(big(x),big(y))
    @test atanh(x) ≈ atanh(big(x))
    @test xatanh(x) ≈ atanh(big(x))
    @test xcbrt(x) ≈ cbrt(big(x))
    @test xcos(x) ≈ cos(big(x))
    @test xcos_u1(x) ≈ cos(big(x))
    @test xcosh(x) ≈ cosh(big(x))
    @test xexp(x) ≈ exp(big(x))
    @test xexp10(x) ≈ exp10(big(x))
    @test xexp2(x) ≈ exp2(big(x))
    @test xexpm1(x) ≈ expm1(big(x))
    @test xlog(x) ≈ log(big(x))
    @test xlog_u1(x) ≈ log(big(x))
    @test xlog10(x) ≈ log10(big(x))
    @test xlog1p(x) ≈ log1p(big(x))
    @test xsin(x) ≈ sin(big(x))
    @test xsin_u1(x) ≈ sin(big(x))
    @test xsinh(x) ≈ sinh(big(x))
    @test xtan(x) ≈ tan(big(x))
    @test xtan_u1(x) ≈ tan(big(x))
    @test tanh(x) ≈ tanh(big(x))

    # Test special values
    @test isequal(xacos(T(1)), T(0))
    @test isequal(xacosh(T(1)), T(0))
    @test_approx_eq_eps xasin(T(1)) T(pi)/2 eps(T)
    @test_approx_eq_eps xasin_u1(T(1)) T(pi)/2 eps(T)
    @test_approx_eq_eps xatan(T(1)) T(pi)/4 eps(T)
    @test_approx_eq_eps xatan_u1(T(1)) T(pi)/4 eps(T)
    @test_approx_eq_eps xatan2(T(1),T(1)) T(pi)/4 eps(T)
    @test_approx_eq_eps xatan2_u1(T(1),T(1)) T(pi)/4 eps(T)
    @test isequal(xcbrt(T(0)), T(0))
    @test isequal(xcbrt(T(1)), T(1))
    @test isequal(xcbrt(T(1000000000)), T(1000))
    @test isequal(xcos(T(0)), T(1))
    @test isequal(xcos_u1(T(0)), T(1))
    @test_approx_eq_eps xcos(T(pi)/2) T(0) eps(T)
    @test_approx_eq_eps xcos_u1(T(pi)/2) T(0) eps(T)
    @test isequal(xcos(T(pi)), T(-1))
    @test isequal(xcos_u1(T(pi)), T(-1))
    @test_approx_eq_eps xexp(T(1)) T(e) 10*eps(T)
    @test isequal(xexp10(T(1)), T(10))
    @test isequal(xexp2(T(1)), T(2))
    @test isequal(xexpm1(T(0)), T(0))
    @test_approx_eq_eps xexpm1(T(1)) T(e)-1 10*eps(T)
    @test isequal(xlog(T(1)), T(0))
    @test isequal(xlog_u1(T(1)), T(0))
    @test_approx_eq_eps xlog(T(e)) T(1) eps(T)
    @test_approx_eq_eps xlog_u1(T(e)) T(1) eps(T)
    @test isequal(xlog10(T(1)), T(0))
    @test isequal(xlog10(T(10)), T(1))
    @test isequal(xlog1p(T(0)), T(0))
    @test_approx_eq_eps xlog1p(T(e)-1) T(1) eps(T)
    @test isequal(xsin(T(0)), T(0))
    @test isequal(xsin_u1(T(0)), T(0))
    @test isequal(xsin(T(pi)/2), T(1))
    @test isequal(xsin_u1(T(pi)/2), T(1))
    @test_approx_eq_eps xsin(T(pi)) T(0) eps(T)
    @test_approx_eq_eps xsin_u1(T(pi)) T(0) eps(T)
    @test isequal(xtan(T(0)), T(0))
    @test isequal(xtan_u1(T(0)), T(0))
    @test_approx_eq_eps xtan(T(pi)/4) T(1) eps(T)
    @test_approx_eq_eps xtan_u1(T(pi)/4) T(1) eps(T)

    # Test inverses
    @test xacos(xcos(x)) ≈ x
    @test xacos_u1(xcos_u1(x)) ≈ x
    @test xacosh(xcosh(x)) ≈ x
    @test xasin(xsin(x)) ≈ x
    @test xasin_u1(xsin_u1(x)) ≈ x
    @test xcbrt(x)^3 ≈ x
    @test xcbrt_u1(x)^3 ≈ x
    @test xcbrt_u1(x^3) ≈ x
    @test xasinh(xsinh(x)) ≈ x
    @test xatan(xtan(x)) ≈ x
    @test xatan_u1(xtan_u1(x)) ≈ x
    @test xatan2(x,y) ≈ xatan(x/y)
    @test xatan2_u1(x,y) ≈ xatan_u1(x/y)
    @test xatanh(xtanh(x)) ≈ x
    @test xcos(xacos(x)) ≈ x
    @test xcos_u1(xacos_u1(x)) ≈ x
    @test xcosh(xacosh(1+x)) ≈ 1+x
    @test xexp(log(x)) ≈ x
    @test xexp10(log10(x)) ≈ x
    @test xexp2(log2(x)) ≈ x
    @test xexpm1(log1p(x)) ≈ x
    @test xlog(xexp(x)) ≈ x
    @test xlog10(xexp10(x)) ≈ x
    @test xlog1p(xexpm1(x)) ≈ x
    @test xsin(xasin(x)) ≈ x
    @test xsin_u1(xasin_u1(x)) ≈ x
    @test xsinh(xasinh(x)) ≈ x
    @test xtan(xatan(x)) ≈ x
    @test xtan_u1(xatan_u1(x)) ≈ x
    @test xtanh(xatanh(x)) ≈ x

    # Test some properties
    @test xcosh(x) ≈ (xexp(x)+xexp(-x))/2
    @test xcosh(x)^2-xsinh(x)^2 ≈ 1
    @test xsin(x)^2+xcos(x)^2 ≈ 1
    @test xsin_u1(x)^2+xcos_u1(x)^2 ≈ 1
    @test xsinh(x) ≈ (xexp(x)-xexp(-x))/2
    @test xtan(x) ≈ xsin(x)/xcos(x)
    @test xtan_u1(x) ≈ xsin_u1(x)/xcos_u1(x)
    @test xtanh(x) ≈ xsinh(x)/xcosh(x)

    #Edge cases
    @test isinf(xlog(zero(T)))
    @test isnan(xlog(convert(T,NaN)))
    @test isinf(xlog1p(-one(T)))
    @test isnan(xlog1p(convert(T,NaN)))
	@test_throws DomainError xlog(-one(T))
	@test_throws DomainError xlog1p(convert(T,-2.0))

end


text:
module taxis
using sequenceCompare
using DataFrames
using JSON
using Iterators

export Point, SequenceRef, LoadData, GetUniqueCoords, GetTableOrderedSubset, ConstructCoordsDatabase, ConstructCoordsDatabase2, PlotPaths, findClosestTrainingExample, findClosestTrainingExampleForTestSet

    immutable Point{T}
      x::T
      y::T
    end

    type SequenceRef
        TripId::Int64
        SequenceIndex::Int64
        LengthRemaining::Int64
    end

    function LoadData(train_path, test_path)
        println("Begin")

        println("loading csv files")
        taxi_df = readtable(train_path)
        taxi_validation_df = readtable(test_path)

        println("loading coords")
        taxi_df[:COORDS] = [float(hcat(JSON.parse(x)...)) for x in  taxi_df[:POLYLINE]]
        taxi_validation_df[:COORDS] = [float(hcat(JSON.parse(x)...)) for x in taxi_validation_df[:POLYLINE]]

        println("getting coords counts")
        taxi_df[:NUM_COORDS] = [int(length(x)/2)::Int64 for x in taxi_df[:COORDS]]
        taxi_validation_df[:NUM_COORDS] = [int(length(x)/2)::Int64 for x in taxi_validation_df[:COORDS]]

        println("deleting unneeded data rows/columns")
        delete!(taxi_validation_df, :POLYLINE)
        delete!(taxi_df, :POLYLINE)

        taxi_df[:START] = [x[:,1] for x in taxi_df[:COORDS]]
        taxi_validation_df[:START] = [x[:,1] for x in taxi_validation_df[:COORDS]]

        taxi_df[:END] = [x[:,end] for x in taxi_df[:COORDS]]
        taxi_validation_df[:END] = [x[:,end] for x in taxi_validation_df[:COORDS]]

        #These examples are not going to be useful!
        deleterows!(taxi_df, find(taxi_df[:NUM_COORDS] .== 0))

    taxi_df[:COORDS_TEST] = [c[:,1:int(floor((0.3 + 0.7*rand(1)[1])*size(c,2)))] for c in taxi_df[:COORDS]]

        println("done!")

        return taxi_df, taxi_validation_df
    end

    function GetUniqueCoords(c, round_len=4)
        c = round(c, round_len)
        c2 = Set([(c[i,1]::Float64, c[i,2]::Float64) for i in 1:size(c, 1)])
        println("all coords: ", length(c), ", unique coords: ", length(c2))
        return c2
    end

    function GetTableOrderedSubset(df, num_records=100000)
        idxs = sort(shuffle([1:size(df,1)])[1:num_records])
        return df[idxs,:]
    end

    # Constructing a coords database

    function ConstructCoordsDatabase(df, round_len=5)
        print(length(df[:TRIP_ID]), ", ", length(df[:COORDS]))
        return ConstructCoordsDatabase2(df[:TRIP_ID], df[:COORDS], round_len)
    end

    function ConstructCoordsDatabase2(TripIds, Coords, round_len=5)
        #println("num trip_ids/coords: ", length(TripIds), " ", length(Coords))
        #print("coords:", Coords)
        D = Dict{(Float64,Float64),Array{SequenceRef,1}}()
        num_coords = size(TripIds,1)
        for i=1:num_coords
            trip_id = TripIds[i]

            coords = Coords[i]
            if length(coords) == 0
                continue
            end

            coords = round(Coords[i],round_len)

            num_coords = size(coords,2)
            #println("num_coords: ", num_coords)
            for j=1:num_coords
                coord_pair = (coords[1,j], coords[2,j])
                #println("coord_pair: ", coord_pair, ", j/num_coords: ", j, "/", num_coords, ", thing: ", [SequenceRef(trip_id, j, num_coords-j)])
                if !haskey(D, coord_pair)
                    D[coord_pair] = [SequenceRef(trip_id, j, num_coords-j)]
                else
                    push!(D[coord_pair], SequenceRef(trip_id, j, num_coords-j))
                end
            end
        end

        return D
    end

    # Finding similar paths using Dynamic Path Warping

    function findClosestTrainingExample(all_train_coords, test_path, w=1)
        num_paths = length(all_train_coords)
        best_dist = 9999.0
        best_path = all_train_coords[1]
        for k=1:num_paths
            train_path = all_train_coords[k]

            if size(train_path, 2) < size(test_path, 2)
                continue
            end

            if w == 1
                dist = sequenceCompare.DTWDistance(train_path, test_path)
            else
                dist = sequenceCompare.DTWDistance(train_path, test_path, w)
            end
            #println("k=", k, " dist=", dist)
            if dist < best_dist
                #println(k, ", old best: ", best_dist, " new best: ", dist)
                best_dist = dist
                best_path = all_train_coords[k]
            end
        end

        return best_path
    end

    function findClosestTrainingExampleForTestSet(all_train_paths, all_test_paths, w=1)
        num_train_paths = length(all_train_paths)
        num_test_paths = length(all_test_paths)
        closest_examples = cell(num_test_paths)
        for k=1:num_test_paths
            if k % 20 == 0
                println(k, "/", num_test_paths, " for ", num_train_paths, " train path examples")
            end
            test_path = all_test_paths[k]
            closest_training_example = findClosestTrainingExample(all_train_paths, test_path, w)
            closest_examples[k] = closest_training_example
        end

        return closest_examples
    end

    function pFindClosestTrainingExampleForTestSet(all_train_paths, all_test_paths)
        getClosestExample = p -> findClosestTrainingExample(all_train_paths, p)
        return pmap(getClosestExample, all_test_paths)
    end

    # Didn't want to throw these away..
    function GetAverageDistanceFromCoordsSequence(coords)
        #Outline of the basic algorithm
        #for coord_pair c_1..c_n
        # estimate <- for i=1:n average(15*(i-1) + average(for c_j in D[c+i] * exp(-|L_j-L_i|) * L_j)
        coords = round(Coords[i],round_len)
    end

    function GetAverageDistancesByCoordsDict(TripIds, Coords, round_len=5)
        distaces = Array{Float64,0}()
        num_coords = size(TripIds,1)
        for i=1:num_coords
            trip_id = TripIds[i]

            coords = Coords[i]
            if length(coords) == 0
                continue
            end

            coords = round(Coords[i],round_len)

            num_coords = size(coords,2)
            #println("num_coords: ", num_coords)
            for j=1:num_coords
                coord_pair = (coords[1,j], coords[2,j])
                #println("coord_pair: ", coord_pair, ", j/num_coords: ", j, "/", num_coords, ", thing: ", [SequenceRef(trip_id, j, num_coords-j)])
                if !haskey(D, coord_pair)
                    D[coord_pair] = [SequenceRef(trip_id, j, num_coords-j)]
                else
                    push!(D[coord_pair], SequenceRef(trip_id, j, num_coords-j))
                end
            end
        end

        return D
    end

# ******************** Visualizations ******************************

    function PlotPaths(coords, show_limit=40)
        num_paths = min(length(coords), show_limit)
        num_coords = length(coords)
        coordsToPlot = coords[randperm(num_coords)[1:num_paths]]
        layers = [layer(x=round(c'[:,1],3),y=round(c'[:,2],3),Geom.point, Geom.path, Theme(default_color=RGB(rand(3)...))) for c in coordsToPlot]
        Gadfly.plot(layers...)
    end

end


text:
export AbstractStage,
    PreExperimentStage,
    PostExperimentStage,
    PreEpisodeStage,
    PostEpisodeStage,
    PreActStage,
    PostActStage,
    PRE_EXPERIMENT_STAGE,
    POST_EXPERIMENT_STAGE,
    PRE_EPISODE_STAGE,
    POST_EPISODE_STAGE,
    PRE_ACT_STAGE,
    POST_ACT_STAGE

#####
# Stage
#####

abstract type AbstractStage end

struct PreExperimentStage <: AbstractStage end
const PRE_EXPERIMENT_STAGE = PreExperimentStage()

struct PostExperimentStage <: AbstractStage end
const POST_EXPERIMENT_STAGE = PostExperimentStage()

struct PreEpisodeStage <: AbstractStage end
const PRE_EPISODE_STAGE = PreEpisodeStage()

struct PostEpisodeStage <: AbstractStage end
const POST_EPISODE_STAGE = PostEpisodeStage()

struct PreActStage <: AbstractStage end
const PRE_ACT_STAGE = PreActStage()

struct PostActStage <: AbstractStage end
const POST_ACT_STAGE = PostActStage()

(p::AbstractPolicy)(::AbstractStage, ::AbstractEnv) = nothing
(p::AbstractPolicy)(::AbstractStage, ::AbstractEnv, action) = nothing


text:
reload("OBC")

data = rand(0:1000,48,18000)

normfac = vec(mapslices(x->quantile(x,0.75), data, 2))

macro test()
    quote
        for i=1:1000
        #cls = MPM.mpm_classifier(rand(0:100,4,3), rand(0:100,4,3); burn=50, thin=50, d1=rand(4), d2=rand(4), usepriors=true)
        #MPM.sample(cls,1)
    end
end

@test

###############################
##### Copy of the above #######
###############################

reload("OBC")

normfac = vec(mapslices(x->quantile(x,0.75), data, 2))

@test

###############################
# Another copy. Probably unnecessary,
# but I've sometimes had one or two reps work
# before before segfaulting
###############################

reload("OBC")

normfac = vec(mapslices(x->quantile(x,0.75), data, 2))

@test


text:
using PerfChecker
using BenchmarkTools

using PatternFolds

target = PatternFolds

function bench() # 0.2.x
    # Intervals
    itv = Interval{Open,Closed}(0.0, 1.0)
    i = IntervalsFold(itv, 2.0, 1000)

    unfold(i)
    collect(i)
    reverse(collect(i))

    # rand(i, 1000)

    # Vectors
    vf = make_vector_fold([0, 1], 2, 1000)
    # @info "Checking VectorFold" vf pattern(vf) gap(vf) folds(vf) length(vf)

    unfold(vf)
    collect(vf)
    reverse(collect(vf))

    rand(vf, 1000)

    return nothing
end

# function bench() # 0.1.1-0.1.5
#     # Intervals
#     i = IntervalsFold(Interval((0.0, true), (1.0, false)), 2.0, 1000)

#     unfold(i)
#     collect(i)
#     reverse(collect(i))

#     # rand(i, 1000)

#     # Vectors
#     vf = VectorFold([0, 1], 2, 1000)
#     # @info "Checking VectorFold" vf pattern(vf) gap(vf) folds(vf) length(vf)

#     unfold(vf)
#     collect(vf)
#     reverse(collect(vf))

#     map(_ -> rand(vf), 1:1000)

#     return nothing
# end

t = @benchmark bench() evals = 1 samples = 1000 seconds = 3600

# Actual call to PerfChecker
store_benchmark(t, target; path=@__DIR__)


text:
module PrivateDynamicData

export gaussianMechConstant, gaussianMechConstant2,
       mean_dp, laplaceMech, gaussianMech, truncatedLaplaceMech,
       staticInputBlock_DPKF_ss, dfactor, evaluateKFperf,
       staticInputBlock_DPLQG_ss, evaluateLQGperf

include("utils.jl")
include("dpkf.jl")  # Differentially private Kalman filtering

#= Not exported
dpkf_ok: boolean, true if Mosek has been installed (for Kalman filtering)
=#

end # module


text:
function JSServe.jsrender(session::Session, scene::Scene)
    three, canvas = WGLMakie.three_display(session, scene)
    Makie.push_screen!(scene, three)
    return canvas
end

function JSServe.jsrender(session::Session, scene::Makie.FigureLike)
    return JSServe.jsrender(session, Makie.get_scene(scene))
end

const WEB_MIMES = (MIME"text/html", MIME"application/vnd.webio.application+html",
                   MIME"application/prs.juno.plotpane+html", MIME"juliavscode/html")

for M in WEB_MIMES
    @eval begin
        function Makie.backend_show(::WGLBackend, io::IO, m::$M, scene::Scene)
            three = nothing
            inline_display = App() do session::Session
                three, canvas = three_display(session, scene)
                Makie.push_screen!(scene, three)
                return canvas
            end
            Base.show(io, m, inline_display)
            return three
        end
    end
end

function scene2image(scene::Scene)
    if !JSServe.has_html_display()
        error("""
        There is no Display that can show HTML.
        If in the REPL and you have a browser,
        you can always run `JSServe.browser_display()` to show plots in the browser
        """)
    end
    three = nothing
    session = nothing
    app = App() do s::Session
        session = s
        three, canvas = three_display(s, scene)
        return canvas
    end
    # display in current
    display(app)
    done = Base.timedwait(()-> isready(session.js_fully_loaded), 30.0)
    if done == :timed_out
        error("JS Session not ready after 30s waiting, possibly errored while displaying")
    end
    return Makie.colorbuffer(three)
end

function Makie.backend_show(::WGLBackend, io::IO, m::MIME"image/png",
                                       scene::Scene)
    img = scene2image(scene)
    return FileIO.save(FileIO.Stream(FileIO.format"PNG", io), img)
end

function Makie.backend_show(::WGLBackend, io::IO, m::MIME"image/jpeg",
                                       scene::Scene)
    img = scene2image(scene)
    return FileIO.save(FileIO.Stream(FileIO.format"JPEG", io), img)
end

function Makie.backend_showable(::WGLBackend, ::T, scene::Scene) where {T<:MIME}
    return T in WEB_MIMES
end

struct WebDisplay <: Makie.AbstractScreen
    three::Base.RefValue{Any}
    display::Any
end

function Makie.backend_display(::WGLBackend, scene::Scene)
    # Reference to three object which gets set once we serve this to a browser
    three_ref = Base.RefValue{Any}(nothing)
    app = App() do s, request
        three, canvas = three_display(s, scene)
        three_ref[] = three
        return canvas
    end
    actual_display = display(app)
    return WebDisplay(three_ref, actual_display)
end

function Base.delete!(td::WebDisplay, scene::Scene, plot::AbstractPlot)
    delete!(get_three(td), scene, plot)
end

function session2image(sessionlike)
    s = JSServe.session(sessionlike)
    to_data = js"document.querySelector('canvas').toDataURL()"
    picture_base64 = JSServe.evaljs_value(s, to_data; time_out=100)
    picture_base64 = replace(picture_base64, "data:image/png;base64," => "")
    bytes = JSServe.Base64.base64decode(picture_base64)
    return ImageMagick.load_(bytes)
end

function Makie.colorbuffer(screen::ThreeDisplay)
    return session2image(screen)
end

function get_three(screen::WebDisplay; timeout = 30)
    # WebDisplay is not guaranteed to get displayed in the browser, so we wait a while
    # to see if anything gets displayed!
    tstart = time()
    while time() - tstart < timeout
        if screen.three[] !== nothing
            three = screen.three[]
            session = JSServe.session(three)
            if isready(session.js_fully_loaded)
                # Error on js during init! We can't continue like this :'(
                if session.init_error[] !== nothing
                    throw(session.init_error[])
                end
                return three
            end
        end
        yield()
    end
    return nothing
end

function Makie.colorbuffer(screen::WebDisplay)
    return session2image(get_three(screen))
end

function Base.insert!(td::WebDisplay, scene::Scene, plot::AbstractPlot)
    disp = get_three(td)
    disp === nothing && error("Plot needs to be displayed to insert additional plots")
    insert!(disp, scene, plot)
end


text:
export dubinsCC, dubinsCC_length, dubinsCC_waypoints

dubinsCC_length(q0::StaticVector, qf::StaticVector; κ_max=1, σ_max=1) = dubinsCC(q0, qf, κ_max=κ_max, σ_max=σ_max).cost
function dubinsCC_waypoints(q0::StaticVector, qf::StaticVector, dt_or_N; v=1, κ_max=1, σ_max=1)
    waypoints(SimpleCarDynamics{0,1}(), SE2κState(q0), dubinsCC(q0, qf, v=v, κ_max=κ_max, σ_max=σ_max).controls, dt_or_N)
end

function dubinsCC(q0::StaticVector{4,T0}, qf::StaticVector{4,Tf}; v::V=1, κ_max::K=1, σ_max::S=1) where {T0,Tf,V,K,S}
    @assert (abs(q0.κ) <= κ_max && abs(qf.κ) <= κ_max) "$q0, $qf endpoint curvatures exceed bound of $κ_max"

    D = SimpleCarDynamics{0,1}()
    T = promote_type(T0, Tf, V, K, S)
    cmin, ctrl = T(Inf), zeros(SVector{11,VelocityCurvRateStep{T}})

    t0 = abs(q0.κ)/σ_max
    tf = abs(qf.κ)/σ_max
    c0_towards = StepControl{T}( t0, VelocityCurvRateControl{T}(v, -flipsign(σ_max, q0.κ)))
    c0_away    = StepControl{T}(-t0, VelocityCurvRateControl{T}(v,  flipsign(σ_max, q0.κ)))
    cf_towards = StepControl{T}( tf, VelocityCurvRateControl{T}(v, -flipsign(σ_max, qf.κ)))
    cf_away    = StepControl{T}(-tf, VelocityCurvRateControl{T}(v,  flipsign(σ_max, qf.κ)))

    # towards 0, towards 0
    q0_0 = SE2State(propagate(D, q0, c0_towards))
    qf_0 = SE2State(propagate(D, qf, cf_towards))
    cnew, ctrl_0 = dubinsCC(q0_0, qf_0, true, false, true, true, qf.κ >= 0, qf.κ <= 0, v=v, κ_max=κ_max, σ_max=σ_max)
    cnew = cnew + (t0 - tf)
    if cnew < cmin
        cmin = cnew
        ctrl_0 = setindex(ctrl_0, StepControl(ctrl_0[9].t - tf, ctrl_0[9].u), 9)
        ctrl = [SVector(c0_towards); ctrl_0; SVector(zero(VelocityCurvRateStep{T}))]
    end

    # towards 0, away from 0
    q0_0 = SE2State(propagate(D, q0, c0_towards))
    qf_0 = SE2State(propagate(D, qf, cf_away))
    cnew, ctrl_0 = dubinsCC(q0_0, qf_0, true, true, true, true, true, true, v=v, κ_max=κ_max, σ_max=σ_max)
    cnew = cnew + (t0 + tf)
    if cnew < cmin
        cmin = cnew
        ctrl = [SVector(c0_towards); ctrl_0; SVector(StepControl(tf, cf_away.u))]
    end

    # away from 0, towards 0
    q0_0 = SE2State(propagate(D, q0, c0_away))
    qf_0 = SE2State(propagate(D, qf, cf_towards))
    cnew, ctrl_0 = dubinsCC(q0_0, qf_0, false, false, q0.κ >= 0, q0.κ <= 0, qf.κ >= 0, qf.κ <= 0, v=v, κ_max=κ_max, σ_max=σ_max)
    cnew = cnew + (-t0 - tf)
    if cnew < cmin
        cmin = cnew
        ctrl_0 = setindex(ctrl_0, StepControl(ctrl_0[1].t - t0, ctrl_0[1].u), 1)
        ctrl_0 = setindex(ctrl_0, StepControl(ctrl_0[9].t - tf, ctrl_0[9].u), 9)
        ctrl = [SVector(zero(VelocityCurvRateStep{T})); ctrl_0; SVector(zero(VelocityCurvRateStep{T}))]
    end

    # away from 0, away from 0
    q0_0 = SE2State(propagate(D, q0, c0_away))
    qf_0 = SE2State(propagate(D, qf, cf_away))
    cnew, ctrl_0 = dubinsCC(q0_0, qf_0, false, true, q0.κ >= 0, q0.κ <= 0, true, true, v=v, κ_max=κ_max, σ_max=σ_max)
    cnew = cnew + (-t0 + tf)
    if cnew < cmin
        cmin = cnew
        ctrl_0 = setindex(ctrl_0, StepControl(ctrl_0[1].t - t0, ctrl_0[1].u), 1)
        ctrl = [SVector(zero(VelocityCurvRateStep{T})); ctrl_0; SVector(StepControl(tf, cf_away.u))]
    end

    (cost=cmin, controls=ctrl)
end

function dubinsCC((x0, y0, θ0)::StaticVector{3,T0},
                  (xf, yf, θf)::StaticVector{3,Tf},
                  allow_short_turn_1=true,
                  allow_short_turn_3=true,
                  LXX=true,
                  RXX=true,
                  XXL=true,
                  XXR=true;
                  v::V=1,
                  κ_max::K=1,
                  σ_max::S=1) where {T0,Tf,V,K,S}
    T = promote_type(T0, Tf, V, K, S)
    θ_lim, r, γ = CC_steering_constants(T(κ_max), T(σ_max/v))
    @inline turn_length(β) = CC_turn_length(β, T(κ_max), T(σ_max/v), θ_lim, r, γ)
    @inline turn_control(β) = CC_turn_control(β, T(κ_max), T(σ_max/v), θ_lim, r, γ)

    dx = (xf - x0)/r
    dy = (yf - y0)/r
    d = hypot(dx, dy)
    θ = atan(dy, dx)
    a = θ0 - θ
    b = θf - θ

    sap, cap = sincos(a + γ)
    sbp, cbp = sincos(b + γ)
    sam, cam = sincos(a - γ)
    sbm, cbm = sincos(b - γ)
    sγ, cγ = sincos(γ)
    cmin = T(Inf)
    ctrl = zeros(SVector{9,VelocityCurvRateStep{T}})

    ### 1. LSL: a-γ, b+γ
    if LXX && XXL
        ca, sa, cb, sb = cam, sam, cbp, sbp
        tmp = 2 + d*d - 2*(ca*cb + sa*sb - d*(sa - sb))
        if tmp >= 4*sγ*sγ    # TODO: M-style connections
            θ = atan(cb - ca, d + sa - sb)
            t = mod2piF(-a + θ)
            p = sqrt(tmp) - 2*sγ
            q = mod2piF(b - θ)
            p = r*p
            if (allow_short_turn_1 || t >= θ_lim) && (allow_short_turn_3 || q >= θ_lim)
                c = turn_length(t) + p + turn_length(q)
                if c < cmin
                    cmin = c
                    ctrl = [
                        turn_control(t);
                        SVector(zero(VelocityCurvRateStep{T}),
                                StepControl(p, VelocityCurvRateControl(T(1), T(0))),
                                zero(VelocityCurvRateStep{T}));
                        turn_control(q)
                    ]
                end
            end
        end
    end

    ### 2. RSR: a+γ, b-γ
    if RXX && XXR
        ca, sa, cb, sb = cap, sap, cbm, sbm
        tmp = 2 + d*d - 2*(ca*cb + sa*sb - d*(sb - sa))
        if tmp >= 4*sγ*sγ    # TODO: M-style connections
            θ = atan(ca - cb, d - sa + sb)
            t = mod2piF(a - θ)
            p = sqrt(tmp) - 2*sγ
            q = mod2piF(-b + θ)
            p = r*p
            if (allow_short_turn_1 || t >= θ_lim) && (allow_short_turn_3 || q >= θ_lim)
                c = turn_length(t) + p + turn_length(q)
                if c < cmin
                    cmin = c
                    ctrl = [
                        turn_control(-t);
                        SVector(zero(VelocityCurvRateStep{T}),
                                StepControl(p, VelocityCurvRateControl(T(1), T(0))),
                                zero(VelocityCurvRateStep{T}));
                        turn_control(-q)
                    ]
                end
            end
        end
    end

    ### 3. RSL: a+γ, b+γ
    if RXX && XXL
        ca, sa, cb, sb = cap, sap, cbp, sbp
        tmp = d*d - 2 + 2*(ca*cb + sa*sb - d*(sa + sb))
        if tmp >= 0
            p = sqrt(tmp + 4*sγ*sγ) - 2*sγ
            x = sqrt(tmp)
            θ = atan(ca + cb, d - sa - sb) - atan(T(2), x)
            ε = acot(x/2) - acot((p/2 + sγ)/cγ)
            t = mod2piF(a - θ - ε)
            q = mod2piF(b - θ - ε)
            p = r*p
            if (allow_short_turn_1 || t >= θ_lim) && (allow_short_turn_3 || q >= θ_lim)
                c = turn_length(t) + p + turn_length(q)
                if c < cmin
                    cmin = c
                    ctrl = [
                        turn_control(-t);
                        SVector(zero(VelocityCurvRateStep{T}),
                                StepControl(p, VelocityCurvRateControl(T(1), T(0))),
                                zero(VelocityCurvRateStep{T}));
                        turn_control(q)
                    ]
                end
            end
        end
    end

    ### 4. LSR: a-γ, b-γ
    if LXX && XXR
        ca, sa, cb, sb = cam, sam, cbm, sbm
        tmp = -2 + d*d + 2*(ca*cb + sa*sb + d*(sa + sb))
        if tmp >= 0
            p = sqrt(tmp + 4*sγ*sγ) - 2*sγ
            x = sqrt(tmp)
            θ = atan(-ca - cb, d + sa + sb) - atan(-T(2), x)
            ε = acot(x/2) - acot((p/2 + sγ)/cγ)
            t = mod2piF(-a + θ - ε)
            q = mod2piF(-b + θ - ε)
            p = r*p
            if (allow_short_turn_1 || t >= θ_lim) && (allow_short_turn_3 || q >= θ_lim)
                c = turn_length(t) + p + turn_length(q)
                if c < cmin
                    cmin = c
                    ctrl = [
                        turn_control(t);
                        SVector(zero(VelocityCurvRateStep{T}),
                                StepControl(p, VelocityCurvRateControl(T(1), T(0))),
                                zero(VelocityCurvRateStep{T}));
                        turn_control(-q)
                    ]
                end
            end
        end
    end

    ### 5. RLR: a+γ, b-γ
    if RXX && XXR
        ca, sa, cb, sb = cap, sap, cbm, sbm
        tmp = (6 - d*d + 2*(ca*cb + sa*sb + d*(sa - sb)))/8
        if abs(tmp) < 1
            p_1 = mod2piF(acos(tmp) - 2*γ)
            p_2 = mod2piF(-p_1 - 4*γ)
            t_2 = mod2piF(a - atan(ca - cb, d - sa + sb) + p_2/2)
            q_2 = mod2piF(a - b - t_2 + p_2)
            t_1 = mod2piF(t_2 + T(pi) - p_2 - 2*γ)
            q_1 = mod2piF(q_2 + T(pi) - p_2 - 2*γ)
            if (allow_short_turn_1 || t_1 >= θ_lim) && (allow_short_turn_3 || q_1 >= θ_lim)
                c_1 = turn_length(t_1) + turn_length(p_1) + turn_length(q_1)
                if c_1 < cmin
                    cmin = c_1
                    ctrl = [
                        turn_control(-t_1);
                        turn_control(p_1);
                        turn_control(-q_1)
                    ]
                end
            end
            if (allow_short_turn_1 || t_2 >= θ_lim) && (allow_short_turn_3 || q_2 >= θ_lim)
                c_2 = turn_length(t_2) + turn_length(p_2) + turn_length(q_2)
                if c_2 < cmin
                    cmin = c_2
                    ctrl = [
                        turn_control(-t_2);
                        turn_control(p_2);
                        turn_control(-q_2)
                    ]
                end
            end
        end
    end

    ### 6. LRL: a-γ, b+γ
    if LXX && XXL
        ca, sa, cb, sb = cam, sam, cbp, sbp
        tmp = (6 - d*d + 2*(ca*cb + sa*sb - d*(sa - sb)))/8
        if abs(tmp) < 1
            p_1 = mod2piF(acos(tmp) - 2*γ)
            p_2 = mod2piF(-p_1 - 4*γ)
            t_2 = mod2piF(-a + atan(-ca + cb, d + sa - sb) + p_2/2)
            q_2 = mod2piF(b - a - t_2 + p_2)
            t_1 = mod2piF(t_2 + T(pi) - p_2 - 2*γ)
            q_1 = mod2piF(q_2 + T(pi) - p_2 - 2*γ)
            if (allow_short_turn_1 || t_1 >= θ_lim) && (allow_short_turn_3 || q_1 >= θ_lim)
                c_1 = turn_length(t_1) + turn_length(p_1) + turn_length(q_1)
                if c_1 < cmin
                    cmin = c_1
                    ctrl = [
                        turn_control(t_1);
                        turn_control(-p_1);
                        turn_control(q_1)
                    ]
                end
            end
            if (allow_short_turn_1 || t_2 >= θ_lim) && (allow_short_turn_3 || q_2 >= θ_lim)
                c_2 = turn_length(t_2) + turn_length(p_2) + turn_length(q_2)
                if c_2 < cmin
                    cmin = c_2
                    ctrl = [
                        turn_control(t_2);
                        turn_control(-p_2);
                        turn_control(q_2)
                    ]
                end
            end
        end
    end

    ctrl = scalespeed.(ctrl, v)
    (cost=cmin/v, controls=ctrl)
end

@inline function CC_steering_constants(κ_max::K, σ_max::S) where {K,S}
    T  = promote_type(K, S)
    xi = sqrt(T(pi)/σ_max)*fresnelC(κ_max/sqrt(T(pi)*σ_max))
    yi = sqrt(T(pi)/σ_max)*fresnelS(κ_max/sqrt(T(pi)*σ_max))
    θi = κ_max*κ_max/(2*σ_max)
    s, c = sincos(θi)
    xΩ = xi - s/κ_max
    yΩ = yi + c/κ_max
    (θ=2*θi, r=hypot(xΩ, yΩ), γ=atan(xΩ/yΩ))
end

@inline function CC_turn_length(β::T, κ_max, σ_max, θ_lim, r, γ) where {T}    # θ_lim = κ_max^2/σ_max
    b = abs(β)
    if b < sqrt(eps(T))
        2*r*sin(γ) + b*r*cos(γ)
    elseif b < θ_lim
        bp = sqrt(b/T(pi))
        σ = (cos(b/2)*fresnelC(bp) + sin(b/2)*fresnelS(bp))/r/sin(b/2 + γ)
        σ = T(pi)*σ*σ
        2*sqrt(b/σ)
    else
        (b + θ_lim)/κ_max
    end
end

@inline function CC_turn_control(β::T, κ_max, σ_max, θ_lim, r, γ) where {T}    # θ_lim = κ_max^2/σ_max
    b = abs(β)
    if b < sqrt(eps(T))
        sγ, cγ = sincos(γ)
        SVector(
            StepControl(r*sγ + b*r*cγ/2, VelocityCurvRateControl(T(1), b/(r*r*sγ*sγ))),
            zero(VelocityCurvRateStep{T}),
            StepControl(r*sγ + b*r*cγ/2, VelocityCurvRateControl(T(1), -b/(r*r*sγ*sγ)))
        )
    elseif b < θ_lim
        bp = sqrt(b/T(pi))
        σ = (cos(b/2)*fresnelC(bp) + sin(b/2)*fresnelS(bp))/r/sin(b/2 + γ)
        σ = T(pi)*σ*σ
        SVector(
            StepControl(sqrt(b/σ), VelocityCurvRateControl(T(1), flipsign(σ, β))),
            zero(VelocityCurvRateStep{T}),
            StepControl(sqrt(b/σ), VelocityCurvRateControl(T(1), -flipsign(σ, β)))
        )
    else
        SVector(
            StepControl(κ_max/σ_max, VelocityCurvRateControl(T(1), flipsign(σ_max, β))),
            StepControl((b - θ_lim)/κ_max, VelocityCurvRateControl(T(1), T(0))),
            StepControl(κ_max/σ_max, VelocityCurvRateControl(T(1), -flipsign(σ_max, β)))
        )
    end
end


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

name = "Tcl"
version = v"8.6.9"

# Collection of sources required to build Tcl
sources = [
    ArchiveSource("https://downloads.sourceforge.net/sourceforge/tcl/tcl$(version)-src.tar.gz",
                  "ad0cd2de2c87b9ba8086b43957a0de3eb2eb565c7159d5f53ccbba3feb915f4e"),
]

# Bash recipe for building across all platforms
script = raw"""
if [[ "${target}" == *-mingw* ]]; then
    cd $WORKSPACE/srcdir/tcl*/win/
    # `make install` calls `tclsh` on Windows
    apk add tcl
else
    cd $WORKSPACE/srcdir/tcl*/unix/
fi

FLAGS=()
if [[ "${target}" == x86_64-* ]]; then
    FLAGS+=(--enable-64bit)
fi
./configure --prefix=${prefix} --host=${target} "${FLAGS[@]}"
make -j${nproc}
make install
# Tk needs private headers
make install-private-headers

# Install license file
install_license $WORKSPACE/srcdir/tcl*/license.terms
"""

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

# The products that we will ensure are always built
products = [
    LibraryProduct(["libtcl8.6", "libtcl8", "tcl86"], :libtcl),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "Zlib_jll",
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: mortenpi/Yggdrasil

max_issues_repo_name: mortenpi/Yggdrasil

max_forks_repo_name: mortenpi/Yggdrasil


text:
struct_nodes = findall("//type[@category='struct']", xroot)
ignored_types = [
    "VkBaseInStructure",
    "VkBaseOutStructure",
]

vk_structs_graph_base = vulkan_types |> @filter(_.category ∈ ["struct", "handle"] && _.name ∉ ignored_types) |> DataFrame
vk_structs_graph = vk_structs_graph_base.name

function remove_pointer(type)
    ptr_match = match(r"Ptr{(.*)}", type)
    isnothing(ptr_match) && return type
    first(ptr_match.captures)
end

function adj_list()
    adjl = Edge{Int}[]
    for (i, row) ∈ enumerate(eachrow(vk_structs_graph_base))
        name = row.name
        name ∈ ignored_types && continue
        if !isnothing(row.alias)
            push!(adjl, Edge(i, findfirst(x -> x == row.alias, vk_structs_graph)))
        elseif row.category == "handle"
            handle_info = vulkan_handles[findfirst(x -> x == name, vulkan_handles.name), :]
            !isnothing(handle_info.parent) && push!(adjl, Edge(i, findfirst(x -> x == handle_info.parent, vk_structs_graph)))
        else
            types = remove_pointer.(vulkan_fields.type[findall(x -> x == name, vulkan_fields.parent)])
            append!(adjl, Edge.(filter(!isnothing, indexin(unique(types), vk_structs_graph)), i))
        end
    end

    adjl
end

check_acyclic(g) = !is_cyclic(g) || error("Cycles: $(getindex.(Ref(vk_structs_graph), simplecycles_hawick_james(g)))")

adj = adj_list()

deps_g = SimpleDiGraph(adj)

# add verts to fill the graph with isolated types (avoid bound errors for isolated types)
nv(deps_g) < length(vk_structs_graph) && add_vertices!(deps_g, length(vk_structs_graph) - nv(deps_g))
check_acyclic(deps_g)

macro check_in_graph(type) :($(esc(type)) ∉ vk_structs_graph && error(string($(esc(type))) * " not in graph")) end

macro graph_index(type)
    quote
        @check_in_graph $(esc(type))
        $(esc(type)) = findfirst(vk_structs_graph .== $(esc(type)))
    end
end

macro graph_index(id, type)
    quote
        @check_in_graph $(esc(type))
        $(esc(id)) = findfirst(vk_structs_graph .== $(esc(type)))
    end
end

"""
    depends_on(a, b)

Returns whether b depends on a

# Examples
```
julia> depends_on("VkPhysicalDeviceProperties", "VkPhysicalDeviceLimits")
false
julia> depends_on("VkPhysicalDeviceLimits", "VkPhysicalDeviceProperties")
true
```
"""
function depends_on(a, b)
    @graph_index a
    @graph_index b
    a ∈ inneighbors(deps_g, b)
end

"""
    vk_dependencies(a)

# Examples
```
julia> vk_dependencies("VkPhysicalDeviceProperties")
2-element Array{String,1}:
 "VkPhysicalDeviceSparseProperties"
 "VkPhysicalDeviceLimits"
```
"""
function vk_dependencies(a)
    @graph_index a
    getindex.(Ref(vk_structs_graph), inneighbors(deps_g, a))
end

@assert vk_dependencies("VkPhysicalDeviceProperties") == ["VkPhysicalDeviceSparseProperties", "VkPhysicalDeviceLimits"]


text:
@system ScapeGrowth(Stage, LeafAppearance, FloralInitiation) begin
    #HACK: can't mixin ScapeRemoval/FlowerAppearance due to cyclic dependency
    scape_removed ~ hold
    flower_appeared ~ hold

    #HACK: use phyllochron
    SR_max(LTAR_max): maximum_scaping_rate ~ track(u"d^-1")

    scape(r=SR_max, β=BF.ΔT) => r*β ~ accumulate(when=scaping)

    scapeable(leaf_appeared & floral_initiated) ~ flag
    scaped(scape_removed | flower_appeared) ~ flag
    scaping(scapeable & !scaped) ~ flag
end

@system ScapeAppearance(Stage, ScapeGrowth) begin
    scape_appearable(scapeable) ~ flag
    scape_appearance_threshold => 3.0 ~ preserve(parameter)
    scape_appeared(scape, t=scape_appearance_threshold) => (scape >= t) ~ flag
    scape_appearing(scape_appearable & !scape_appeared) ~ flag

    # def finish(self):
    #     print(f"* Scape Tip Visible: time = {self.time}, leaves = {self.pheno.leaves_appeared} / {self.pheno.leaves_initiated}")
end

@system ScapeRemoval(Stage, ScapeGrowth, ScapeAppearance) begin
    #FIXME handling default (non-removal) value?
    scape_removal_date => nothing ~ preserve::datetime(optional, parameter)

    scape_removeable(scape_appeared) ~ flag
    scape_removed(scape_removal_date, scape_removeable, t=calendar.time) => begin
        isnothing(scape_removal_date) ? false : scape_removeable && (t >= scape_removal_date)
    end ~ flag
    scape_removing(scape_appeared & !scape_removed) ~ flag

    # def finish(self):
    #     print(f"* Scape Removed and Bulb Maturing: time = {self.time}")
end

@system FlowerAppearance(Stage, ScapeGrowth) begin
    flower_appearable(scapeable) ~ flag
    flower_appearance_threshold => 5.0 ~ preserve(parameter)
    flower_appeared(scape, t=flower_appearance_threshold, scape_removed) => (scape >= t && !scape_removed) ~ flag
    flower_appearing(flower_appearable & !flower_appeared) ~ flag

    # def finish(self):
    #     print(f"* Inflorescence Visible and Flowering: time = {self.time}")
end

@system BulbilAppearance(Stage, ScapeGrowth) begin
    bulbil_appearable(scapeable) ~ flag
    bulbil_appearance_threshold => 5.5 ~ preserve(parameter)
    bulbil_appeared(scape, t=bulbil_appearance_threshold, scape_removed) => (scape >= t && !scape_removed) ~ flag
    bulbil_appearing(bulbil_appearable & !bulbil_appeared) ~ flag

    # def finish(self):
    #     print(f"* Bulbil and Bulb Maturing: time = {self.time}")
end


text:
function ec_radar_init()
    opt = ec_init()
    opt["tooltip"]["trigger"] = "axis"
    opt["calculable"] = true
    opt["polar"] = Array(Any,0)

    return opt
end

function ec_radar{S<:String, T<:String, N<:Real}(legendlist::Vector{S}, data::Array{N}, labels::Vector{T}, bounds::Vector{N})
    opt = ec_radar_init()
    opt["legend"]["data"] = legendlist

    m, n = size(data)
    circle_label = Dict("indicator" => Array(Any, m))

    for i = 1:m
        circle_label["indicator"][i] = Dict("text" => labels[i], "max" => bounds[i])
    end

    push!(opt["polar"], circle_label)

    data_dict = Dict(
        "name" => "",
        "type" => "radar")

    data_dict["data"] = Array(Any, n)

    for i = 1:n
        data_dict["data"][i] = Dict(
                                "name" => legendlist[i],
                                "value" => data[:,i]
                                )
    end

    push!(opt["series"], data_dict)

    return opt
end


text:
import Images

save_clamped(filename::AbstractString, img::AbstractMatrix{<:RGB}) =
    Images.save(filename, map(clamp01, img))

export save_clamped


text:
distance(p1::Complex, p2::Complex) = hypot(p2 - p1)
distance((p1, p2)::Tuple{Complex, Complex}) = distance(p1, p2)
function distance(p1::AbstractPoint2D, p2::AbstractPoint2D)
    z1 = complex(p1._x, p1._y)
    z2 = complex(p2._x, p2._y)
    return distance(z1, z2)
end

"""
    longedge(edge, tolerance, geom2fcn)

Return true if `edge` has length greater than `tolerance`.
"""
function longedge(edge::DelaunayEdge, tolerance, g2f::Geometry2Function)
    return distance(g2f(edge)) > tolerance
end

"""
    fcn2geom(x, xmin, xmax)

Linearly scale `x`, which exists ∈ (`xmin`, `xmax`) into `RootsAndPoles` domain.
"""
function fcn2geom(x, xmin, xmax)
    return (MAXCOORD - MINCOORD)*(x - xmin)/(xmax - xmin) + MINCOORD
end

# TODO: See pull request #50 on VoronoiDelaunay.jl which handles the space
# mapping automatically.

"""
    fcn2geom(z, rmin, rmax, imin, imax)

Linearly scale complex `z` within real range (`rmin`, `rmax`) and imaginary range (`imin`,
`imax`) into `RootsAndPoles` domain.
"""
function fcn2geom(z, rmin, rmax, imin, imax)
    zr = fcn2geom(real(z), rmin, rmax)
    zi = fcn2geom(imag(z), imin, imax)
    return complex(zr, zi)
end

"""
    geom2fcn(g, xmin, xmax)

Linearly scale `g`, which exists in the `RootsAndPoles` domain into the range of the
original function, bounded by (`xmin`, `xmax`).
"""
function geom2fcn(g, xmin, xmax)
    return (xmax - xmin)*(g - MINCOORD)/(MAXCOORD - MINCOORD) + xmin
end

"""
    geom2fcn(z, rmin, rmax, imin, imax)

Scale complex value `z` from the `RootsAndPoles` geometry domain to the comlex function""
domain with real range (`rmin`, `rmax`) and imaginary range (`imin`, `imax`).
"""
geom2fcn(z, rmin, rmax, imin, imax) = geom2fcn(real(z), imag(z), rmin, rmax, imin, imax)

"""
    geom2fcn(zr, zi, rmin, rmax, imin, imax)

Scale real and complex components `zr`, `zi` from the `RootsAndPoles` geometry domain to the
function domain.
"""
geom2fcn(zr, zi, rmin, rmax, imin, imax) = complex(geom2fcn(zr, rmin, rmax),
                                                   geom2fcn(zi, imin, imax))

"""
    geom2fcn(pt::AbstractPoint2D, rmin, rmax, imin, imax)

Scale `pt` from the geometry domain to the function domain.
"""
geom2fcn(pt::AbstractPoint2D, rmin, rmax, imin, imax) = geom2fcn(getx(pt), gety(pt), rmin,
                                                                 rmax, imin, imax)

"""
    geom2fcn(edge::DelaunayEdge, rmin, rmax, imin, imax)

Scale both points of `edge` from the geometry domain to the function domain.
"""
function geom2fcn(edge::DelaunayEdge{P}, rmin, rmax, imin, imax) where P <: AbstractPoint2D
    return (geom2fcn(geta(edge), rmin, rmax, imin, imax),
        geom2fcn(getb(edge), rmin, rmax, imin, imax))
end

"""
    getplotdata(tess, quadrants, phasediffs, g2f)

Return complex coordinates of every tesselation edge and a vector with the corresponding
index specifiying its phase quadrant.

This function is useful for plotting the final tesselation using results from calling
`grpf` with the `PlotData()` argument.

The edge colors are (where Q means "quadrant"):

| Color index |         Meaning         |  Quadrant phase   |
|:-----------:|:-----------------------:|:-----------------:|
|      1      |           Q 1           | 0 ≤ arg f < π/2   |
|      2      |           Q 2           | π/2 ≤ arg f < π   |
|      3      |           Q 3           | π ≤ arg f < 3π/2  |
|      4      |           Q 4           | 3π/2 ≤ arg f < 2π |
|      5      | phase change (boundary) |         -         |
|      6      | candidate edge          |         -         |

# Example

Although not self contained, this example shows how `getplotdata` can be used:

```julia
roots, poles, quads, phasediffs, tess, g2f = grpf(f, origcoords, PlotData());
z, edgecolors = getplotdata(tess, quadrants, phasediffs, g2f);

using Plots
plot(real(z), imag(z), group=edgecolors)
```
"""
function getplotdata(tess, quadrants, phasediffs, g2f)
    edges = [e for e in delaunayedges(tess)]

    edgecolors = fill(5, length(edges))
    for ii in eachindex(edges)
        if phasediffs[ii] == 2  # candidate edges
            edgecolors[ii] = 6
        elseif phasediffs[ii] == 0
            edgecolors[ii] = quadrants[getindex(geta(edges[ii]))]
        end
    end

    edgecolors = repeat(edgecolors, inner=3)  # x, y is xₐ, xᵦ, NaN, repeat

    x, y = getplotxy(edges)
    z = g2f.(x, y)

    I = sortperm(edgecolors)

    return z[I], edgecolors[I]
end


text:
# This is voxel_gravity.jl
#
# This function allows the computation of bouguer anomalies over a regular grid
#
# Author: Arne Spang 06/2021

using Printf            # do we need this?
using Statistics        # do we need that?

export voxGrav

"""
    voxGrav(X::Array{Float64, 3}, Y::Array{Float64, 3}, Z::Array{Float64, 3}, RHO::Array{Float64, 3};
    refMod="AVG", lengthUnit="m", rhoTol=1e-9, Topo=[], outName="Bouguer", printing=true)

    Computes Bouguer anomalies and gradients

    Required arguments:
    X,Y,Z:       3D matrices with the coordinates of the grid (X should vary in the first dimension, Y in the second, Z (vertical) in the thrid)
    RHO:         3D matrix with the densitiy at each grid point [kg/m^3]

    Optional arguments:
    refMod:      1D vector with the reference density for each depth
                 Alternatively, the strings "NE", "SE", "SW", "NW", "AVG" can be used.
                 In that case, one of the corners of `RHO` is used as reference model.
                 In case of "AVG" the reference model is the average of each depth slice.
    lengthUnit:  The unit of the coordinates and topography file. Either "m" or "km"
    rhoTol:      density differences smaller than rhoTol will be ignored [kg/m^3]
    Topo:        2D matrix with the topography of the surface (only relevant for the paraview output)
    outName:     name of the paraview output (do not include file type)
    printing:    activate printing of additional information [true or false]
"""
function voxGrav(X::Array{Float64, 3}, Y::Array{Float64, 3}, Z::Array{Float64, 3}, RHO::Array{Float64, 3};
                 refMod="AVG", lengthUnit="m", rhoTol=1e-9, Topo=[], outName="Bouguer", printing=true)

    ## check input
    X, Y, Z, RHO, RefMod, rhoTol, Topo, orient = checkInput(X, Y, Z, RHO, refMod, lengthUnit, rhoTol, Topo, outName, printing)

    ################ precompute things ################
    # define constants
    G      = 6.67408e-11

    # get coordinate vectors
    x_vec  = X[:,1,1]
    y_vec  = Y[1,:,1]
    z_vec  = Z[1,1,:]

    # cut everything above sea level
    ind    = findall(x->x<=0, z_vec)
    X      = X[:,:,ind]
    Y      = Y[:,:,ind]
    Z      = Z[:,:,ind]
    RHO    = RHO[:,:,ind]
    RefMod = RefMod[ind]

    # check dimensions
    nx     = size(X,1);
    ny     = size(X,2);
    nz     = size(X,3);

    # substract reference model
    for i = 1 : nz
        RHO[:,:,i] .= RHO[:,:,i] .- RefMod[i]
    end

    # interpolate density grid to cell centers
    DRHO = RHO[1:end-1,1:end-1,1:end-1] .+ RHO[2:end,1:end-1,1:end-1] .+ RHO[2:end,2:end,1:end-1] .+ RHO[1:end-1,2:end,1:end-1] +
            RHO[1:end-1,1:end-1,2:end]   .+ RHO[2:end,1:end-1,2:end]   .+ RHO[2:end,2:end,2:end]   .+ RHO[1:end-1,2:end,2:end]
    DRHO = DRHO ./ 8

    # voxel volume
    dx     = X[2,1,1] - X[1,1,1]
    dy     = Y[1,2,1] - Y[1,1,1]
    dz     = Z[1,1,2] - Z[1,1,1]
    dV     = abs(dx*dy*dz)

    # dV * G is a constants
    VG_vox = dV * G

    # coordinate vector of cell center grid
    xCells = x_vec[1:end-1] .+ dx/2
    yCells = y_vec[1:end-1] .+ dy/2
    zCells = z_vec[1:end-1] .+ dz/2

    # precompute distances
    if printing
        @printf "Precomputing Distances:"
        @time d_cube = precompDist(x_vec, y_vec, xCells, yCells, zCells, nx, ny, nz)
        @printf "\n"
    else
        d_cube = precompDist(x_vec, y_vec, xCells, yCells, zCells, nx, ny, nz)
    end
    ###################################################

    ############# compute bouguer anomaly #############
    if printing
        @printf "Computing Bouguer anomaly:"
        @time dg = computeBoug(nx,ny,nz,DRHO,d_cube,VG_vox,zCells,rhoTol)
        @printf "\n"
    else
        dg = computeBoug(nx,ny,nz,DRHO,d_cube,VG_vox,zCells,rhoTol)
    end
    ###################################################

    ############ compute bouguer gradients ############
    if printing
        @printf "Computing Bouguer gradients:"
        @time gradX, gradY = computeBougGrads(nx,ny,dg)
        @printf "\n"
    else
        gradX, gradY = computeBougGrads(nx,ny,dg)
    end
    ###################################################

    ############## write info and output ##############
    numRel = length(findall(x->abs(x)>rhoTol, DRHO))
    frac   = numRel/length(DRHO)
    coords = cat(X[:,:,1],Y[:,:,1],Topo,dims=4)
    coords = permutedims(coords,[4,1,2,3])
    vtkfile = vtk_grid(outName,coords)

    if printing
        @printf "%.3f %% of the domain contained anomalous densities. If this is more than expected, adjust rhoTol for faster computation.\n" 100*frac
        @printf "Writing output...\n"
        vtkfile["Bouguer Anomaly [mGal]"]   = dg
        vtkfile["Boug_Gradient_X [mGal/m]"] = gradX
        vtkfile["Boug_Gradient_Y [mGal/m]"] = gradY
        outfiles = vtk_save(vtkfile)
        @printf "Wrote output to: %s\n\n" outfiles
    else
        vtkfile["Bouguer Anomaly [mGal]"]   = dg
        vtkfile["Boug_Gradient_X [mGal/m]"] = gradX
        vtkfile["Boug_Gradient_Y [mGal/m]"] = gradY
        outfiles = vtk_save(vtkfile)
    end
    ###################################################

    if orient == 1
        return dg, gradX, gradY
    else
        return permutedims(dg, [2,1]), permutedims(gradX, [2,1]), permutedims(gradY, [2,1])
    end
end

function checkInput(X, Y, Z, RHO, refMod, lengthUnit, rhoTol, Topo, outName, printing)
    # orientation
    if X[1,1,1] ≠ X[2,1,1] && X[1,1,1] == X[1,2,1] && X[1,1,1] == X[1,1,2]
        orientation = 1
    elseif X[1,1,1] == X[2,1,1] && X[1,1,1] ≠ X[1,2,1] && X[1,1,1] == X[1,1,2]
        orientation = 2
        X   = permutedims(X,   [2,1,3])
        Y   = permutedims(Y,   [2,1,3])
        Z   = permutedims(Z,   [2,1,3])
        RHO = permutedims(RHO, [2,1,3])
    else
        error("Coordinate orientation looks wrong!")
    end

    # dimensions
    if !(size(X) == size(Y) && size(X) == size(Z) && size(X) == size(RHO))
        error("X, Y, Z, RHO must be 3D matrices of the same size.")
    end
    nz = size(X,3)

    # check if grid is regular
    dx = diff(X,dims=1)[:]
    dy = diff(Y,dims=2)[:]
    dz = diff(Z,dims=3)[:]
    tol = 1e-12
    if !(all(a->a<dx[1]+tol && a > dx[1]-tol,dx) && all(a->a<dy[1]+tol && a > dy[1]-tol,dy) && all(a->a<dz[1]+tol && a > dz[1]-tol,dz))
        error("Non-regular grids are not supported yet")
    end

    # Topo
    if !isempty(Topo)
        if !(size(Topo,1) == size(X,1) && size(Topo,2) == size(X,2))
            if printing
                @printf "Topo input dimensions do not fit X, Y, Z, RHO. Using a flat topography.\n"
            end
            Topo = zeros(size(X,1),size(X,2))
        end
    else
        Topo = zeros(size(X,1),size(X,2))
    end

    # tolerance
    if printing
        @printf "Using density tolerance: %.3e\n" rhoTol
    end

    # length unit
    if lengthUnit == "m" || lengthUnit == "M"
        if printing
            @printf "Assuming coordinates to be in meters.\n\n"
        end
    elseif lengthUnit == "km" || lengthUnit == "KM"
        if printing
            @printf "Converting coordinates to meters.\n\n"
        end
        X     = X    .* 1000;
        Y     = Y    .* 1000;
        Z     = Z    .* 1000;
        Topo  = Topo .* 1000;
    else
        error("lengthUnit should be \"m\" or \"km\".")
    end

    # reference model
    if !(typeof(refMod) == String)
        if length(refMod) == nz
            RefMod = refMod
        else
            error("RefMod must have the same length as the third dimension of X, Y, Z, RHO.")
        end
    else
        if refMod == "NE"
            RefMod = RHO[end,end,:]
        elseif refMod == "SE"
            RefMod = RHO[end,1,:]
        elseif refMod == "SW"
            RefMod = RHO[1,1,:]
        elseif refMod == "NW"
            RefMod = RHO[1,end,:]
        elseif refMode == "AVG"
            RefMod = !mean([1.,1.,1.],RHO)
        else
            error("RefMod should be NE, SE, SW, NW, AVG or a vector with one value for each depth.")
        end
    end

    # outName
    if !(typeof(outName) == String)
        error("outName must be a string.")
    end

    return X, Y, Z, RHO, RefMod, rhoTol, Topo, orientation

end

function precompDist(x_vec, y_vec, xCells, yCells, zCells, nx, ny, nz)
    d_square = zeros(nx-1,ny-1,nz-1)

    for iX = 1 : nx - 1
        for iY = 1 : ny - 1
            for iZ = 1 : nz - 1
                d_square[iX,iY,iZ] = (xCells[iX]-x_vec[1])^2 + (yCells[iY]-y_vec[1])^2 + zCells[iZ]^2
            end
        end
    end

    d_cube = d_square .^1.5

    return d_cube
end

function computeBoug(nx,ny,nz,DRHO,d_cube,VG_vox,zCells,rhoTol)
    mGal   = 1e5

    dg     = zeros(nx,ny)

    for jX = 1 : nx-1
        for jY = 1 : ny-1
            for jZ = 1 : nz-1
                if (DRHO[jX,jY,jZ] > rhoTol || DRHO[jX,jY,jZ] < -rhoTol)
                    for iX = 1 : nx
                        for iY = 1 : ny
                            d_indX = jX-iX
                            if d_indX < 0; d_indX = abs(d_indX); else d_indX = d_indX + 1; end
                            d_indY = jY-iY;
                            if d_indY < 0; d_indY = abs(d_indY); else d_indY = d_indY + 1; end
                            d_c       = d_cube[d_indX,d_indY,jZ];
                            dg[iX,iY] = dg[iX,iY] + VG_vox * DRHO[jX,jY,jZ] * -zCells[jZ]/d_c;
                        end
                    end
                end
            end
        end
    end

    dg = dg * mGal

    return dg
end

function computeBougGrads(nx,ny,dg)
    gradX  = zeros(nx,ny)
    gradY  = zeros(nx,ny)

    for iY = 1 : ny
        itp  = interpolate(dg[:,iY], BSpline(Quadratic(Reflect(OnCell()))))
        for iX = 1 : nx
            grad         = Interpolations.gradient(itp,iX)
            gradX[iX,iY] = grad[1]
        end
    end

    for iX = 1 : nx
        itp  = interpolate(dg[iX,:], BSpline(Quadratic(Reflect(OnCell()))))
        for iY = 1 : ny
            grad         = Interpolations.gradient(itp,iY)
            gradY[iX,iY] = grad[1]
        end
    end

    return gradX, gradY
end