text:
using RayTracer

@testset "creating_a_ray" begin
    origin = RayTracer.point(1,2,3)
    direction = RayTracer.vector(4,5,6)
    r = RayTracer.ray(origin, direction)
    @test r.origin ≈ origin
    @test r.direction ≈ direction
end

@testset "computing_point_from_distance" begin
    origin = RayTracer.point(2,3,4)
    direction = RayTracer.vector(1,0,0)
    r = RayTracer.ray(origin, direction)
    @test RayTracer.position(r,0) ≈ RayTracer.point(2,3,4)
    @test RayTracer.position(r,1) ≈ RayTracer.point(3,3,4)
    @test RayTracer.position(r,-1) ≈ RayTracer.point(1,3,4)
    @test RayTracer.position(r,2.5) ≈ RayTracer.point(4.5,3,4)
end

@testset "ray_intersecting_two_points_shpere" begin
    origin = RayTracer.point(0,0,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].t ≈ 4
    @test xs[2].t ≈ 6
end

@testset "ray_intersecting_tangent_point" begin
    origin = RayTracer.point(0,1,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].t ≈ 5
    @test xs[2].t ≈ 5
end

@testset "ray_missing_intersection_sphere" begin
    origin = RayTracer.point(0,2,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 0
end

@testset "ray_intersecting_center_shpere" begin
    origin = RayTracer.point(0,0,0)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].t ≈ -1
    @test xs[2].t ≈ 1
end

@testset "ray_intersecting_front_sphere" begin
    origin = RayTracer.point(0,0,5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].t ≈ -6
    @test xs[2].t ≈ -4
end

@testset "intersection_creation" begin
    s = RayTracer.sphere()
    i = RayTracer.intersection(3,s)
    @test i.t ≈ 3
    @test i.object === s
end

@testset "aggregating_intersections" begin
    s = RayTracer.sphere()
    i1 = RayTracer.intersection(1,s)
    i2 = RayTracer.intersection(2,s)
    xs = [i1,i2]
    @test length(xs) == 2
    @test xs[1].t ≈ 1
    @test xs[2].t ≈ 2
end

@testset "intersect_has_obhect_on_intersection" begin
    origin = RayTracer.point(0,0,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].object === s
    @test xs[2].object === s
end

@testset "hit_when_all_t_are_positive" begin
    s = RayTracer.sphere()
    i1 = RayTracer.intersection(1,s)
    i2 = RayTracer.intersection(2,s)
    xs = [i2,i1]
    i = RayTracer.hit(xs)
    @test i === i1
end

@testset "hit_when_one_t_is_positive" begin
    s = RayTracer.sphere()
    i1 = RayTracer.intersection(-1,s)
    i2 = RayTracer.intersection(1,s)
    xs = [i2,i1]
    i = RayTracer.hit(xs)
    @test i === i2
end

@testset "hit_when_all_t_are_negative" begin
    s = RayTracer.sphere()
    i1 = RayTracer.intersection(-2,s)
    i2 = RayTracer.intersection(-1,s)
    xs = [i2,i1]
    i = RayTracer.hit(xs)
    @test i === nothing
end

@testset "hit_lowest_positive_t" begin
    s = RayTracer.sphere()
    i1 = RayTracer.intersection(5,s)
    i2 = RayTracer.intersection(7,s)
    i3 = RayTracer.intersection(-3,s)
    i4 = RayTracer.intersection(2,s)
    xs = [i1,i2,i3,i4]
    i = RayTracer.hit(xs)
    @test i === i4
end

@testset "translating_ray" begin
    origin = RayTracer.point(1,2,3)
    direction = RayTracer.vector(0,1,0)
    r = RayTracer.ray(origin, direction)
    m = RayTracer.translation(3,4,5)
    r2 = RayTracer.transform(r,m)
    @test r2.origin ≈ RayTracer.point(4,6,8)
    @test r2.direction ≈ RayTracer.vector(0,1,0)
end

@testset "scaling_ray" begin
    origin = RayTracer.point(1,2,3)
    direction = RayTracer.vector(0,1,0)
    r = RayTracer.ray(origin, direction)
    m = RayTracer.scaling(2,3,4)
    r2 = RayTracer.transform(r,m)
    @test r2.origin ≈ RayTracer.point(2,6,12)
    @test r2.direction ≈ RayTracer.vector(0,3,0)
end

@testset "sphere_default_transformation" begin
    s = RayTracer.sphere()
    @test s.transform ≈ Matrix(I,4,4)
end

@testset "changing_sphere_transformation" begin
    s = RayTracer.sphere()
    t = RayTracer.translation(2,3,4)
    RayTracer.set_transform!(s,t)
    @test s.transform ≈ t
end

@testset "intersecting_scaled_sphere_with_ray" begin
    origin = RayTracer.point(0,0,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    RayTracer.set_transform!(s,RayTracer.scaling(2,2,2))
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 2
    @test xs[1].t ≈ 3
    @test xs[2].t ≈ 7
end

@testset "intersecting_translated_sphere_with_ray" begin
    origin = RayTracer.point(0,0,-5)
    direction = RayTracer.vector(0,0,1)
    r = RayTracer.ray(origin, direction)
    s = RayTracer.sphere()
    RayTracer.set_transform!(s,RayTracer.translation(5,0,0))
    xs = RayTracer.intersect(s,r)
    @test length(xs) == 0
end


text:
module ThriftUtilsTests

using Thrift
using Test

struct _enum_TestEnum
    BOOLEAN::Int32
    INT32::Int32
    INT64::Int32
    INT96::Int32
    FLOAT::Int32
    DOUBLE::Int32
    BYTE_ARRAY::Int32
    FIXED_LEN_BYTE_ARRAY::Int32
end
const TestEnum = _enum_TestEnum(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

function test_enum()
    println("    enum...")
    @test enumstr(TestEnum, TestEnum.BOOLEAN) == "BOOLEAN"
    @test_throws ErrorException enumstr(TestEnum, Int32(11))
end

function test_container_check()
    println("    iscontainer/isplain...")
    for T in (Bool, UInt8, Float64, Int16, Int32, Int64, Vector{UInt8}, String)
        @test Thrift.isplain(T)
        @test !Thrift.iscontainer(T)
    end
    for T in (Dict, Set, Any)
        @test !Thrift.isplain(T)
        @test Thrift.iscontainer(T)
    end
end

mutable struct TestMetaAllTypes <: Thrift.TMsg
  bool_val::Bool
  byte_val::UInt8
  i16_val::Int16
  i32_val::Int32
  i64_val::Int64
  double_val::Float64
  string_val::String
  TestMetaAllTypes() = (o=new(); fillunset(o); o)
end # mutable struct AllTypes

function test_meta()
    @test !Thrift.isplain(TestMetaAllTypes)
    @test Thrift.iscontainer(TestMetaAllTypes)

    types = TestMetaAllTypes()
    @test !isfilled(types)
    @test !isfilled(types, :bool_val)

    set_field!(types, :bool_val, true)
    @test isfilled(types, :bool_val)
    @test !isfilled(types)

    set_field!(types, :byte_val,    UInt8(1))
    set_field!(types, :i16_val,     Int16(1))
    set_field!(types, :i32_val,     Int32(1))
    set_field!(types, :i64_val,     Int64(1))
    set_field!(types, :double_val,  1.1)
    set_field!(types, :string_val,  "1")
    @test isfilled(types)
    @test isinitialized(types)

    types2 = TestMetaAllTypes()
    @test !isfilled(types2)
    copy!(types2, types)
    @test isfilled(types2)

    nothing
end

println("Testing utils functions...")
test_enum()
test_container_check()
test_meta()
println("passed.")

end

ext: jl


text:
using Xyce
using Test

@testset "Xyce" begin
@testset "VRC.sp tran" begin
	res = Xyce.simulate("VRC.sp")
	@test haskey(res, "TIME")
	@test haskey(res, "V(INP)")
	@test haskey(res, "V(OUT)")
	time = res["TIME"]
	inp = res["V(INP)"]
	out = res["V(OUT)"]
	@test length(time) > 100
	@test time[1] == 0.0
	@test time[end] == 0.002
	@test maximum(inp) ≈ 1.0 atol=0.01
	@test minimum(inp) ≈ -1.0 atol=0.01
	@test inp[1] == 0
	@test inp[end] ≈ 0 atol=1e-12
	@test out[1] == 0
	@test maximum(out)/maximum(inp) ≈ 0.862 atol=0.001
	i1 = argmax(inp)
	@test out[i1]/inp[i1] ≈ 0.710 atol=0.001
end
@testset "V1D1.sp tran" begin
	res = Xyce.simulate("V1D1.sp")
	@test haskey(res, "TIME")
	@test haskey(res, "V(NET1)")
	@test haskey(res, "I(V1)")
	time = res["TIME"]
	v1 = res["V(NET1)"]
	i1 = res["I(V1)"]
	@test time[1] == 0.0
	@test time[end] == 0.001
	@test maximum(v1) ≈ 20 atol=0.01
	@test minimum(v1) ≈ -20 atol=0.01
	@test v1[1] == 0
	@test v1[end] ≈ 0 atol=1e-12
	@test i1[1] ≈ 0 atol=1e-20
	@test minimum(i1) < -29
	@test maximum(i1) ≈ 0 atol=1e-4
	@test length(v1) > 75
	idx = argmin(i1)
	@test v1[idx] > 19.9
end
end


text:
using Combinatorics

function part1(input)
    D = parse_input(input)
    return foldl(max, (total_happiness(D, order) for order in permutations(1:8)))
end

function part2(input)
    D = parse_input(input)
    return foldl(max, (total_happiness2(D, order) for order in permutations(1:8)))
end

function parse_input(input)
    D = fill(0, 8, 8)
    lines = eachline(input)
    for i = 1:8
        for j = 1:7
            words = split(first(lines))
            happiness = parse(Int, words[4])
            if words[3] == "lose"
                happiness *= -1
            end
            D[i, j + (j >= i)] += happiness
            D[j + (j >= i), i] += happiness
        end
    end
    return D
end

total_happiness(D, order) = sum(D[order[i], order[mod1(i + 1, 8)]] for i = 1:8)
total_happiness2(D, order) = sum(D[order[i], order[i + 1]] for i = 1:7)


text:
const mseed_support_data = "# Mini-SEED Support
SeisIO supports mini-SEED, the \"data-only\" extension of the SEED (Standard for
the Exchange of Earthquake Data) file format.

## Supported Blockette Types

| Blockette                                     | Key in `:misc`  |
|:----------------------------------------------|:----------------|
| [100] Sample Rate Blockette                   |                 |
| [201] Murdock Event Detection Blockette       | seed_event ⁽¹⁾  |
| [300] Step Calibration Blockette              | seed_calib ⁽¹⁾  |
| [310] Sine Calibration Blockette              | seed_calib ⁽¹⁾  |
| [320] Pseudo-random Calibration Blockette     | seed_calib ⁽¹⁾  |
| [390] Generic Calibration Blockette           | seed_calib ⁽¹⁾  |
| [500] Timing Blockette                        | seed_timing     |
| [1000] Data Only SEED Blockette               |                 |
| [1001] Data Extension Blockette               |                 |
| [2000] Variable Length Opaque Data Blockette  | seed_opaque ⁽²⁾ |

Notes on the Table
1. Stored in `:misc` in String arrays; each blockette gets a single String in the named key, separated by a newline character (\\n).
2. Stored in `:misc[\"seed_opaque\"]`, which contains raw (UInt8) byte
vectors of all data in each packet.

## Supported Data Encodings

| Format  | Data Encoding                                               |
|---------|:------------------------------------------------------------|
| 0       | ASCII [Saved to `:misc[\"seed_ascii\"]`, not `:x`]          |
| 1       | Int16 unencoded                                             |
| 3       | Int32 unencoded                                             |
| 4       | Float32 unencoded                                           |
| 5       | Float64 unencoded [Converted to Float32]                    |
| 10      | Steim-1                                                     |
| 11      | Steim-2                                                     |
| 12      | GEOSCOPE multiplexed, 16-bit gain ranged, 3-bit exponent    |
| 15      | GEOSCOPE multiplexed, 16-bit gain ranged, 4-bit exponent    |
| 16      | CDSN, 16-bit gain ranged                                    |
| 30      | SRO                                                         |
| 32      | DWWSSN gain ranged                                          |

### Unsupported Data Encodings
These have never been encountered by SeisIO. If support is needed, please send example files.

| Format  | Data Encoding                                               |
|---------|:------------------------------------------------------------|
| 2       | Int24 unencoded                                             |
| 12      | GEOSCOPE multiplexed, 24-bit integer                        |
| 15      | US National Network                                         |
| 17      | Graefenberg, 16-bit gain ranged                             |
| 18      | IPG - Strasbourg, 16-bit gain ranged                        |
| 19      | Steim-3                                                     |
| 31      | HGLP                                                        |
| 33      | RSTN 16-bit gain ranged                                     |
"

const dataless_support_data = "# Dataless SEED Support
SeisIO supports reading of dataless SEED meta-data files, a popular extension of the SEED (Standard for the Exchange of Earthquake Data) file format.

## Supported Blockette Types

| Blockette                                             |
|:------------------------------------------------------|
| [10] Volume Identifier Blockette                      |
| [11] Volume Station Header Index Blockette            |
| [12] Volume Time Span Index Blockette                 |
| [31] Comment Description Blockette                    |
| [33] Generic Abbreviation Blockette                   |
| [34] Units Abbreviations Blockette                    |
| [41] FIR Dictionary Blockette                         |
| [43] Response (Poles & Zeros) Dictionary Blockette    |
| [44] Response (Coefficients) Dictionary Blockette     |
| [45] Response List Dictionary Blockette               |
| [47] Decimation Dictionary Blockette                  |
| [48] Channel Sensitivity/Gain Dictionary Blockette    |
| [50] Station Identifier Blockette                     |
| [52] Channel Identifier Blockette                     |
| [53] Response (Poles & Zeros) Blockette               |
| [54] Response (Coefficients) Blockette                |
| [57] Decimation Blockette                             |
| [58] Channel Sensitivity/Gain Blockette               |
| [59] Channel Comment Blockette                        |
| [60] Response Reference Blockette                     |
| [61] FIR Response Blockette                           |

## Unsupportable Blockette Types
These blockettes will probably never be supported as
their information lies outside the scope of SeisIO. At
high verbosity (v > 2), their information is dumped to
stdout.

| Blockette                                             |
|:------------------------------------------------------|
| [30] Data Format Dictionary Blockette                 |
| [32] Cited Source Dictionary Blockette                |
| [51] Station Comment Blockette                        |

"

const broken_resp_data = "The following is a list of breaking SEED RESP issues that we've encountered in real data. Files with these issues don't read correctly into any known program (e.g., ObsPy, SeisIO, SEED C libraries).

| Network | Station(s)    | Problem(s)
| :----   | :----         | :----
| CN      | (broadbands)  | B058F05-06 contain units; should be B053F05-06
"

"mseed_support() shall info. dump mini-SEED blockette support to stdout."
mseed_support
function mseed_support()
  show(stdout, MIME("text/plain"), Markdown.parse(mseed_support_data))
  return nothing
end

"dataless_support() shall info. dump dataless SEED blockette support to stdout."
function dataless_support()
  show(stdout, MIME("text/plain"), Markdown.parse(dataless_support_data))
  return nothing
end

"seed_support() shall info. dump ALL SEED support info to stdout."
function seed_support()
  show(stdout, MIME("text/plain"), Markdown.parse(mseed_support_data))
  show(stdout, MIME("text/plain"), Markdown.parse(dataless_support_data))
  show(stdout, MIME("text/plain"), Markdown.parse(broken_resp_data))
  return nothing
end

"resp_wont_read() shall info. dump to stdout about broken resp files."
function resp_wont_read()
  show(stdout, MIME("text/plain"), Markdown.parse(broken_resp_data))
  return nothing
end


text:
using Revise
	using DiffEqOperators, ForwardDiff, DifferentialEquations
	using BifurcationKit, LinearAlgebra, Plots, SparseArrays, Parameters, Setfield
	const BK = BifurcationKit

norminf(x) = norm(x, Inf)

function Laplacian2D(Nx, Ny, lx, ly, bc = :Dirichlet)
	hx = 2lx/Nx
	hy = 2ly/Ny
	D2x = CenteredDifference(2, 2, hx, Nx)
	D2y = CenteredDifference(2, 2, hy, Ny)
	if bc == :Neumann
		Qx = Neumann0BC(hx)
		Qy = Neumann0BC(hy)
	elseif  bc == :Dirichlet
		Qx = Dirichlet0BC(typeof(hx))
		Qy = Dirichlet0BC(typeof(hy))
	end
	D2xsp = sparse(D2x * Qx)[1]
	D2ysp = sparse(D2y * Qy)[1]

	A = kron(sparse(I, Ny, Ny), D2xsp) + kron(D2ysp, sparse(I, Nx, Nx))
	return A, D2x
end

function NL!(f, u, p, t = 0.)
	@unpack r, μ, ν, c3, c5 = p
	n = div(length(u), 2)
	u1 = @view u[1:n]
	u2 = @view u[n+1:2n]

	ua = u1.^2 .+ u2.^2

	f1 = @view f[1:n]
	f2 = @view f[n+1:2n]

	@. f1 .= r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	@. f2 .= r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

NL(u, p) = NL!(similar(u), u, p)

function Fcgl!(f, u, p, t = 0.)
	mul!(f, p.Δ, u)
	f .= f .+ NL(u, p)
end

Fcgl(u, p, t = 0.) = Fcgl!(similar(u), u, p, t)

# computation of the first derivative
# d1Fcgl(x, p, dx) = ForwardDiff.derivative(t -> Fcgl(x .+ t .* dx, p), 0.)

d1NL(x, p, dx) = ForwardDiff.derivative(t -> NL(x .+ t .* dx, p), 0.)

function dFcgl(x, p, dx)
	f = similar(dx)
	mul!(f, p.Δ, dx)

	nl = d1NL(x, p, dx)
	f .= f .+ nl
end

function Jcgl(u, p, t = 0.)
	@unpack r, μ, ν, c3, c5, Δ = p

	n = div(length(u), 2)
	u1 = @view u[1:n]
	u2 = @view u[n+1:2n]

	ua = u1.^2 .+ u2.^2

	f1u = zero(u1)
	f2u = zero(u1)
	f1v = zero(u1)
	f2v = zero(u1)

	@. f1u =  r - 2 * u1 * (c3 * u1 - μ * u2) - c3 * ua - 4 * c5 * ua * u1^2 - c5 * ua^2
	@. f1v = -ν - 2 * u2 * (c3 * u1 - μ * u2)  + μ * ua - 4 * c5 * ua * u1 * u2
	@. f2u =  ν - 2 * u1 * (c3 * u2 + μ * u1)  - μ * ua - 4 * c5 * ua * u1 * u2
	@. f2v =  r - 2 * u2 * (c3 * u2 + μ * u1) - c3 * ua - 4 * c5 * ua * u2 ^2 - c5 * ua^2

	jacdiag = vcat(f1u, f2v)

	Δ + spdiagm(0 => jacdiag, n => f1v, -n => f2u)
end

####################################################################################################
Nx = 41*1
	Ny = 21*1
	n = Nx*Ny
	lx = pi
	ly = pi/2

	Δ = Laplacian2D(Nx, Ny, lx, ly)[1]
	par_cgl = (r = 0.5, μ = 0.1, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ))
	sol0 = 0.1rand(2Nx, Ny)
	sol0_f = vec(sol0)

