text:
using ACE1docs
using Documenter

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

makedocs(;
    modules=[ACE1docs],
    authors="Christoph Ortner <christophortner0@gmail.com> and contributors",
    repo="https://github.com/ACEsuit/ACE1docs.jl/blob/{commit}{path}#{line}",
    sitename="ACE1docs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://ACEsuit.github.io/ACE1docs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Getting Started" => Any[
            "gettingstarted/installation.md",
            "gettingstarted/developing.md",
            "gettingstarted/pkg.md",
        ],
        "ACE" => Any[
            "ACE/aceintro.md",
            "ACE/datatypes.md",
            "ACE/create_ACE.md",
        ],
        "IPFitting" => Any[
            "IPFitting/IPFitting.md",
            "IPFitting/File IO.md",
            "IPFitting/Atomic Configurations in Julia.md",
            "IPFitting/Solvers.md",
            "IPFitting/Manipulating potentials.md",
        ],
        "Using ACE potentials" => Any[
            "Using_ACE/python_ase.md",
            "Using_ACE/lammps.md",
            "Using_ACE/openmm.md",
        ],
        "Tutorials" => Any[
            "tutorials/developing.md",
            "tutorials/createsavefit.md",
            "tutorials/lammps.md",
        ],

    ],
)

deploydocs(;
    repo="github.com/ACEsuit/ACE1docs.jl",
    devbranch="main",
)


text:
using ReachabilityAnalysis, SparseArrays

LazySets.set_rtol(Float64, 1e-13)
LazySets.set_ztol(Float64, 1e-13)

const eₓ = sparsevec([2], [1.], 4)
const x0 = 0.05

# model without parameter variation
function embrake_no_pv(; Tsample=1.E-4, ζ=1e-6, x0=0.05)
    # model's constants
    L = 1.e-3
    KP = 10000.
    KI = 1000.
    R = 0.5
    K = 0.02
    drot = 0.1
    i = 113.1167

    # state variables: [I, x, xe, xc]
    A = Matrix([-(R+K^2/drot)/L 0      KP/L  KI/L;
                          K/i/drot      0      0      0;
                          0             0      0      0;
                          0             0      0      0])

    EMbrake = @system(x' = Ax)

    # initial conditions
    I₀  = Singleton([0.0])
    x₀  = Singleton([0.0])
    xe₀ = Singleton([0.0])
    xc₀ = Singleton([0.0])
    X₀ = I₀ × x₀ × xe₀ × xc₀

    # reset map
    Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1., 1., -1., -Tsample, 1.], 4, 4)
    br = sparsevec([3, 4], [x0, Tsample*x0], 4)
    reset_map(X) = Ar * X + br

    # hybrid system with clocked linear dynamics
    ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
    return IVP(ha, X₀)
end

# model with parameter variation changing only 1 coefficient, corresponds to the Flow* settings in [SO15]
function embrake_pv_1(; Tsample=1.E-4, ζ=1e-6, Δ=3.0, x0=0.05)
    # model's constants
    L = 1.e-3
    KP = 10000.
    KI = 1000.
    R = 0.5
    K = 0.02
    drot = 0.1
    i = 113.1167
    p = 504. + (-Δ .. Δ)

    # state variables: [I, x, xe, xc]
    A = IntervalMatrix([-p            0      KP/L   KI/L;
                        K/i/drot      0      0      0;
                        0             0      0      0;
                        0             0      0      0])

    EMbrake = @system(x' = Ax)

    # initial conditions
    I₀  = Singleton([0.0])
    x₀  = Singleton([0.0])
    xe₀ = Singleton([0.0])
    xc₀ = Singleton([0.0])
    X₀ = I₀ × x₀ × xe₀ × xc₀

    # reset map
    Ar = sparse([1, 2, 3, 4, 4], [1, 2, 2, 2, 4], [1., 1., -1., -Tsample, 1.], 4, 4)
    br = sparsevec([3, 4], [x0, Tsample*x0], 4)
    reset_map(X) = Ar * X + br

    # hybrid system with clocked linear dynamics
    ha = HACLD1(EMbrake, reset_map, Tsample, ζ)
    return IVP(ha, X₀)
end

# maximum time by which the solution is below the x0 threshold
function max_time(sol)
    sfpos = []
    times = []

    for fp in sol
        for (j, X) in enumerate(fp)
            push!(sfpos, ρ(eₓ, X))
            push!(times, tspan(fp, j))
        end
    end
    idx = findfirst(x -> x > x0, sfpos)

    if isnothing(idx)
        return tend(sol)
    else
        return sup(times[idx-1])
    end
end

ext: jl


text:
using ZChop
using Test

@testset "basic" begin
    @test zchop(1.0) == 1.0
    @test zchop(1e-15) == 0.0
    @test zchop(1) == 1
    @test zchop(1, 3) == 0
end

@testset "complex" begin
    z1 = complex(1, 1e-15)
    zz1 = complex(1, 0.0)
    @test zchop(z1) == zz1

    z2 = complex(1e-15, 1)
    zz2 = complex(0.0, 1)
    @test zchop(z2) == zz2

    za = [z1, z2]
    zza = [zz1, zz2]
    @test zchop(za) == zza
    t = (za, za)
    @test zchop(t) == (zza, zza)
    @test zchop([za za]) == [zza zza]
end

@testset "big" begin
    @test zchop(BigFloat[1e-15, 1]) == BigFloat[0, 1]
    if Int != Int32 @test zchop(1//10^15) == 0//1 end
    @test zchop(BigInt(1)) == BigInt(1)
    @test zchop(BigInt(1), 3) == BigInt(0)
end

@testset "misc types" begin
    r = r"cat"
    @test zchop(r) == r
    @test zchop(Int) == Int
    @test zchop([true, false]) == [true, false]
    @test zchop('x') == 'x'
    @test zchop(stdout) == stdout
    @test zchop(Base.pi) == float(Base.pi)
    @test zchop(Base.pi, 4) == 0.0
    @test zchop(:a) == :a
    @test zchop(:(1+1)) == :(1+1)
    @test zchop(:(1+1e-16)) == :(1+0.0)
end

@testset "collections" begin
    @test zchop([1e-16, "cat"]) == [0.0, "cat"]
    m = [[1.0, 2.0, 1e-16] [3.0, 4.0, 5.0]]
    @test zchop(m) == [[1.0, 2.0, 0.0] [3.0, 4.0, 5.0]]
    @test zchop(:("cat" + [1e-16, 2.0])) == :("cat" + [0.0, 2.0])
    b = collect(1:5) * 1e-15
    c = zchop(b)
    zchop!(b)
    @test c == b
    d = (m[:, 1]...,)
    @test zchop(d) == (1.0, 2.0, 0.0)
    @test isa(zchop(d), Tuple)
    @test zchop((x for x in d)) == zchop(collect(d))
end

struct NumberSubtype <: Number
end

@testset "Number" begin
    x = NumberSubtype()
    @test x === zchop(x)
end

@testset "error" begin
    @test zchop((1e-15, 1, 1)) == (0.0, 1, 1)
    @test_throws MethodError zchop()
end


text:
function check_prec(LS::LoadingSequence, task::Task, precedences::Dict)
    for req in precedences[task.p]
        if !(req in LS.filled)
            print("trying to add position, but still missing")
            print(setdiff(precedences[task.p],LS.filled))
            return(false)
        end
    end
    return(true)
end

function add_task(LS::LoadingSequence, task::Task, precedences::Dict, qc::Array{1, QuayCrane})
    if !(task.c in LS.loaded) && !(task.p in LS.filled)
        if check_precedences == true
            LS.len = LS.len + 1
            LS.tasks_left = LS.P - LS.len
            push!(LS.order, task)
            push!(LS.filled, task.p)
            push!(LS.loaded,task.c)
        else
            return(false)
        end
    else
        return("already done")
    end
end

function subset_pos(set, dict_by_position, pos)
    sub_set=[]
    for task in dict_by_position[pos]
        push!(sub_set, task.c)
    end
    return(sub_set)
end

function subset_bay(J, Q, delta, q)
    sub_set=[(q-1)*(delta+1)+1:J-(Q-q)*(delta+1)]
    return(sub_set)
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
# canonical data version: 1.2.0

using Test

include("reverse-string.jl")

@testset "an empty string" begin
    @test myreverse("") == ""
end

@testset "a word" begin
    @test myreverse("robot") == "tobor"
end

@testset "a capitalized word" begin
    @test myreverse("Ramen") == "nemaR"
end

@testset "a sentence with punctuation" begin
    @test myreverse("I'm hungry!") == "!yrgnuh m'I"
end

@testset "a palindrome" begin
    @test myreverse("racecar") == "racecar"
end

@testset "an even-sized word" begin
    @test myreverse("drawer") == "reward"
end

@testset "reversing a string twice" begin
    @test myreverse(myreverse("gift")) == "gift"
end

@testset "emoji" begin
    @test myreverse("hi 🐱") == "🐱 ih"
end

if @isdefined(TEST_GRAPHEMES)
    @eval @testset "graphemes" begin
        @test myreverse("as⃝df̅") == "f̅ds⃝a"
        @test myreverse("hi 👋🏾") == "👋🏾 ih"
    end
end


text:
"""

$(TYPEDEF)

Structure to contain a list of all the distances smaller than
the cutoff.

$(TYPEDFIELDS)

"""
struct CutoffDistances

    # To store all distances smaller than the cutoff
    nd::Vector{Int} # Number of small distances (vector size = 1), just to be mutable
    d::Vector{Float64} # All distances smaller than the cutoff
    iat::Vector{Int} # atom of the solute
    jat::Vector{Int} # atom of the solvent
    imol::Vector{Int} # molecule of the solute
    jmol::Vector{Int} # molecule of the solvent

    # Size of the arrays
    maxdim::Vector{Int}

end

# Generator
CutoffDistances(natoms::Int) = CutoffDistances(
    zeros(Int, 1), # nd
    zeros(Float64, natoms), # d
    zeros(Int, natoms), # iat
    zeros(Int, natoms), # jat
    zeros(Int, natoms), # imol
    zeros(Int, natoms), # jmol
    [natoms],
) # maxdim

# Function that zeroes all the values in this structure
function reset!(c::CutoffDistances)
    c.nd[1] = 0
    @. c.d = 0.0
    @. c.iat = 0
    @. c.jat = 0
    @. c.imol = 0
    @. c.jmol = 0
    return nothing
end

function reset!(c::Vector{CutoffDistances})
    for i = 1:size(c, 1)
        reset!(c[i])
    end
    return nothing
end


text:
include("run_cutest.jl")

function perturb_cons(nr::CUTEstModel,scale::Float64)
    srand(1)

    lcon = deepcopy(nr.meta.lcon);
    ucon = deepcopy(nr.meta.ucon);
    lvar = deepcopy(nr.meta.lvar);
    uvar = deepcopy(nr.meta.uvar);

    f(x) = obj(nr, x)
    #shift = scale * (rand(length(lcon)) - 0.5)
    shift = scale * ones(length(lcon)) #(rand(length(lcon)) - 0.5)
    c(x) = cons(nr, x) #
    g(x) = grad(nr, x)
    J(x) = jac(nr, x)
    H(x; obj_weight=1.0, y=zeros) = hess(nr, x, obj_weight=obj_weight, y=y)
    Hcoord(x; obj_weight=1.0, y=zeros) = hess_coord(nr, x, obj_weight=obj_weight, y=y)
    Jcoord(x) = jac_coord(nr, x)

    snlp = SimpleNLPModel(f, nr.meta.x0;
    lvar = lvar, uvar = uvar,
    lcon = lcon - shift, ucon = ucon - shift,
    name = "Generic",
    g=g,
    c=c,
    J=J,
    H=H,
    Jcoord=Jcoord,
    Hcoord=Hcoord)

    return snlp
end

function run_infeas(test_name, problem_list, my_par, solve_func!)
    summary = Dict{String, problem_summary2}()

    if_mkdir("../results/$test_name")
    if_mkdir("../results/$test_name/log")
    if_mkdir("../results/$test_name/jld")

    for problem_name in problem_list
        nr = CUTEstModel(problem_name)
        snlp = perturb_cons(nr, 1.0)
        summary[problem_name] = solve_func!(snlp, problem_name, test_name, my_par)
        finalize(nr)

        write_summary("../results/$(test_name)/summary.txt", summary)
        save("../results/$(test_name)/summary.jld","summary",summary, "pars", my_par)
    end
end

problem_list = default_list()

my_par = OnePhase.Class_parameters()
#run_infeas("one_phase/infeas-Dec26", problem_list, my_par, one_phase_solve)
#run_infeas("ipopt/infeas-Dec26", problem_list, my_par, ipopt_solve)
run_infeas("one_phase/Jan4/infeas", problem_list, my_par, one_phase_run_and_store)


text:
using LaTeXStrings
using PyPlot
using Printf

function dsdemo7()
    println("Invariant Set for MOD2 (find2dPIS)")

    A = [1 0; 1 1]
    B = [1 -1; 1 -2]
    C = [0 1]
    D = [0 0]
    ABCD = [A B; C D]
    order = 2

    fig = figure(1)
    fig.clf()
    fig.canvas.set_window_title("Working Invariant Set")

    u = 1 / pi
    s = find2dPIS(u, ABCD, true)

    N = 10000
    skip = 100
    _, x, = simulateDSM(u * ones(N+skip), ABCD, 2)
    x = x[:,1+skip:N+skip]
    nv = size(s, 2)
    splus, eplus, sminus, eminus = dssplit2d(u, ABCD, s)
    Buv = B * [u, 1]
    s1 = A*splus + Buv*ones(1,size(splus, 2))
    Buv = B * [u, -1]
    s2 = A*sminus + Buv*ones(1,size(sminus, 2))
    ns = [s1 s2]
    out = outconvex2d(ns, s)

    fig = figure(2)
    fig.clf()
    fig.canvas.set_window_title("Final Invariant Set")

    grid()
    dotplot(x, "k", ".")
    polyplot(s, "b")
    polyplot(s1, "m")
    polyplot(s2, "c")
    outi = (out .!= 0)[:]
    dotplot(ns[:,outi], "r", "s")
    str = @sprintf("Final Object: %d image vertices outside", sum(outi))
    title(str)
    xlabel(L"x_1")
    ylabel(L"y_1")

    @printf("%d points from the %d simulated states are outside.\n", sum(outconvex2d(x, s)), N)
    @printf("%d image points are outside.\n", sum(out))
    @printf("The returned polygon has %d vertices.\n", size(s, 2))

    return nothing
end

max_stars_repo_name: Mahmoud-Kharsa/DeltaSigma

max_issues_repo_name: Mahmoud-Kharsa/DeltaSigma

max_forks_repo_name: Mahmoud-Kharsa/DeltaSigma


text:
import  Base: eltype, eps, length, ndims, real, size, *, \,
        A_mul_B!, Ac_mul_B, Ac_mul_B!

export  A_mul_B

# Improve readability of iterative methods
\(f::Function, b) = f(b)
*(f::Function, b) = f(b)

#### Type-handling
Adivtype(A, b) = typeof(one(eltype(b))/one(eltype(A)))
Amultype(A, x) = typeof(one(eltype(A))*one(eltype(x)))
if VERSION < v"0.4.0-dev+6068"
    real{T<:Real}(::Type{Complex{T}}) = T
    real{T<:Real}(::Type{T}) = T
end
eps{T<:Real}(::Type{Complex{T}}) = eps(T)

function randx(A, b)
    T = Adivtype(A, b)
    x = initrand!(Array(T, size(A, 2)))
end

function zerox(A, b)
    T = Adivtype(A, b)
    x = zeros(T, size(A, 2))
end

#### Numerics
function update!(x, α::Number, p::AbstractVector)
    for i = 1:length(x)
        x[i] += α*p[i]
    end
    x
end

function initrand!(v::Vector)
    _randn!(v)
    nv = norm(v)
    for i = 1:length(v)
        v[i] /= nv
    end
    v
end
_randn!(v::Array{Float64}) = randn!(v)
_randn!(v) = copy!(v, randn(length(v)))

#### Errors
export PosSemidefException

type PosSemidefException <: Exception
    msg :: AbstractString
    PosSemidefException(msg::AbstractString="Matrix was not positive semidefinite") = new(msg)
end

#### Support for linear-operators-defined-as-functions
#
# Suppose you have a function implementing multiplication of b by A.
# This function can have any syntax, but for the purposes of
# illustration let's suppose it's defined as
#   mulbyA!(output, b, Adata)
# Where Adata might be some parameters that your function needs.
# You can represent it as a linear operator using
#   A = MatrixFcn{T}(m, n, (output, b) -> mulbyA!(output, b, Adata))
# where T is the "element type" of Adata.
# Note that there are a couple of requirements:
#   - mulbyA! stores the result in the pre-allocated output
#   - mulbyA! should also return output as its sole return value

# If the algorithm also needs multiplication by A', use MatrixCFcn
# instead and initialize it with both functions.

export MatrixFcn, MatrixCFcn

abstract AbstractMatrixFcn{T}

type MatrixFcn{T} <: AbstractMatrixFcn{T}
    m::Int
    n::Int
    mul::Function
end

type MatrixCFcn{T} <: AbstractMatrixFcn{T}
    m::Int
    n::Int
    mul::Function
    mulc::Function
end

MatrixFcn(A::AbstractMatrix) = MatrixFcn{eltype(A)}(size(A, 1), size(A, 2), (output, x)-> A_mul_B!(output, A, x))
MatrixCFcn(A::AbstractMatrix) = MatrixFcn{eltype(A)}(size(A, 1), size(A, 2), (output, x)->A_mul_B!(output, A, x), (output, b) -> Ac_mul_B(output, A, b))

eltype{T}(op::AbstractMatrixFcn{T}) = T

ndims(op::AbstractMatrixFcn) = 2

size(op::AbstractMatrixFcn, dim::Integer) = (dim == 1) ? op.m :
                                            (dim == 2) ? op.n : 1
size(op::AbstractMatrixFcn) = (op.m, op.n)

length(op::AbstractMatrixFcn) = op.m*op.n

*(op::AbstractMatrixFcn, b) = A_mul_B(op, b)

A_mul_B{R,S}(op::AbstractMatrixFcn{R}, b::AbstractVector{S}) = A_mul_B!(Array(promote_type(R,S), op.m), op, b)

A_mul_B!(output, op::AbstractMatrixFcn, b) = op.mul(output, b)

Ac_mul_B{R,S}(op::MatrixCFcn{R}, b::AbstractVector{S}) = Ac_mul_B!(Array(promote_type(R,S), op.n), op, b)

Ac_mul_B!(output, op::AbstractMatrixFcn, b) = op.mulc(output, b)


text:
using Mosek
using MosekTools
using LinearAlgebra
using MinimaxEstimation
using Plots
using Random

Random.seed!(10)
n_steps = 20
F = [1.1 -.5 0.1
    1 0 0
    0 1 0]
B = [-1;2;3]
H = [1 0 0]
x0 = [0; 0; 0]
Q = I(3)
R = 1
P = I(3)
gamma = 2

kf1 = KalmanFilter(x0, F, H,B, P, Q, R)
kf2 = KalmanFilter(x0, -F, H,B, P, Q, R)
mini = MinimaxMMAE(copy([kf1, kf2]), gamma, Mosek.Optimizer)
bayesian = BayesianMMAE(copy([kf1, kf2]))
v = randn((1,n_steps))
w = randn((3, n_steps))
u(k) = sin(k/5)
y = zeros(n_steps+1)
z = zeros(n_steps+1)
y[1] = 0
z[1] = 0

zhat_mini = zeros(n_steps+1)
vals_mini = zeros(n_steps+1)
zhat_bayesian = zeros(n_steps+1)

zhat_mini[1]= predict(mini, true)[1][1]
zhat_bayesian[1] = predict(bayesian)[1]
x = x0
for k=1:n_steps
    update!(mini,y[k], u(k))
    update!(bayesian, y[k], u(k))

    zhat_mini[k+1] = predict(mini, true)[1][1]
    zhat_bayesian[k+1] = predict(bayesian)[1]
    global x = F*x + B*u(k) + w[:,k]
    y[k+1] = (H*x)[1] + v[k]
    z[k+1] = (H*x)[1]
    k +=1
end
print("||z - zhat_mini|| = ", norm(z - zhat_mini), "\n")
print("||z - zhat_bayesian|| = ", norm(z - zhat_bayesian), "\n")
print("||z_{n+1} - y_n|| = ", norm(z[2:end]-y[1:end-1]), "\n")
plot(
    [z zhat_mini zhat_bayesian],
    legend=:bottomright,
    labels = ["z" "zhat_minimax" "zhat_bayesian"],
    linewidth=3,
    markershape = :o,
    xlabel = "time-step",
    ylabel = "z")

savefig("experiment.png")


text:
using ProbabilisticCircuits
using DocumenterTools

#DocumenterTools.genkeys(ProbabilisticCircuits)

DocumenterTools.genkeys(user="khosravipasha", repo="git@github.com:Juice-jl/ProbabilisticCircuits.jl.git")


text:
CachedOperator(::Type{Matrix},op::Operator;padding::Bool=false) =
    CachedOperator(op,Array{eltype(op)}(0,0),padding)

# Grow cached operator

function resizedata!(B::CachedOperator{T,Matrix{T}},n::Integer,m::Integer) where T<:Number
    if n > size(B,1) || m > size(B,2)
        throw(ArgumentError("Cannot resize beyound size of operator"))
    end

    # this does nothing if already in dimensions
    N,M=size(B.data)
    if n > N && m > M
        B.data = unsafe_resize!(B.data,n,m)
    elseif n > N
        B.data = unsafe_resize!(B.data,n,:)
    elseif m > M
        B.data = unsafe_resize!(B.data,:,m)
    end

    if n ≤ B.datasize[1] && m ≤ B.datasize[2]
        # do nothing
        B
    elseif n ≤ B.datasize[1]
        kr,jr=1:B.datasize[1],B.datasize[2]+1:m
        B.data[kr,jr] = B.op[kr,jr]
        B.datasize = (B.datasize[1],m)
        B
    elseif m ≤ B.datasize[2]
        kr,jr=B.datasize[1]+1:n,1:B.datasize[2]
        B.data[kr,jr] = B.op[kr,jr]
        B.datasize = (n,B.datasize[2])
        B
    else
        # resize rows then columns
        resizedata!(resizedata!(B,n,B.datasize[2]),n,m)
    end
end

# Apply Householder

function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
                      B::AbstractVector{T},tolerance,maxlength) where {RR,T}
    A = parent(Ac)
    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    H=A.QR.H
    M=size(H,1)
    m=length(B)
    Y=pad(B,m+M+10)

    k=1
    yp=view(Y,1:M)
    while (k ≤ m+M || norm(yp) > tolerance )
        if k > maxlength
            @warn "Maximum length $maxlength reached."
            break
        end

        if k+M-1>length(Y)
            pad!(Y,2*(k+M))
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,k+M+50)
            H=A.QR.H
        end

        wp=view(H,:,k)
        yp=view(Y,k:k+M-1)

        dt=dot(wp,yp)
        LinearAlgebra.axpy!(-2*dt,wp,yp)
        k+=1
    end
    resize!(Y,k)  # chop off zeros
end

# BLAS apply Q

function mulpars(Ac::Adjoint{T,<:QROperatorQ{QROperator{RR,Matrix{T},T},T}},
                       B::AbstractVector{T},
                       tolerance,maxlength) where {RR,T<:BlasFloat}

    A = parent(Ac)

    A_dim = size(A,1)

    if length(B) > A.QR.ncols
        # upper triangularize extra columns to prepare for \
        resizedata!(A.QR,:,length(B)+size(A.QR.H,1)+10)
    end

    if size(A.QR.H,1) == 1  # diagonal scaling, avoid growing
        diagv=view(A.QR.H,1,1:length(B))
        ret=similar(B)
        @simd for k=1:length(ret)
            @inbounds ret[k]=(1-2A.QR.H[1,k]^2)*B[k]
        end
        return ret
    end

    H=A.QR.H
    h=pointer(H)
    M=size(H,1)
    st=stride(H,2)

    sz=sizeof(T)

    m=length(B)
    Y=pad(B,m+M+10)
    y=pointer(Y)

    k=1
    yp=y
    while (k ≤ min(m+M,A_dim) || BLAS.nrm2(M,yp,1) > tolerance )
        if k > maxlength
            @warn "Maximum length $maxlength reached."
            break
        end

        if k+M-1>length(Y)
            pad!(Y,2*(k+M))
            y=pointer(Y)
        end
        if k > A.QR.ncols
            # upper triangularize extra columns to prepare for \
            resizedata!(A.QR,:,2*(k+M))
            H=A.QR.H
            h=pointer(H)
        end

        wp=h+sz*st*(k-1)
        yp=y+sz*(k-1)

        dt = BandedMatrices.dot(M,wp,1,yp,1)
        BLAS.axpy!(M,-2*dt,wp,1,yp,1)
        k+=1
    end
    resize!(Y,min(k,A_dim))  # chop off zeros
end


text:
import Iterators: product

using NURBS.Utils
using NURBS.Bases

info("Testing NURBS.Bases")

# BSplineBasis inner constructor
# ========================================================================

basis = BSplineBasis(linspace(0, 4, 5), 3)
@test basis.knots == [0, 0, 0, 1, 2, 3, 4, 4, 4]
@test basis.order == 3
@test basis.deriv == 0

basis = BSplineBasis([0, 0, 0, 1, 2, 3, 4, 4, 4], 3, 1; extend=false)
@test basis.knots == [0, 0, 0, 1, 2, 3, 4, 4, 4]
@test basis.order == 3
@test basis.deriv == 1

@test_throws ArgumentError BSplineBasis([2, 2, 1, 0, 0], 2, 0; extend=false)
@test_throws ArgumentError BSplineBasis([0, 0, 1, 2, 2], 3, 0; extend=false)
@test_throws ArgumentError BSplineBasis([0, 0, 1, 2, 3], 2, 0; extend=false)
@test_throws ArgumentError BSplineBasis([0, 1, 2, 3, 3], 2, 0; extend=false)
@test_throws ArgumentError BSplineBasis([0, 1], 2, 2)

# BSplineBasis outer constructors
# ========================================================================

basis = BSplineBasis(0, 1, 4, 4)
@test basis.knots == [0, 0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1, 1]
@test basis.order == 4
@test basis.deriv == 0

# BSplineBasis length, size and domain
# ========================================================================

basis = BSplineBasis(0, 1, 4, 4)
@test length(basis) == 7
@test size(basis) == (7,)
@test domain(basis) == Interval(0, 1)

basis = BSplineBasis(0, 2, 5, 3)
@test length(basis) == 7
@test size(basis) == (7,)
@test domain(basis) == Interval(0, 2)

# BSplineBasis evaluation
# ========================================================================

# Single point, no derivatives
function testeval_snn(b, x, vals, idxs)
    (rvals, ridxs) = b(x)
    @test ridxs == idxs
    @test_approx_eq rvals vals

    coeffs = rand(length(b))
    @test_approx_eq dot(coeffs[ridxs], vals) b(x, coeffs)

    for (v, i) in zip(vals, idxs)
        @test_approx_eq b[i](x) v
    end

    for i in 1:length(b)
        if i ∉ idxs
            @test_approx_eq b[i](x) 0.0
        end
    end
end

basis = BSplineBasis(0, 1, 2, 2)
testeval_snn(basis, 0, [1, 0], 1:2)
testeval_snn(basis, 0.2, [0.6, 0.4], 1:2)
testeval_snn(basis, 0.4, [0.2, 0.8], 1:2)
testeval_snn(basis, 0.5, [1, 0], 2:3)
testeval_snn(basis, 0.75, [0.5, 0.5], 2:3)
testeval_snn(basis, 1, [0, 1], 2:3)

basis = BSplineBasis(0, 1, 1, 3)
testeval_snn(basis, 0, [1, 0, 0], 1:3)
testeval_snn(basis, 0.3, [0.49, 0.42, 0.09], 1:3)
testeval_snn(basis, 0.8, [0.04, 0.32, 0.64], 1:3)

basis = BSplineBasis(0, 2, 2, 4)
testeval_snn(basis, 0.2, [0.512, 0.434, 0.052, 0.002], 1:4)
testeval_snn(basis, 0.7, [0.027, 0.49525, 0.392, 0.08575], 1:4)
testeval_snn(basis, 1.0, [0.25, 0.5, 0.25, 0], 2:5)
testeval_snn(basis, 1.5, [0.03125, 0.25, 0.59375, 0.125], 2:5)
testeval_snn(basis, 2.0, [0, 0, 0, 1], 2:5)

basis = BSplineBasis([0, 1, 3], 4)
testeval_snn(basis, 0.5, [0.125, 0.680555555556, 0.18055555555555564, 0.013888888888888902], 1:4)
testeval_snn(basis, 1.4, [0.227555555555557, 0.483555555556, 0.280888888889, 0.008], 2:5)
testeval_snn(basis, 2.1, [0.0405, 0.26325, 0.529875, 0.166375], 2:5)

# Multiple points, no derivatives
function testeval_mnn(b, xs)
    res = b(xs)
    for (x, (vals, idxs)) in zip(xs, res)
        @test typeof(vals) <: Vector
        testeval_snn(b, x, vals, idxs)
    end

    coeffs = rand(length(b), 3)
    resc = b(xs, coeffs)
    for (i, (vals, idxs)) in enumerate(res)
        @test_approx_eq vals'*coeffs[idxs,:] resc[i,:]
    end
end

basis = BSplineBasis(0, 2, 2, 2)
pts = [0.0, 0.4, 0.8, 1.0, 1.5, 2.0]
testeval_mnn(basis, pts)
@test_approx_eq basis[1](pts) [1.0; 0.6; 0.2; zeros(3)]
@test_approx_eq basis[2](pts) [0.0, 0.4, 0.8, 1.0, 0.5, 0.0]
@test_approx_eq basis[3](pts) [zeros(4); 0.5; 1.0]