####################################################################################################
eigls = EigArpack(1.0, :LM)
# eigls = eig_MF_KrylovKit(tol = 1e-8, dim = 60, x₀ = rand(ComplexF64, Nx*Ny), verbose = 1)
opt_newton = NewtonPar(tol = 1e-9, verbose = true, eigsolver = eigls, maxIter = 20)
opts_br = ContinuationPar(dsmax = 0.02, ds = 0.01, pMax = 2., detectBifurcation = 3, nev = 15, newtonOptions = (@set opt_newton.verbose = false), nInversion = 4)

	br, = @time continuation(Fcgl, Jcgl, vec(sol0), par_cgl, (@lens _.r), opts_br, verbosity = 0)

plot(br)
####################################################################################################
# Look for periodic orbits
f1 = DiffEqArrayOperator(par_cgl.Δ)
f2 = NL!
prob_sp = SplitODEProblem(f1, f2, sol0_f, (0.0, 120.0), @set par_cgl.r = 1.2; reltol = 1e-8, dt = 0.1)
prob = ODEProblem(Fcgl, sol0_f, (0.0, 120.0), (@set par_cgl.r = 1.2))#, jac = Jcgl, jac_prototype = Jcgl(sol0_f, par_cgl))
####################################################################################################
# sol = @time solve(prob, Vern9(); abstol=1e-14, reltol=1e-14)
sol = @time solve(prob_sp, ETDRK2(krylov=true); abstol=1e-14, reltol=1e-14, dt = 0.1) #1.78s
# sol = @time solve(prob, LawsonEuler(krylov=true, m=50); abstol=1e-14, reltol=1e-14, dt = 0.1)
# sol = @time solve(prob_sp, CNAB2(linsolve=LinSolveGMRES()); abstol=1e-14, reltol=1e-14, dt = 0.03)

plot(sol.t, [norm(v[1:Nx*Ny], Inf) for v in sol.u], xlims=(115, 120))

# plotting the solution as a movie
for ii = 1:20:length(sol.t)
	# heatmap(reshape(sol[1:Nx*Ny,ii],Nx,Ny),title="$(sol.t[ii])") |> display
end

####################################################################################################
# this encodes the functional for the Shooting problem
probSh = ShootingProblem(
	# we pass the ODEProblem encoding the flow and the time stepper
	prob_sp, ETDRK2(krylov = true),
	[sol[:, end]], abstol = 1e-10, reltol = 1e-8)

initpo = vcat(sol(116.), 4.9) |> vec
	probSh(initpo, @set par_cgl.r = 1.2) |> norminf

ls = GMRESIterativeSolvers(reltol = 1e-4, N = 2n + 1, maxiter = 50, verbose = false)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 25, linsolver = ls)
outpo, = @time newton(probSh, initpo, (@set par_cgl.r = 1.2), optn; normN = norminf)

heatmap(reshape(outpo[1:Nx*Ny], Nx, Ny), color = :viridis)

eig = EigKrylovKit(tol = 1e-7, x₀ = rand(2Nx*Ny), verbose = 2, dim = 40)
	opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= -0.01, pMax = 2.5, maxSteps = 32, newtonOptions = (@set optn.eigsolver = eig), nev = 15, precisionStability = 1e-3, detectBifurcation = 0, plotEveryStep = 1)
br_po, upo, = @time continuation(probSh, outpo, (@set par_cgl.r = 1.2), (@lens _.r),
		opts_po_cont;
		verbosity = 3,
		plot = true,
		# tangentAlgo = BorderedPred(),
		linearAlgo = MatrixFreeBLS(@set ls.N = probSh.M*2n+2),
		# callbackN = cb_ss,
		plotSolution = (x, p; kwargs...) -> heatmap!(reshape(x[1:Nx*Ny], Nx, Ny); color=:viridis, kwargs...),
		recordFromSolution = (u, p; k...) -> BK.getAmplitude(probSh, u, (@set par_cgl.r = p); ratio = 2), normC = norminf)

####################################################################################################
# automatic branch switching
jet  = BK.getJet(Fcgl, Jcgl)

ls = GMRESIterativeSolvers(reltol = 1e-4, maxiter = 50, verbose = false)
	optn = NewtonPar(verbose = true, tol = 1e-9,  maxIter = 25, linsolver = ls)
eig = EigKrylovKit(tol = 1e-7, x₀ = rand(2Nx*Ny), verbose = 2, dim = 40)
	opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.02, ds= 0.01, pMax = 2.5, maxSteps = 32, newtonOptions = (@set optn.eigsolver = eig), nev = 15, precisionStability = 1e-3, detectBifurcation = 0, plotEveryStep = 1)

Mt = 1 # number of time sections
br_po, = continuation(
	jet...,	br, 2,
	# arguments for continuation
	opts_po_cont,
	# probSh;
	ShootingProblem(Mt, prob_sp, ETDRK2(krylov = true); abstol = 1e-10, reltol = 1e-8) ;
	verbosity = 3, plot = true, ampfactor = 1.5, δp = 0.01,
	# callbackN = (state; kwargs...) -> (println("--> amplitude = ", BK.amplitude(state.x, n, M; ratio = 2));true),
	linearAlgo = MatrixFreeBLS(@set ls.N = Mt*2n+2),
	finaliseSolution = (z, tau, step, contResult; k...) ->begin
		BK.haseigenvalues(contResult) && Base.display(contResult.eig[end].eigenvals)
		return true
	end,
	recordFromSolution = (u, p; k...) -> BK.getAmplitude(p.prob, u, (@set par_cgl.r = p.p); ratio = 2),
	normC = norminf)

#ShootingProblem(1, par_cgl, prob_sp, ETDRK2(krylov = true)) ;


text:
using Distributed
import LineSearches: InitialStatic

dataDir = joinpath(dirname(pathof(gcSolver)), "..", "data")

""" Creates full vector of unknown values to be fit """
function getUnkVec()
    #kfwd, k4, k5, k16, k17, k22, k23, k27, endo, aendo, sort, krec, kdeg, k34, k35, k36, k37, k38, k39
    p = 0.1ones(18)

    p[1] = 0.001 # means of prior distributions from gc-cytokines paper
    p[8] = 1.0 #27
    p[9] = 1.0 # initial Treg stat
    p[10] = 0.5 # initial Thelp stat
    p[11] = 0.3 # initial NK stat
    p[12] = 0.3 # initial CD8 stat
    p[13:18] .= 0.001

    return p
end

""" Takes in full unkvec and constructs it into full fit parameters vector - TODO move this. """
function fitParams(ILs, unkVec::Vector{T}, recAbundances, CellType::String) where {T}
    kfbnd = 0.60
    farhatVec = getFarhatVec()
    paramvec = zeros(T, Nparams)
    paramvec[1:3] = ILs
    paramvec[4] = unkVec[1] #kfwd
    paramvec[5] = kfbnd * 10.0 #k1rev
    paramvec[6] = kfbnd * 144.0 #k2rev
    paramvec[7] = unkVec[2] #k4rev
    paramvec[8] = unkVec[3] #k5rev
    paramvec[9] = 12.0 * paramvec[8] / 1.5 #k10
    paramvec[10] = 63.0 * paramvec[8] / 1.5 #k11
    paramvec[11] = kfbnd * 0.065 #k13
    paramvec[12] = kfbnd * 438.0 #k14
    paramvec[13:16] = unkVec[4:7] #k16, k17, k22, k23
    paramvec[17] = kfbnd * 59.0 #k25
    paramvec[18] = unkVec[8]  #k27
    paramvec[19] = 5.0 #endo
    ke = 0.09
    aendo = 2.5
    sortF = 0.125
    krec = 0.009
    kdeg = 0.0068
    paramvec[20] = ke
    paramvec[21] = aendo
    paramvec[22] = sortF
    paramvec[23] = krec
    paramvec[24] = kdeg
    paramvec[25:29] = recAbundances * ke / (1.0 + (krec * (1.0 - sortF) / kdeg / sortF))
    if CellType == "Treg"
        paramvec[30] = unkVec[9] #initial stat
    elseif CellType == "Thelper"
        paramvec[30] = unkVec[10] #initial stat
    elseif CellType == "NK"
        paramvec[30] = unkVec[11] #initial stat
    elseif CellType == "CD8"
        paramvec[30] = unkVec[12] #initial stat
    else
        @assert false
    end
    paramvec[31:36] = unkVec[13:18]

    return paramvec
end

""" Adjusts the binding activity of ligand to match that of mutein """
function mutAffAdjust(p::Vector{T}, dfRow) where {T}
    p[5] = dfRow.IL2RaKD[1] * 0.6
    bgAdjust = (dfRow.IL2RBGKD[1] * 0.6) / p[8]
    p[6:10] .*= bgAdjust
    return p
end

""" Calculates squared error for a given unkVec. """
function resids(x::Vector{T})::T where {T}
    @assert all(x .>= 0.0)
    df = importData(true)
    #df = df[df.Ligand .!= "IL15", :]

    sort!(df, :Time)
    df.Time *= 60.0

    df.MeanPredict = similar(df.Mean, T)
    FutureDict = Dict()
    cost = 0.0

    for ligand in unique(df.Ligand)
        # Put the highest dose first so we catch a solving error early
        for dose in reverse(sort(unique(df.Dose)))
            if ligand == "IL15"
                ligVec = [0.0, dose, 0.0]
            else
                ligVec = [dose, 0.0, 0.0]
            end

            for cell in unique(df.Cell)
                idxx = findfirst(df.Cell .== cell)
                recpE = 10.0 .^ Vector{Float64}(df[idxx, [:IL2Ra, :IL2Rb, :gc, :IL15Ra, :IL7Ra]])
                vector = vec(fitParams(ligVec, x, recpE, cell))
                if ligand !== "IL15" && ligand !== "IL2"
                    vector = mutAffAdjust(vector, df[findfirst(df.Ligand .== ligand), [:IL2RaKD, :IL2RBGKD]])
                end
                idxs = (df.Dose .== dose) .& (df.Ligand .== ligand) .& (df.Cell .== cell)

                tpss = df[idxs, :Time]
                # Make sure duplicate times are not considered duplicates
                tpss += range(0.0, 0.01; length = length(tpss))

                # Regularize for exploding values
                cost += sum(softplus.(vector .- 1.0e6))
                cost += sum(softplus.(x .^ -1 .- 1.0e4))
                #cost += sum(softplus.(x[2] .^ -1))

                FutureDict[(dose, ligand, cell)] = @spawnat :any runCkine(tpss, vector; pSTAT5 = true)
            end
        end
    end

    for (key, val) in FutureDict
        idxs = (df.Dose .== key[1]) .& (df.Ligand .== key[2]) .& (df.Cell .== key[3])

        df[idxs, :MeanPredict] = fetch(val)
    end

    @assert all(df.MeanPredict .>= -0.01)

    predictsFrame = copy(df)
    delete!(predictsFrame, 1:size(predictsFrame, 1))
    for date in unique(df.Date)
        dateFilt = filter(row -> string(row["Date"]) .== date, df)
        correction = dateFilt.MeanPredict \ dateFilt.Mean
        #cost += softplus.(correction .- 1.0e3)
        dateFilt.MeanPredict .*= correction
        append!(predictsFrame, dateFilt)
    end
    cost += norm(predictsFrame.MeanPredict - predictsFrame.Mean)

    return cost
end

""" Gets inital unkowns, optimizes them, and returns parameters of best fit"""
function runFit(; itern = 1000000)
    x₀ = invsoftplus.(getUnkVec())

    opts = Optim.Options(iterations = itern, show_trace = true, extended_trace = true, allow_f_increases = true)
    lsi = InitialStatic(; alpha = 0.0001)
    fit = optimize((x) -> resids(softplus.(x)), x₀, LBFGS(; m = 100, alphaguess = lsi), opts, autodiff = :forward)

    @show fit

    return fit.minimizer
end

export getExpression, getUnkVec, fitParams, runCkine

""" Writes a CSV for conversion factors between predictions and pSTAT values """
function getDateConvDict()
    fitVec = importFit()
    fitVec = convert(Vector{Float64}, fitVec[!, :Fit])
    x = softplus.(fitVec)

    df = importData(true)

    sort!(df, :Time)
    df.Time *= 60.0

    df.MeanPredict = similar(df.Mean, Float64)
    FutureDict = Dict()

    for ligand in unique(df.Ligand)
        # Put the highest dose first so we catch a solving error early
        for dose in reverse(sort(unique(df.Dose)))
            if ligand == "IL15"
                ligVec = [0.0, dose, 0.0]
            else
                ligVec = [dose, 0.0, 0.0]
            end

            for cell in unique(df.Cell)
                idxx = findfirst(df.Cell .== cell)
                recpE = 10.0 .^ Vector{Float64}(df[idxx, [:IL2Ra, :IL2Rb, :gc, :IL15Ra, :IL7Ra]])
                vector = vec(fitParams(ligVec, x, recpE, cell))

                if ligand != "IL15" && ligand !== "IL2"
                    vector = mutAffAdjust(vector, df[findfirst(df.Ligand .== ligand), [:IL2RaKD, :IL2RBGKD]])
                end
                idxs = (df.Dose .== dose) .& (df.Ligand .== ligand) .& (df.Cell .== cell)

                tpss = df[idxs, :Time]
                # Make sure duplicate times are not considered duplicates
                tpss += range(0.0, 0.001; length = length(tpss))

                # Regularize for exploding values

                FutureDict[(dose, ligand, cell)] = @spawnat :any runCkine(tpss, vector; pSTAT5 = true)
            end
        end
    end

    for (key, val) in FutureDict
        idxs = (df.Dose .== key[1]) .& (df.Ligand .== key[2]) .& (df.Cell .== key[3])

        df[idxs, :MeanPredict] = fetch(val)
    end

    @assert all(df.MeanPredict .>= -0.01)

    dateConvDF = DataFrame(Date = [], Conv = [])
    for date in unique(df.Date)
        dateFilt = filter(row -> string(row["Date"]) .== date, df)
        convFact = dateFilt.MeanPredict \ dateFilt.Mean
        miniDF = DataFrame(Date = date, Conv = convFact)
        dateFilt.MeanPredict .*= dateFilt.MeanPredict \ dateFilt.Mean
        append!(dateConvDF, miniDF)
    end

    CSV.write(joinpath(dataDir, "DateConvFrame.csv"), dateConvDF)
end

""" Writes results of Farhat et al. model fitting to unknown vector """
function getFarhatVec()
    #kfwd, k4, k5, k16, k17, k22, k23, k27, endo, aendo, sort, krec, kdeg, k34, k35, k36, k37, k38, k39
    p = 0.1ones(23)
    p[1] = 0.003 #0.001 # means of prior distributions from gc-cytokines paper
    p[2] = 12.25
    p[3] = 0.08
    p[4] = 28
    p[5] = 0.04
    p[6] = 4.57
    p[7] = 25.12
    p[8] = 1.0
    """
    p[9] = 0.09
    p[10] = 2.5
    p[11] = 0.125
    p[12] = 0.009
    p[13] = 0.0068
    """
    p[9] = 2.0 # initial Treg stat
    p[10] = 0.5 # initial Thelp stat
    p[11] = 0.2 # initial NK stat
    p[12] = 0.2 # initial CD8 stat
    p[13:18] .= 0.001 # pSTAT Rates
    return p
end

"""Use this to get string of parameter names"""
function getParamNames()
    return [
        "kfwd"
        "k4"
        "k5"
        "k16"
        "k17"
        "k22"
        "k23"
        "k27"
        "pSTAT Treg"
        "pSTAT Thelp"
        "pSTAT NK"
        "pSTAT CD8"
        "k34"
        "k35"
        "k36"
        "k37"
        "k38"
        "k39"
    ]
end


text:
include("test_gage_array.jl")
# include("test_gage_array_of_arrays.jl")


text:
using Test
@isdefined(TestHelpers) || include("../../TestHelpers.jl")

TestHelpers.@timed_testset "trained_weights/" begin
    include("trained_weights/MNIST.n1.jl")
    include("trained_weights/MNIST.WK17a_linf0.1_authors.jl")
end


text:
using MultiThreadedCaches

using Test

@info "Testing with:" Threads.nthreads()

@testset "basics" begin
    cache = MultiThreadedCache{Int,Int}()
    init_cache!(cache)

    @test length(cache.thread_caches) == Threads.nthreads()

    get!(cache, 1) do
        return 10
    end

    # The second get!() shouldn't have any effect, and the result should still be 10.
    @test get!(cache, 1) do
        return 100
    end == 10

    # Even on different threads, the cache should still have a value for key `1`.
    @sync for _ in 1:100
        Threads.@spawn begin
            @test get!(cache, 1) do
                return 100
            end == 10
        end
    end

    # Internals test:
    @test length(cache.base_cache_futures) == 0
end
@testset "KV types" begin
    cache = MultiThreadedCache{Int,String}()
    init_cache!(cache)

    @test get!(()->"hi", cache, 1) == "hi"
    @test get!(()->"bye", cache, 1) == "hi"

    cache = MultiThreadedCache{Any,Any}()
    init_cache!(cache)

    @test get!(()->"hi", cache, 1) == "hi"
    @test get!(()->2.0, cache, 1) == "hi"
    @test get!(()->3.0, cache, 1.0) == "hi"
end

# Helper function for stress-test: returns true if all elements in iterable `v` are equal.
function all_equal(v)
    length(v) < 2 && return true
    e1 = first(v)
    i = 2
    @inbounds for i=2:length(v)
        v[i] == e1 || return false
    end
    return true
end

@testset "stress-test" begin
    cache = MultiThreadedCache{Int,Int}()
    init_cache!(cache)

    len = 100
    outputs = [Channel{Int}(Inf) for _ in 1:len]

    @sync for _ in 1:100_000
        Threads.@spawn begin
            key = rand(1:len)
            value = get!(cache, key) do
                # Randomize a value for (key,value) the *first time* it's requested
                rand(Int)
            end
            # Every subsequent get!() on the cache should see the exact same value,
            # even if requested in parallel from multiple threads.
            put!(outputs[key], value)
        end
    end

    for ch in outputs
        close(ch)
        @test all_equal(collect(ch))
    end

    # Internals test:
    @test length(cache.base_cache_futures) == 0
end

@testset "constructors" begin
    # Constructing with precomputed cache key/value pairs:
    cache = MultiThreadedCache{Int64, Int64}(Dict(2 => 3, 1 => 2))
    init_cache!(cache)

    # Existing keys use the cached value:
    @test get!(()->10, cache, 1) == 2
    @test get!(()->10, cache, 2) == 3
    # New keys use the newly provided value:
    @test get!(()->10, cache, 3) == 10
end

@testset "exceptions" begin
    cache = MultiThreadedCache{Int64, Int64}()
    init_cache!(cache)

    # Test that if the construction Task throws an exception, all Tasks see that exception:
    exceptions = Channel(Inf)
    @sync for i in 1:1000
        Threads.@spawn begin
            try
                get!(cache, 1) do
                    sleep(1)  # sleep, to give all Tasks time to spin up
                    throw(ArgumentError("$i"))
                end
            catch e
                @test e isa ArgumentError
                put!(exceptions, e)
            end
        end
    end

    # Test that all Tasks saw the same exception thrown.
    close(exceptions)
    @test all_equal(collect(exceptions))

    # Test that after throwing an exception during a get! function, the cache is still
    # robust and working as intended.

    @test get!(()->10, cache, 1) == 10
    @test get!(()->20, cache, 1) == 10

    # Internals test:
    # (despite the exception, all futures are closed, and the base cache only contains
    # expected values.)
    @test length(cache.base_cache_futures) == 0
    @test cache.base_cache == Dict(1=>10)
end

@testset "show" begin
    cache = MultiThreadedCache{Int64, Int64}()
    # Exercise both show functions
    show(cache)
    display(cache)
end

# ---------------------------------------------------------------
# What follows is a benchmark, testing that scaling the number of threads (increasing
# contention on the cache) doesn't negatively affect performance:

@testset "perf-benchmark" begin
    cache = MultiThreadedCache{Int,Int}()
    init_cache!(cache)

    len = 100
    outputs = [Channel{Int}(Inf) for _ in 1:len]

    run_serial(N) = @sync for _ in 1:N
        begin
            key = rand(1:len)
            value = get!(cache, key) do
                # Randomize a value for (key,value) the *first time* it's requested
                rand(Int)
            end
            # Every subsequent get!() on the cache should see the exact same value,
            # even if requested in parallel from multiple threads.
            put!(outputs[key], value)
        end
    end

    run_parallel(N) = @sync for _ in 1:N
        Threads.@spawn begin
            key = rand(1:len)
            value = get!(cache, key) do
                # Randomize a value for (key,value) the *first time* it's requested
                rand(Int)
            end
            # Every subsequent get!() on the cache should see the exact same value,
            # even if requested in parallel from multiple threads.
            put!(outputs[key], value)
        end
    end

    # Measure a baseline against a dict with a lock
    mutex = ReentrantLock()
    dict = Dict{Int,Int}()
    baseline_outputs = [Channel{Int}(Inf) for _ in 1:len]
    run_baseline(N) = @sync for _ in 1:N
        Threads.@spawn begin
            key = rand(1:len)
            Base.@lock mutex begin
                value = get!(dict, key) do
                    rand(Int)
                end
                put!(baseline_outputs[key], value)
            end
        end
    end

    # warmup
    run_serial(10); run_parallel(10); run_baseline(10)

    # measurement
    n = 100_000
    time_serial = @elapsed run_serial(n)
    time_parallel = @elapsed run_parallel(n)
    time_baseline = @elapsed run_baseline(n)

    for ch in outputs
        close(ch)
        @test all_equal(collect(ch))
    end
    for ch in baseline_outputs
        close(ch)
        @test all_equal(collect(ch))
    end

    # Performance test:
    @info "benchmark results" Threads.nthreads() time_serial time_parallel time_baseline
end


text:
module DoublyPeriodicTaylorGreen

using Printf

using Oceananigans, Oceananigans.OutputWriters

# Advected vortex: ψ(x, y, t) = exp(-2t) * cos(x - U*t) * cos(y)
u(x, y, t, U=1) = U + exp(-2t) * cos(x - U*t) * sin(y)
v(x, y, t, U=1) =   - exp(-2t) * sin(x - U*t) * cos(y)

#####
##### x, y
#####

function setup_simulation(; Nx, Δt, stop_iteration, U=1, architecture=CPU(), dir="data")

    grid = RegularCartesianGrid(size=(Nx, Nx, 1), x=(0, 2π), y=(0, 2π), z=(0, 1),
                                topology=(Periodic, Periodic, Bounded))

    model = IncompressibleModel(architecture = architecture,
                                        grid = grid,
                                    coriolis = nothing,
                                    buoyancy = nothing,
                                     tracers = nothing,
                                     closure = ConstantIsotropicDiffusivity(ν=1))

    set!(model, u = (x, y, z) -> u(x, y, 0, U),
                v = (x, y, z) -> v(x, y, 0, U))

    simulation = Simulation(model, Δt=Δt, stop_iteration=stop_iteration, progress_frequency=stop_iteration)

    simulation.output_writers[:fields] = JLD2OutputWriter(model, FieldOutputs(model.velocities);
                                                          dir = dir, force = true,
                                                          prefix = @sprintf("taylor_green_Nx%d_Δt%.1e", Nx, Δt),
                                                          interval = stop_iteration * Δt / 10)

    return simulation
end

function setup_and_run(; setup...)

    simulation = setup_simulation(; setup...)

    println("""
            Running decaying Taylor-Green vortex simulation in x, y with

                Nx = $(setup[:Nx])
                Δt = $(setup[:Δt])

            """)

    @time run!(simulation)
    return nothing
end

end # module


text:
"""
$(SIGNATURES)

Given 3 points and their velocities, calculate the time `t` required to bring the area of that triangle to zero
"""
function find_maximum_t(p1::Vertex3D, p2::Vertex3D, p3::Vertex3D)
    # Three points A, B and C, with coordinates (x, y)
    # The area of a triangle is 1/2 * [ Ax (By - Cy) + Bx (Cy - Ay) + Cx (Ay - By)]
    # where each point of the triangle is where it will be after time t
    # i.e. a point goes from P to P+tV where V is the velocity of that point.
    # Notation is here B -> B + t_vB

    # To make the calculation simpler, everything is translated so that A is at the
    # origin of the plane and its velocity is nil.
    Br = p2.r - p1.r
    Bc = p2.c - p1.c
    Cr = p3.r - p1.r
    Cc = p3.c - p1.c

    t_vBr = p2.vr - p1.vr
    t_vBc = p2.vc - p1.vc
    t_vCr = p3.vr - p1.vr
    t_vCc = p3.vc - p1.vc

    # After this, given that Ar = Ac = t_vAr = t_vAc = 0, the area is nil iff
    # (Br + t_vBr) (Cc + t_vCc ) - (Cr + t_vCr) (Bc + t_vBc) = 0.
    # After expansion and reshuffling to have a quadratic equation where t
    # is the variable, the coefficients of that equation are:
    a = t_vCc * t_vBr - t_vBc * t_vCr
    b = -Bc * t_vCr - Cr * t_vBc + Br * t_vCc + Cc * t_vBr
    c = Br * Cc - Cr * Bc

    # if a = 0, this is just a linear equation.
    if a == 0 && b != 0
        return smallest_positive(-c / b, c / b)
    else
        discriminant = b^2 - 4a * c

        # If there is a solution
        if discriminant >= 0
            d = sqrt(discriminant)
            return smallest_positive((-b - d) / 2a, (-b + d) / 2a)
        end
    end
    # There can be no solution if, after translation, B abd C move in parallel direction.
    # C will never end up on the line AB.
    # Very unlikely with Float64.
    # If nothing found, no move
    return nothing
end

find_maximum_t(p::Tuple{Vertex3D,Vertex3D,Vertex3D}) = find_maximum_t(p[1], p[2], p[3])

"""
$(SIGNATURES)

Create a mash of the solid block to be carved.
Imagining a mid-plane where the heightmap is carved, start with a block whose lower plane is
bottom_distance away and the highr plane is top_distance above.
The distances are measured in meters and converted to pixels when necessary.

The final mesh is a vector of vertices each with an ID number and a position.

"""
function create_solid_as_dict(
    mesh::FaceMesh;
    bottom_distance = Bottom_Offset,
    top_distance = Top_Offset,
)

    height, width = size(mesh)
    bottom_offset = bottom_distance / Meters_Per_Pixel
    top_offset = top_distance / Meters_Per_Pixel

    # To do this, we use several dictionaries before saving
    # Dictionary 1: Triangle to 3 strings each naming a triangle vertex.
    # Dictionary 2: Named vertex to its x, y, z coordinates.
    # Dictionary 3: Named vertex to a unique numerical index.
    triangle_dict = Dict{String,Vector{String}}()
    vertex_dict = Dict{String,Tuple{Float64,Float64,Float64}}()

    mesh_height, mesh_width = size(mesh.corners.ϕ)
    for row ∈ 1:mesh_height-1, col ∈ 1:mesh_width-1

        # Top face, add vertices, top and bottom triangles
        for side ∈ [:top, :bottom]
            side_name = side == :top ? "Top" : "Bottom"
            t = side == :top ? mesh.toptriangles[row, col] : mesh.bottriangles[row, col]

            # Add the vertices 1 by 1
            vertices_names = ["", "", ""]
            for i ∈ 1:3
                v_row, v_col = t[i]
                r = mesh.corners.r[v_row, v_col]
                c = mesh.corners.c[v_row, v_col]
                ϕ = mesh.corners.ϕ[v_row, v_col]
                name = vertex_name("Top", row, col)
                vertices_names[i] = name
                vertex_dict[name] = (r, c, ϕ)
            end

            # Add the triangle
            triangle_dict[triangle_name(side_name, row, col, side)] = vertices_names
        end

        # Bottom face
        v_coord_top = [(row, col), (row + 1, col), (row, col + 1)]
        vertices_names_top = ["", "", ""]

        v_coord_bot = [(row, col + 1), (row + 1, col), (row + 1, col + 1)]
        vertices_names_bot = ["", "", ""]

        # Add the vertices 1 by one
        for i ∈ 1:3
            (v_row, v_col) = v_coord_top[i]
            name = vertex_name("Bottom", v_row, v_col)
            vertices_names_top[i] = name
            vertex_dict[name] =
                (Float64(v_row), Float64(v_col), Float64(-bottom_offset * 0))

            (v_row, v_col) = v_coord_bot[i]
            name = vertex_name("Bottom", v_row, v_col)
            vertices_names_bot[i] = name
            vertex_dict[name] = (Float64(v_row), Float64(v_col), Float64(-bottom_offset))
        end

        # Add the triangles
        triangle_dict[triangle_name("Bottom", row, col, :top)] = vertices_names_top
        triangle_dict[triangle_name("Bottom", row, col, :bottom)] = vertices_names_bot
    end

    ###
    ### Side meshes
    ###
    # Build triangles to create side meshes and close the mesh
    # The mesh is made of a single couple top/bottom triangles joining the bottom face
    # and the top face.
    # All the vertices already exist

    # (Looking from above), west/east sides.
    for (col, face_name) ∈ [(1, "West"), (width, "East")]
        for row = 1:height-1

            v_coord_top = [(:top, row, col), (:bottom, row, col), (:top, row + 1, col)]
            vertices_names_top = ["", "", ""]
            for i ∈ 1:3
                vertices_names_top[i] = vertex_name(
                    (v_coord_top[i][1] == :top ? "Top" : "Bottom"),
                    v_coord_top[i][2],
                    v_coord_top[i][3],
                )
            end

            v_coord_bot =
                [(:top, row + 1, col), (:bottom, row + 1, col), (:bottom, row, col)]
            vertices_names_bot = ["", "", ""]
            for i ∈ 1:3
                vertices_names_bot[i] = vertex_name(
                    (v_coord_bot[i][1] == :top ? "Top" : "Bottom"),
                    v_coord_bot[i][2],
                    v_coord_bot[i][3],
                )
            end

            # Add the triangles
            triangle_dict[triangle_name(face_name, row, col, :top)] = vertices_names_top
            triangle_dict[triangle_name(face_name, row, col, :bottom)] = vertices_names_bot
        end
    end

    # (Looking from above), north/south sides.
    for (row, face_name) ∈ [(1, "North"), (height, "South")]
        for col ∈ 1:width-1

            v_coord_top = [(:top, row, col), (:bottom, row, col), (:top, row, col + 1)]
            vertices_names_top = ["", "", ""]
            for i ∈ 1:3
                vertices_names_top[i] = vertex_name(
                    (v_coord_top[i][1] == :top ? "Top" : "Bottom"),
                    v_coord_top[i][2],
                    v_coord_top[i][3],
                )
            end

            v_coord_bot = [(:top, row, col + 1), (:bottom, row, col + 1), (:top, row, col)]
            vertices_names_bot = ["", "", ""]
            for i ∈ 1:3
                vertices_names_bot[i] = vertex_name(
                    (v_coord_bot[i][1] == :top ? "Top" : "Bottom"),
                    v_coord_bot[i][2],
                    v_coord_bot[i][3],
                )
            end

            # Add the triangles
            triangle_dict[triangle_name(face_name, row, col, :top)] = vertices_names_top
            triangle_dict[triangle_name(face_name, row, col, :bottom)] = vertices_names_bot
        end
    end

    return triangle_dict, vertex_dict
end

"""
$(SIGNATURES)

"""
function create_solid(
    mesh::FaceMesh;
    bottom_distance = Bottom_Offset,
    top_distance = Top_Offset,
)

    triangle_dict, vertex_dict = create_solid_as_dict(
        mesh;
        bottom_distance = Bottom_Offset,
        top_distance = Top_Offset,
    )

    _, _, triangle_index, vertex_index = create_solid(
        triangle_dict,
        vertex_dict;
        bottom_distance = bottom_distance,
        top_distance = top_distance,
    )

    return triangle_dict, vertex_dict, triangle_index, vertex_index
end

"""
$(SIGNATURES)

"""
function create_solid(
    triangle_dict::Dict{String,Vector{String}},
    vertex_dict::Dict{String,Tuple{Float64,Float64,Float64}};
    bottom_distance = Bottom_Offset,
    top_distance = Top_Offset,
)

    index_dict = Dict{String,Int64}()
    index_of_vertex = 0
    for name ∈ keys(vertex_dict)
        index_of_vertex += 1
        index_dict[name] = index_of_vertex
    end

    vertex_index = Vector{Tuple{Float64,Float64,Float64}}(undef, length(vertex_dict))
    # For each vertex using its name and coordinates
    for (name, v) ∈ vertex_dict
        # Determine the index of the vertex
        i = index_dict[name]
        x = v[1]
        y = v[2]
        z = v[3]
        vertex_index[i] = (x, y, z)
    end

    # The connections are listed as a tuple of 3 arrays.
    # First array contains the first point of a triangle, and so on
    triangle_index = [
        (index_dict[v[1]], index_dict[v[2]], index_dict[v[3]]) for
        v ∈ values(triangle_dict)
    ]

    return triangle_dict, vertex_dict, triangle_index, vertex_index

end


text:
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

@testset "OpenEphysLoader.jl" begin
    include("util.jl")
    include("original.jl")
    include("continuous.jl")
    include("test_continuous.jl")
    include("metadata.jl")
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module FlatBuffers

# Compat
if !isdefined(Core, :String)
    typealias String UTF8String