basis = BSplineBasis(0, 3, 4, 3)
pts = [0.0, 0.4, 0.8, 1.0, 1.5, 2.0, 2.2, 2.8, 3.0]
testeval_mnn(basis, pts)
@test_approx_eq basis[1](pts) [1.0; 0.2177777777777774; zeros(7)]
@test_approx_eq basis[3](pts) [0.0, 0.1422222222222222, 0.56222222222222222, 0.722222222222222,
                               0.5, 0.0555555555555555, 0.00222222222222222, 0.0, 0.0]
@test_approx_eq basis[6](pts) [zeros(7); 0.53777777777777777; 1.0]

basis = BSplineBasis([0, 1, 3, 4], 4)
pts = [0.0, 0.8, 1.6, 2.4, 3.2, 4.0]
testeval_mnn(basis, pts)
@test_approx_eq basis[2](pts) [0.0; 0.5795555555555555; 0.15244444444444444; 0.012; zeros(2)]
@test_approx_eq basis[4](pts) [0.0, 0.0426666666666666, 0.29333333333333333, 0.542222222222222,
                               0.3697777777777777, 0.0]

# Single point, derivatives
function testeval_sdn(b, x, dt=1e-6, tol=1e-8)
    dvals, didxs = deriv(b)(x)
    (dvp, didp), (dvn, didn) = b([x+dt, x-dt])
    @test didxs == didp == didn
    nderiv = (dvp - dvn) / 2dt
    @test_approx_eq_eps nderiv dvals tol

    coeffs = rand(length(b), 4)
    @test_approx_eq dvals'*coeffs[didxs,:] deriv(b)(x, coeffs)

    for (v, i) in zip(dvals, didxs)
        @test_approx_eq deriv(b[i])(x) v
    end

    for i in 1:length(b)
        if !(i in didxs)
            @test_approx_eq deriv(b[i])(x) 0.0
        end
    end
end

basis = BSplineBasis(0, 1, 2, 2)
for pt in [0.25, 0.51, 0.9]
    testeval_sdn(basis, pt)
end
basis = deriv(basis)
@test_throws ArgumentError deriv(basis)

basis = BSplineBasis([0, 0.5, 0.7, 0.8], 3)
for i in 1:2
    for pt in [0.1, 0.3, 0.6, 0.73]
        testeval_sdn(basis, pt)
    end
    basis = deriv(basis)
end
@test_throws ArgumentError deriv(basis)

basis = BSplineBasis(0, 2, 4, 4)
for i in 1:3
    for pt in [0.1, 0.7, 1.2, 1.6]
        testeval_sdn(basis, pt)
    end
    basis = deriv(basis)
end
@test_throws ArgumentError deriv(basis)

# Multiple points, derivatives
function testeval_mdn(b, xs, dt=1e-6, tol=1e-8)
    res = deriv(b)(xs)
    resn = b(xs+dt)
    resp = b(xs-dt)
    for ((v, id), (nv, nid), (pv, pid)) in zip(res, resn, resp)
        @test typeof(v) <: Vector
        @test id == nid == pid
        nderiv = (nv - pv) / 2dt
        @test_approx_eq_eps nderiv v tol
    end

    coeffs = rand(length(b))
    resc = deriv(b)(xs, coeffs)
    for (i, (v, id)) in enumerate(res)
        @test_approx_eq dot(v, coeffs[id]) resc[i]
    end
end

basis = BSplineBasis([0, 0.1, 0.8, 1.2], 2)
testeval_mdn(basis, [0.02, 0.3, 0.7, 0.9])

basis = BSplineBasis(6, 12, 3, 3)
for i in 1:2
    testeval_mdn(basis, [6.43, 7.8, 8.3, 9.1, 10.1, 11.8])
    basis = deriv(basis)
end

basis = BSplineBasis(5, 15, 10, 4)
for i in 1:3
    testeval_mdn(basis, collect(linspace(5.5, 14.5, 10)))
    basis = deriv(basis)
end

# BSplineBasis indexing
# ========================================================================

basis = BSplineBasis(0, 4, 4, 2)
@test_throws BoundsError basis[0]
@test_throws BoundsError basis[6]
@test basis[1].basis == basis
@test basis[2].index == 2
@test basis[5].deriv == 0
@test domain(basis[1]) == Interval(0, 1)
@test domain(basis[2]) == Interval(0, 2)
@test domain(basis[5]) == Interval(3, 4)
@test deriv(basis[1]).basis == basis
@test deriv(basis[2]).index == 2
@test deriv(basis[3]).deriv == 1
@test_throws ArgumentError deriv(deriv(basis[2]))

basis = BSplineBasis([0, 4, 5, 6], 4)
@test_throws BoundsError basis[7]
@test basis[6].basis == basis
@test basis[2].deriv == 0
@test basis[1].index == 1
@test domain(basis[1]) == Interval(0, 4)
@test domain(basis[2]) == Interval(0, 5)
@test domain(basis[3]) == Interval(0, 6)
@test domain(basis[4]) == Interval(0, 6)
@test domain(basis[5]) == Interval(4, 6)
@test domain(basis[6]) == Interval(5, 6)
@test deriv(basis[1]).basis == basis
@test deriv(basis[4]).index == 4
@test deriv(deriv(basis[4])).deriv == 2
@test_throws ArgumentError deriv(deriv(deriv(deriv(basis[3]))))

basis = deriv(basis)
@test deriv(basis[2]).deriv == 2
@test deriv(basis[2]).basis.deriv == 1
@test_throws ArgumentError deriv(deriv(deriv(basis[1])))


text:
type Bam_pileup
    b::Ptr{Record}
    qpos::Cint
    indel::Cint
    level::Cint
    flag::Cuint
end
typealias Bam_plp Ptr{Void}
typealias Bam_plp_auto_f Ptr{Void}
typealias Bam_mplp Ptr{Void}

#todo ret status signal handle
function bam_plp_init(func::Function,data::Ptr{Void})
    ccall((:bam_plp_init,libhts),Bam_pileup,(Function,Ptr{Void}),func,data)
end

function bam_plp_destroy(iter::Bam_plp)
    ccall((:bam_plp_destroy,libhts),Void,(Ptr{Bam_plt},),iter)
end

function bam_plp_push(iter::Bam_plp,b::Ptr{Record})
    ccall((:bam_plp_push, libhts), Cint, (Bam_plp, Ptr{Record}),iter,b)
end

function bam_plp_next(iter::Bam_plp, _tid::Ptr{Cint}, _pos::Ptr{Cint}, _n_plp::Ptr{Cint})
    ccall((:bam_plp_next, libhts), Ptr{Bam_pileup}, (Bam_plp,Ptr{Cint},Ptr{Cint},Ptr{Cint}),iter,_tid,_pos,_n_plp)
end
function bam_plp_auto(iter::Bam_plp, _tid::Ptr{Cint}, _pos::Ptr{Cint}, _n_plp::Ptr{Cint})
    ccall((:bam_plp_auto, libhts), Ptr{Bam_pileup}, (Bam_plp,Ptr{Cint},Ptr{Cint},Ptr{Cint}),iter,_tid,_pos,_n_plp)
end
function bam_plp_set_maxcnt(iter::Bam_plp,maxcnt::Cint)
    ccall((:bam_plp_maxcnt, libhts),Void,(Bam_plp,Cint),iter,maxcnt)
end
function bam_mplp_init(n::Cint,func::Bam_plp_auto_f,data::Ptr{Ptr{Void}})
    ccall((:bam_mplp_init, libhts),Bam_mplp,(Cint,Bam_plp_auto_f,Ptr{Ptr{Void}}),n,func,data)
end
function bam_mplp_init_overlaps(iter::Bam_mplp)
    ccall((:bam_mplp_init_overlaps, libhts),Void,(Bam_mplp,), iter)
end
function bam_mplp_destroy(iter::Bam_mplp)
    ccall((:bam_mplp_destroy, libhts),Void,(Bam_mplp,), iter)
end
function bam_mplp_set_maxcnt(iter::Bam_mplp,maxcnt::Cint)
    ccall((:bam_mplp_set_maxcnt, libhts),Void,(Bam_mplp,Cint),iter,maxcnt)
end
function bam_mplp_auto(iter::Bam_mplp, _tid::Ptr{Cint}, _pos::Ptr{Cint}, _n_plp::Ptr{Cint}, plp::Ptr{Ptr{Bam_pileup}})
    ccall((:bam_mplp_auto, libhts),Cint,(Bam_plp,Ptr{Cint},Ptr{Cint},Ptr{Cint}),iter,_tid,_pos,_n_plp, plp)
end


text:
using CAGD;

model = CAGD.Model([
    2.0 1.0 3.0 2.0 1.0 -3.
    2.0 4.0 5.0 1.0 9.0 1.0
    1.0 4.0 1.0 4.0 3.0 4.0
]);

CAGD.addModelCells!(model, 1, [[1,2],[2,3],[1,3],[1,4],[2,4],[1,5],[2,5],[1,6],[2,6]]);
CAGD.addModelCells!(model, 2, [[1,2,3],[1,4,5],[1,6,7],[1,8,9]]);