end

if !isdefined(Base, :view)
    view = sub
end

if VERSION < v"0.5.0-dev+977"
    export foreach

    foreach(f) = (f(); nothing)
    foreach(f, itr) = (for x in itr; f(x); end; nothing)
    foreach(f, itrs...) = (for z in zip(itrs...); f(z...); end; nothing)
end

if VERSION < v"0.5.0-dev+4631"
    unsafe_wrap{A<:Array}(::Type{A}, ptr, len) = pointer_to_array(ptr, len)
end

# utils
immutable UndefinedType end
const Undefined = UndefinedType()

getfieldvalue{T}(obj::T, i) = isdefined(obj, i) ? getfield(obj, i) : Undefined

"""
     Scalar
A Union of the Julia types `T <: Number` that are allowed in FlatBuffers schema
"""
typealias Scalar Union{UndefinedType, Bool,
                        Int8, Int16, Int32, Int64,
                        UInt8, UInt16, UInt32, UInt64,
                        Float32, Float64}

isstruct(T) = !T.mutable && isleaftype(T)
isbitstype(T) = nfields(T) == 0

default(T, TT, sym) = default(TT)
default(::Type{UndefinedType}) = Undefined
default{T<:Scalar}(::Type{T}) = zero(T)
default{T<:AbstractString}(::Type{T}) = ""
default{T<:Enum}(::Type{T}) = enumtype(T)(T(0))
default{T}(::Type{Vector{T}}) = T[]
# fallback that recursively builds a default; for structs/tables
default{T}(::Type{T}) = isa(T, Union) ? nothing : T(map(TT->TT == T ? TT() : default(TT),T.types)...)

function typeorder end

enumtype{T<:Enum}(::Type{T}) = UInt8

# Types
"""
    Table

The object containing the flatbuffer and positional information specific to the table.
The `vtable` containing the offsets for specific members precedes `pos`.
The actual values in the table follow `pos` offset and size of the vtable.

- `bytes::Vector{UInt8}`: the flatbuffer itself
- `pos::Int`:  the base position in `bytes` of the table
"""
type Table{T}
    bytes::Vector{UInt8}
    pos::Int
end

"""
Builder is a state machine for creating FlatBuffer objects.
Use a Builder to construct object(s) starting from leaf nodes.

A Builder constructs byte buffers in a last-first manner for simplicity and
performance.
"""
type Builder{T}
    bytes::Vector{UInt8}
    minalign::Int
    vtable::Vector{Int}
	objectend::Int
	vtables::Vector{Int}
	head::Int
	nested::Bool
	finished::Bool
end

function Base.show{T}(io::IO, x::Union{Builder{T},Table{T}})
    println(io, "FlatBuffers.$(typeof(x)): ")
    buffer = typeof(x) <: Table ? x.bytes : x.bytes[x.head+1:end]
    if isempty(buffer)
        print(io, " (empty flatbuffer)")
    else
        pos = Int(typeof(x) <: Table ? x.pos : get(buffer, 0, Int32))
        # print vtable offset
        syms = T.name.names
        maxpad = max(length(" vtable rel. start pos: "), maximum(map(x->length(string(x)), syms)))
        stringify(buf, x, y, msg) = replace(string(rpad(string(lpad("$(x): ", 6, ' '),lpad(msg, maxpad, ' ')),maxpad+6,' '),string(map(z->lpad(string(Int(z)), 4, ' '),buf[x:y]))[9:end-1]),'"',"")
        println(io, stringify(buffer, 1, 4, " root position: "))
        vtaboff = get(buffer, pos, Int32)
        vtabstart = pos - vtaboff + 5

        println(io, stringify(buffer, 5, 6, " vtable size: "))
        println(io, stringify(buffer, 7, 8, " data size: "))
        i = vtabstart
        x = 1
        for y = 1:length(syms)
            println(io, stringify(buffer, i, i+1, "$(syms[x]): "))
            i += 2
            x += 1
        end
        # print rel pos. of vtable
        println(io, stringify(buffer, i, i+3, " vtable rel. start pos: "))
        i += 4
        # now we're pointing at data
        while i < length(buffer)
            println(io, stringify(buffer, i, i+3, " "))
            i += 4
        end
    end
end

include("internals.jl")

function Table{T}(::Type{T}, buffer::Vector{UInt8}, pos::Integer)
    return Table{T}(buffer, pos)
end

Table{T}(b::Builder{T}) = Table(T, b.bytes[b.head+1:end], get(b, b.head, Int32))

function untilindex(func, itr)
    for (i,x) in enumerate(itr)
        func(x) && return i
    end
    return 0
end

getvalue(t, o, ::Type{Void}) = nothing
getvalue{T<:Scalar}(t, o, ::Type{T}) = get(t, t.pos + o, T)
getvalue{T<:Enum}(t, o, ::Type{T}) = T(get(t, t.pos + o, enumtype(T)))
function getvalue{T<:AbstractString}(t, o, ::Type{T})
    o += get(t, t.pos + o, Int32)
    strlen = get(t, t.pos + o, Int32)
    o += t.pos + sizeof(Int32)
    return String(t.bytes[o + 1:o + strlen])
end
function getvalue(t, o, ::Type{Vector{UInt8}})
    o += get(t, t.pos + o, Int32)
    len = get(t, t.pos + o, Int32)
    o += t.pos + sizeof(Int32)
    return t.bytes[o + 1:o + len] #TODO: maybe not make copy here?
end

getarray{T<:Scalar}(t, vp, len, ::Type{T}) = unsafe_wrap(Array, convert(Ptr{T}, pointer(view(t.bytes, (vp + 1):length(t.bytes)))), len)
getarray{T<:Enum}(t, vp, len, ::Type{T}) = convert(Vector{T}, unsafe_wrap(Array, convert(Ptr{enumtype(T)}, pointer(view(t.bytes, (vp + 1):length(t.bytes)))), len))
function getarray{T<:Union{AbstractString,Vector{UInt8}}}(t, vp, len, ::Type{T})
    A = Vector{T}(len)
    for i = 1:len
        A[i] = getvalue(t, vp - t.pos, T)
        vp += sizeof(Int32)
    end
    return A
end
function getarray{T}(t, vp, len, ::Type{T})
    if isstruct(T)
        return unsafe_wrap(Array, convert(Ptr{T}, pointer(view(t.bytes, (vp + 1):length(t.bytes)))), len)
    else
        A = Vector{T}(len)
        for i = 1:len
            A[i] = getvalue(t, vp - t.pos, T)
            vp += sizeof(Int32)
        end
        return A
    end
end

function getvalue{T}(t, o, ::Type{Vector{T}})
    vl = vectorlen(t, o)
    vp = vector(t, o)
    return getarray(t, vp, vl, T)
end

# fallback which recursively calls read
function getvalue{T}(t, o, ::Type{T})
    if isstruct(T)
        if any(x-> x <: Enum, T.types)
            args = []
            o = t.pos + o + 1
            for typ in T.types
                val = unsafe_load(convert(Ptr{typ <: Enum ? enumtype(typ) : typ}, pointer(view(t.bytes, o:length(t.bytes)))))
                push!(args, val)
                o += sizeof(typ <: Enum ? enumtype(typ) : typ)
            end
            return T(args...)
        else
            return unsafe_load(convert(Ptr{T}, pointer(view(t.bytes, (t.pos + o + 1):length(t.bytes)))))
        end
    else
        newt = Table{T}(t.bytes, t.pos + o + get(t, t.pos + o, Int32))
        return FlatBuffers.read(newt, T)
    end
end

"""
`FlatBuffers.read` parses a `T` at `t.pos` in Table `t`.
Will recurse as necessary for nested types (Arrays, Tables, etc.)
"""
function FlatBuffers.read{T1,T}(t::Table{T1}, ::Type{T}=T1)
    args = []
    numfields = length(T.types)
    for i = 1:numfields
        TT = T.types[i]
        # if it's a Union type, use the previous arg to figure out the true type that was serialized
        if isa(TT, Union)
            TT = typeorder(TT, args[end])
        end
        o = offset(t, 4 + ((i - 1) * 2))
        if o == 0
            push!(args, default(T, TT, T.name.names[i]))
        else
            push!(args, getvalue(t, o, TT))
        end
    end
    return T(args...)
end

FlatBuffers.read{T}(::Type{T}, buffer::Vector{UInt8}, pos::Integer) = FlatBuffers.read(Table(T, buffer, pos))
FlatBuffers.read{T}(b::Builder{T}) = FlatBuffers.read(Table(T, b.bytes[b.head+1:end], get(b, b.head, Int32)))

function Builder{T}(::Type{T}=Any, size=0)
    objectend = 0
    vtables = zeros(Int, 0)
    head = size
    nested = false
    bytes = zeros(UInt8, size)
    minalign = 1
    vtable = zeros(Int, 0)
    finished = false
    b = Builder{T}(bytes, minalign, vtable, objectend,
                vtables, head, nested, finished)
    return b
end

# build!
"`alignment` looks for the largest scalar member of `T` that represents a flatbuffer Struct"
function alignment{T}(::Type{T})
    largest = 0
    for typ in T.types
        largest = isbitstype(typ) ? max(largest,sizeof(typ)) : alignment(typ)
    end
    return largest
end

"""
`buildvector!` is for building vectors with all kinds of element types,
even building its elements recursively if needed (Array of Arrays, Array of tables, etc.).
"""
function buildvector! end

# empty vector
function buildvector!(b, A::Vector{Void}, len)
    startvector(b, 1, 0, 1)
    return endvector(b, 0)
end
# scalar type vector
function buildvector!{T<:Scalar}(b, A::Vector{T}, len)
    startvector(b, sizeof(T), len, sizeof(T))
    foreach(x->prepend!(b, A[x]), len:-1:1)
    return endvector(b, len)
end
function buildvector!{T<:Enum}(b, A::Vector{T}, len)
    startvector(b, sizeof(enumtype(T)), len, sizeof(enumtype(T)))
    foreach(x->prepend!(b, enumtype(T)(A[x])), len:-1:1)
    return endvector(b, len)
end

function putoffsetvector!(b, offsets, len)
    startvector(b, 4, len, 4) #TODO: elsize/alignment correct here?
    foreach(x->prependoffset!(b, offsets[x]), len:-1:1)
    return endvector(b, len)
end
# byte vector vector
function buildvector!(b, A::Vector{Vector{UInt8}}, len)
    offsets = map(x->createbytevector(b, A[x]), 1:len)
    return putoffsetvector!(b, offsets, len)
end
# string vector
function buildvector!{T<:AbstractString}(b, A::Vector{T}, len)
    offsets = map(x->createstring(b, A[x]), 1:len)
    return putoffsetvector!(b, offsets, len)
end
# array vector
function buildvector!{T}(b, A::Vector{Vector{T}}, len)
    offsets = map(x->buildbuffer!(b, A[x]), 1:len)
    return putoffsetvector!(b, offsets, len)
end
# struct or table/object vector
function buildvector!{T}(b, A::Vector{T}, len)
    if isstruct(T)
        # struct
        startvector(b, sizeof(T), len, alignment(T)) #TODO: forced elsize/alignment correct here?
        foreach(x->buildbuffer!(b, A[x]), len:-1:1)
        return endvector(b, len)
    else # table/object
        offsets = map(x->buildbuffer!(b, A[x]), 1:len)
        return putoffsetvector!(b, offsets, len)
    end
end

"""
`getoffset` checks if a given field argument needs to be built
offset (Arrays, Strings, other tables) or can be inlined (Scalar or Struct types).
`getoffset` has a recursive nature in that it will build offset types
down to their last leaf scalar types before returning the highest-level offset.
"""
function getoffset end

getoffset(b, arg::Void) = 0
getoffset{T<:Scalar}(b, arg::T) = 0
getoffset{T<:Enum}(b, arg::T) = 0
getoffset(b, arg::Vector{UInt8}) = createbytevector(b, arg)
getoffset(b, arg::AbstractString) = createstring(b, arg)
getoffset{T}(b, arg::Vector{T}) = buildbuffer!(b, arg)
# structs or table/object
getoffset{T}(b, arg::T) = isstruct(T) ? 0 : buildbuffer!(b, arg)

"""
`putslot!` is one of the final steps in building a flatbuffer.
It puts the final value in the "data" section of the flatbuffer,
whether that be an actual value (for `Scalar`, Struct types) or an offset
to the actual data (Arrays, Strings, other tables)
"""
function putslot! end

putslot!{T<:Scalar}(b, i, arg::T, off) = prependslot!(b, i, arg, default(T))
putslot!{T<:Enum}(b, i, arg::T, off) = prependslot!(b, i, enumtype(T)(arg), default(T))
putslot!(b, i, arg::AbstractString, off) = prependoffsetslot!(b, i, off, 0)
putslot!{T}(b, i, arg::Vector{T}, off) = prependoffsetslot!(b, i, off, 0)
# structs or table/object
putslot!{T}(b, i, arg::T, off) =
    isstruct(T) ? prependstructslot!(b, i, buildbuffer!(b, arg), 0) : prependoffsetslot!(b, i, off, 0)

function buildbuffer!{T1,T}(b::Builder{T1}, arg::T)
    if T <: Array
        # array of things
        n = buildvector!(b, arg, length(arg))
    elseif isstruct(T)
        # build a struct type with provided `arg`
        all(isstruct, T.types) || throw(ArgumentError("can't seralize flatbuffer, $T is not a pure struct"))
        align = alignment(T)
        prep!(b, align, 2align)
        for i = length(T.types):-1:1
            typ = T.types[i]
            if typ <: Enum
                prepend!(b, enumtype(typ)(getfield(arg,i)))
            elseif isbitstype(typ)
                prepend!(b, getfield(arg,i))
            else
                buildbuffer!(b, getfield(arg,i))
            end
        end
        n = offset(b)
    else
        # build a table type
        # check for string/array/table types
        numfields = length(T.types)
        offsets = [getoffset(b, getfieldvalue(arg,i)) for i = 1:numfields]

        # all nested have been written, with offsets in `offsets[]`
        startobject(b, numfields)
        foreach(i->putslot!(b, i, getfieldvalue(arg,i), offsets[i]), 1:numfields)
        n = endobject(b)
    end
    return n
end

function build!(b, arg)
    n = buildbuffer!(b, arg)
    finish!(b, n)
    return b
end

build!{T}(arg::T) = build!(Builder(T), arg)

include("macros.jl")

end # module


text:
const CLVSolver = StagedSolver{<: CLVProblem, <: CLVSolution}

"""
    CLVSolver(prob::CLVProblem; <keyword arguments>)
    CLVSolver(prob::LEProblem; <keyword arguments>)

The preferred and equivalent method to get a solver for a `CLVProblem`
is `init(prob::CLVProblem)`.  Note that `CLVSolver(prob::LEProblem)`
is equivalent to `init(CLVProblem(prob))`.

# Arguments
- `record::Vector{Symbol}`: Variables to be saved.  A subset of
  `[:G, :C, :D, :x]`.
"""
function CLVSolver(prob::CLVProblem;
                   record::Vector{Symbol} = Symbol[],
                   forward_relaxer::Type = ForwardRelaxer,
                   forward_dynamics::Type = ForwardDynamics,
                   backward_relaxer::Type = BackwardRelaxer,
                   backward_dynamics::Type = if :D in record
                       BackwardDynamicsWithD
                   else
                       BackwardDynamics
                   end,
                   )
    stage_types = [
        # PhaseRelaxer???
        forward_relaxer,
        forward_dynamics,
        backward_relaxer,
        backward_dynamics,
    ]
    return CLVSolver(prob, stage_types, record)
end

function CLVSolver(prob::CLVProblem, stage_types::AbstractVector,
                   record::Vector{Symbol})
    sol = CLVSolution(prob, record)
    args = (prob, sol)  # additional arguments to each `stage_types`
    # FIXME: Currently sol is used as "cache" (args[2]) as well.
    return StagedSolver(prob, sol, stage_types, args)
end

CLVSolver(prob::LEProblem; kwargs...) = CLVSolver(CLVProblem(prob); kwargs...)

init(prob::CLVProblem; kwargs...) = CLVSolver(prob; kwargs...)

"""
    forward_dynamics!(solver::CLVSolver) :: ForwardDynamics

Solve the CLV problem up to the forward dynamics stage and return an
iterator to step through the forward dynamics.
See also: [`backward_dynamics!`](@ref), [`goto!`](@ref).
"""
forward_dynamics!(solver::CLVSolver; kwargs...) =
    goto!(solver, ForwardDynamics; kwargs...)

"""
    indexed_forward_dynamics!(solver::CLVSolver)
    indexed_forward_dynamics!(stage::ForwardDynamics)

Just a short-hand for `enumerate(forward_dynamics!(solver))`.
It's for symmetry with [`indexed_backward_dynamics!`](@ref).
"""
indexed_forward_dynamics!(solver::CLVSolver; kwargs...) =
    indexed_forward_dynamics!(forward_dynamics!(solver; kwargs...))
indexed_forward_dynamics!(fitr::ForwardDynamics) = enumerate(fitr)

"""
    backward_dynamics!(solver::CLVSolver) :: BackwardDynamics

Solve the CLV problem up to the (final) backward dynamics stage and
return an iterator to step through the backward dynamics.
See also [`forward_dynamics!`](@ref), [`goto!`](@ref).

Note that finishing iteration of the returned iterator does not
finalize all the solver stages (namely, recording to `solver.sol`, if
non-default `backward_dynamics` is used).  In this case,
`solve!(solver)` has to be called after the iteration.

### Example
```julia
angles = [acos(abs(dot(C[:, 1], C[:, 2]))) * 2 / π for C
          in backward_dynamics!(solver)]
```
"""
backward_dynamics!(solver::CLVSolver; kwargs...) =
    goto!(solver, BackwardDynamics; kwargs...)

"""
    indexed_backward_dynamics!(solver::CLVSolver)
    indexed_backward_dynamics!(stage::BackwardDynamics)

It is equivalent to `zip(some_counter, backward_dynamics!(solver))`
where `some_counter` is an iterator over integers such that the same
indices returned by [`indexed_forward_dynamics!`](@ref) indicate that
those events are at the same time point of the backward and forward
passes.  This is useful when combining matrices [`CLV.G`](@ref) and
[`CLV.C`](@ref) to obtain the CLV in the (original) tangent space.

For such example, see:
[Covariant Lyapunov vectors on the Lorenz attractor](@ref)
in the online manual.

Note that `indexed_backward_dynamics!(solver::CLVSolver)` is
equivalent to
```julia
backward = backward_dynamics!(solver)
indexed_backward_dynamics!(backward)
```
Separately calling [`backward_dynamics!`](@ref) is useful when the
quantities other than [`CLV.G`](@ref) (e.g., [`CLV.R`](@ref)) are
required.

See also [`indexed_forward_dynamics!`](@ref).
"""
indexed_backward_dynamics!(solver::CLVSolver;
                           until::Int = 1,
                           kwargs...) =
    indexed_backward_dynamics!(backward_dynamics!(solver; kwargs...);
                               until = until)

function indexed_backward_dynamics!(bitr::BackwardDynamics;
                                    until::Int = 1)
    return zip(length(bitr)-1:-1:until, bitr)
end


text:
@doc raw"""
Sellmeier formula:
```math
n^2 - 1 = c_1 + \sum \limits_{i=1}^{N} \frac{c_{2i} \lambda^2}{\lambda^2 - c_{2i+1}^2}
```
"""
struct Sellmeier{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Sellmeier{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] * λ^2 / (λ^2 - c[i + 1]^2)
    end
    return sqrt(rhs + 1)
end

@doc raw"""
Sellmeier-2 formula:
```math
n^2 - 1 = c_1 + \sum \limits_{i=1}^{N} \frac{c_{2i} \lambda^2}{\lambda^2 - c_{2i+1}}
```
"""
struct Sellmeier2{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Sellmeier2{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] * λ^2 / (λ^2 - c[i + 1])
    end
    return sqrt(rhs + 1)
end

@doc raw"""
Polynomial formula:
```math
n^2 = c_1 + \sum \limits_{i=1}^{N}  c_{2i} \lambda^{c_{2i+1}}
```
"""
struct Polynomial{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Polynomial{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] * λ^c[i + 1]
    end
    return sqrt(rhs)
end

@doc raw"""
RefractiveIndex.INFO formula:
```math
n^2 = c_1
    + \frac{c_{2} \lambda^{c_{3}}}{\lambda^2 - c_{4}^{c_{5}}}
    + \frac{c_{6} \lambda^{c_{7}}}{\lambda^2 - c_{8}^{c_{9}}}
    + \sum \limits_{i=5}^{N} c_{2i} \lambda^{c_{2i+1}}
```
"""
struct RIInfo{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::RIInfo{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:4:min(N, 9)
        rhs += (c[i] * λ^c[i + 1]) / (λ^2 - c[i + 2]^c[i + 3])
    end
    for i = 10:2:N
        rhs += c[i] * λ^c[i + 1]
    end
    return sqrt(rhs)
end

@doc raw"""
Cauchy formula:
```math
n = c_1 + \sum \limits_{i=5}^{N} c_{2i} \lambda^{c_{2i+1}}
```
"""
struct Cauchy{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Cauchy{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] * λ^c[i + 1]
    end
    return rhs
end

@doc raw"""
Gases formula:
```math
n - 1 = c_1
    + \sum \limits_{i=1}^{N} \frac{c_{2i}}{c_{2i+1} - \lambda^{-2} }
```
"""
struct Gases{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Gases{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] / (c[i + 1] - 1 / λ^2)
    end
    return rhs + 1
end

@doc raw"""
Herzberger formula:
```math
n = c_1
    + \frac{c_{2}}{\lambda^2 - 0.028}
    + c_{3} \left( \frac{1}{\lambda^2 - 0.028} \right)^2
    + \sum \limits_{i=1}^{N} c_{i+3} \lambda^{2i}
```
"""
struct Herzberger{N} <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{N,Float64}
end

function (f::Herzberger{N})(λ) where {N}
    c = f.coeffs
    rhs = c[1]
    rhs += c[2] / (λ^2 - 0.028)
    rhs += c[3] * (1 / (λ^2 - 0.028))^2
    for i = 4:N
        pow = 2 * (i - 3)
        rhs += c[i] * λ^pow
    end
    return rhs
end

@doc raw"""
Retro formula:
```math
\frac{n^2 - 1}{n^2 +2} = c_1
    + \frac{c_{2} \lambda^2}{\lambda^2 - c_{3}}
    +  c_{4} \lambda^2
```
"""
struct Retro <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{4,Float64}
end

function (f::Retro)(λ)
    c = f.coeffs
    rhs = c[1] + c[2] * λ^2 / (λ^2 - c[3]) + c[4] * λ^2
    return sqrt((-2rhs - 1) / (rhs - 1))
end

@doc raw"""
Exotic formula:
```math
n^2 = c_1
    + \frac{c_{2}}{\lambda^2 - c_{3}}
    + \frac{c_{4}(\lambda - c_{5})}{(\lambda - c_{5})^2 + c_{6}}
```
"""
struct Exotic <: Formula
    λrange::NTuple{2,Float64}
    coeffs::NTuple{6,Float64}
end

function (f::Exotic)(λ)
    c = f.coeffs
    rhs = c[1] + c[2] / (λ^2 - c[3]) + c[4] * (λ - c[5]) / ((λ - c[5])^2 + c[6])
    return sqrt(rhs)
end


text:
using EGraphs
using Test
#=
Write invariant checking functions

every node in EClass find_root! is the EClass

=#

in_right_class(e::EGraph) = all(   all( find_root!(e,e.memo[n]) == k for n in cls.nodes)  for (k,cls) in e.classes )
memo_canonical(e::EGraph) = all(   all( find_root!(e,a) == id for a in n.args)  for (n,id) in e.memo )

@testset "Atomic equivalence" begin
    # Write your tests here.
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    @test length(e.classes) == 1
    a = constant!(e, :a)
    @test length(e.classes) == 1
    b = constant!(e, :b)
    @test length(e.classes) == 2
    c = constant!(e, :c)
    @test length(e.classes) == 3
    union!(e, a, b)
    @test length(e.classes) == 2
    union!(e, a, c)
    @test length(e.classes) == 1

end

@testset "Simple congruence" begin
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    f = term!(e, :f)
    apply(n, f, x) = n == 0 ? x : apply(n-1,f,f(x))

    union!(e, a, apply(6,f,a))
    union!(e, a, apply(9,f,a))
    println(e.memo)
    # @test memo_canonical(e) # yeah. the parents still need to be fixed right? And that's ok?
    rebuild!(e)
    println(e.memo)
    println(e.classes)
    println(find_root!(e, Id(6)))
    @test in_right_class(e)
    @test length(e.classes) == 3
    @test length(e.memo) == 4

    union!(e, a, apply(11,f,a))
    rebuild!(e)
    @test in_right_class(e)
    @test memo_canonical(e)
    @test length(e.classes) == 1
    @test length(e.memo) == 2

end

EMPTY_DICT2 = Base.ImmutableDict{PatVar, Id}(PatVar(:____),  Id(-1))

@testset "Simple match" begin
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    f = term!(e, :f)
    p = PatVar

    for sub in EGraphs.ematch(e , PatVar(:x) , a,  EMPTY_DICT2 )
        #println(sub)
    end
    #f(f(f(f(a))))
    for sub in EGraphs.ematch(e , PatTerm(:f, [PatVar(:x)]) , f(f(f(a))),  EMPTY_DICT2 )
        #println(sub)
    end
    union!(e, a, f(a))
    rebuild!(e)
    @test length(e.classes) == 1
    #@test length(e.memo) == 2
    println(e.memo) # what? This dict has two entries for the same thing? Could this be my mutation problem rearing it's head?
    @test length(e.classes[Id(1)].parents) == 1
    f(f(f(a)))
    println(e.memo)
    println(e.classes) # The class has duplicates in it's node array
    for sub in EGraphs.ematch(e , PatTerm(:f, [PatVar(:x)]) , f(f(f(a))),  EMPTY_DICT2 )
        println(sub)
    end

end

@testset "Assoc1" begin
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    b = constant!(e, :b)
    c = constant!(e, :c)
    f = term!(e, :f)

    f(f(a,b),c)
    graphplot(e)
    x = PatVar(:x)
    y = PatVar(:y)
    z = PatVar(:z)
    pf = (x,y) -> PatTerm(:f, [x,y])
    r = Rule(  pf(pf(x,y),z) ,     pf(x,pf(y,z)) )