#test = eval_ord_faces(model, 1)
τ = 1;
faces = model.T[2][:, τ].nzind
edge = model.T[1][τ, :].nzind
FVs = [setdiff((model.T[1]'*model.T[2][f, :]).nzind, edge) for f in faces]
Gface, Gfaceidx = CAGD.getModelCellGeometry(model, 2, faces[1], true)
M = CAGD.build_projection_matrix(Gface, y0 = true)
PG = (M * [model.G; ones(1, size(model, 0, 2))])[2:3, :]
angles = [[mod(atan(PG[2, v], PG[1, v]), 2π) for v in fv] for fv in FVs]
angles = [[α for α in a if !(α ≈ 0)] for a in angles]
angles = [[faces[i], angles[i][1]] for i = 1 : length(faces)]
sort!(angles, lt = (a, b) -> a[2] < b[2])
map(a->Int(a[1]), angles)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
JULIA="${JULIA:-julia --color=yes --startup-file=no}"
export JULIA_PROJECT="$(dirname ${BASH_SOURCE[0]})"
exec ${JULIA} "${BASH_SOURCE[0]}" "$@"
=#

printstyled(let
    message = " Warnings below (if any) are fine. "
    margin = (displaysize(stdout)[2] - length(message)) ÷ 2
    ("=" ^ margin) * message * ("=" ^ margin)
end, color=:blue)
println()
flush(stdout)
using Documenter
using DiffEqBase: AbstractODEProblem
using Setfield: Lens
import Plots
import OrdinaryDiffEq
flush(stdout)
flush(stderr)
printstyled("=" ^ displaysize(stdout)[2], color=:blue)
println()

using Bifurcations

let strict = get(ENV, "DOCS_STRICT", "yes") == "yes"
    @info """
    `makedocs` with:
    strict = $strict
    """
    makedocs(
        # modules = [Bifurcations],
        pages = [
            "Home" => "index.md",
            "API" => "api.md",
            "Examples" => [
                hide("examples.md"),
                "examples/calcium.md",
                "examples/van_der_pol.md",
                "examples/morris_lecar.md",
            ],
            hide("Internals" => "internals.md"),
        ],
        repo = "https://github.com/tkf/Bifurcations.jl/blob/{commit}{path}#L{line}",
        sitename = "Bifurcations.jl",
        authors = "Takafumi Arakaki",
        strict = strict,
    )
end


text:
using ImageFiltering, ImageFeatures, OffsetArrays, CoordinateTransformations, PaddedViews, Distances, StaticArrays, Kalman, Plots, Gtk, Images, ImageView, GtkReactive, HDF5, ProgressMeter, Interpolations
gr()
default(show=false)
const Mark = SVector{2, Float64}
const sz = 100
const intensity_radius = 2
struct Track
    coordinates::Vector{Mark}
    color::RGB{Float64}
    id::Int
end
img = rand(100,100)
formatlabel(x) = 2x+1
rs = [Track([Mark(rand(), rand()) for j in 1:100], rand(RGB{Float64}), rand(1:10)) for i in 1:5]
path = homedir()
function saveoverview(img::Matrix{Float64}, formatlabel::Function, rs::Vector{Track}, path::String)
    img .= adjust_gamma(img, 2.4)
    mM = quantile(vec(img), [.5, .995])
    img .= imadjustintensity(img, mM)
    heatmap(flipdim(img, 1), aspect_ratio = 1, yformatter = formatlabel, xformatter = formatlabel, xlabel = "X (mm)", ylabel = "Y (mm)", color=:inferno, yflip=true, colorbar=false, legend=false, size=(sz,sz), dpi=50/3)
    for r in rs
        x = last.(r.coordinates)
        y = first.(r.coordinates)
        plot!(x, y, color=r.color, linestyle = :dot, linewidth = 5, linewidth = intensity_radius, annotations=(x[1],y[1],text(string(r.id), :left, r.color, 25))) # works
    end
    png(joinpath(path, "roots"))
end
saveoverview(img, formatlabel, rs, path)


text:
# Adjoint test for F and J
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: May 2020
#

using Distributed

parsed_args = parse_commandline()

nlayer = parsed_args["nlayer"]
tti = parsed_args["tti"]
fs =  parsed_args["fs"]

# # Set parallel if specified
nw = parsed_args["parallel"]
if nw > 1 && nworkers() < nw
    addprocs(nw-nworkers() + 1; exeflags=["--code-coverage=user", "--inline=no", "--check-bounds=yes"])
end

@everywhere using JUDI, LinearAlgebra, Test, Distributed

### Model
model, model0, dm = setup_model(parsed_args["tti"], parsed_args["nlayer"])
q, srcGeometry, recGeometry, info = setup_geom(model; nsrc=nw)
dt = srcGeometry.dt[1]

tol = 5f-4
(tti && fs) && (tol = 5f-3)
###################################################################################################
# Modeling operators
@testset "Adjoint test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin

    opt = Options(sum_padding=true, dt_comp=dt, free_surface=parsed_args["fs"])
    F = judiModeling(model0, srcGeometry, recGeometry; options=opt)

    # Nonlinear modeling
    y = F*q

    # Generate random noise data vector with size of d_hat in the range of F
    wave_rand = (.5f0 .+ rand(Float32, size(q.data[1]))).*q.data[1]
    q_rand = judiVector(srcGeometry, wave_rand)

    # Forward-adjoint
    d_hat = F*q_rand
    q_hat = adjoint(F)*y

    # Result F
    a = dot(y, d_hat)
    b = dot(q_rand, q_hat)
    @printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", a, b, (a - b)/(a + b))
    @test isapprox(a/(a+b), b/(a+b), atol=tol, rtol=0)

    # Linearized modeling
    J = judiJacobian(F, q)

    ld_hat = J*dm
    dm_hat = J'*y

    c = dot(ld_hat, y)
    @show norm(dm), norm(dm_hat)
    d = dot(dm_hat, dm)
    @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, (c - d)/(c + d))
    @test isapprox(c/(c+d), d/(c+d), atol=tol, rtol=0)

end
###################################################################################################
# Extended source modeling
if parsed_args["tti"] &&  parsed_args["fs"]
    # FS + tti leads to slightly worst (still fairly ok) accuracy
    tol = 5f-3
end

@testset "Extended source adjoint test with $(nlayer) layers and tti $(tti) and freesurface $(fs)" begin

    opt = Options(sum_padding=true, dt_comp=dt, free_surface=parsed_args["fs"])
    F = judiModeling(info, model0, srcGeometry, recGeometry; options=opt)
    # Nonlinear modeling
    y = F*q

    Pr = judiProjection(info, recGeometry)
    Fw = judiModeling(info, model0; options=opt)
    Pw = judiLRWF(info, q.data[1])
    Fw = Pr*Fw*adjoint(Pw)

    # Extended source weights
    w = .5f0 .+ rand(model0.n...)
    parsed_args["fs"] ? w[:, 1:2] .= 0f0 : nothing
    w = judiWeights(w; nsrc=nw)

    # Forward-Adjoint computation
    dw_hat = Fw*w
    w_hat = adjoint(Fw)*y

    # Result F
    a = dot(y, dw_hat)
    b = dot(w, w_hat)
    @printf(" <F x, y> : %2.5e, <x, F' y> : %2.5e, relative error : %2.5e \n", a, b, (a - b)/(a + b))
    @test isapprox(a/(a+b), b/(a+b), atol=tol, rtol=0)

    # Linearized modeling
    Jw = judiJacobian(Fw, w)

    ddw_hat = Jw*dm
    dmw_hat = adjoint(Jw)*y

    c = dot(y, ddw_hat)
    d = dot(dm, dmw_hat)
    @printf(" <J x, y> : %2.5e, <x, J' y> : %2.5e, relative error : %2.5e \n", c, d, (c - d)/(c + d))
    @test isapprox(c/(c+d), d/(c+d), atol=tol, rtol=0)

end

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
# This file "resolves" the variable names used by the World Ocean Database (WOD).
# The idea is that the WOD has a naming convention (plus some exceptions)
# that are taken care of by the code below.

# Fallback for using direct download
function fallback_download(remotepath, localdir)
    @assert(isdir(localdir))
    filename = basename(remotepath)  # only works for URLs with filename as last part of name
    localpath = joinpath(localdir, filename)
    Downloads.download(remotepath, localpath)
    return localpath
end

"""
    register_WOA(product_year, tracer, period, resolution)

Registers a `datadep` for the variable `tracer` averaged over `period` at resolution `resolution`.
"""
function register_WOA(product_year, tracer, period, resolution)
    register(DataDep(
        my_DataDeps_name(product_year, tracer, period, resolution),
        string(citation(product_year, tracer)),
        url_WOA(product_year, tracer, period, resolution),
        sha2_256,
        fetch_method = fallback_download
    ))
    return nothing
end

#============================================================
DataDeps registering name
============================================================#
my_DataDeps_name(product_year, tracer, period, resolution) = string(
    "WOA",
    my_product_year(product_year), "_",
    my_averaging_period(period), "_",
    WOA_path_varname(tracer), "_",
    surface_grid_size(resolution)
)

#============================================================
WOA product_year of the data product
============================================================#
my_product_year(product_year) = @match product_year begin
    2005 || 05 || "2005" || "05" => "05"
    2009 || 09 || "2009" || "09" => "09"
    2013 || 13 || "2013" || "13" => "13"
    2018 || 18 || "2018" || "18" => "18"
    _ => error("Cannot register WOA data from year $product_year")
end

#============================================================
Variable names
============================================================#
incorrect_varname(tracer) = """
"$tracer" is an incorrect variable name.
Use one of these `String`s:
    - "t" for Temperature
    - "s" for Salinity
    - "I" for Density
    - "C" for Conductivity
    - "o" for Dissolved Oxygen
    - "O" for Percent Oxygen Saturation
    - "A" for Apparent Oxygen Utilization
    - "i" for Silicate
    - "p" for Phosphate
    - "n" for Nitrate
"""
WOA_path_varname(tracer) = @match my_varname(tracer) begin
    "Temp"  => "temperature"
    "Salt"  => "salinity"
    "Dens"  => "density"
    "O2"    => "oxygen"
    "O2sat" => "o2sat"
    "AOU"   => "AOU"
    "DSi"   => "silicate"
    "DIP"   => "phosphate"
    "DIN"   => "nitrate"
    "Cond"  => "conductivity"
end
WOA_filename_varname(tracer) = @match my_varname(tracer) begin
    "Temp"  => "t"
    "Salt"  => "s"
    "Dens"  => "I"
    "O2"    => "o"
    "O2sat" => "O"
    "AOU"   => "A"
    "DSi"   => "i"
    "DIP"   => "p"
    "DIN"   => "n"
    "Cond"  => "C"
end
my_varname(tracer) = @match tracer begin
    "t" || "T" || "Temperature" || "temperature" || "Temp" || "temp"                   => "Temp"
    "s" || "Salinity" || "salinity" || "Salt" || "salt"                                => "Salt"
    "I" || "Density" || "density" || "Dens" || "dens" || "σ"                           => "Dens"
    "o" || "O2" || "O₂" || "Oxygen" || "oxygen"  || "Dissolved oxygen"                 => "O2"
    "O" || "o2sat" || "O₂sat" || "O2sat" || "O2Sat" || "oxygen saturation" || "Oxygen saturation" => "O2sat"
    "A" || "AOU" || "Apparent oxygen utilization"                                      => "AOU"
    "i" || "silicate" || "DSi" || "Silicic Acid" || "Si(OH)4" || "SiOH4" || "sioh4"    => "DSi"
    "p" || "phosphate" || "PO4" || "Phosphate" || "DIP" || "po4" || "PO₄"              => "DIP"
    "n" || "nitrate" || "NO3" || "Nitrate" || "DIN" || "no3" || "NO₃"                  => "DIN"
    "C" || "conductivity" || "Conductivity" || "Cond" || "cond"                        => "Cond"
    _ => error(incorrect_varname(tracer))
end
WOA_varname(tracer, field) = string(WOA_filename_varname(tracer), "_", my_field_type_code(field))

#============================================================
Averaging period
============================================================#
incorrect_averagingperiod(period) = """
"$period" is an incorrect averaging period.
The averaging period must fit the World Ocean Atlas naming convention.
That is, it must be one of these:
    - "00" for annual statistics, all data used
    - "13" to "16" for seasonal statistics:
        - "13" for winter (first three months of the year - Jan–Mar)
        - "14" for spring (Apr–Jun)
        - "15" for summer (Jul–Sep)
        - "16" for autumn (Oct–Dec)
    - "01" to "12" for monthly statistics (starting with "01" = January, to "12" = December)
"""
my_averaging_period(period::String) = @match lowercase(period) begin
    "00" ||  "0" ||    "annual"          =>    "Annual"
    "13" || "13" ||    "winter"          =>    "Winter"
    "14" || "14" ||    "spring"          =>    "Spring"
    "15" || "15" ||    "summer"          =>    "Summer"
    "16" || "16" ||    "autumn"          =>    "Autumn"
    "01" ||  "1" ||   "january" || "jan" =>   "January"
    "02" ||  "2" ||  "february" || "feb" =>  "February"
    "03" ||  "3" ||     "march" || "mar" =>     "March"
    "04" ||  "4" ||     "april" || "apr" =>     "April"
    "05" ||  "5" ||       "may" || "may" =>       "May"
    "06" ||  "6" ||      "june" || "jun" =>      "June"
    "07" ||  "7" ||      "july" || "jul" =>      "July"
    "08" ||  "8" ||    "august" || "aug" =>    "August"
    "09" ||  "9" || "september" || "sep" => "September"
    "10" || "10" ||   "october" || "oct" =>   "October"
    "11" || "11" ||  "november" || "nov" =>  "November"
    "12" || "12" ||  "december" || "dec" =>  "December"
    _ => error(incorrect_averagingperiod(period))
end
my_averaging_period(period::Int) = my_averaging_period(string(period))
WOA_averaging_period(period) = @match my_averaging_period(period) begin
    "Annual"    => "00"
    "Winter"    => "13"
    "Spring"    => "14"
    "Summer"    => "15"
    "Autumn"    => "16"
    "January"   => "01"
    "February"  => "02"
    "March"     => "03"
    "April"     => "04"
    "May"       => "05"
    "June"      => "06"
    "July"      => "07"
    "August"    => "08"
    "September" => "09"
    "October"   => "10"
    "November"  => "11"
    "December"  => "12"
end
seasonal_annual_monthly(period) = @match WOA_averaging_period(period) begin
    "00" => "annual"
    "13" || "14" || "15" || "16" => "seasonal"
    _ => "monthly"
end

#============================================================
Field type
============================================================#
incorrect_field_type_code(field) = """
"$field" is an incorrect "field type code".
The "field type code" must be must be one of these:
    - "an" for the objectively analyzed climatology
    - "mn" for the statistical mean
    - "sd" for the standard deviation
You need to edit the `my_function` function to add more!
"""
my_field_type_code(field) = @match lowercase(field) begin
    "an" || "objectively analyzed climatology"  => "an"
    "mn" || "mean" || "statistical mean"        => "mn"
    "sd" || "std"  || "standard deviation"      => "sd"
    "dd" || "number of observations"            => "dd"
    _ => error(incorrect_varfunc(field))
end

#============================================================
Resolution names
============================================================#
incorrect_resolution(resolution) = """
"$resolution" is an incorrect resolution.
Only these resolutions are available:
    - "1" for 1°×1° resolution
    - "5" for 5°×5° resolution
    - "0.25" for 0.25°×0.25° resolution
You need to edit the `myresolution` function to add more!
"""
WOA_path_resolution(resolution) = @match my_resolution(resolution) begin
    "0.25°" => "0.25"
    "1°"    => "1.00"
    "5°"    => "5deg"
end
WOA_filename_resolution(resolution) = @match my_resolution(resolution) begin
    "0.25°" => "04"
    "1°"    => "01"
    "5°"    => "5d"
end
surface_grid_size(resolution) = @match my_resolution(resolution) begin
    "0.25°" => "1440x720"
    "1°"    => "360x180"
    "5°"    => "73x36"
end
my_resolution(resolution) = @match resolution begin
    "0.25°×0.25°" || "04" || "0.25d" || "0.25" || "0.25°" => "0.25°"
       "1°×1°"    || "01" || "1d"    || "1"    || "1°"    => "1°"
       "5°×5°"    || "05" || "5d"    || "5"    || "5°"    => "5°"
    _ => error(incorrect_resolution(resolution))
end
WOA09_file_resolution(resolution) = @match my_resolution(resolution) begin
    "1°" => "1deg"
    "5°" => "5deg"
    _ => error("No such resolution for WOA09 data")
end

#============================================================
Decade names
============================================================#
WOA_decade(tracer) = @match my_varname(tracer) begin
    "Temp" || "Salt" || "Dens" || "Cond" => "decav"
    "O2" || "O2sat" || "AOU" || "DSi" || "DIP" || "DIN" => "all"
    _ => error(incorrect_varname(tracer))
end

WOA_v2(tracer) = @match my_varname(tracer) begin
    "Salt" => "v2"
    "Dens" || "Temp" || "Cond" || "O2" || "O2sat" || "AOU" || "DSi" || "DIP" || "DIN" => ""
    _ => error(incorrect_varname(tracer))
end

#============================================================
NetCDF file name
============================================================#
function WOA_NetCDF_filename(product_year, tracer, period, resolution)
    return string("woa",
                  my_product_year(product_year), "_",
                  WOA_decade(tracer), "_",
                  WOA_filename_varname(tracer),
                  WOA_averaging_period(period), "_",
                  WOA_filename_resolution(resolution),
                  WOA_v2(tracer), ".nc")
end

#============================================================
URLs
============================================================#
"""
    url_WOA(product_year, tracer, period, resolution)

Returns the URL (`String`) for the NetCDF file of the World Ocean Atlas.
This URL `String` typically looks like
    ```
    "https://data.nodc.noaa.gov/woa/WOA13/DATAv2/phosphate/netcdf/all/1.00/woa13_all_p00_01.nc"
    ```
"""
url_WOA(product_year, tracer, period, resolution) = string("https://data.nodc.noaa.gov/woa/WOA",
                                   my_product_year(product_year), "/",
                                   url_DATA(product_year), "/",
                                   WOA_path_varname(tracer), "/netcdf/",
                                   WOA_decade(tracer), "/",
                                   WOA_path_resolution(resolution), "/",
                                   WOA_NetCDF_filename(product_year, tracer, period, resolution))
url_WOA_THREDDS_18(tracer, period, resolution) = string("https://data.nodc.noaa.gov/thredds/dodsC/ncei/woa/",
                                        WOA_path_varname(tracer), "/",
                                        WOA_decade(tracer), "/",
                                        WOA_path_resolution(resolution), "/",
                                        WOA_NetCDF_filename(2018, tracer, period, resolution))
url_WOA_THREDDS_09(tracer, period, resolution) = string("https://data.nodc.noaa.gov/thredds/dodsC/woa09/",
                                        WOA_path_varname(tracer), "_",
                                        seasonal_annual_monthly(period), "_",
                                        WOA09_file_resolution(resolution), ".nc")
url_WOA_THREDDS(product_year, tracer, period, resolution) = @match my_product_year(product_year) begin
    "18" => url_WOA_THREDDS_18(tracer, period, resolution)
    "09" => url_WOA_THREDDS_09(tracer, period, resolution)
    _ => "No THREDDS links for year $product_year"
end
url_DATA(product_year) = @match my_product_year(product_year) begin
    "09" => "DATA"
    "13" => "DATAv2"
    "18" => "DATA"
end

"""
    varname(tracer)

Returns the World Ocean Data variable name that "matches" `tracer`.
"""
varname(tracer::String) = @match lowercase(tracer) begin
    "lat" || "latitude"   => "lat"
    "lon" || "longitude"  => "lon"
    "depth" || "depths"   => "depth"
end


text:
# *****************************************************************************
# Written by Ritchie Lee, ritchie.lee@sv.cmu.edu
# *****************************************************************************
# Copyright ã 2015, United States Government, as represented by the
# Administrator of the National Aeronautics and Space Administration. All
# rights reserved.  The Reinforcement Learning Encounter Simulator (RLES)
# platform is 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.
# _____________________________________________________________________________
# Reinforcement Learning Encounter Simulator (RLES) includes the following
# third party software. The SISLES.jl package is licensed under the MIT Expat
# License: Copyright (c) 2014: Youngjun Kim.
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED
# "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
# NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR
# PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
# HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN
# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
# *****************************************************************************

import Compat.ASCIIString
using DerivationTrees

#note: tds columns are not used, but left in the data
const COLNAMES = [
  "RA_1",
  "vert_rate_1",
  "alt_diff_1",
  "psi_1",
  "intr_sr_1",
  "intr_chi_1",
  "intr_vrc0_1",
  "intr_vrc1_1",
  "intr_vrc2_1",
  "cc0_1",
  "cc1_1",
  "cc2_1",
  "vc0_1",
  "vc1_1",
  "vc2_1",
  "ua0_1",
  "ua1_1",
  "ua2_1",
  "da0_1",
  "da1_1",
  "da2_1",
  "target_rate_1",
  "crossing_1",
  "alarm_1",
  "alert_1",
  "intr_out_vrc0_1",
  "intr_out_vrc1_1",
  "intr_out_vrc2_1",
  "intr_out_tds_1",
  "response_none_1",
  "response_stay_1",
  "response_follow_1",
  "response_timer_1",
  "response_h_d_1",
  "response_psi_d_1",
  "v_1",
  "h_1",
  "RA_2",
  "vert_rate_2",
  "alt_diff_2",
  "psi_2",
  "intr_sr_2",
  "intr_chi_2",
  "intr_vrc0_2",
  "intr_vrc1_2",
  "intr_vrc2_2",
  "cc0_2",
  "cc1_2",
  "cc2_2",
  "vc0_2",
  "vc1_2",
  "vc2_2",
  "ua0_2",
  "ua1_2",
  "ua2_2",
  "da0_2",
  "da1_2",
  "da2_2",
  "target_rate_2",
  "crossing_2",
  "alarm_2",
  "alert_2",
  "intr_out_vrc0_2",
  "intr_out_vrc1_2",
  "intr_out_vrc2_2",
  "intr_out_tds_2",
  "response_none_2",
  "response_stay_2",
  "response_follow_2",
  "response_timer_2",
  "response_h_d_2",
  "response_psi_d_2",
  "v_2",
  "h_2",
  "converging",
  "abs_alt_diff",
  "horizontal_range"]

function get_map(num_aircraft::Int64 = 2)
  D = Dict{ASCIIString,ASCIIString}()
  for i = 1:num_aircraft
    D["RA_$i"] = "RA is active on aircraft $i"
    D["vert_rate_$i"] = "aircraft $i's vertical rate"
    D["alt_diff_$i"] = "altitude difference relative to aircraft $i"
    D["psi_$i"] = "aircraft $i's heading angle"
    D["intr_sr_$i"] = "aircraft $i's intruder slant range"
    D["intr_chi_$i"] = "aircraft $i's intruder bearing"
    D["intr_vrc0_$i"] = "aircraft $i's input vrc bit 0 is set"
    D["intr_vrc1_$i"] = "aircraft $i's input vrc bit 1 is set"
    D["intr_vrc2_$i"] = "aircraft $i's input vrc bit 2 is set"
    D["cc0_$i"] = "aircraft $i's cc bit 0 is set"
    D["cc1_$i"] = "aircraft $i's cc bit 1 is set"
    D["cc2_$i"] = "aircraft $i's cc bit 2 is set"
    D["vc0_$i"] = "aircraft $i's vc bit 0 is set"
    D["vc1_$i"] = "aircraft $i's vc bit 1 is set"
    D["vc2_$i"] = "aircraft $i's vc bit 2 is set"
    D["ua0_$i"] = "aircraft $i's ua bit 0 is set"
    D["ua1_$i"] = "aircraft $i's ua bit 1 is set"
    D["ua2_$i"] = "aircraft $i's ua bit 2 is set"
    D["da0_$i"] = "aircraft $i's da bit 0 is set"
    D["da1_$i"] = "aircraft $i's da bit 1 is set"
    D["da2_$i"] = "aircraft $i's da bit 2 is set"
    D["target_rate_$i"] = "aircraft $i's RA target rate"
    D["crossing_$i"] = "crossing RA is issued to aircraft $i"
    D["alarm_$i"] = "RA alarm occurs on aircraft $i"
    D["alert_$i"] = "RA alert occurs aircraft $i"
    D["intr_out_vrc0_$i"] = "aircraft $i's output vrc bit 0 is set"
    D["intr_out_vrc1_$i"] = "aircraft $i's output vrc bit 1 is set"
    D["intr_out_vrc2_$i"] = "aircraft $i's output vrc bit 2 is set"
    D["intr_out_tds_$i"] = "aircraft $i's output tds"
    D["response_none_$i"] = "pilot $i is flying intended trajectory"
    D["response_stay_$i"] = "pilot $i is responding to previous RA"
    D["response_follow_$i"] = "pilot $i is responding to current RA"
    D["response_timer_$i"] = "number of seconds remaining in pilot $i's response delay"
    D["response_h_d_$i"] = "pilot $i's commanded vertical rate"
    D["response_psi_d_$i"] = "pilot $i's commanded turn rate"
    D["v_$i"] = "aircraft $i's velocity"
    D["h_$i"] = "aircraft $i's altitude"
  end
  D["converging"] = "aircraft are converging"
  D["abs_alt_diff"] = "absolute altitude difference"
  D["horizontal_range"] = "horizontal range"
  return D
end

const U_LTE = "<=" #"\u2264"
const U_GTE = ">=" #"\u2265"
const U_IMPLIES = "=>"

function get_format_pretty{T<:AbstractString}(colnames::Vector{T})
  fmt = Format()

  fmt["always"] = (cmd, args) -> "G($(args[1]))"
  fmt["eventually"] = (cmd, args) -> "F($(args[1]))"

  bin_infix(cmd, args, insym) = "($(args[1]) $insym $(args[2]))"
  fmt["and"] = (cmd, args) -> bin_infix(cmd, args, "&")
  fmt["or"] = (cmd, args) -> bin_infix(cmd, args, "|")
  fmt["not"] = (cmd, args) -> "!($(args[1]))"
  fmt["implies"] = (cmd, args) -> "$(args[1]) => $(args[2])"

  count_op(cmd, args, op) = "(count($(args[1])) $op $(args[2]))"
  fmt["count.1.1.1.1"] = (cmd, args) -> count_op(cmd, args, "<")
  fmt["count.1.1.1.2"] = (cmd, args) -> count_op(cmd, args, "<=")
  fmt["count.1.1.2"] = (cmd, args) -> count_op(cmd, args, ">")
  fmt["count.1.2"] = (cmd, args) -> count_op(cmd, args, ">=")
  fmt["count.2"] = (cmd, args) -> count_op(cmd, args, "==")

  bin_infix_eq(cmd, args) = bin_infix(cmd, args, ".==")
  bin_infix_lt(cmd, args) = bin_infix(cmd, args, ".<")
  bin_infix_lte(cmd, args) = bin_infix(cmd, args, ".<=")
  fmt["vrate_eq.1"] = bin_infix_eq
  fmt["vrate_eq.2"] = bin_infix_eq
  fmt["altdiff_eq"] = bin_infix_eq
  fmt["chi_angle_eq.1"] = bin_infix_eq
  fmt["chi_angle_eq.2"] = bin_infix_eq
  fmt["psi_angle_eq"] = bin_infix_eq
  fmt["sr_eq"] = bin_infix_eq
  fmt["timer_eq.1"] = bin_infix_eq
  fmt["timer_eq.2"] = bin_infix_eq
  fmt["psid_eq.1"] = bin_infix_eq
  fmt["psid_eq.2"] = bin_infix_eq
  fmt["v_eq.1"] = bin_infix_eq
  fmt["v_eq.2"] = bin_infix_eq
  fmt["alt_eq.1"] = bin_infix_eq
  fmt["alt_eq.2"] = bin_infix_eq
  fmt["abs_altdiff_eq"] = bin_infix_eq
  fmt["vrate_lt.1"] = bin_infix_lt
  fmt["vrate_lt.2"] = bin_infix_lt
  fmt["altdiff_lt"] = bin_infix_lt
  fmt["chi_angle_lt.1"] = bin_infix_lt
  fmt["chi_angle_lt.2"] = bin_infix_lt
  fmt["psi_angle_lt"] = bin_infix_lt
  fmt["sr_lt"] = bin_infix_lt
  fmt["timer_lt.1"] = bin_infix_lt
  fmt["timer_lt.2"] = bin_infix_lt
  fmt["psid_lt.1"] = bin_infix_lt
  fmt["psid_lt.2"] = bin_infix_lt
  fmt["v_lt.1"] = bin_infix_lt
  fmt["v_lt.2"] = bin_infix_lt
  fmt["alt_lt.1"] = bin_infix_lt
  fmt["alt_lt.2"] = bin_infix_lt
  fmt["abs_altdiff_lt"] = bin_infix_lt
  fmt["vrate_lte.1"] = bin_infix_lte
  fmt["vrate_lte.2"] = bin_infix_lte
  fmt["altdiff_lte"] = bin_infix_lte
  fmt["chi_angle_lte.1"] = bin_infix_lte
  fmt["chi_angle_lte.2"] = bin_infix_lte
  fmt["psi_angle_lte"] = bin_infix_lte
  fmt["sr_lte"] = bin_infix_lte
  fmt["timer_lte.1"] = bin_infix_lte
  fmt["timer_lte.2"] = bin_infix_lte
  fmt["psid_lte.1"] = bin_infix_lte
  fmt["psid_lte.2"] = bin_infix_lte
  fmt["v_lte.1"] = bin_infix_lte
  fmt["v_lte.2"] = bin_infix_lte
  fmt["alt_lte.1"] = bin_infix_lte
  fmt["alt_lte.2"] = bin_infix_lte
  fmt["abs_altdiff_lte"] = bin_infix_lte

  abs_bin_infix_op(cmd, args, op) = "(abs($(args[1])) $op $(args[2]))"
  abs_bin_infix_eq(cmd, args) = abs_bin_infix_op(cmd, args, ".==")
  abs_bin_infix_lt(cmd, args) = abs_bin_infix_op(cmd, args, ".<")
  abs_bin_infix_lte(cmd, args) = abs_bin_infix_op(cmd, args, ".<=")
  fmt["vrate_abseq.1"] = abs_bin_infix_eq
  fmt["vrate_abseq.2"] = abs_bin_infix_eq
  fmt["altdiff_abseq"] = abs_bin_infix_eq
  fmt["chi_angle_abseq.1"] = abs_bin_infix_eq
  fmt["chi_angle_abseq.2"] = abs_bin_infix_eq
  fmt["psi_angle_abseq"] = abs_bin_infix_eq
  fmt["sr_abseq"] = abs_bin_infix_eq
  fmt["timer_abseq.1"] = abs_bin_infix_eq
  fmt["timer_abseq.2"] = abs_bin_infix_eq
  fmt["psid_abseq.1"] = abs_bin_infix_eq
  fmt["psid_abseq.2"] = abs_bin_infix_eq
  fmt["v_abseq.1"] = abs_bin_infix_eq
  fmt["v_abseq.2"] = abs_bin_infix_eq
  fmt["alt_abseq.1"] = abs_bin_infix_eq
  fmt["alt_abseq.2"] = abs_bin_infix_eq
  fmt["vrate_abslt.1"] = abs_bin_infix_lt
  fmt["vrate_abslt.2"] = abs_bin_infix_lt
  fmt["altdiff_abslt"] = abs_bin_infix_lt
  fmt["chi_angle_abslt.1"] = abs_bin_infix_lt
  fmt["chi_angle_abslt.2"] = abs_bin_infix_lt
  fmt["psi_angle_abslt"] = abs_bin_infix_lt
  fmt["sr_abslt"] = abs_bin_infix_lt
  fmt["timer_abslt.1"] = abs_bin_infix_lt
  fmt["timer_abslt.2"] = abs_bin_infix_lt
  fmt["psid_abslt.1"] = abs_bin_infix_lt
  fmt["psid_abslt.2"] = abs_bin_infix_lt
  fmt["v_abslt.1"] = abs_bin_infix_lt
  fmt["v_abslt.2"] = abs_bin_infix_lt
  fmt["alt_abslt.1"] = abs_bin_infix_lt
  fmt["alt_abslt.2"] = abs_bin_infix_lt
  fmt["vrate_abslte.1"] = abs_bin_infix_lte
  fmt["vrate_abslte.2"] = abs_bin_infix_lte
  fmt["altdiff_abslte"] = abs_bin_infix_lte
  fmt["chi_angle_abslte.1"] = abs_bin_infix_lte
  fmt["chi_angle_abslte.2"] = abs_bin_infix_lte
  fmt["psi_angle_abslte"] = abs_bin_infix_lte
  fmt["sr_abslte"] = abs_bin_infix_lte
  fmt["timer_abslte.1"] = abs_bin_infix_lte
  fmt["timer_abslte.2"] = abs_bin_infix_lte
  fmt["psid_abslte.1"] = abs_bin_infix_lte
  fmt["psid_abslte.2"] = abs_bin_infix_lte
  fmt["v_abslte.1"] = abs_bin_infix_lte
  fmt["v_abslte.2"] = abs_bin_infix_lte
  fmt["alt_abslte.1"] = abs_bin_infix_lte
  fmt["alt_abslte.2"] = abs_bin_infix_lte

  diff_infix_op(cmd, args, op) = "($(args[1]) - $(args[2]) $op $(args[3]))"
  diff_infix_eq(cmd, args) = diff_infix_op(cmd, args, ".==")
  diff_infix_lt(cmd, args) = diff_infix_op(cmd, args, ".<")
  diff_infix_lte(cmd, args) = diff_infix_op(cmd, args, ".<=")
  fmt["vrate_diff_eq"] = diff_infix_eq
  fmt["chi_angle_diff_eq"] = diff_infix_eq
  fmt["psi_angle_diff_eq"] = diff_infix_eq
  fmt["timer_diff_eq"] = diff_infix_eq
  fmt["psid_diff_eq"] = diff_infix_eq
  fmt["v_diff_eq"] = diff_infix_eq
  fmt["vrate_diff_lt"] = diff_infix_lt
  fmt["chi_angle_diff_lt"] = diff_infix_lt
  fmt["psi_angle_diff_lt"] = diff_infix_lt
  fmt["timer_diff_lt"] = diff_infix_lt
  fmt["psid_diff_lt"] = diff_infix_lt
  fmt["v_diff_lt"] = diff_infix_lt
  fmt["vrate_diff_lte"] = diff_infix_lte
  fmt["chi_angle_diff_lte"] = diff_infix_lte
  fmt["psi_angle_diff_lte"] = diff_infix_lte
  fmt["timer_diff_lte"] = diff_infix_lte
  fmt["psid_diff_lte"] = diff_infix_lte
  fmt["v_diff_lte"] = diff_infix_lte

  abs_diff_infix_op(cmd, args, op) = "(abs($(args[1]) - $(args[2])) $op $(args[3]))"
  abs_diff_infix_eq(cmd, args) = abs_diff_infix_op(cmd, args, ".==")
  abs_diff_infix_lt(cmd, args) = abs_diff_infix_op(cmd, args, ".<")
  abs_diff_infix_lte(cmd, args) = abs_diff_infix_op(cmd, args, ".<=")
  fmt["vrate_absdiff_eq"] = abs_diff_infix_eq
  fmt["chi_angle_absdiff_eq"] = abs_diff_infix_eq
  fmt["psi_angle_absdiff_eq"] = abs_diff_infix_eq
  fmt["timer_absdiff_eq"] = abs_diff_infix_eq
  fmt["psid_absdiff_eq"] = abs_diff_infix_eq
  fmt["v_absdiff_eq"] = abs_diff_infix_eq
  fmt["vrate_absdiff_lt"] = abs_diff_infix_lt
  fmt["chi_angle_absdiff_lt"] = abs_diff_infix_lt
  fmt["psi_angle_absdiff_lt"] = abs_diff_infix_lt
  fmt["timer_absdiff_lt"] = abs_diff_infix_lt
  fmt["psid_absdiff_lt"] = abs_diff_infix_lt
  fmt["v_absdiff_lt"] = abs_diff_infix_lt
  fmt["vrate_absdiff_lte"] = abs_diff_infix_lte
  fmt["chi_angle_absdiff_lte"] = abs_diff_infix_lte
  fmt["psi_angle_absdiff_lte"] = abs_diff_infix_lte
  fmt["timer_absdiff_lte"] = abs_diff_infix_lte
  fmt["psid_absdiff_lte"] = abs_diff_infix_lte
  fmt["v_absdiff_lte"] = abs_diff_infix_lte

  feat(cmd, args) = "$(colnames[parse(Int, args[1])])"
  fmt["bin_feat"] = feat
  fmt["vrate_feat"] = feat
  fmt["altdiff_feat"] = feat
  fmt["abs_altdiff_feat"] = feat
  fmt["angle_feat"] = feat
  fmt["psi_angle_feat"] = feat
  fmt["chi_angle_feat"] = feat
  fmt["sr_feat"] = feat
  fmt["timer_feat"] = feat
  fmt["psid_feat"] = feat
  fmt["v_feat"] = feat
  fmt["alt_feat"] = feat

  valunit(cmd, args, unit) = "$(args[1]) $unit"
  fmt["vrate_val"] = (cmd, args) -> valunit(cmd, args, "ft/s")
  fmt["altdiff_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["abs_altdiff_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["angle_val"] = (cmd, args) -> valunit(cmd, args, "deg")
  fmt["sr_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["timer_val"] = (cmd, args) -> valunit(cmd, args, "s")
  fmt["psid_val"] = (cmd, args) -> valunit(cmd, args, "deg/s")
  fmt["v_val"] = (cmd, args) -> valunit(cmd, args, "ft/s")
  fmt["alt_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["timestep"] = (cmd, args) -> valunit(cmd, args, "s")

  bin_sign(cmd, args) = "(sign($(args[1])) .== sign($(args[2])))"
  fmt["vrate_sign"] = bin_sign
  fmt["chi_angle_sign"] = bin_sign
  fmt["psid_sign"] = bin_sign

  return fmt
end

function get_format_natural{T<:AbstractString}(colnames::Vector{T})
  fmt = Format()

  fmt["always"] = (cmd, args) -> "for all time, $(args[1])"
  fmt["eventually"] = (cmd, args) -> "at some point, $(args[1])"

  bin_infix(cmd, args, insym) = "[$(args[1]) $insym $(args[2])]"
  fmt["and"] = (cmd, args) -> bin_infix(cmd, args, "and")
  fmt["or"] = (cmd, args) -> bin_infix(cmd, args, "or")
  fmt["not"] = (cmd, args) -> "[it is not true that $(args[1])]"
  fmt["implies"] = (cmd, args) -> "whenever $(args[1]), it is also true that $(args[2])"

  count_op(cmd, args, op) = "[the number of times $(args[1]) $op $(args[2])]"
  fmt["count.1.1.1.1"] = (cmd, args) -> count_op(cmd, args, "is less than")
  fmt["count.1.1.1.2"] = (cmd, args) -> count_op(cmd, args, "is less than or equal to")
  fmt["count.1.1.2"] = (cmd, args) -> count_op(cmd, args, "is greater than")
  fmt["count.1.2"] = (cmd, args) -> count_op(cmd, args, "is greater than or equal to")
  fmt["count.2"] = (cmd, args) -> count_op(cmd, args, "equals")

  bin_infix_eq(cmd, args) = bin_infix(cmd, args, "equals")
  bin_infix_lt(cmd, args) = bin_infix(cmd, args, "is less than")
  bin_infix_lte(cmd, args) = bin_infix(cmd, args, "is less than or equal to")
  fmt["vrate_eq.1"] = bin_infix_eq
  fmt["vrate_eq.2"] = bin_infix_eq
  fmt["altdiff_eq"] = bin_infix_eq
  fmt["chi_angle_eq.1"] = bin_infix_eq
  fmt["chi_angle_eq.2"] = bin_infix_eq
  fmt["psi_angle_eq"] = bin_infix_eq
  fmt["sr_eq"] = bin_infix_eq
  fmt["timer_eq.1"] = bin_infix_eq
  fmt["timer_eq.2"] = bin_infix_eq
  fmt["psid_eq.1"] = bin_infix_eq
  fmt["psid_eq.2"] = bin_infix_eq
  fmt["v_eq.1"] = bin_infix_eq
  fmt["v_eq.2"] = bin_infix_eq
  fmt["alt_eq.1"] = bin_infix_eq
  fmt["alt_eq.2"] = bin_infix_eq
  fmt["abs_altdiff_eq"] = bin_infix_eq
  fmt["vrate_lt.1"] = bin_infix_lt
  fmt["vrate_lt.2"] = bin_infix_lt
  fmt["altdiff_lt"] = bin_infix_lt
  fmt["chi_angle_lt.1"] = bin_infix_lt
  fmt["chi_angle_lt.2"] = bin_infix_lt
  fmt["psi_angle_lt"] = bin_infix_lt
  fmt["sr_lt"] = bin_infix_lt
  fmt["timer_lt.1"] = bin_infix_lt
  fmt["timer_lt.2"] = bin_infix_lt
  fmt["psid_lt.1"] = bin_infix_lt
  fmt["psid_lt.2"] = bin_infix_lt
  fmt["v_lt.1"] = bin_infix_lt
  fmt["v_lt.2"] = bin_infix_lt
  fmt["alt_lt.1"] = bin_infix_lt
  fmt["alt_lt.2"] = bin_infix_lt
  fmt["abs_altdiff_lt"] = bin_infix_lt
  fmt["vrate_lte.1"] = bin_infix_lte
  fmt["vrate_lte.2"] = bin_infix_lte
  fmt["altdiff_lte"] = bin_infix_lte
  fmt["chi_angle_lte.1"] = bin_infix_lte
  fmt["chi_angle_lte.2"] = bin_infix_lte
  fmt["psi_angle_lte"] = bin_infix_lte
  fmt["sr_lte"] = bin_infix_lte
  fmt["timer_lte.1"] = bin_infix_lte
  fmt["timer_lte.2"] = bin_infix_lte
  fmt["psid_lte.1"] = bin_infix_lte
  fmt["psid_lte.2"] = bin_infix_lte
  fmt["v_lte.1"] = bin_infix_lte
  fmt["v_lte.2"] = bin_infix_lte
  fmt["alt_lte.1"] = bin_infix_lte
  fmt["alt_lte.2"] = bin_infix_lte
  fmt["abs_altdiff_lte"] = bin_infix_lte

  abs_bin_infix_op(cmd, args, op) = "[the absolute value of $(args[1]) $op $(args[2]))]"
  abs_bin_infix_eq(cmd, args) = abs_bin_infix_op(cmd, args, "equals")
  abs_bin_infix_lt(cmd, args) = abs_bin_infix_op(cmd, args, "is less than")
  abs_bin_infix_lte(cmd, args) = abs_bin_infix_op(cmd, args, "is less than or equal to")
  fmt["vrate_abseq.1"] = abs_bin_infix_eq
  fmt["vrate_abseq.2"] = abs_bin_infix_eq
  fmt["altdiff_abseq"] = abs_bin_infix_eq
  fmt["chi_angle_abseq.1"] = abs_bin_infix_eq
  fmt["chi_angle_abseq.2"] = abs_bin_infix_eq
  fmt["psi_angle_abseq"] = abs_bin_infix_eq
  fmt["sr_abseq"] = abs_bin_infix_eq
  fmt["timer_abseq.1"] = abs_bin_infix_eq
  fmt["timer_abseq.2"] = abs_bin_infix_eq
  fmt["psid_abseq.1"] = abs_bin_infix_eq
  fmt["psid_abseq.2"] = abs_bin_infix_eq
  fmt["v_abseq.1"] = abs_bin_infix_eq
  fmt["v_abseq.2"] = abs_bin_infix_eq
  fmt["alt_abseq.1"] = abs_bin_infix_eq
  fmt["alt_abseq.2"] = abs_bin_infix_eq
  fmt["vrate_abslt.1"] = abs_bin_infix_lt
  fmt["vrate_abslt.2"] = abs_bin_infix_lt
  fmt["altdiff_abslt"] = abs_bin_infix_lt
  fmt["chi_angle_abslt.1"] = abs_bin_infix_lt
  fmt["chi_angle_abslt.2"] = abs_bin_infix_lt
  fmt["psi_angle_abslt"] = abs_bin_infix_lt
  fmt["sr_abslt"] = abs_bin_infix_lt
  fmt["timer_abslt.1"] = abs_bin_infix_lt
  fmt["timer_abslt.2"] = abs_bin_infix_lt
  fmt["psid_abslt.1"] = abs_bin_infix_lt
  fmt["psid_abslt.2"] = abs_bin_infix_lt
  fmt["v_abslt.1"] = abs_bin_infix_lt
  fmt["v_abslt.2"] = abs_bin_infix_lt
  fmt["alt_abslt.1"] = abs_bin_infix_lt
  fmt["alt_abslt.2"] = abs_bin_infix_lt
  fmt["vrate_abslte.1"] = abs_bin_infix_lte
  fmt["vrate_abslte.2"] = abs_bin_infix_lte
  fmt["altdiff_abslte"] = abs_bin_infix_lte
  fmt["chi_angle_abslte.1"] = abs_bin_infix_lte
  fmt["chi_angle_abslte.2"] = abs_bin_infix_lte
  fmt["psi_angle_abslte"] = abs_bin_infix_lte
  fmt["sr_abslte"] = abs_bin_infix_lte
  fmt["timer_abslte.1"] = abs_bin_infix_lte
  fmt["timer_abslte.2"] = abs_bin_infix_lte
  fmt["psid_abslte.1"] = abs_bin_infix_lte
  fmt["psid_abslte.2"] = abs_bin_infix_lte
  fmt["v_abslte.1"] = abs_bin_infix_lte
  fmt["v_abslte.2"] = abs_bin_infix_lte
  fmt["alt_abslte.1"] = abs_bin_infix_lte
  fmt["alt_abslte.2"] = abs_bin_infix_lte

  diff_infix_op(cmd, args, op) = "[the difference between $(args[1]) and $(args[2]) $op $(args[3])]"
  diff_infix_eq(cmd, args) = diff_infix_op(cmd, args, "equals")
  diff_infix_lt(cmd, args) = diff_infix_op(cmd, args, "is less than")
  diff_infix_lte(cmd, args) = diff_infix_op(cmd, args, "is less than or equal to")
  fmt["vrate_diff_eq"] = diff_infix_eq
  fmt["chi_angle_diff_eq"] = diff_infix_eq
  fmt["psi_angle_diff_eq"] = diff_infix_eq
  fmt["timer_diff_eq"] = diff_infix_eq
  fmt["psid_diff_eq"] = diff_infix_eq
  fmt["v_diff_eq"] = diff_infix_eq
  fmt["vrate_diff_lt"] = diff_infix_lt
  fmt["chi_angle_diff_lt"] = diff_infix_lt
  fmt["psi_angle_diff_lt"] = diff_infix_lt
  fmt["timer_diff_lt"] = diff_infix_lt
  fmt["psid_diff_lt"] = diff_infix_lt
  fmt["v_diff_lt"] = diff_infix_lt
  fmt["vrate_diff_lte"] = diff_infix_lte
  fmt["chi_angle_diff_lte"] = diff_infix_lte
  fmt["psi_angle_diff_lte"] = diff_infix_lte
  fmt["timer_diff_lte"] = diff_infix_lte
  fmt["psid_diff_lte"] = diff_infix_lte
  fmt["v_diff_lte"] = diff_infix_lte

  abs_diff_infix_op(cmd, args, op) = "[the absolute difference between $(args[1]) and $(args[2]) $op $(args[3])]"
  abs_diff_infix_eq(cmd, args) = abs_diff_infix_op(cmd, args, "equals")
  abs_diff_infix_lt(cmd, args) = abs_diff_infix_op(cmd, args, "is less than")
  abs_diff_infix_lte(cmd, args) = abs_diff_infix_op(cmd, args, "is less than or equal to")
  fmt["vrate_absdiff_eq"] = abs_diff_infix_eq
  fmt["chi_angle_absdiff_eq"] = abs_diff_infix_eq
  fmt["psi_angle_absdiff_eq"] = abs_diff_infix_eq
  fmt["timer_absdiff_eq"] = abs_diff_infix_eq
  fmt["psid_absdiff_eq"] = abs_diff_infix_eq
  fmt["v_absdiff_eq"] = abs_diff_infix_eq
  fmt["vrate_absdiff_lt"] = abs_diff_infix_lt
  fmt["chi_angle_absdiff_lt"] = abs_diff_infix_lt
  fmt["psi_angle_absdiff_lt"] = abs_diff_infix_lt
  fmt["timer_absdiff_lt"] = abs_diff_infix_lt
  fmt["psid_absdiff_lt"] = abs_diff_infix_lt
  fmt["v_absdiff_lt"] = abs_diff_infix_lt
  fmt["vrate_absdiff_lte"] = abs_diff_infix_lte
  fmt["chi_angle_absdiff_lte"] = abs_diff_infix_lte
  fmt["psi_angle_absdiff_lte"] = abs_diff_infix_lte
  fmt["timer_absdiff_lte"] = abs_diff_infix_lte
  fmt["psid_absdiff_lte"] = abs_diff_infix_lte
  fmt["v_absdiff_lte"] = abs_diff_infix_lte

  feat(cmd, args) = "[$(colnames[parse(Int, args[1])])]"
  fmt["bin_feat"] = feat
  fmt["vrate_feat"] = feat
  fmt["altdiff_feat"] = feat
  fmt["abs_altdiff_feat"] = feat
  fmt["angle_feat"] = feat
  fmt["psi_angle_feat"] = feat
  fmt["chi_angle_feat"] = feat
  fmt["sr_feat"] = feat
  fmt["timer_feat"] = feat
  fmt["psid_feat"] = feat
  fmt["v_feat"] = feat
  fmt["alt_feat"] = feat

  valunit(cmd, args, unit) = "$(args[1]) $unit"
  fmt["vrate_val"] = (cmd, args) -> valunit(cmd, args, "ft/s")
  fmt["altdiff_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["abs_altdiff_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["angle_val"] = (cmd, args) -> valunit(cmd, args, "deg")
  fmt["sr_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["timer_val"] = (cmd, args) -> valunit(cmd, args, "s")
  fmt["psid_val"] = (cmd, args) -> valunit(cmd, args, "deg/s")
  fmt["v_val"] = (cmd, args) -> valunit(cmd, args, "ft/s")
  fmt["alt_val"] = (cmd, args) -> valunit(cmd, args, "ft")
  fmt["timestep"] = (cmd, args) -> valunit(cmd, args, "s")

  bin_sign(cmd, args) = "[the sign of $(args[1]) is equal to the sign of $(args[2])]"
  fmt["vrate_sign"] = bin_sign
  fmt["chi_angle_sign"] = bin_sign
  fmt["psid_sign"] = bin_sign

  return fmt
end

const DESCRIP_MAP = get_map()
const COLNAMES_FULL = map(x -> DESCRIP_MAP[x], COLNAMES)
const FMT_PRETTY = get_format_pretty(COLNAMES)
const FMT_NATURAL = get_format_natural(COLNAMES_FULL)


text:
using Test

using PotentialFlow
using LinearAlgebra: norm, ⋅

@test isempty(detect_ambiguities(Vortex))

#=
# Copied verbatim from ssfrr/TestSetExtensions.jl
# The package is not 0.7/1.0 compatible yet, but
# this function is too useful
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(f->endswith(f, ".jl") && f!= basename(rootfile), tests)
        else
            tests = map(f->string(f, ".jl"), tests)
        end
        println();
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

@testset "All tests" begin
    @includetests ARGS
end
=#

const GROUP = get(ENV, "GROUP", "All")

if GROUP == "All" || GROUP == "Autodiff"
  include("forwarddiff.jl")
end

if GROUP == "All" || GROUP == "Bodies"
  include("Bodies.jl")
  include("utils/circle_plane.jl")
end

if GROUP == "All" || GROUP == "Plates"
  include("chebyshev.jl")
  include("Plates.jl")
end

if GROUP == "All" || GROUP == "Sheets"
  include("Sheets.jl")
end

if GROUP == "All" || GROUP == "TimeMarching"
  include("TimeMarching.jl")
end

if GROUP == "All" || GROUP == "Dispatch"
  include("dispatch.jl")
end

if GROUP == "All" || GROUP == "Utils"
  include("Utils.jl")
end

if GROUP == "All" || GROUP == "Notebooks"
  include("notebooks.jl")
end


text:
using PGFPlotsX, SpectralKit, Colors

function plot_family(family, ns, order, title;
                     xmin = domain_extrema(family)[1], xmax = domain_extrema(family)[2],
                     rel_offsets = [0.01, 0.01])
    axis = @pgf Axis({ xmin = xmin, xmax = xmax, xlabel = "y", title = title,
                       legend_pos = "outer north east"})
    xgrid = range(xmin, xmax; length = 400)
    mi, ma = domain_extrema(family)
    colors = distinguishable_colors(length(ns); lchoices = 0:60)
    offsets = [1, -1] .* rel_offsets .* (xmax - xmin)
    for (n, col) in zip(ns, colors)
        y = evaluate.(family, n, xgrid, order)
        @pgf push!(axis,
                   Plot({ no_marks, thick, color = col }, Table(xgrid, y)),
                   LegendEntry(string(n)))
    end
    ex = [mi, ma]
    for (i, (n, col)) in enumerate(zip(ns, colors))
        y = evaluate.(family, n, ex, order)
        @pgf push!(axis, Plot({ only_marks, color = col },
                              Table([xmin, xmax] .+ offsets .* i, y)))
    end
    axis
end

###
### Chebyshev
###

p = plot_family(Chebyshev(), 0:5, Val(0), "Chebyshev polynomials (value)")
pgfsave("chebyshev.png", p; dpi = 600)

p = plot_family(Chebyshev(), 0:5, Val(1), "Chebyshev polynomials (derivative)";
                rel_offsets = [0, 0])
pgfsave("chebyshev_deriv.png", p; dpi = 600)

###
### Rational Chebysev [0,∞)
###

F = ChebyshevSemiInf(0.0, 1.0)
p = plot_family(F, 0:5, Val(0), raw"Rational Chebyshev functions $[0,\infty)$ (value)";
                xmax = 1)
pgfsave("semiinf.png", p; dpi = 600)

p = plot_family(F, 0:5, Val(0), raw"Rational Chebyshev functions $[0,\infty)$ (value)";
                xmax = 200)
pgfsave("semiinf_birdseye.png", p; dpi = 600)

p = plot_family(F, 0:5, Val(1), raw"Rational Chebyshev functions $[0,\infty)$ (derivative)";
                xmax = 1, rel_offsets = [0, 0.01])
pgfsave("semiinf_deriv.png", p; dpi = 600)

###
### Rational Chebysev (-∞,∞)
###

F = ChebyshevInf(0.0, 1.0)
p = plot_family(F, 0:5, Val(0),
                raw"Rational Chebyshev functions $(-\infty,\infty)$ (value)";
                xmin = -1, xmax = 1)
pgfsave("inf.png", p; dpi = 600)

p = plot_family(F, 0:5, Val(0),
                raw"Rational Chebyshev functions $(-\infty,\infty)$ (value)";
                xmin = -30, xmax = 30)
pgfsave("inf_birdseye.png", p; dpi = 600)

p = plot_family(F, 0:5, Val(1),
                raw"Rational Chebyshev functions $(-\infty,\infty)$ (derivative)";
                xmin = -5, xmax = 5, rel_offsets = [0.01, 0.01])
pgfsave("inf_deriv.png", p; dpi = 600)


text:
module MonteCarlo

using LinearAlgebra: AbstractMatrix
using Reexport
# Loading the RNG will fail if Random is nto exported
@reexport using MonteCarloObservable, Random
import MonteCarloObservable.AbstractObservable
using Parameters, Requires
using TimerOutputs, LoopVectorization, StructArrays
using Printf, SparseArrays, LinearAlgebra, Dates, Statistics, Random, Distributed

import JLD, JLD2
using CodecZlib

# Because we fully load all data directly for JLD we lose access to the path
# This is supposed to keep track of path information
struct FileWrapper{T}
    file::T
    path::String
end

Base.getindex(fw::FileWrapper, k) = getindex(fw.file, k)
Base.setindex!(fw::FileWrapper, k, v) = setindex!(fw.file, k, v)
Base.haskey(fw::FileWrapper, k) = haskey(fw.file, k)
Base.write(fw::FileWrapper, x) = write(fw.file, x)
Base.write(fw::FileWrapper, k, x) = write(fw.file, k, x)

# To allow switching between JLD and JLD2:
const UnknownType = Union{JLD.UnsupportedType, JLD2.UnknownType}
const JLDFile = Union{FileWrapper{<: JLD.JldFile}, FileWrapper{<: JLD2.JLDFile}, JLD.JldFile, JLD2.JLDFile}

include("helpers.jl")
export enable_benchmarks, disable_benchmarks, print_timer, reset_timer!
include("linalg/general.jl")
include("linalg/UDT.jl")
# include("linalg/complex.jl") # TODO
include("linalg/blockdiagonal.jl")
include("flavors/abstract.jl")
include("models/abstract.jl")
include("lattices/abstract.jl")

include("configurations.jl")
export Discarder, ConfigRecorder, BufferedConfigRecorder, RelativePath, AbsolutePath
include("Measurements.jl")
export measurements, observables

include("lattices/lattice_iterators.jl")
include("lattices/square.jl")
include("lattices/chain.jl")
include("lattices/cubic.jl")
include("lattices/honeycomb.jl")
include("lattices/triangular.jl")
include("lattices/ALPS.jl")
include("lattices/deprecated.jl")
# export directions, RawMask, DistanceMask # maybe getorder?
export AbstractLattice, Chain, SquareLattice, CubicLattice, TriangularLattice, ALPSLattice
export EachSite, EachSiteAndFlavor, OnSite, EachSitePair, EachSitePairByDistance,
        EachLocalQuadByDistance, EachLocalQuadBySyncedDistance,
        Sum, ApplySymmetries
export neighbors, directions

include("flavors/MC/MC.jl")
include("flavors/DQMC/main.jl")
export Greens, GreensAt, CombinedGreensIterator
export GreensMatrix, dagger
export boson_energy_measurement, greens_measurement, occupation, magnetization
export charge_density, charge_density_correlation, charge_density_susceptibility
export spin_density, spin_density_correlation, spin_density_susceptibility
export pairing, pairing_correlation, pairing_susceptibility
export current_current_susceptibility, superfluid_density
export noninteracting_energy, interacting_energy, total_energy

export EmptyScheduler, SimpleScheduler, AdaptiveScheduler
export Adaptive, NoUpdate, LocalSweep
export GlobalFlip, GlobalShuffle, SpatialShuffle, TemporalShuffle,
        Denoise, DenoiseFlip, StaggeredDenoise
export ReplicaExchange, ReplicaPull, connect, disconnect
# export mask, uniform_fourier, structure_factor, SymmetryWrapped, swave, eswave

include("models/Ising/IsingModel.jl")
include("models/HubbardModel/HubbardModel.jl")
export IsingEnergyMeasurement, IsingMagnetizationMeasurement

include("FileIO.jl")
export save, load, resume!

export reset!
export run!, resume!, replay!
export Model, IsingModel, HubbardModel, HubbardModelAttractive, HubbardModelRepulsive
export MonteCarloFlavor, MC, DQMC
export greens, greens!, lattice, model, parameters

# For extending
export AbstractMeasurement, Model

function __init__()
    # @require LatPhysBase="eec5c15a-e8bd-11e8-0d23-6799ca40c963" include("lattices/LatPhys.jl")
    @require LatticePhysics = "53011200-ee7a-11e8-39f1-5f3e57afe4fd" include("lattices/LatPhys.jl")
    @require LightXML = "9c8b4983-aa76-5018-a973-4c85ecc9e179" include("lattices/ALPS.jl")
    @require MPI = "da04e1cc-30fd-572f-bb4f-1f8673147195" begin
        include("mpi.jl")
        export mpi_queue

        include("flavors/DQMC/updates/mpi_updates.jl")
        export MPIReplicaExchange, MPIReplicaPull
    end
end

end # module


text:
module AppendedTriangulationsTests

using Test
using Gridap.ReferenceFEs
using Gridap.Arrays
using Gridap.Geometry
using Gridap.Fields
using Gridap.Integration
using LinearAlgebra: ⋅
using Gridap.CellData

domain = (0,1,0,1)
partition = (10,10)
model = CartesianDiscreteModel(domain,partition)

ncells = num_cells(model)
nin = ceil(Int,2*ncells/3)
cell_to_mask = fill(false,ncells)
cell_to_mask[1:nin] .= true

grid = get_grid(model)

trian_in = RestrictedTriangulation(grid,cell_to_mask)
trian_out = RestrictedTriangulation(grid,collect(Bool, .! cell_to_mask))

trian = lazy_append(trian_out,trian_in)
test_triangulation(trian)

order = 1
quad = CellQuadrature(trian,2*order)
quad_in = CellQuadrature(trian_in,2*order)
quad_out = CellQuadrature(trian_out,2*order)

q = get_coordinates(quad)
w = get_weights(quad)
@test isa(q,AppendedArray)
@test isa(w,AppendedArray)

# Append triangulations of different cell type

domain = (0,1,0,1)
partition = (10,10)
grid1 = CartesianGrid(domain,partition)

domain = (1,2,0,1)
partition = (10,10)
grid2 = simplexify(CartesianGrid(domain,partition))

trian = lazy_append(grid1,grid2)
test_triangulation(trian)

end # module


text:
######################################################################
# smc_multi.jl
#
# Like smc.jl, but for multi-agent trajectories
#
# In theory, this shouldn't be so different from single version
######################################################################

function smc_trajectory(em::ErgodicManagerR2, vtm::VTM; verbose::Bool=true, umax::Float64=1.0)

	num_agents = length(vtm)
	N = vtm[1].N
	dt = vtm[1].h
	#ud = VVF(tm.N)
	#xd = VVF(tm.N+1)
	xds = VVVF(num_agents)
	uds = VVVF(num_agents)
	for j = 1:num_agents
		xds[j] = VVF(N+1)
		uds[j] = VVF(N)
		xds[j][1] = deepcopy(vtm[j].x0)
	end

	ck = zeros(em.K+1, em.K+1)
	for n = 1:N

		# update ck
		update_ck(em, xds, n, ck)

		# compute B
		for j = 1:num_agents
			Bx, By = compute_Bj(em, xds[j], n, dt, ck, num_agents)

			# normalize
			den = sqrt(Bx*Bx + By*By)
			ux = -umax * Bx / den
			uy = -umax * By / den

			uds[j][n] = [ux, uy]
			xds[j][n+1] = integrate(vtm[j], xds[j][n], uds[j][n])
		end
	end

	return vvvf2vvf(xds, uds)
end

function update_ck(em::ErgodicManagerR2, xds::VVVF, n::Int, ck::MF)
	num_agents = length(xds)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	for k1 = 0:em.K
		for k2 = 0:em.K
			hk = em.hk[k1+1, k2+1]
			ck[k1+1, k2+1] *= (n-1)
			for j = 1:num_agents
				dxn = xds[j][n][1] - x_min(em)
				dyn = xds[j][n][2] - y_min(em)

				cx = cos(k1*pi*dxn/Lx)
				cy = cos(k2*pi*dyn/Ly)

				# compute ck
				ck[k1+1, k2+1] += cx*cy/(hk*num_agents)
			end
			ck[k1+1, k2+1] /= n
		end
	end
end

# TODO: this can be done a lot quicker
function compute_Bj(em::ErgodicManagerR2, xd::VVF,n::Int,h::Float64,ck::MF, num_agents::Int)
	Lx = em.domain.lengths[1]
	Ly = em.domain.lengths[2]
	dxn = xd[n][1] - x_min(em)
	dyn = xd[n][2] - y_min(em)
	Bx = By = 0.0
	for k1 = 0:em.K
		for k2 = 0:em.K
			hk = em.hk[k1+1, k2+1]
			cx = cos(k1*pi*dxn/Lx)
			cy = cos(k2*pi*dyn/Ly)

			# compute ck
			#fk = ( (n-1)*ck[k1+1,k2+1] + cx*cy/hk ) / n
			#ck[k1+1, k2+1] = fk

			# t * Lambda_k * (c_k - phi_k)
			de = ck[k1+1, k2+1] - em.phik[k1+1,k2+1]
			LS = em.Lambda[k1+1,k2+1] * h*n * num_agents * de

			# multiplying by gradient
			Bx += -LS * k1*pi * sin(k1*pi*dxn/Lx) * cy / (hk*Lx)
			By += -LS * k2*pi * cx * sin(k2*pi*dyn/Ly) / (hk*Ly)
		end
	end
	return Bx, By
end


text:
include("euler/euler.jl")
using .Primes: prime_numbers

compute(n::Int)::Int = prime_numbers(20n)[n]

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
"""
    ElementVectorSampler

A sampler that makes it easy to draw signed element vectors (of a specified dimension)
uniformly at random.

Notice that, w.p. `1/2`, each element vector is multiplied by `-1`. This is useful for
implementing random walks.

An `element vector` `e_i` is a vector in which the `i`th element is a `1`, and all other
elements are `0`.

##### Fields
- `dimension::Int` : the dimension (i.e. length) of element vectors.
"""
struct ElementVectorSampler <: Sampleable{Multivariate, Discrete}
    dimension::Int
end

"""
    EVecSampler

Alias for `ElementVectorSampler`
"""
EVecSampler = ElementVectorSampler

"""
    Base.show(io::IO, elementvectorsampler::ElementVectorSampler) -> nothing

Show an `ElementVectorSampler` in a human-friendly manner.
"""
function Base.show(io::IO, elementvectorsampler::ElementVectorSampler)
    println(io, typeof(elementvectorsampler))
    println(io, " dimension: ", elementvectorsampler.dimension)
    return nothing
end

"""
    Base.rand(rng::AbstractRNG, elementvectorsampler::ElementVectorSampler) -> Vector{Int}

Return a vector with a specified `length` and with exactly one element equal to `1`.
"""
function Base.rand(rng::AbstractRNG, elementvectorsampler::ElementVectorSampler)
    vector = zeros(Int, elementvectorsampler.dimension)
    # pick an index to replace with 1 or -1 (uses multiple dispatch)
    vector[rand(1:elementvectorsampler.dimension)] = 1*rand((-1,1))
    return vector
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
# submarine cablesb

Pkg.add("Downloads")

using GLMakie, CairoMakie, Downloads
using GeoMakie, GeoMakie.GeoJSON, GeoMakie.GeoInterface
CairoMakie.activate!()
let
    # data from
    # https://github.com/telegeography/www.submarinecablemap.com
    urlPoints = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/landing-point/landing-point-geo.json"
    urlCables = "https://raw.githubusercontent.com/telegeography/www.submarinecablemap.com/master/web/public/api/v3/cable/cable-geo.json"

    landPoints = Downloads.download(urlPoints, IOBuffer())
    landCables = Downloads.download(urlCables, IOBuffer())

    land_geoPoints = GeoJSON.read(seekstart(landPoints))
    land_geoCables = GeoJSON.read(seekstart(landCables))

    toPoints = GeoMakie.geo2basic(land_geoPoints)
    #toLines = GeoMakie.geo2basic(land_geoCables) # this should probably be supported.
    feat = GeoInterface.features(land_geoCables)
    toLines = GeoInterface.coordinates.(GeoInterface.geometry.(feat))
    splitLines = []
    for i in 1:length(toLines)
        for j in 1:length(toLines[i])
            push!(splitLines, Point2f.(toLines[i][j]))
        end
    end

    function plotCables()
        fig, ax, = CairoMakie.scatter(toPoints; markersize= 5, color = 1:length(toPoints),
            colormap = :plasma, figure = (;resolution = (1200,800), fontsize = 24))
        [CairoMakie.lines!(ax, splitLines[i]; linewidth=0.85) for i in 1:length(splitLines)]
        CairoMakie.lines!(GeoMakie.coastlines(); color = :white, linewidth = 0.35)
        limits!(ax,-185,185,-95,95)
        ax.xticks = -180:60:180
        ax.yticks = -90:30:90
        ax.xtickformat = "{:d}ᵒ"
        ax.ytickformat = "{:d}ᵒ"
        fig
    end
    fig = with_theme(plotCables, theme_dark())
    display(fig)
end

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
include("ncmonomial.jl")
include("ncoperators.jl")
include("nccomparison.jl")


text:
using MLDatasets: MNIST
using Plots
using Flux: onehotbatch, onecold, onehot
using Random
include("src/NeuralNet.jl")

### Preparando dados

xtrain, ytrain = MNIST.traindata();
xtest, ytest = MNIST.testdata();

ytrain_OneHot = onehotbatch(ytrain,0:9);
ytest_OneHot = onehotbatch(ytest,0:9);

Gray.(xtrain[:,:,1])

xtrain_vec = zeros(28*28,60000)
xtest_vec = zeros(28*28,10000)

for i in 1:size(xtrain)[3]
    xtrain_vec[:,i] = vec(xtrain[:,:,i])
end

for i in 1:size(xtest)[3]
    xtest_vec[:,i] = vec(xtest[:,:,i])
end

nn = NeuralNet.NeuralNetwork([28*28,30,10],"sigmoid");

train_data = NeuralNet.train!(nn, xtrain_vec, ytrain_OneHot, batch_size=10, epochs=30, η=3.0)
# train_data = NeuralNet.train!(nn, xtrain_vec, ytrain_OneHot, batch_size=10, epochs=30, η=3.0)

function getPlot_accuracy(train_data)
    the_plot = plot(train_data.accuracy*100,
             label="Accuracy",
             xlabel="Epochs",
             ylabel="Accuracy as %",
             title="Development of accuracy at each iteration for training set");
    return the_plot
end

function getPlot_cost(train_data)
    the_plot = plot(train_data.costs,
             label="Cost Function",
             xlabel="Epochs",
             ylabel="Cost",
             color="red",
             title="Development of cost at each iteration for training set");
    return the_plot
end

function plot_nn(train_data)
    Plots.plot(getPlot_accuracy(train_data), getPlot_cost(train_data), layout = (2, 1), size = (800, 600))
end

plot_nn(train_data)

acc = 0
for i in 1:10000
    x = onecold(NeuralNet.predict(nn,xtest_vec[:,i])).-1
    y = ytest[i]
    if x[1] ≈ y
        acc += 1
    end
end
acc /= 10000
## \approx 94.98%


text:
using Distributed
using JLD2, FileIO
using Logging
#if nprocs() == 1
#    addprocs(7)
#end
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)
#Pkg.instantiate()
@everywhere using LadderDGA

#io = open("lDGA.log","w+")
io = stdout
#logger = ConsoleLogger(io, Logging.Info,
#                  meta_formatter=Logging.default_metafmt,
#                  show_limited=true, right_justify=0)
logger = SimpleLogger(io)
global_logger(logger)

function run_sim(; cfg_file=nothing, res_prefix="", res_postfix="", save_results=true)
    @warn "assuming linear, continuous nu grid for chi/trilex"
    if cfg_file === nothing
        print("specify location of config file: ")
        cfg_file = readline()
    end
    mP, sP, env, kGrids, kGridLoc, freqRed_map, freqList, freqList_min, parents, ops, nFermi, nBose, shift, base, offset = readConfig(cfg_file)#

    last_λsp = nothing
    last_λspch = nothing
    for kIteration in 1:length(kGrids)
        kG = kGrids[kIteration]
        @info " ===== Iteration $(kIteration)/$(length(kGrids)) with $(kG.Nk) k points. ===== "
        νGrid, sumHelper_f, impQ_sp, impQ_ch, GImp_fft, GLoc_fft, Σ_loc, FUpDo, gImp, gLoc = setup_LDGA(kG, freqList, mP, sP, env);

        @info "Calculating local quantities: "
        flush(io)
        bubbleLoc = calc_bubble(νGrid, GImp_fft, kGridLoc, mP, sP)
        locQ_sp = calc_χ_trilex(impQ_sp.Γ, bubbleLoc, kGridLoc, νGrid, sumHelper_f, mP.U, mP, sP);
        locQ_ch = calc_χ_trilex(impQ_ch.Γ, bubbleLoc, kGridLoc, νGrid, sumHelper_f, -mP.U, mP, sP);

        Σ_ladderLoc = calc_Σ(locQ_sp, locQ_ch, bubbleLoc, GImp_fft, FUpDo,
                             kGridLoc, sumHelper_f, mP, sP)
        Σ_ladderLoc = Σ_ladderLoc .+ mP.n * mP.U/2.0;

        @info "Calculating bubble: "
        flush(io)
        bubble = calc_bubble(νGrid, GLoc_fft, kG, mP, sP);

        @info "Calculating χ and γ: "
        flush(io)
        nlQ_sp = calc_χ_trilex(impQ_sp.Γ, bubble, kG, νGrid, sumHelper_f, mP.U, mP, sP);
        nlQ_ch = calc_χ_trilex(impQ_ch.Γ, bubble, kG, νGrid, sumHelper_f, -mP.U, mP, sP);

        @info "Calculating λ correction: "
        flush(io)
        λ_sp, λ_spch  = λ_correction!(impQ_sp, impQ_ch, FUpDo, Σ_loc, Σ_ladderLoc, nlQ_sp, nlQ_ch, bubble, GLoc_fft, kG, νGrid, mP, sP, init_sp=last_λsp, init_spch=last_λspch)
        last_λsp = λ_sp
        last_λspch = λ_spch

        @info "Calculating Σ ladder: "
        flush(io)
        Σ_ladder = calc_Σ(nlQ_sp, nlQ_ch, bubble, GLoc_fft, FUpDo, kG, sumHelper_f, mP, sP)
        Σ_ladder_corrected = Σ_ladder .- Σ_ladderLoc .+ Σ_loc[1:size(Σ_ladder,1)]
        @info "Done."
        flush(io)

        if save_results
        fname = res_prefix*"lDGA_b$(mP.β)_U$(mP.U)_k$(kG.Ns)_"*String(sP.tc_type_f)*"_lambda"*String(sP.λc_type)*res_postfix*".jld2"
        @info "Writing to $(fname)"
        save(fname, "λ_sp", λ_sp, "λ_spch", λ_spch, "bubbleLoc", bubbleLoc, "locQ_sp", locQ_sp, "locQ_ch", locQ_ch, "bubble", bubble, "nlQ_ch", nlQ_ch, "nlQ_sp", nlQ_sp, "Σ_ladder", Σ_ladder, "Σ_ladderLoc", Σ_ladderLoc)
        end
    end
end

function run2(cfg_file)
     λ_sp, λ_new_sp, λ_new_ch, _, _, _, _, nlQ_ch, nlQ_sp, Σ_ladder, _ = run_sim(cfg_file=cfg_file)
    return λ_sp, λ_new_sp, λ_new_ch, nlQ_ch, nlQ_sp, Σ_ladder
end


text:
push!(LOAD_PATH, pwd()*"/../../src/")
using Distributions
using PyPlot
using BVB

true_d = Gamma(2.0, 3.0)
ll(x) = logpdf(true_d, x)
x = linspace(-5,30,500)

init_approx_d = Normal(0, 20)

max_iter = 300

# include("run_laplacian_gradient_boosting.jl")
include("run_laplacian_gradient_boosting_exclusive.jl")
show()

println("Done.")
readline(STDIN)

max_stars_repo_name: richardkwo/boosting-VB

max_issues_repo_name: richardkwo/boosting-VB

max_forks_repo_name: richardkwo/boosting-VB


text:
sparse_options = (; conic_solver = sparse_hypatia)

insts = OrderedDict()

insts["test"] = [
    ((5, 3, 3, true),),
    ((5, 3, 3, false),),
    ((15, 7, 5, true), sparse_options),
    ((15, 7, 5, false), sparse_options),
    ((50, 15, 10, true), sparse_options),
    ((50, 15, 10, false), sparse_options),
]

function matrixregression_insts(use_nat::Bool, max_p::Int)
    return [((2 * p, p, 5, use_nat), sparse_options) for p in vcat(10, 10:10:max_p)]
end

insts["nat"] = matrixregression_insts(true, 120)
insts["ext"] = matrixregression_insts(false, 60)

return (MatrixRegression, insts)


text:
module AWSQueue

using ...SimpleTasks.Types

import AWS, AWS.SQS
import SimpleTasks.Services.Queue

export AWSQueueService

type AWSQueueService <: QueueService
    env::AWS.AWSEnv
    name::ASCIIString
    url::ASCIIString
end

AWSQueueService(env::AWS.AWSEnv, name::ASCIIString) =
        AWSQueueService(env, name, get_queue_url(env, name))

function Base.string(queue::AWSQueueService)
    return "$(queue.name) from $(queue.url)"
end

"""
    get_queue_url(env::AWS.AWSEnv, queue_name::AbstractString)

Find the correct url with our aws environment and bucket name
"""
function get_queue_url(env::AWS.AWSEnv, queue_name::AbstractString)
    response = SQS.CreateQueue(env, queueName = queue_name)

    # try creating a queue, does not overwrite existing queue so it's safe
    if response.http_code != 200
        error("Unable to create/get queue: $queue_name, response:
            ($(response.http_code))")
    end

    return response.obj.queueUrl
end

"""
    Queue.pop_message(queue::AWSQueueService)

Pop a message of the aws queue.
"""
function Queue.pop_message(queue::AWSQueueService)
    receive_response = SQS.ReceiveMessage(queue.env; queueUrl = queue.url)

    if receive_response.http_code != 200
        error("Unable to retrieve message from $queue")
    end

    if length(receive_response.obj.messageSet) < 1
        return ""
    end

    receiptHandle=receive_response.obj.messageSet[1].receiptHandle

    delete_response = SQS.DeleteMessage(queue.env; queueUrl = queue.url,
        receiptHandle = receiptHandle)

    return strip(receive_response.obj.messageSet[1].body)
end

"""
    to_SQS_attribute(attribute::Array{Tuple{AbstractString,
        Union{AbstractString, Number, Array{Uint8, 1}}}, 1})

Convert the attribute pair into AWS.SQSMessageAttributeType
"""
function to_SQS_attribute(attribute::Pair{AbstractString,
        Union{AbstractString, Number}})
    return AWS.SQS.MessageAttributeType(name = attribute[1],
        value = AWS.SQS.MessageAttributeValueType(
            dataType = typeof(attribute[2]), stringValue = attribute[2]))
end
function to_SQS_attribute(attribute::Pair{AbstractString, Array{UInt8, 1}})
    return AWS.SQS.MessageAttributeType(name = attribute[1],
        value = AWS.SQS.MessageAttributeValueType(
            dataType = "Binary", stringValue = attribute[2]))
end

"""
Push a message onto the aws queue.
Inputs:
    attributes - should be a Dict of attributes that can either bet a String,
        a Number, or a Byte array for Binary data
"""
function Queue.push_message(queue::AWSQueueService;
        attributes::Dict{AbstractString,
            Union{AbstractString, Number, Array{UInt8, 1}}} = Dict{AbstractString,
            Union{AbstractString, Number, Array{UInt8, 1}}}(),
        message_body::AbstractString = "")
    message_attributes = nothing
    if !isempty(attributes)
        message_attributes = map(to_SQS_attribute, attributes)
    end
    send_response = AWS.SQS.SendMessage(queue.env; queueUrl = queue.url, messageBody =
        message_body, messageAttributeSet=message_attributes)

    if send_response.http_code != 200
        error("Unable to send message $message_body to $queue")
    end

    return true
end

end # module AWSQueue

max_forks_repo_name: seung-lab/SimpleTasks


text:
"""
Information Gain.

The uncertainty of the starting node, minus the weighted impurity of
two child nodes.
"""
function info_gain(left, right, current_uncertainty::Float64)::Float64
    p = Float64(length(left)) / (length(left) + length(right))
    current_uncertainty - p * gini(left) - (1 - p) * gini(right)
end