    @test !in_same_set(e,     f(f(a,b),c),     f(a, f(b,c)))
    rewrite!(e, r)
    @test in_same_set(e,     f(f(a,b),c),     f(a, f(b,c)))
    rewrite!(e, r) # shouldn't do anything extra
    graphplot(e)

end

@testset "Assoc2" begin
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    b = constant!(e, :b)
    c = constant!(e, :c)
    d = constant!(e, :d)
    #h = constant!(e, :h)
    f = term!(e, :f)

    t1 = f(f(f(a,b),c),d)
    t2 = f(a,f(b,f(c,d)))
    graphplot(e)
    x = PatVar(:x)
    y = PatVar(:y)
    z = PatVar(:z)
    pf = (x,y) -> PatTerm(:f, [x,y])
    r = Rule(  pf(pf(x,y),z) ,     pf(x,pf(y,z)) )

    @test !in_same_set(e,    t1, t2)
    rewrite!(e, r)
    graphplot(e)
    @test ! in_same_set(e,     t1, t2)
    rewrite!(e, r)
    @test in_same_set(e,     t1, t2)
    graphplot(e)

end

#=

So we should also write the extraction functions to get good terms out.

=#

@testset "Cat1" begin
    e = EGraphs.EGraph()
    a = constant!(e, :a)
    #b = constant!(e, :b)
    #c = constant!(e, :c)
    f = constant!(e, :f)

    id = term!(e, :id) # allow optionally to specify the arity
    comp = term!(e, :∘)
    swap = term!(e, :σ)

    pid = ob -> PatTerm(:id, [ob])
    pcomp = (f,g) -> PatTerm(:∘, [f,g])
    x = PatVar(:x)
    pf = PatVar(:f)

    r1 = Rule( pcomp( pf, pid(x) ), pf)
    r2 = Rule( pcomp( pid(x), pf ), pf)
    r3 = Rule( pf, pcomp( pf, pid(x) ))
    comp(f,id(a))
    graphplot(e)
    rewrite!(e, r1)
    graphplot(e)
    rewrite!(e, r3)
    graphplot(e)

end


text:
# ---
# title: 835. Image Overlap
# id: problem835
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/image-overlap/description/>
# hidden: true
# ---
#
# You are given two images `img1` and `img2` both of size `n x n`, represented
# as binary, square matrices of the same size. (A binary matrix has only 0s and
# 1s as values.)
#
# We translate one image however we choose (sliding it left, right, up, or down
# any number of units), and place it on top of the other image.  After, the
# _overlap_ of this translation is the number of positions that have a 1 in both
# images.
#
# (Note also that a translation does **not** include any kind of rotation.)
#
# What is the largest possible overlap?
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/09/09/overlap1.jpg)
#
#
#
#     Input: img1 = [[1,1,0],[0,1,0],[0,1,0]], img2 = [[0,0,0],[0,1,1],[0,0,1]]
#     Output: 3
#     Explanation: We slide img1 to right by 1 unit and down by 1 unit.
#     ![](https://assets.leetcode.com/uploads/2020/09/09/overlap_step1.jpg)
#     The number of positions that have a 1 in both images is 3. (Shown in red)
#     ![](https://assets.leetcode.com/uploads/2020/09/09/overlap_step2.jpg)
#
#
# **Example 2:**
#
#
#
#     Input: img1 = [[1]], img2 = [[1]]
#     Output: 1
#
#
# **Example 3:**
#
#
#
#     Input: img1 = [[0]], img2 = [[0]]
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `n == img1.length`
#   * `n == img1[i].length`
#   * `n == img2.length `
#   * `n == img2[i].length`
#   * `1 <= n <= 30`
#   * `img1[i][j]` is `0` or `1`.
#   * `img2[i][j]` is `0` or `1`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
# https://raw.githubusercontent.com/sswatson/cheatsheets/master/plotsjl-cheatsheet.pdf

using Logging
logger = SimpleLogger(stdout, Logging.Info)
global_logger(logger)

using Plots
plotly()

using Polynomials

pt = plot(title="Learn how to plot", size=(600,600))
f = Polynomial([0,0,1])
plot!(pt, f, -2, 2)

gui(pt)

nothing


text:
using LinearAlgebra
using Random

Random.seed!(0)

function gradient_fd(f,x) #compute gradient using finite differences
	gradfd = zero(x)
	xeps1 = zero(x)
	xeps2 = zero(x)
	delta = sqrt(eps())
	for i in eachindex(gradfd)
		xeps1 .= x
		xeps2 .= x
		xeps1[i] -= delta
		xeps2[i] += delta
		gradfd[i] = (f(xeps2)-f(xeps1))/(2*delta)
	end
	return gradfd
end

stuff = [
  Dict( "f"      => ElasticNet(2.0, 3.0),
        "x"      => [-2., -1., 0., 1., 2., 3.],
        "∇f(x)"  => 2*[-1., -1., 0., 1., 1., 1.] + 3*[-2., -1., 0., 1., 2., 3.],
      ),
  Dict( "f"      => NormL1(2.0),
        "x"      => [-2., -1., 0., 1., 2., 3.],
        "∇f(x)"  => 2*[-1., -1., 0., 1., 1., 1.],
      ),
  Dict( "f"      => NormL2(2.0),
        "x"      => [-1., 0., 2.],
        "∇f(x)"  => 2/sqrt(5)*[-1., 0., 2.],
      ),
  Dict( "f"      => NormL2(2.0),
        "x"      => [-1., 0., 2.],
        "∇f(x)"  => 2/sqrt(5)*[-1., 0., 2.],
      ),
  Dict( "f"      => NormL2(2.0),
        "x"      => [0., 0., 0.],
        "∇f(x)"  => [0., 0., 0.],
      ),
  # Dict( "f"      => NormL21(),
  #       "x"      => ,
  #       "∇f(x)"  => ,
  #     ),
  Dict( "f"      => NormLinf(2.),
        "x"      => [-2., -1., 0., 1., 2., 3.],
        "∇f(x)"  => 2*[0., 0., 0., 0., 0., 1.],
      ),
  Dict( "f"      => NormLinf(2.),
        "x"      => [-4., -1., 0., 1., 2., 3.],
        "∇f(x)"  => 2*[-1., 0., 0., 0., 0., 0.],
      ),
  # Dict( "f"      => NuclearNorm(),
  #       "x"      => ,
  #       "∇f(x)"  => ,
  #     ),
  Dict( "f"      => SqrNormL2(2.0),
        "x"      => [-2., -1., 0., 1., 2., 3.],
        "∇f(x)"  => 2*[-2., -1., 0., 1., 2., 3.],
      ),
  Dict( "f"      => HingeLoss([1., 2., 1., 2., 1.], 2.0),
        "x"      => [-2., -1., 0., 2., 3.],
        "∇f(x)"  => -2*[1., 2., 1.,  0., 0.],
      ),
  Dict( "f"      => HuberLoss(2., 3.),
        "x"      => [-1., 0.5],
        "∇f(x)"  => 3*[-1., 0.5],
      ),
  Dict( "f"      => HuberLoss(2., 3.),
        "x"      => [-2., 1.5],
        "∇f(x)"  => [-4.8, 3.6], # 3*2*[-2., 1.5]/norm([-2., 1.5]),
      ),
  Dict( "f"      => Linear([-1.,2.,3.]),
        "x"      => [1., 5., 3.14],
        "∇f(x)"  => [-1.,2.,3.],
      ),
  Dict( "f"      => LogBarrier(2.0, 1.0, 1.5),
        "x"      => [1.0, 2.0, 3.0],
        "∇f(x)"  => -1.5*2.0./(2.0*[1.0, 2.0, 3.0].+1.0),# -μ*a*/(a*x+b)
      ),
  Dict( "f"      => LeastSquares([1.0 2.0; 3.0 4.0], [0.1, 0.2], 2.0, iterative=false),
        "x"      => [-1., 2.],
        "∇f(x)"  => [34.6, 50.], #λ*(A'A*x-A'b),
      ),
  Dict( "f"      => LeastSquares([1.0 2.0; 3.0 4.0], [0.1, 0.2], 2.0, iterative=true),
        "x"      => [-1., 2.],
        "∇f(x)"  => [34.6, 50.], #λ*(A'A*x-A'b),
      ),
  Dict( "f"      => Quadratic([2. -1.; -1. 2.], [0.1, 0.2], iterative=false),
        "x"      => [3., 4.],
        "∇f(x)"  => [2.1, 5.2], #[2. -1.; -1. 2.]*[3., 4.]+[0.1, 0.2],
      ),
  Dict( "f"      => Quadratic([2. -1.; -1. 2.], [0.1, 0.2], iterative=true),
        "x"      => [3., 4.],
        "∇f(x)"  => [2.1, 5.2], #[2. -1.; -1. 2.]*[3., 4.]+[0.1, 0.2],
      ),
  Dict( "f"      => SumPositive(),
        "x"      => [-1., 1., 2.],
        "∇f(x)"  => [0., 1., 1.],
      ),
  # Dict( "f"      => Maximum(2.0),
  #       "x"      => [-4., 2., 3.],
  #       "∇f(x)"  => 2*[0., 0., 1.],
  #     ),
  Dict( "f"      => DistL2(IndZero(),2.0),
        "x"      => [1., -2],
        "∇f(x)"  => 2*[1., -2]./norm([1., -2]),
      ),
  Dict( "f"      => DistL2(IndBallL2(1.0),2.0),
        "x"      => [2., -2],
        "∇f(x)"  => [sqrt(2), -sqrt(2)],
      ),
  Dict( "f"      => SqrDistL2(IndZero(),2.0),
        "x"      => [1., -2],
        "∇f(x)"  => 2*[1., -2],
      ),
  Dict( "f"      => SqrDistL2(IndBallL2(1.0),2.0),
        "x"      => [2., -2],
        "∇f(x)"  => 2.0*([2., -2]-[1/sqrt(2), -1/(sqrt(2))]),
      ),
  Dict( "f"      => LogisticLoss([1.0, -1.0, 1.0, -1.0, 1.0], 1.5),
        "x"      => [-1.0, -2.0, 3.0, 2.0, 1.0],
        "∇f(x)"  => [-1.0965878679450072, 0.17880438303317633, -0.07113880976635019, 1.3211956169668235, -0.4034121320549927]
      ),
  Dict( "f"      => SqrHingeLoss([1.0, -1.0, 1.0, -1.0, 1.0], 1.5),
       "x"      => randn(MersenneTwister(0),Float64,5),
       "∇f(x)"  => gradient_fd(SqrHingeLoss([1.0, -1.0, 1.0, -1.0, 1.0], 1.5),randn(MersenneTwister(0),Float64,5))
      ),
  Dict( "f"      => CrossEntropy([1.0, 0., 1.0,  0., 1.0]),
       "x"      => rand(MersenneTwister(0),Float64,5),
       "∇f(x)"  => gradient_fd(CrossEntropy([1.0, 0., 1.0,  0., 1.0]),rand(MersenneTwister(0),Float64,5))
      ),
  Dict( "f"      => CrossEntropy([true, false, true,  false, true]),
       "x"      => rand(MersenneTwister(0),Float64,5),
       "∇f(x)"  => gradient_fd(CrossEntropy([true, false, true,  false, true]),rand(MersenneTwister(0),Float64,5))
      ),
  Dict( "f"      => CrossEntropy([true, false, true,  false, true]),
       "x"      => rand(MersenneTwister(0),Float64,1,5),
       "∇f(x)"  => gradient_fd(CrossEntropy([true, false, true,  false, true]),rand(MersenneTwister(0),Float64,1,5))
      ),
]

for i = 1:length(stuff)

  f = stuff[i]["f"]
  x = stuff[i]["x"]

  ref_∇f = stuff[i]["∇f(x)"]

  ref_fx = f(x)
  ∇f = similar(x)
  fx = gradient!(∇f, f, x)
  @test fx ≈ ref_fx
  @test ∇f ≈ ref_∇f

  for j = 1:11
    #For initial point x and 10 other random points
    fx = gradient!(∇f, f, x)
    for k = 1:10
      # Test conditions in different directions
      if ProximalOperators.is_convex(f)
        # Test ∇f is subgradient
	if typeof(f) <: CrossEntropy
		d = x.*(rand(Float64, size(x)).-1)./2 # assures 0 <= x+d <= 1
	else
		d = randn(Float64, size(x))
	end
        @test f(x+d) ≥ fx + dot(d, ∇f) - TOL_ASSERT
      else
        # Assume smooth function
        d = randn(Float64, size(x))
        d ./= norm(d).*sqrt(TOL_ASSERT)
        @test f(x+d) ≈ fx + dot(d, ∇f) atol=TOL_ASSERT
      end
    end
    x = rand(Float64, size(x))
  end
end


text:
const τ = 1.7724538509055159

function derivatives1!( β, Xh, xh, h, z, zr, g, dg )
    for t ∈ eachindex( xh )
        v = Xh .- xh[t]
        u = v[ -z[t] .≤ v .≤ zr[t] ]
        left = right = zero( eltype( β ) )
        for i ∈ eachindex( u )
            left += dg( u[i], 1, z[t], zr[t] )
            right += g( u[i], 1, z[t], zr[t] )
        end
        β[t,1] = - left / ( right * h )
    end
    return nothing
end

function derivativesmulti!( β, S, Xh, xh, h, z, zr, g, dg )
    mult = [ h^(s-1) / factorial(s-1) for s ∈ 1:S ]
    for t ∈ eachindex( xh )
        v = Xh .- xh[t]
        u = v[ -z[t] .≤ v .≤ zr[t] ]
        Tp = eltype( β )
        left = zeros( Tp, S ); right = zeros( Tp, S, S )
        for j ∈ 1:S
            for i ∈ eachindex( u )
                left[j] += dg( u[i], j, z[t], zr[t] )
                Rij = g( u[i], j, z[t], zr[t] )
                for s ∈ 1:S
                    right[j,s] += Rij * u[i]^(s-1) * mult[s]
                end
            end
        end
        β[t,:] = - ( right \ left ) / h
    end
    return nothing
end

function common( X, x, S, h, minx, maxx, anal )
    sanity_checks_basic( X, x, S, h, minx, maxx )
    Xh = ( X .- minx ) / h;  xh = ( x .- minx ) / h;  maxxh = (maxx - minx ) / h
    Tp = eltype( X )
    z = min.(xh, one( Tp ) );  zr = min.(maxxh .- xh, one( Tp ) )
    anal && sanity_checks_g( g, S, z, zl, zr )
    β = Array{Tp}(undef, length(x), S)
    return ( β, Xh, xh, z, zr, S )
end

"""
    β = logdensityderivatives( X, x, h; various optional arguments )

    Computes estimates of the derivatives of the log density of the data
    in the vector X evaluated at the elements in the vector x, using the
    bandwidth h. Returns a matrix β with rows corresponding to the
    evaluation points and columns to the derivatives.  This is a matrix
    even if S = 1.

# Mandatory arguments
- `X::Vector{T}`: vector of data
- `x::Vector{T}`: vector of evaluation points
- `h::T`: bandwidth

T is a floating point type

# Optional arguments
- `g::Function=g`: the function g to be used (see Pinkse and Schurter)
- `dg::Function=dg`: its derivative
- `S::Int=1`: the number of derivatives to be computed
- `minx::T=0.0`: the left hand side of the support
- `maxx::T=Inf`: the right hand side of the support
- `anal::Bool=true`: whether to run many sanity checks
"""
function logdensityderivatives(
    X::Vector{Tp},                    # data
    x::Vector{Tp},                    # evaluation points
    h::Tp;                            # bandwidth
    g = g,                            # g-function in paper
    dg = dg,                          # its derivative
    S = 1,                            # degree
    minx = 0.0,                       # bottom of support
    maxx = Inf,                       # top of support
    anal = true                       # whether to be fussy about checking arguments
    ) where {Tp<:Real}

    ( β, Xh, xh, z, zr, S ) = common( X, x, S, h, minx, maxx, anal )
    if S == 1
        derivatives1!( β, Xh, xh, h, z, zr, g, dg )
    else
        derivativesmulti!( β, S, Xh, xh, h, z, zr, g, dg )
    end
    return β
end

denfunc( x, c ) = 0.75 * ((c == 0.0) ? x - x^3 / 3.0 : exp( c * x ) * ( c*c * (1.0-x*x) + 2.0 * c * x - 2.0 ) / c^3)
function denfunc(x, a, b )
    b==0.0 && return denfunc( x, a )
    if b > 0.0
        sqrtb = sqrt(b)
        γ = 0.5 * (a + 2.0*b*x) / sqrtb
        return 0.75 * (- 0.125 * exp(-0.25*a*a/b) * ( τ * (a*a - 2.0 * b * (2.0*b+1.0) ) * erfi(γ) -
            2 * sqrtb * exp(γ*γ) * (a-2.0*b * x) ) / sqrtb^5)
    end
    b = -b
    sqrtb = sqrt(b)
    0.75 * (exp( x * (a-b*x)) * ( a + 2.0 * b *x ) * 0.25 / (b*b) - τ  * ( a * a + 2.0 * (1.0 - 2.0 * b ) * b) * exp( a*a *0.25/b) * erf( ( b*x -0.5*a ) / sqrtb ) * 0.125 / sqrtb^5)
end

"""
    (logf, β) = logdensity( X, x, h; various optional arguments )

    Computes estimates of the log density and its derivatives of the data
    in the vector X evaluated at the elements in the vector x, using the
    bandwidth h. Returns a tuple ( logf, β ) where logf is a vector whose
    elements correspond to the evaluation points and β a matrix with rows
    corresponding to the evaluation points and columns to the derivatives;
    this is a matrix even if S = 1.

# Mandatory arguments
- `X::Vector{T}`: vector of data
- `x::Vector{T}`: vector of evaluation points
- `h::T`: bandwidth

T is a floating point type

# Optional arguments
- `g::Function=g`: the function g to be used (see Pinkse and Schurter)
- `dg::Function=dg`: its derivative
- `S::Int=1`: the number of derivatives to be computed
- `minx::T=0.0`: the left hand side of the support
- `maxx::T=Inf`: the right hand side of the support
- `logf::Bool=true`: whether the log density itself should be computed
- `mz::Function=LogDensity.epanechnikov`: the kernel to be used for the log density itself
- `anal::Bool=true`: whether to run many sanity checks

Other choices for the kernel are:
LogDensity.quartic
LogDensity.gaussian
LogDensity.triweight
LogDensity.uniform
LogDensity.cosinus
"""
function logdensity(
    X::Vector{Tp},                    # data
    x::Vector{Tp},                    # evaluation points
    h::Tp;                            # bandwidth
    g = g,                            # g-function in paper
    dg = dg,                          # its derivative
    S = 1,                            # degree
    minx = 0.0,                       # bottom of support
    maxx = Inf,                       # top of support
    logf = true,                      # whether to provide log density function estimates
    mz::Function = epanechnikov,      # kernel used to provide estimates of f
    anal = true                       # whether to be fussy about checking arguments
    ) where {Tp<:Real}

    ( β, Xh, xh, z, zr, S ) = common( X, x, S, h, minx, maxx, anal )
    if S == 1
        derivatives1!( β, Xh, xh, h, z, zr, g, dg )
    else
        derivativesmulti!( β, S, Xh, xh, h, z, zr, g, dg )
    end

    logf || return ( nothing, β )

    N = length( X )
    num = zeros( Tp, length( xh ) )
    for t ∈ eachindex( xh )
        v = Xh .- xh[t]
        u = v[ -z[t] .≤ v .≤ zr[t] ]
        num[t] = sum( mz( u[i] ) for i ∈ eachindex( u ) ) / ( N * h )
    end

    if mz == epanechnikov && S ≤ 2
        if S == 1
            for t ∈ eachindex( x )
                num[t] /= ( denfunc( zr[t], β[t,1] * h ) - denfunc( -z[t], β[t,1] * h ) )
            end
        else
            for t ∈ eachindex( x )
                num[t] /= ( denfunc( zr[t],  β[t,1] * h, β[t,2] * h^2 * 0.5 ) - denfunc( -z[t], β[t,1] * h, β[t,2] * h^2 * 0.5 ) )
            end
        end
    else
        fac = [ factorial( s ) for s ∈ 1 : S ]
        for t ∈ eachindex( x )
            function numint( xx )
                mz(xx) * exp( sum( β[t,s] * (xx * h)^s / fac[s] for s = 1:S ) )
            end
            denom = try quadgk(numint, -z[t], zr[t], order = 15)[1]
                    catch e
                        jassert( false. "integration failure in denominator" )
                    end
            num[t] /= denom
        end
    end
    return ( log.( num ), β )
end


text:
function instantiate_fplane_1(T)
    coriolis = FPlane(T, f=π)
    return coriolis.f == T(π)
end

function instantiate_fplane_2(T)
    coriolis = FPlane(T, rotation_rate=2, latitude=30)
    return coriolis.f == T(2)
end

function instantiate_betaplane_1(T)
    coriolis = BetaPlane(T, f₀=π, β=2π)
    return coriolis.f₀ == T(π)
end

function instantiate_betaplane_2(T)
    coriolis = BetaPlane(T, latitude=70, radius=2π, rotation_rate=3π)
    return coriolis.f₀ == T(6π * sind(70))
end

@testset "Coriolis" begin
    println("Testing Coriolis...")

    @testset "Coriolis" begin
        for T in float_types
            @test instantiate_fplane_1(T)
            @test instantiate_fplane_2(T)
            @test instantiate_betaplane_1(T)
            @test instantiate_betaplane_2(T)
            # Test that FPlane throws an ArgumentError
            @test_throws ArgumentError FPlane(T, rotation_rate=7e-5)
            @test_throws ArgumentError FPlane(T, latitude=40)
            @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5)
            @test_throws ArgumentError FPlane(T, f=1, latitude=40)
            @test_throws ArgumentError FPlane(T, f=1, rotation_rate=7e-5, latitude=40)
            # Non-exhaustively test that BetaPlane throws an ArgumentError
            @test_throws ArgumentError BetaPlane(T, f₀=1)
            @test_throws ArgumentError BetaPlane(T, β=1)
            @test_throws ArgumentError BetaPlane(T, latitude=70)
            @test_throws ArgumentError BetaPlane(T, f₀=1e-4, β=1e-11, latitude=70)
        end
    end
end


text:
# Utility function which returns a Dual number of the index matches the target.
@inline function context_getindex(ctx, am, addr)
    if (addr, ) == ctx.target
        Dual(getindex(am, addr), 1.0)
    else
        getindex(am, addr)
    end
end

function forward(addr, params, cl::DynamicCallSite, seed)
    ctx = ForwardMode(addr, params, cl, seed)
    ret = ctx(cl.fn, cl.args...)
    ret, ctx.weight
end

# ------------ Choice sites ------------ #

@inline function (ctx::ForwardModeContext)(call::typeof(trace),
                                           addr::A,
                                           d::Distribution{K}) where {A <: Address, K}
    visit!(ctx, addr)
    v = context_getindex(ctx, ctx.map, addr)
    ctx.weight += logpdf(d, v)
    v
end

# ------------ Learnable ------------ #

@inline function (ctx::ForwardModeContext)(fn::typeof(learnable), addr::Address)
    visit!(ctx, addr)
    haskey(ctx.params, addr) && return context_getindex(ctx, ctx.params, addr)
    error("Parameter not provided at address $addr.")
end

# ------------ Call sites ------------ #

@inline function (ctx::ForwardModeContext)(c::typeof(trace),
                                           addr::A,
                                           call::Function,
                                           args...) where A <: Address
    visit!(ctx, addr)
    ps = get_sub(ctx.params, addr)
    ret, w = forward(ctx.target[2 : end], ps, get_sub(ctx.map, addr), Dual(1.0, 0.0))
    ctx.weight += w
    return ret
end

# ------------ Convenience ------------ #

function get_choice_gradient(addr::T, cl::C) where {T <: Tuple, C <: CallSite}
    fn = seed -> begin
        ret, w = forward(addr, Empty(), cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    cl[addr], d.partials.values[1]
end

function get_choice_gradient(addr::T, ps::P, cl::C) where {P <: AddressMap, T <: Tuple, C <: CallSite}
    fn = seed -> begin
        ret, w = forward(addr, ps, cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    cl[addr], d.partials.values[1]
end

function get_learnable_gradient(addr::T, ps::P, cl::C) where {P <: AddressMap, T <: Tuple, C <: CallSite}
    fn = seed -> begin
        ret, w = forward(addr, ps, cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    ps[addr], d.partials.values[1]
end

function get_choice_gradient(addr::T, cl::DynamicCallSite) where T <: Tuple
    fn = seed -> begin
        ret, w = forward(addr, Empty(), cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    cl[addr], d.partials.values[1]
end

function get_choice_gradient(addr::T, ps::P, cl::DynamicCallSite) where {P <: AddressMap, T <: Tuple}
    fn = seed -> begin
        ret, w = forward(addr, ps, cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    cl[addr], d.partials.values[1]
end

function get_learnable_gradient(addr::T, ps::P, cl::DynamicCallSite) where {P <: AddressMap, T <: Tuple}
    fn = seed -> begin
        ret, w = forward(addr, ps, cl, seed)
        w
    end
    d = fn(Dual(1.0, 0.0))
    ps[addr], d.partials.values[1]
end

# ------------ Documentation ------------ #


text:
module Reddit

import Requests
import TimeSeries

export countdailyposts

###### include ##################

include("functions.jl")

end


text:
# in file 'MyPrecompiledModule.jl' (e.g. in the same directory as myscript.jl)
__precompile__()
module MyPrecompiledModule
  export exportedfun;
  using ArgParse;
  function innerfun()
    println("Hello from MyPrecompiledModule.innerfun");
  end

  function main()
    for x in ARGS print(x) end
    s = ArgParseSettings()
    @add_arg_table s begin
        "-k"
            help = "low freeuncy kudge"
            arg_type = Int
            default = 2
        "-m"
            help = "low freeuncy kudge"
            arg_type = Int
            default = 1
        "-v"
            help = "set veborse"
            action = :store_true
    end

    args = parse_args(ARGS, s)

    for (arg,val) in args
            println("  $arg  =>  $val")
    end
  end

  function exportedfun()
    innerfun()
    print("Hello from MyPrecompiledModule.exportedfun");
  end
end

max_stars_repo_name: timm/misc

max_issues_repo_name: timm/misc

max_forks_repo_name: timm/misc