text:
module ElectronicCalculationInterface
abstract type coordinate_type end

struct real_space{T} <= coordinate_type where {T<=Real}
    x::T
    y::T
    z::T
end

struct momentum_space{T} <= coordinate_type where {T<=Real}
    kx::T
    ky::T
    kz::T
end

"""
    is_electronic_structure_solution(x)
Returns `true` if `x` represents the solution to some electronic structure calculation.
If true, `calculation_type` must return the method used (ex hartree_fock, or dft)
If true, `basis_set` must return the basis set used, where 'is_basis_set(x.basis_set)' is true
If true, `orbital_type` must return the type of orbitals used where is_orbital_set(x.orbital_set) is true
If true, `density_matrix` must return (see page 139 of sabo)
If true, then `coefficient_matrix` must return (see 136 of sabo)
"""
is_electronic_structure(x) = istemperature(typeof(x))
is_electronic_structure(x::Type{T}) where {T} = false

export is_electronic_structure

"""
    is_basis_set(x)
Returns `true` if `x` represents a set of functions (called basis functions), where the
basis functions can be used to represent an electronic wave function for a given method.
From wikipedia, the use of the basis set is equivalent to the use of an approximate resolution of the identity.
The basis set must be able to expand orbitals |ψ_i>, where

|ψ_i> ≈ Σ_j c_{i j} |phi_j > where expansion coefficients c_{i j} are given by the user to provide the
best fit for a given orbital

If true, `coordinate_type` must be defined for `x` as either `real_space` or `momentum_space`.
If true, `x` must impliment iterable interface (see julia docs)
if true, `index_type` must be defined for `x` appropriately
if true, `overlap(i::T,j::T) where {T <= x.index_type}` must be defined for `x` appropriately
If true, `x[idx::index_type]' returns a function `f(coords::T) where {T<=x.coordinate_type} appropriately
"""
is_basis_set(x) = is_basis_set(typeof(x))
is_basis_set(x::Type{T}) where {T} = false

export is_basis_set

end

abstract type CalculationResult{T} end

struct DFT_Result{T,P} <: CalculationResult{T}
    orbitals::AbstractVector
    number_of_iterations::Int
    input_parameters:P
end

electron_density(x::DFT_Result) = mapreduce(x->abs2(x),+, x.orbitals)

ElectronicCalculationInterface.is_electronic_structure(t::Type{<=DFT_Result}) = true


text:
module TwoBody

using LinearAlgebra: norm, ×, ⋅
using StaticArrays: SVector
using ReferenceFrameRotations: angle_to_dcm

using Roots: find_zero, Order5, Order2

using ..Util: angle, azimuth, normalize_angle

export
    EccentricAnomaly,
    Elliptic,
    Hyperbolic,
    MeanAnomaly,
    Parabolic,
    TrueAnomaly,
    cartesian,
    eccentric_anomaly,
    elliptic,
    hyperbolic,
    isprograde,
    isretrograde,
    keplerian,
    mean_anomaly,
    parabolic,
    perifocal,
    semilatus,
    transform,
    true_anomaly

include("stumpff.jl")
include("kepler.jl")

semilatus(a, ecc) = ecc ≈ 0.0 ? a : a * (1.0 - ecc^2)

function cartesian(a, ecc, i, Ω, ω, ν, μ)
    p = semilatus(a, ecc)
    r_pqw, v_pqw = perifocal(p, ecc, ν, μ)
    M = angle_to_dcm(-ω, -i, -Ω,:ZXZ)

    M * r_pqw, M * v_pqw
end

function perifocal(p, ecc, ν, μ)
    sinν, cosν = sincos(ν)
    sqrt_μop = sqrt(μ/p)

    r_pqw = [cosν, sinν, zero(p)] .* p ./ (1 + ecc * cosν)
    v_pqw = [-sinν, ecc + cosν, zero(sqrt_μop)] .* sqrt_μop

    r_pqw, v_pqw
end

iscircular(ecc, tol=1e-8) = isapprox(ecc, 0.0, atol=tol)
isequatorial(inc, tol=1e-8) = isapprox(abs(inc), 0.0, atol=tol)

function keplerian(pos, vel, µ, tol=1e-8)
    r = norm(pos)
    v = norm(vel)
    # Angular momentum
    h = pos × vel
    hm = norm(h)
    k = SVector(0.0, 0.0, 1.0)
    # Node vector
    node = k × h
    # Eccentricity vector
    e = ((v^2 - µ / r) * pos - (pos ⋅ vel) * vel) / µ
    ecc = norm(e)
    inc = angle(h, k)

    equatorial = isequatorial(inc, tol)
    circular = iscircular(ecc, tol)

    if circular
        # Semi-latus rectum
        a = hm^2 / µ
    else
        ξ = v^2 / 2 - µ / r
        a = -µ / (2ξ)
    end

    if equatorial && !circular
        Ω = 0.0
        # Longitude of pericenter
        ω = azimuth(e)
        ν = atan(h ⋅ (e × pos) / hm, pos ⋅ e)
    elseif !equatorial && circular
        Ω = azimuth(node)
        ω = 0.0
        # Argument of latitude
        ν = atan((pos ⋅ (h × node)) / hm, pos ⋅ node)
    elseif equatorial && circular
        Ω = 0.0
        ω = 0.0
        # True longitude
        ν = azimuth(pos)
    else
        if a > 0
            # Elliptic
            E_se = pos ⋅ vel / sqrt(μ * a)
            E_ce = r * v^2 / μ - 1
            ν = transform(elliptic, eccentric_anomaly, true_anomaly,
                          atan(E_se, E_ce), ecc)
        else
            # Hyperbolic
            E_sh = pos ⋅ vel / sqrt(-μ * a)
            E_ch = r * v^2 / μ - 1
            ν = transform(hyperbolic, eccentric_anomaly, true_anomaly,
                          log((E_ch + E_sh) / (E_ch - E_sh)) / 2, ecc)
        end
        Ω = azimuth(node)
        px = pos ⋅ node
        py = pos ⋅ (h × node) / hm
        ω = atan(py, px) - ν
    end
    Ω = mod2pi(Ω)
    ω = mod2pi(ω)
    ν = normalize_angle(ν, 0.0)

    a, ecc, inc, Ω, ω, ν
end

isprograde(inc) = inc < π/2 && !(inc ≈ π/2)
isretrograde(inc) = inc > π/2 && !(inc ≈ π/2)
ispolar(inc) = inc ≈ π/2

abstract type OrbitType end

struct Elliptic <: OrbitType end
struct Hyperbolic <: OrbitType end
struct Parabolic <: OrbitType end

const elliptic = Elliptic()
const hyperbolic = Hyperbolic()
const parabolic = Parabolic()

abstract type Anomaly end

struct MeanAnomaly <: Anomaly end
struct TrueAnomaly <: Anomaly end
struct EccentricAnomaly <: Anomaly end

const mean_anomaly = MeanAnomaly()
const true_anomaly = TrueAnomaly()
const eccentric_anomaly = EccentricAnomaly()

kepler(::Elliptic, E, M, ecc) = E - ecc * sin(E) - M
kepler(::Hyperbolic, E, M, ecc) = -E + ecc * sinh(E) - M
kepler(::Parabolic, E, M, ecc) = mean_parabolic(E, ecc) - M

function mean_parabolic(E, ecc, tolerance=1e-16)
    x = (ecc - 1.0) / (ecc + 1.0) * E^2
    small = false
    S = 0.0
    k = 0
    while !small
        term = (ecc - 1.0 / (2.0k + 2.0)) * x^k
        small = abs(term) < tolerance
        S += term
        k += 1
    end
    sqrt(2.0 / (1.0 + ecc)) * E + sqrt(2.0 / (1.0 + ecc)^3) * E^3 * S
end

"""
    transform(from, to, a, ecc)

Transform anomaly `a` `from` one anomaly type `to` another for an orbit with eccentricity `ecc`.

### Arguments ###

- `from`, `to`: Anomaly types
    - `true_anomaly`
    - `eccentric_anomaly`
    - `mean_anomaly`
- `a`: Current value of the anomaly
- `ecc`: Eccentricity of the orbit

### Output ###

Returns the transformed anomaly.

### References ###

- Farnocchia, Davide, Davide Bracali Cioci, and Andrea Milani.
    "Robust resolution of Kepler’s equation in all eccentricity regimes."
    Celestial Mechanics and Dynamical Astronomy 116, no. 1 (2013): 21-34.
"""
function transform(from, to, a, ecc)
    ecc ≈ 1.0 && return transform(parabolic, from, to, a, 1.0)
    ecc > 1.0 && return transform(hyperbolic, from, to, a, ecc)

    transform(elliptic, from, to, a, ecc)
end

transform(::Parabolic, ::EccentricAnomaly, ::TrueAnomaly, a, _) = 2.0 * atan(a)
transform(::Parabolic, ::TrueAnomaly, ::EccentricAnomaly, a, _) = tan(a / 2.0)

function transform(::Elliptic, ::TrueAnomaly, ::EccentricAnomaly, ν, ecc)
    β = ecc / (1 + sqrt(1 - ecc^2))
    ν - 2 * atan(β * sin(ν) / (1 + β * cos(ν)))
end
function transform(::Elliptic, ::EccentricAnomaly, ::TrueAnomaly, E, ecc)
    β = ecc / (1 + sqrt((1 - ecc) * (1 + ecc)))
    E + 2 * atan(β * sin(E) / (1 - β * cos(E)))
end

function transform(::Hyperbolic, ::TrueAnomaly, ::EccentricAnomaly, a, ecc)
    log((sqrt(ecc + 1) + sqrt(ecc - 1) * tan(a / 2)) /
        (sqrt(ecc + 1) - sqrt(ecc - 1) * tan(a / 2)))
end
function transform(::Hyperbolic, ::EccentricAnomaly, ::TrueAnomaly, a, ecc)
    2 * atan((exp(a) * sqrt(ecc + 1) - sqrt(ecc + 1)) /
             (exp(a) * sqrt(ecc - 1) + sqrt(ecc - 1)))
end

function transform(::Elliptic, ::MeanAnomaly, ::EccentricAnomaly, M, ecc)
    find_zero(E -> kepler(elliptic, E, M, ecc), M, Order5())
end
transform(::Elliptic, ::EccentricAnomaly, ::MeanAnomaly, E, ecc) = kepler(elliptic, E, 0.0, ecc)

function transform(::Hyperbolic, ::MeanAnomaly, ::EccentricAnomaly, M, ecc)
    find_zero(E -> kepler(hyperbolic, E, M, ecc), asinh(M / ecc), Order2())
end
transform(::Hyperbolic, ::EccentricAnomaly, ::MeanAnomaly, E, ecc) = kepler(hyperbolic, E, 0.0, ecc)

function transform(::Parabolic, ::MeanAnomaly, ::EccentricAnomaly, M, ecc)
    B = 3.0 * M / 2.0
    A = (B + sqrt(1.0 + B^2))^(2.0 / 3.0)
    guess = 2 * A * B / (1 + A + A^2)
    find_zero(E -> kepler(parabolic, E, M, ecc), guess, Order5())
end
transform(::Parabolic, ::EccentricAnomaly, ::MeanAnomaly, E, ecc) = kepler(parabolic, E, 0.0, ecc)

function transform(::OrbitType, ::MeanAnomaly, ::TrueAnomaly, M, ecc, delta=1e-2)
    if ecc > 1 + delta
        typ = hyperbolic
    elseif ecc < 1 - delta
        typ = elliptic
    else
        typ = parabolic
    end

    E = transform(typ, mean_anomaly, eccentric_anomaly, M, ecc)
    transform(typ, eccentric_anomaly, true_anomaly, E, ecc)
end
function transform(::OrbitType, ::TrueAnomaly, ::MeanAnomaly, ν, ecc, delta=1e-2)
    if ecc > 1 + delta
        typ = hyperbolic
    elseif ecc < 1 - delta
        typ = elliptic
    else
        typ = parabolic
    end

    E = transform(typ, true_anomaly, eccentric_anomaly, ν, ecc)
    transform(typ, eccentric_anomaly, mean_anomaly, E, ecc)
end

end

ext: jl


text:
# Running these tests requires StaticArrays

using Test
using StaticArrays
using StaticNumbers

@testset "mmatrix test" begin

    m = Matrix{Float64}(undef, static(3), static(2))
    @test m isa MMatrix
    @test eltype(m) == Float64
    @test size(m) === (3, 2)

    m = Matrix(undef, static(3), static(2))
    @test m isa MMatrix
    @test eltype(m) == Any
    @test size(m) === (3, 2)

end


text:
# module to hold optimization interface

module OptimizationInterface

using PDESolver  # evalFunctional, calcFunctionalDeriv
using SummationByParts
using ODLCommonTools
using Utils
using LinearSolvers
using NonlinearSolvers  # need jacobian calculation function  # unneded?
using PdePumiInterface  # need write vtk files
using MPI

# steady_adjoint.jl
export calcAdjoint

# error_estimate.jl
export doHAdaptation, solveAdaptive, calcErrorEstimate, AdaptOpts

include("steady_adjoint.jl")
include("error_estimate.jl")
include("target_els.jl")

end  # end module


text:
using CFLPLib
using Test

@testset "LoadSymbol" begin
    data = loadCFLP(:cap41)
    @test data.name == "cap41"
    @test data.facilities[3].capacity == 5000
    @test data.facilities[5].fixed_cost == 7500.0
    @test data.customers[7].demand == 2370
    @test length(data.facilities) == 16
    @test length(data.customers) == 50
    @test data.lb == -Inf
    @test data.ub == Inf
    @test_nowarn println(data)
    @test_nowarn println(data.facilities[2])
    @test_nowarn println(data.customers[6])
end

@testset "ErrorSymbol" begin
    data = loadCFLP(:notaninstance)
    @test data === nothing
end

@testset "LoadString" begin
    data = loadCFLP(joinpath(pkgdir(CFLPLib), "test/data/cap41.txt"))
    @test data.name == "cap41"
    @test data.facilities[3].capacity == 5000
    @test data.facilities[5].fixed_cost == 7500.0
    @test data.customers[7].demand == 2370
    @test length(data.facilities) == 16
    @test length(data.customers) == 50
    @test data.lb == -Inf
    @test data.ub == Inf
    @test_nowarn println(data)
    @test_nowarn println(data.facilities[2])
    @test_nowarn println(data.customers[6])
end

@testset "ErrorString" begin
    data = loadCFLP("notaninstance")
    @test data === nothing
end


text:
#=

    Data on the primitive binary polynomials (a) and the corresponding
    starting values m, for Sobol sequences in up to 21201 dimensions,
    taken from:
        S. Joe and F. Y. Kuo, Constructing Sobol sequences with better two-dimensional projections,
        SIAM J. Sci. Comput. 30, 2635-2654 (2008).

    Figures used were derived from files provided on Joe Kuo's website
    http://web.maths.unsw.edu.au/~fkuo/sobol/ on 16/10/2010
    and are provided under the license below.

    Copyright (c) 2008, Frances Y. Kuo and Stephen Joe
    All rights reserved.

    Redistribution and use in source and binary forms, with or without
    modification, are permitted provided that the following conditions are met:

        * Redistributions of source code must retain the above copyright
          notice, this list of conditions and the following disclaimer.

        * Redistributions in binary form must reproduce the above copyright
          notice, this list of conditions and the following disclaimer in the
          documentation and/or other materials provided with the distribution.

        * Neither the names of the copyright holders nor the names of the
          University of New South Wales and the University of Waikato
          and its contributors may be used to endorse or promote products derived
          from this software without specific prior written permission.

    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ``AS IS'' AND ANY
    EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
    WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
    DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS BE LIABLE FOR ANY
    DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
    (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
    LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
    ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
    (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
    SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using DelimitedFiles

const sobol_a_file = joinpath(@__DIR__, "sobol_a.csv")
const sobol_minit_file = joinpath(@__DIR__, "sobol_minit.csv")

# we need to re-compile if these files change:
include_dependency(sobol_a_file)
include_dependency(sobol_minit_file)

#successive primitive binary-coefficient polynomials p(z)
#   = a_0 + a_1 z + a_2 z^2 + ... a_31 z^31, where a_i is the
#     i-th bit of sobol_a[j] for the j-th polynomial.
const sobol_a = vec(readdlm(sobol_a_file, ',', UInt32))

# starting direction #'s m[i] = sobol_minit[i][j] for i=0..d of the
# degree-d primitive polynomial sobol_a[j].
const sobol_minit = readdlm(sobol_minit_file, ',', UInt32)


text:
#
# tabulated.jl -
#
# Implement unidimensional interpolation operator using precomputed indices and
# coefficients.
#
#------------------------------------------------------------------------------
#
# This file is part of the LinearInterpolators package licensed under the MIT
# "Expat" License.
#
# Copyright (C) 2016-2018, Éric Thiébaut.
#

module TabulatedInterpolators

export
    TabulatedInterpolator

using InterpolationKernels
using LazyAlgebra
using LazyAlgebra.Foundations
import LazyAlgebra: vcreate, apply!, apply

using ..LinearInterpolators
using ..LinearInterpolators: limits, getcoefs
import ..LinearInterpolators.Meta

struct TabulatedInterpolator{T<:AbstractFloat,S,D} <: LinearMapping
    nrows::Int     # length of output dimension
    ncols::Int     # length of input dimension
    J::Matrix{Int} # indices along input dimension
    W::Matrix{T}   # interpolation weights
    d::D           # dimension to interpolate or nothing
end

# Interpolator can be used as a function.
(A::TabulatedInterpolator)(x::AbstractArray) = apply(A, x)

nrows(A::TabulatedInterpolator) = A.nrows
ncols(A::TabulatedInterpolator) = A.ncols
width(A::TabulatedInterpolator{T,S,D}) where {T,S,D} = S
coefficients(A::TabulatedInterpolator) = A.W
columns(A::TabulatedInterpolator) = A.J

"""
    TabulatedInterpolator([T,] [d = nothing,] ker, pos, nrows, ncols)

yields a linear map to interpolate with kernel `ker` along a dimension of
length `ncols` to produce a dimension of length `nrows`.  The function `pos(i)`
for `i ∈ 1:nrows` gives the positions to interpolate in fractional index unit
along a dimension of length `ncols`.  Argument `d` is the rank of the dimension
along which to interpolate when the operator is applied to a multidimensional
array.  If it is unspecifed when the operator is created, it will have to be
specified each time the operator is applied (see below).

The positions to interpolate can also be specified by a vector `x` as in:

    TabulatedInterpolator([T,] [d = nothing,] ker, x, ncols)

to produce an interpolator whose output dimension is `nrows = length(x)`.  Here
`x` can be an abstract vector.

Optional argument `T` is the floating-point type of the coefficients of the
linear map.  By default, it is given by the promotion of the element type of
the arguments `ker` and, if specified, `x`.

The tabulated operator, say `A`, can be applied to an argument `x`:

    apply!(α, P::Operations, A, [d,] x, scratch, β, y) -> y

to overwrite `y` with `α*P(A)⋅x + β*y`.  If `x` and `y` are multi-dimensional,
the dimension `d` to interpolate must be specified.

"""
function TabulatedInterpolator(::Type{T},
                               d::D,
                               ker::Kernel{T,S},
                               pos::Function,
                               nrows::Integer,
                               ncols::Integer) where {T<:AbstractFloat,
                                                      D<:Union{Nothing,Int},S}
    nrows ≥ 1 || error("number of rows too small")
    ncols ≥ 1 || error("number of columns too small")
    J, W = __maketables(ker, pos, convert(Int, nrows), convert(Int, ncols))
    return TabulatedInterpolator{T,S,D}(nrows, ncols, J, W, d)
end

function TabulatedInterpolator(::Type{T},
                               d::D,
                               ker::Kernel{T,S},
                               x::AbstractVector{<:Real},
                               ncols::Integer) where {T<:AbstractFloat,
                                                      D<:Union{Nothing,Int},S}
    nrows = length(x)
    nrows ≥ 1 || error("number of positions too small")
    ncols ≥ 1 || error("number of columns too small")
    J, W = __maketables(ker, x, nrows, convert(Int, ncols))
    return TabulatedInterpolator{T,S,D}(nrows, ncols, J, W, d)
end

# This one forces the kernel to have the requested precision.
function TabulatedInterpolator(::Type{T},
                               d::Union{Nothing,Int},
                               ker::Kernel,
                               args...) where {T<:AbstractFloat}
    TabulatedInterpolator(T, d, T(ker), args...)
end

function TabulatedInterpolator(d::Union{Nothing,Int},
                               ker::Kernel{T},
                               pos::Function,
                               nrows::Integer,
                               ncols::Integer) where {T}
    TabulatedInterpolator(float(T), d, ker, pos, nrows, ncols)
end

function TabulatedInterpolator(d::Union{Nothing,Int},
                               ker::Kernel{Tk},
                               x::AbstractVector{Tx},
                               ncols::Integer) where {Tk,Tx}
    TabulatedInterpolator(float(promote_type(Tk, Tx)), d, ker, x, ncols)
end

TabulatedInterpolator(T::DataType, ker::Kernel, args...) =
    TabulatedInterpolator(T, nothing, ker, args...)

TabulatedInterpolator(ker::Kernel, args...) =
    TabulatedInterpolator(nothing, ker, args...)

TabulatedInterpolator(T::DataType, d::Integer, ker::Kernel, args...) =
    TabulatedInterpolator(T, convert(Int, d), ker, args...)

TabulatedInterpolator(d::Integer, ker::Kernel, args...) =
    TabulatedInterpolator(convert(Int, d), ker, args...)

@generated function __maketables(ker::Kernel{T,S},
                                 X::AbstractVector,
                                 nrows::Int,
                                 ncols::Int) where {T,S}
    J_ = [Symbol(:j_,s) for s in 1:S]
    W_ = [Symbol(:w_,s) for s in 1:S]
    code = (Meta.generate_getcoefs(J_, W_, :ker, :lim, :x),
            [:( J[$s,i] = $(J_[s]) ) for s in 1:S]...,
            [:( W[$s,i] = $(W_[s]) ) for s in 1:S]...)
    quote
        J = Array{Int}(undef, S, nrows)
        W = Array{T}(undef, S, nrows)
        lim = limits(ker, ncols)
        # The following assertion is to cope with
        # https://github.com/JuliaLang/julia/issues/40276
        if size(J,1) == S && size(W,1) == S
            @inbounds for i in 1:nrows
                x = convert(T, X[i])
                $(code...)
            end
        else
            throw_unexpected_stride()
        end
        return J, W
    end
end

@generated function __maketables(ker::Kernel{T,S},
                                 pos::Function,
                                 nrows::Int,
                                 ncols::Int) where {T,S}
    J_ = [Symbol(:j_,s) for s in 1:S]
    W_ = [Symbol(:w_,s) for s in 1:S]
    code = (Meta.generate_getcoefs(J_, W_, :ker, :lim, :x),
            [:( J[$s,i] = $(J_[s]) ) for s in 1:S]...,
            [:( W[$s,i] = $(W_[s]) ) for s in 1:S]...)
    quote
        J = Array{Int}(undef, S, nrows)
        W = Array{T}(undef, S, nrows)
        lim = limits(ker, ncols)
        # The following assertion is to cope with
        # https://github.com/JuliaLang/julia/issues/40276
        if size(J,1) == S && size(W,1) == S
            @inbounds for i in 1:nrows
                x = convert(T, pos(i))
                $(code...)
            end
        else
            throw_unexpected_stride()
        end
        return J, W
    end
end

throw_unexpected_stride() =
    throw(AssertionError("unexpected stride"))

# Source and destination arrays may have elements which are either reals or
# complexes.
const RorC{T<:AbstractFloat} = Union{T,Complex{T}}

@inline __errdims(msg::String) = throw(DimensionMismatch(msg))

@inline function __check(A::TabulatedInterpolator{T,S}) where {T,S}
    if !(size(A.J) == size(A.W) == (S, A.nrows))
        error("corrupted interpolation table")
    end
end

@inline function __fullcheck(A::TabulatedInterpolator{T,S}) where {T,S}
    # The loop below is a bit faster than using `extrema`.
    __check(A)
    J, ncols = A.J, A.ncols
    @inbounds for k in eachindex(J)
        if !(1 ≤ J[k] ≤ ncols)
            error("out of bound index in interpolation table")
        end
    end
end

# For vectors, the dimension to interpolate is 1.
function vcreate(::Type{P},
                 A::TabulatedInterpolator{R,S,Nothing},
                 src::AbstractVector{T},
                 scratch::Bool=false) where {P<:Operations,
                                             R<:AbstractFloat,S,
                                             T<:RorC{R}}
    __vcreate(P, A, 1, src)
end

function vcreate(::Type{P},
                 A::TabulatedInterpolator{R,S,Nothing},
                 src::AbstractArray{T,N},
                 scratch::Bool=false) where {P<:Operations,
                                             R<:AbstractFloat,S,
                                             T<:RorC{R},N}
    error("dimension to interpolate must be provided")
end

function vcreate(::Type{P},
                 A::TabulatedInterpolator{R,S,Nothing},
                 d::Integer,
                 src::AbstractArray{T,N},
                 scratch::Bool=false) where {P<:Operations,
                                             R<:AbstractFloat,S,
                                             T<:RorC{R},N}
    __vcreate(P, A, convert(Int, d), src)
end

function apply!(alpha::Real,
                ::Type{P},
                A::TabulatedInterpolator{R,S,Nothing},
                d::Integer,
                src::AbstractArray{T,N},
                scratch::Bool,
                beta::Real,
                dst::AbstractArray{T,N}) where {P<:Operations,
                                                R<:AbstractFloat,S,
                                                T<:RorC{R},N}
    __apply!(alpha, P, A, convert(Int, d), src, beta, dst)
end

# For vectors, the dimension to interpolate is 1.
function apply!(alpha::Real,
                ::Type{P},
                A::TabulatedInterpolator{R,S,Nothing},
                src::AbstractVector{T},
                scratch::Bool,
                beta::Real,
                dst::AbstractVector{T}) where {P<:Operations,
                                               R<:AbstractFloat,S,
                                               T<:RorC{R}}
    __apply!(alpha, P, A, 1, src, beta, dst)
end

function apply!(alpha::Real,
                ::Type{P},
                A::TabulatedInterpolator{R,S,Nothing},
                src::AbstractArray{T,N},
                scratch::Bool,
                beta::Real,
                dst::AbstractArray{T,N}) where {P<:Operations,
                                                R<:AbstractFloat,S,
                                                T<:RorC{R},N}
    error("dimension to interpolate must be provided")
end

function vcreate(::Type{P},
                 A::TabulatedInterpolator{R,S,Int},
                 src::AbstractArray{T,N},
                 scratch::Bool=false) where {P<:Operations,
                                             R<:AbstractFloat,S,
                                             T<:RorC{R},N}
    __vcreate(P, A, A.d, src)
end

function apply!(alpha::Real,
                ::Type{P},
                A::TabulatedInterpolator{R,S,Int},
                src::AbstractArray{T,N},
                scratch::Bool,
                beta::Real,
                dst::AbstractArray{T,N}) where {P<:Operations,
                                                R<:AbstractFloat,S,
                                                T<:RorC{R},N}
    __apply!(alpha, P, A, A.d, src, beta, dst)
end

function __vcreate(::Type{Direct},
                   A::TabulatedInterpolator{R,S,D},
                   d::Int,
                   src::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                   T<:RorC{R},N}
    1 ≤ d ≤ N || error("out of range dimension $d")
    nrows = A.nrows
    srcdims = size(src)
    dstdims = ntuple(i -> (i == d ? nrows : srcdims[i]), Val(N))
    return Array{T}(undef, dstdims)
end

function __vcreate(::Type{Adjoint},
                   A::TabulatedInterpolator{R,S,D},
                   d::Int,
                   src::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                   T<:RorC{R},N}
    1 ≤ d ≤ N || error("out of range dimension $d")
    ncols = A.ncols
    srcdims = size(src)
    dstdims = ntuple(i -> (i == d ? ncols : srcdims[i]), Val(N))
    return Array{T}(undef, dstdims)
end

function  __apply!(alpha::Real,
                   ::Type{Direct},
                   A::TabulatedInterpolator{R,S,D},
                   d::Int,
                   src::AbstractArray{T,N},
                   beta::Real,
                   dst::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                   T<:RorC{R},N}
    # Check arguments.
    __fullcheck(A)
    1 ≤ d ≤ N || error("out of range dimension $d")
    nrows, ncols = A.nrows, A.ncols
    srcdims, dstdims = size(src), size(dst)
    dstdims[d] == nrows ||
        __errdims("dimension $d of destination must be $nrows")
    srcdims[d] == ncols ||
        __errdims("dimension $d of source must be $ncols")
    predims = srcdims[1:d-1]
    postdims = srcdims[d+1:end]
    (dstdims[1:d-1] == predims && dstdims[d+1:end] == postdims) ||
        __errdims("incompatible destination and source dimensions")

    # Get rid of the case alpha = 0 and apply direct interpolation.
    if alpha == 0
        vscale!(dst, beta)
    else
        __direct!(convert(R, alpha), A, src,
                  CartesianIndices(predims), nrows, CartesianIndices(postdims),
                  convert(R, beta), dst)
    end
    return dst
end

function __apply!(alpha::Real,
                  ::Type{Adjoint},
                  A::TabulatedInterpolator{R,S,D},
                  d::Int,
                  src::AbstractArray{T,N},
                  beta::Real,
                  dst::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                  T<:RorC{R},N}
    # Check arguments.
    __fullcheck(A)
    1 ≤ d ≤ N || error("out of range dimension $d")
    nrows, ncols = A.nrows, A.ncols
    srcdims, dstdims = size(src), size(dst)
    dstdims[d] == ncols ||
        __errdims("dimension $d of destination must be $ncols")
    srcdims[d] == nrows ||
        __errdims("dimension $d of source must be $nrows")
    predims = srcdims[1:d-1]
    postdims = srcdims[d+1:end]
    (dstdims[1:d-1] == predims && dstdims[d+1:end] == postdims) ||
        __errdims("incompatible destination and source dimensions")

    # Scale destination by beta, and apply adjoint interpolation unless
    # alpha = 0.
    vscale!(dst, beta)
    if alpha != 0
        __adjoint!(convert(R, alpha), A, src, CartesianIndices(predims), nrows,
                   CartesianIndices(postdims), dst)
    end
    return dst
end

function __direct!(α::R,
                   A::TabulatedInterpolator{R,S,D},
                   src::AbstractArray{T,N},
                   I_pre::CartesianIndices,
                   nrows::Int,
                   I_post::CartesianIndices,
                   β::R,
                   dst::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                   T<:RorC{R},N}
    J, W = A.J, A.W
    # The following assertion is to cope with
    # https://github.com/JuliaLang/julia/issues/40276
    if size(J,1) == S && size(W,1) == S
        # We already know that α != 0.
        if α == 1 && β == 0
            @inbounds for i_post in I_post, i_pre in I_pre, i in 1:nrows
                a = zero(T)
                @simd for s in 1:S
                    j, w = J[s,i], W[s,i]
                    a += src[i_pre,j,i_post]*w
                end
                dst[i_pre,i,i_post] = a
            end
        elseif β == 0
            @inbounds for i_post in I_post, i_pre in I_pre, i in 1:nrows
                a = zero(T)
                @simd for s in 1:S
                    j, w = J[s,i], W[s,i]
                    a += src[i_pre,j,i_post]*w
                end
                dst[i_pre,i,i_post] = α*a
            end
        else
            @inbounds for i_post in I_post, i_pre in I_pre, i in 1:nrows
                a = zero(T)
                @simd for s in 1:S
                    j, w = J[s,i], W[s,i]
                    a += src[i_pre,j,i_post]*w
                end
                dst[i_pre,i,i_post] = α*a + β*dst[i_pre,i,i_post]
            end
        end
    else
        throw_unexpected_stride()
    end
end

function __adjoint!(α::R,
                    A::TabulatedInterpolator{R,S,D},
                    src::AbstractArray{T,N},
                    I_pre::CartesianIndices,
                    nrows::Int,
                    I_post::CartesianIndices,
                    dst::AbstractArray{T,N}) where {R<:AbstractFloat,S,D,
                                                    T<:RorC{R},N}
    J, W = A.J, A.W
    # The following assertion is to cope with
    # https://github.com/JuliaLang/julia/issues/40276
    if size(J,1) == S && size(W,1) == S
        # We already know that α != 0.
        if α == 1
            @inbounds for i_post in I_post, i_pre in I_pre, i in 1:nrows
                x = src[i_pre,i,i_post]
                if x != zero(T)
                    @simd for s in 1:S
                        j, w = J[s,i], W[s,i]
                        dst[i_pre,j,i_post] += w*x
                    end
                end
            end
        else
            @inbounds for i_post in I_post, i_pre in I_pre, i in 1:nrows
                x = α*src[i_pre,i,i_post]
                if x != zero(T)
                    @simd for s in 1:S
                        j, w = J[s,i], W[s,i]
                        dst[i_pre,j,i_post] += w*x
                    end
                end
            end
        end
    else
        throw_unexpected_stride()
    end
end

end # module


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

@testset "Test Formats for Spatial Reference Systems" begin
    @testset "PROJ4 Format" begin
        proj4326 = "+proj=longlat +datum=WGS84 +no_defs"
        proj26912 = "+proj=utm +zone=12 +datum=NAD83 +units=m +no_defs"
        AG.importPROJ4(proj4326) do spatialref
            spatialref2 = AG.importPROJ4(proj26912)
            @test AG.toPROJ4(spatialref2) == proj26912
            AG.importPROJ4!(spatialref2, AG.toPROJ4(spatialref))
            @test AG.toPROJ4(spatialref2) == proj4326
        end
    end

    @testset "WKT Format" begin
        wkt4326 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
        wkt26912 = "PROJCS[\"NAD83 / UTM zone 12N\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],TOWGS84[0,0,0,0,0,0,0],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4269\"]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-111],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"26912\"]]"
        AG.importWKT(wkt4326) do spatialref
            spatialref2 = AG.importWKT(wkt26912)
            @test AG.toWKT(spatialref2) == wkt26912
            AG.importWKT!(spatialref2, AG.toWKT(spatialref))
            @test AG.toWKT(spatialref2) == wkt4326
        end
    end

    @testset "ESRI Format" begin
        esri4326 = "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433],AXIS[\"Longitude\",EAST],AXIS[\"Latitude\",NORTH]]"
        esri26912 = "PROJCS[\"NAD83 / UTM zone 12N\",GEOGCS[\"NAD83\",DATUM[\"North_American_Datum_1983\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"6269\"]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]],PROJECTION[\"Transverse_Mercator\"],PARAMETER[\"latitude_of_origin\",0],PARAMETER[\"central_meridian\",-111],PARAMETER[\"scale_factor\",0.9996],PARAMETER[\"false_easting\",500000],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH]]"
        AG.importESRI(esri4326) do spatialref
            spatialref2 = AG.importESRI(esri26912)
            @test AG.toWKT(spatialref2) == esri26912
            AG.importESRI!(spatialref2, AG.toWKT(spatialref))
            @test AG.toWKT(spatialref2) == esri4326
        end
    end

    @testset "XML Format" begin
        xml4326 = """<gml:GeographicCRS gml:id="ogrcrs1">
          <gml:srsName>GCS_WGS_1984</gml:srsName>
          <gml:usesEllipsoidalCS>
            <gml:EllipsoidalCS gml:id="ogrcrs2">
              <gml:csName>ellipsoidal</gml:csName>
              <gml:csID>
                <gml:name codeSpace="urn:ogc:def:cs:EPSG::">6402</gml:name>
              </gml:csID>
              <gml:usesAxis>
                <gml:CoordinateSystemAxis gml:id="ogrcrs3" gml:uom="urn:ogc:def:uom:EPSG::9102">
                  <gml:name>Geodetic latitude</gml:name>
                  <gml:axisID>
                    <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9901</gml:name>
                  </gml:axisID>
                  <gml:axisAbbrev>Lat</gml:axisAbbrev>
                  <gml:axisDirection>north</gml:axisDirection>
                </gml:CoordinateSystemAxis>
              </gml:usesAxis>
              <gml:usesAxis>
                <gml:CoordinateSystemAxis gml:id="ogrcrs4" gml:uom="urn:ogc:def:uom:EPSG::9102">
                  <gml:name>Geodetic longitude</gml:name>
                  <gml:axisID>
                    <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9902</gml:name>
                  </gml:axisID>
                  <gml:axisAbbrev>Lon</gml:axisAbbrev>
                  <gml:axisDirection>east</gml:axisDirection>
                </gml:CoordinateSystemAxis>
              </gml:usesAxis>
            </gml:EllipsoidalCS>
          </gml:usesEllipsoidalCS>
          <gml:usesGeodeticDatum>
            <gml:GeodeticDatum gml:id="ogrcrs5">
              <gml:datumName>D_WGS_1984</gml:datumName>
              <gml:usesPrimeMeridian>
                <gml:PrimeMeridian gml:id="ogrcrs6">
                  <gml:meridianName>Greenwich</gml:meridianName>
                  <gml:greenwichLongitude>
                    <gml:angle uom="urn:ogc:def:uom:EPSG::9102">0</gml:angle>
                  </gml:greenwichLongitude>
                </gml:PrimeMeridian>
              </gml:usesPrimeMeridian>
              <gml:usesEllipsoid>
                <gml:Ellipsoid gml:id="ogrcrs7">
                  <gml:ellipsoidName>WGS_1984</gml:ellipsoidName>
                  <gml:semiMajorAxis uom="urn:ogc:def:uom:EPSG::9001">6378137</gml:semiMajorAxis>
                  <gml:secondDefiningParameter>
                    <gml:inverseFlattening uom="urn:ogc:def:uom:EPSG::9201">298.257223563</gml:inverseFlattening>
                  </gml:secondDefiningParameter>
                </gml:Ellipsoid>
              </gml:usesEllipsoid>
            </gml:GeodeticDatum>
          </gml:usesGeodeticDatum>
        </gml:GeographicCRS>
        """
        xml26912 = """<gml:ProjectedCRS gml:id="ogrcrs8">
          <gml:srsName>NAD_1983_UTM_Zone_12N</gml:srsName>
          <gml:baseCRS>
            <gml:GeographicCRS gml:id="ogrcrs9">
              <gml:srsName>GCS_North_American_1983</gml:srsName>
              <gml:usesEllipsoidalCS>
                <gml:EllipsoidalCS gml:id="ogrcrs10">
                  <gml:csName>ellipsoidal</gml:csName>
                  <gml:csID>
                    <gml:name codeSpace="urn:ogc:def:cs:EPSG::">6402</gml:name>
                  </gml:csID>
                  <gml:usesAxis>
                    <gml:CoordinateSystemAxis gml:id="ogrcrs11" gml:uom="urn:ogc:def:uom:EPSG::9102">
                      <gml:name>Geodetic latitude</gml:name>
                      <gml:axisID>
                        <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9901</gml:name>
                      </gml:axisID>
                      <gml:axisAbbrev>Lat</gml:axisAbbrev>
                      <gml:axisDirection>north</gml:axisDirection>
                    </gml:CoordinateSystemAxis>
                  </gml:usesAxis>
                  <gml:usesAxis>
                    <gml:CoordinateSystemAxis gml:id="ogrcrs12" gml:uom="urn:ogc:def:uom:EPSG::9102">
                      <gml:name>Geodetic longitude</gml:name>
                      <gml:axisID>
                        <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9902</gml:name>
                      </gml:axisID>
                      <gml:axisAbbrev>Lon</gml:axisAbbrev>
                      <gml:axisDirection>east</gml:axisDirection>
                    </gml:CoordinateSystemAxis>
                  </gml:usesAxis>
                </gml:EllipsoidalCS>
              </gml:usesEllipsoidalCS>
              <gml:usesGeodeticDatum>
                <gml:GeodeticDatum gml:id="ogrcrs13">
                  <gml:datumName>D_North_American_1983</gml:datumName>
                  <gml:usesPrimeMeridian>
                    <gml:PrimeMeridian gml:id="ogrcrs14">
                      <gml:meridianName>Greenwich</gml:meridianName>
                      <gml:greenwichLongitude>
                        <gml:angle uom="urn:ogc:def:uom:EPSG::9102">0</gml:angle>
                      </gml:greenwichLongitude>
                    </gml:PrimeMeridian>
                  </gml:usesPrimeMeridian>
                  <gml:usesEllipsoid>
                    <gml:Ellipsoid gml:id="ogrcrs15">
                      <gml:ellipsoidName>GRS_1980</gml:ellipsoidName>
                      <gml:semiMajorAxis uom="urn:ogc:def:uom:EPSG::9001">6378137</gml:semiMajorAxis>
                      <gml:secondDefiningParameter>
                        <gml:inverseFlattening uom="urn:ogc:def:uom:EPSG::9201">298.257222101</gml:inverseFlattening>
                      </gml:secondDefiningParameter>
                    </gml:Ellipsoid>
                  </gml:usesEllipsoid>
                </gml:GeodeticDatum>
              </gml:usesGeodeticDatum>
            </gml:GeographicCRS>
          </gml:baseCRS>
          <gml:definedByConversion>
            <gml:Conversion gml:id="ogrcrs16">
              <gml:coordinateOperationName>Transverse_Mercator</gml:coordinateOperationName>
              <gml:usesMethod xlink:href="urn:ogc:def:method:EPSG::9807" />
              <gml:usesValue>
                <gml:value uom="urn:ogc:def:uom:EPSG::9102">0</gml:value>
                <gml:valueOfParameter xlink:href="urn:ogc:def:parameter:EPSG::8801" />
              </gml:usesValue>
              <gml:usesValue>
                <gml:value uom="urn:ogc:def:uom:EPSG::9102">-111</gml:value>
                <gml:valueOfParameter xlink:href="urn:ogc:def:parameter:EPSG::8802" />
              </gml:usesValue>
              <gml:usesValue>
                <gml:value uom="urn:ogc:def:uom:EPSG::9001">0.9996</gml:value>
                <gml:valueOfParameter xlink:href="urn:ogc:def:parameter:EPSG::8805" />
              </gml:usesValue>
              <gml:usesValue>
                <gml:value uom="urn:ogc:def:uom:EPSG::9001">500000</gml:value>
                <gml:valueOfParameter xlink:href="urn:ogc:def:parameter:EPSG::8806" />
              </gml:usesValue>
              <gml:usesValue>
                <gml:value uom="urn:ogc:def:uom:EPSG::9001">0</gml:value>
                <gml:valueOfParameter xlink:href="urn:ogc:def:parameter:EPSG::8807" />
              </gml:usesValue>
            </gml:Conversion>
          </gml:definedByConversion>
          <gml:usesCartesianCS>
            <gml:CartesianCS gml:id="ogrcrs17">
              <gml:csName>Cartesian</gml:csName>
              <gml:csID>
                <gml:name codeSpace="urn:ogc:def:cs:EPSG::">4400</gml:name>
              </gml:csID>
              <gml:usesAxis>
                <gml:CoordinateSystemAxis gml:id="ogrcrs18" gml:uom="urn:ogc:def:uom:EPSG::9001">
                  <gml:name>Easting</gml:name>
                  <gml:axisID>
                    <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9906</gml:name>
                  </gml:axisID>
                  <gml:axisAbbrev>E</gml:axisAbbrev>
                  <gml:axisDirection>east</gml:axisDirection>
                </gml:CoordinateSystemAxis>
              </gml:usesAxis>
              <gml:usesAxis>
                <gml:CoordinateSystemAxis gml:id="ogrcrs19" gml:uom="urn:ogc:def:uom:EPSG::9001">
                  <gml:name>Northing</gml:name>
                  <gml:axisID>
                    <gml:name codeSpace="urn:ogc:def:axis:EPSG::">9907</gml:name>
                  </gml:axisID>
                  <gml:axisAbbrev>N</gml:axisAbbrev>
                  <gml:axisDirection>north</gml:axisDirection>
                </gml:CoordinateSystemAxis>
              </gml:usesAxis>
            </gml:CartesianCS>
          </gml:usesCartesianCS>
        </gml:ProjectedCRS>
        """
        AG.importXML(xml4326) do spatialref
            spatialref2 = AG.importXML(xml26912)
            @test startswith(AG.toXML(spatialref2), "<gml:ProjectedCRS")
            AG.importXML!(spatialref2, xml4326)
            @test startswith(AG.toXML(spatialref2), "<gml:GeographicCRS")
        end

        # ERROR: GDALError (CE_Failure, code 6):
        #   GDAL/OGR not compiled with libcurl support, remote requests not supported.
        # url4326 = "http://spatialreference.org/ref/epsg/4326/ogcwkt/"
        # proj4326 = "+proj=longlat +datum=WGS84 +no_defs "
        # urlsample = "http://spatialreference.org/ref/epsg/2039/esriwkt/"
        # projsample = "+proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +units=m +no_defs "
        # AG.importURL(url4326) do spatialref
        #     spatialref2 = AG.importURL(urlsample)
        #     @test AG.toWKT(spatialref2) == projsample
        #     AG.importURL!(spatialref2, url4326)
        #     @test AG.toWKT(spatialref2) == proj4326
        # end
    end
end

@testset "Cloning NULL SRS" begin
    @test sprint(print, AG.clone(AG.ISpatialRef())) == "NULL Spatial Reference System"
    AG.clone(AG.ISpatialRef()) do spatialref
        @test sprint(print, spatialref) == "NULL Spatial Reference System"
    end
end

@testset "Getting and Setting Attribute Values" begin
    AG.importEPSG(4326) do spatialref
        @test AG.toWKT(spatialref) == "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
        @test sprint(print, spatialref) == "Spatial Reference System: +proj=longlat +datum=WGS84 +no_defs"
        AG.setattrvalue!(spatialref, "GEOGCS|AUTHORITY|EPSG") # tests seems to be broken in gdal 3.0
        @test AG.toWKT(spatialref) == "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
        AG.setattrvalue!(spatialref, "GEOGCS|NEWATTRIBUTE", "7031") # tests seems to be broken in gdal 3.0
        @test AG.toWKT(spatialref) == "GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AXIS[\"Latitude\",NORTH],AXIS[\"Longitude\",EAST],AUTHORITY[\"EPSG\",\"4326\"]]"
        @test AG.getattrvalue(spatialref, "AUTHORITY", 0) == "EPSG"
        @test AG.getattrvalue(spatialref, "AUTHORITY", 1) == "4326"
    end
end


text:
module Model

using Compat
import ForwardDiff
using StaticArrays

# parameter types
export Image, SkyPatch, PsfComponent,
       GalaxyComponent, GalaxyPrototype,
       PriorParams, CanonicalParams, BrightnessParams, StarPosParams,
       GalaxyPosParams, GalaxyShapeParams,
       PsfParams, RawPSF, CatalogEntry,
       ParamSet, SkyIntensity

# functions
export align

# constants
export NUM_COLOR_COMPONENTS, NUM_SOURCE_TYPES, NUM_BANDS,
       BAND_LETTERS, galaxy_prototypes, prior,
       shape_standard_alignment,
       brightness_standard_alignment,
       gal_shape_alignment,
       ids_names,
       ids, star_ids, gal_ids, gal_shape_ids, psf_ids, bids

import Base: convert, length, +
import Distributions
import FITSIO, WCS
import WCS.WCSTransform
import ..Log
import ..Celeste: Const, @aliasscope, @unroll_loop
import ..SensitiveFloats
import ..SensitiveFloats: SensitiveFloat
import ..BivariateNormals
using ..BivariateNormals: BivariateNormalDerivatives, BvnComponent,
                          GalaxySigmaDerivs, get_bvn_cov, eval_bvn_pdf!,
                          get_bvn_derivs!, transform_bvn_derivs!,
                          transform_bvn_ux_derivs!

const cfgdir = joinpath(Pkg.dir("Celeste"), "cfg")

const BAND_LETTERS = "ugriz"

# The number of bands (colors).
const NUM_BANDS = length(BAND_LETTERS)

include("model/light_source_model.jl")
include("model/psf_model.jl")
include("model/image_model.jl")
include("model/param_set.jl")
include("model/imaged_sources.jl")
include("model/wcs_utils.jl")
include("model/fsm_util.jl")
include("model/log_prob.jl")

end  # module


text:
## Ziyi Yin converted the Compass velocity model to permeability and porosity

using DrWatson
@quickactivate "CompassTimeLapseCCS"

using Random, PyPlot
using LinearAlgebra, JLD2, MAT, Images, FFTW

n = (637, 82)
d = (25.0f0, 25.0f0)

JLD2.@load "../model/CompassSlice.jld2" v rho d # original Vp in 25m, 6m grid spacing

function find_water_bottom(m::AbstractArray{avDT,2};eps = 1e-4) where {avDT}
  #return the indices of the water bottom of a seismic image
  ### function from JUDI package
  n = size(m)
  idx = zeros(Integer, n[1])
  for j=1:n[1]
      k=1
      while true
          if abs(m[j,k]) > eps
              idx[j] = k
              break
          end
          k += 1
      end
  end
  return idx
end

function down_sample(temp1)
    ## downsample the model to 25m grid spacing w/ low-pass filter
    fft1 = fftshift(fft(temp1))
    mask = zeros(size(fft1))
    taper_length = 8
    half_taper_length = Int(taper_length/2)
    taper = (1f0 .+ sin.((Float32(pi)*(0:taper_length-1))/(taper_length - 1) .- Float32(pi)/2f0))/2f0
    mask[:,171-41:171+41] .= 1f0
    mask[:,171-41-taper_length+1+half_taper_length:171-41+half_taper_length] .= reshape(taper,1,taper_length)
    mask[:,171+41-half_taper_length:171+41+taper_length-1-half_taper_length] .= reshape(reverse(taper),1,taper_length)
    fft2 = mask .* fft1
    temp2 = real.(ifft(ifftshift(fft2)))
    temp2[:,1:35] .= temp1[:,1:35]  # fix water layer
    v = imresize(temp2,n)
    v[v.<=minimum(temp1)] .= minimum(temp1)
    v[v.>=maximum(temp1)] .= maximum(temp1)
    return v
end

v = down_sample(deepcopy(v))

extentx = (n[1]-1)*d[1]
extentz = (n[2]-1)*d[2]

idx_wb = find_water_bottom(v.-minimum(v))[1]

perm = 5f0*(v.-minimum(v))/(maximum(v)-minimum(v))

for i = 1:n[1]
  perm[i,1:idx_wb] .= 1f-4
  unconformaty = findall(v[i,:].>=3500f0)[1]
  perm[i,idx_wb+1:unconformaty] = perm[i,idx_wb+1:unconformaty] .+ 15f0
  perm[i,unconformaty+1:unconformaty+2] = perm[i,unconformaty+1:unconformaty+2]*1f-3
  perm[i,unconformaty+3:end] = perm[i,unconformaty+3:end] .+ 170f0
end

function poro_to_perm(phi)
    # Kozeny-carman equation
  c = (1.527/0.0314)^2f0
  perm = c*phi.^3f0./(1f0.-phi).^2f0
  return perm
end

phi = (perm*(1-0.22)^2/(1.527/0.0314)^2f0).^(1/3f0) # a proxy

JLD2.@save "../model/PermPoro.jld2" perm phi d


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

struct OperationsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Lists all of the available Microsoft.Resources REST API operations.
Param: api_version::String (required)
Return: OperationListResult
"""
function _swaggerinternal_operationsList(_api::OperationsApi, api_version::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", OperationListResult, "/providers/Microsoft.Resources/operations", ["azure_auth"])
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function operationsList(_api::OperationsApi, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_operationsList(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function operationsList(_api::OperationsApi, response_stream::Channel, api_version::String; _mediaType=nothing)
    _ctx = _swaggerinternal_operationsList(_api, api_version; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export operationsList


text:
# Section 4 of
# [PEDJ16] M. Philippe, R. Essick, G. E. Dullerud and R. M. Jungers.
# Stability of discrete-time switching systems with constrained switching sequences.
# Automatica, 72:242-250, 2016
# and running example of
# [LJP16] Legat, B.; Jungers, R. M. & Parrilo, P. A.
# Generating unstable trajectories for Switched Systems via Dual Sum-Of-Squares techniques
# Proceedings of the 19th International Conference on Hybrid Systems: Computation and Control, ACM, 2016, 51-60
# The JSR is 0.9748171979372074

using HybridSystems
using StaticArrays

A = @SMatrix [0.94 0.56; 0.14 0.46]
B = [0; 1]
k1 = -0.49
k2 = 0.27
As = [A + B * [k1 k2], A + B * [0 k2], A + B * [k1 0], A]
@assert eltype(As) <: SMatrix
G = LightAutomaton(4)
add_transition!(G, 1, 2, 3)
add_transition!(G, 2, 1, 2)
add_transition!(G, 1, 3, 1)
add_transition!(G, 3, 1, 2)
add_transition!(G, 2, 3, 1)
add_transition!(G, 3, 2, 3)
add_transition!(G, 3, 3, 1)
add_transition!(G, 3, 4, 4)
add_transition!(G, 4, 3, 1)
hs = discreteswitchedsystem(As, G)


text:
# Autogenerated wrapper script for TDLib_jll for x86_64-linux-musl-cxx03
export libtdjson

using OpenSSL_jll
using Zlib_jll
JLLWrappers.@generate_wrapper_header("TDLib")
JLLWrappers.@declare_library_product(libtdjson, "libtdjson.so.1.8.0")
function __init__()
    JLLWrappers.@generate_init_header(OpenSSL_jll, Zlib_jll)
    JLLWrappers.@init_library_product(
        libtdjson,
        "lib/libtdjson.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
filename = "ff_collapse.d"
infile = open(filename)
ρ      = log(1.0E4)
av     = log(10.0E0)
temp   = log(10.0E0)

outfile = open("new_isochoric_$filename","w")
header1 = "! time    log(Av)    log(n)    log(T)   \n"
header2 = "! (yr)   (mag)      (cm-3)    (K)      \n"
write(outfile,header1)
write(outfile,header2)
for ln in eachline(infile)
#  println(ln)
  newln = ln
  ρ_str = @sprintf("%9.3e",ρ)
  av_str = @sprintf("%9.3e",av)
  temp_str = @sprintf("%9.3e",temp)
  new_time = "$(newln[3:14])"
  new_time = float(new_time)
  time_str = @sprintf("%9.3e",new_time)
  println(time_str,ρ_str,av_str,temp_str)
  outline = "$(time_str)  $(ρ_str) $(av_str) $(temp_str)\n"
  write(outfile,outline)
end

close(infile)
close(outfile)


text:
mutable struct SingFunc
    Core::Array{Expr,1}
    Param::Symbol
end

function exec(SingFun,V)
    SingFunc = deepcopy(SingFun)
    prepend!(SingFunc.Core,[Expr(:(=),SingFunc.Param,V)])
    println(Expr(:block,SingFunc.Core...))
    eval(Expr(:block,SingFunc.Core...))

end # function

function func(SingFunc,Type)
    eval(Expr(:function,
         Expr(:tuple,Expr(:(::),SingFunc.Param,Type)),
         Expr(:block,SingFunc.Core...)))
end


text:
export IsRoleAllowedAsHand
IsRoleAllowedAsHand(arg1 :: CppEnum{Symbol("vr::ETrackedControllerRole"),UInt32}) :: Bool =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :IsRoleAllowedAsHand}}(), arg1)
IsRoleAllowedAsHand(arg1) =
  IsRoleAllowedAsHand(reinterpret4cxx(Cxx.CxxCore.CppEnum{Symbol("vr::ETrackedControllerRole"),UInt32}, arg1))

export ButtonMaskFromId
ButtonMaskFromId(arg1 :: CppEnum{Symbol("vr::EVRButtonId"),UInt32}) :: UInt64 =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :ButtonMaskFromId}}(), arg1)
ButtonMaskFromId(arg1) =
  ButtonMaskFromId(reinterpret4cxx(Cxx.CxxCore.CppEnum{Symbol("vr::EVRButtonId"),UInt32}, arg1))

export VRToken
VRToken() :: CppRef{UInt32,(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRToken}}())

export OpenVRInternal_ModuleContext
OpenVRInternal_ModuleContext() :: CppRef{CppBaseType{Symbol("vr::COpenVRContext")},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :OpenVRInternal_ModuleContext}}())

export VRSystem
VRSystem() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSystem")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRSystem}}())

export VRChaperone
VRChaperone() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRChaperone")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRChaperone}}())

export VRChaperoneSetup
VRChaperoneSetup() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRChaperoneSetup")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRChaperoneSetup}}())

export VRCompositor
VRCompositor() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRCompositor")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRCompositor}}())

export VROverlay
VROverlay() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVROverlay")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VROverlay}}())

export VRScreenshots
VRScreenshots() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRScreenshots")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRScreenshots}}())

export VRRenderModels
VRRenderModels() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRRenderModels")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRRenderModels}}())

export VRApplications
VRApplications() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRApplications")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRApplications}}())

export VRSettings
VRSettings() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSettings")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRSettings}}())

export VRResources
VRResources() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRResources")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRResources}}())

export VRExtendedDisplay
VRExtendedDisplay() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRExtendedDisplay")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRExtendedDisplay}}())

export VRTrackedCamera
VRTrackedCamera() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRTrackedCamera")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRTrackedCamera}}())

export VRDriverManager
VRDriverManager() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRDriverManager")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRDriverManager}}())

export VRInput
VRInput() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRInput")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRInput}}())

export VRIOBuffer
VRIOBuffer() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRIOBuffer")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRIOBuffer}}())

export VRSpatialAnchors
VRSpatialAnchors() :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSpatialAnchors")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VRSpatialAnchors}}())

export VR_Init
VR_Init(
    arg1 :: Ptr{CppEnum{Symbol("vr::EVRInitError"),UInt32}}
  , arg2 :: CppEnum{Symbol("vr::EVRApplicationType"),UInt32}
  , arg3 :: Ptr{UInt8}
  )      :: CppPtr{CxxQualType{CppBaseType{Symbol("vr::IVRSystem")},(false, false, false)},(false, false, false)} =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VR_Init}}(), arg1, arg2, arg3)
VR_Init(arg1, arg2, arg3) =
  VR_Init(
      reinterpret4cxx(Ptr{Cxx.CxxCore.CppEnum{Symbol("vr::EVRInitError"),UInt32}}, arg1)
    , reinterpret4cxx(Cxx.CxxCore.CppEnum{Symbol("vr::EVRApplicationType"),UInt32}, arg2)
    , reinterpret4cxx(Ptr{UInt8}, arg3)
    )

export VR_Shutdown
VR_Shutdown() :: Nothing =
  cppcall(__current_compiler__, CppNNS{(Tuple){:vr, :VR_Shutdown}}())


text:
using ..Continuations: ContinuationOptions
using ..Codim1
using ..Codim1: Codim1Solver, resolve_point
using ..BifurcationsBase: AbstractSpecialPoint, SpecialPoint
using ...Bifurcations: FixedPointBifurcationProblem, DiffEqWrapper
using DiffEqBase: ODEProblem, remake, solve

"""
    LimitCycleProblem(point::AbstractSpecialPoint,
                      solver)

Construct limit cycle continuation (codimension-1 bifurcation) problem
given a bifurcation `point`.
"""
function LimitCycleProblem(point::AbstractSpecialPoint,
                           solver::Codim1Solver;
                           kwargs...)

    cd1_prob = solver.prob :: FixedPointBifurcationProblem
    de_wrapper = cd1_prob.p :: DiffEqWrapper
    de_prob = de_wrapper.de_prob :: ODEProblem

    return LimitCycleProblem(
        resolve_point(point, solver),
        solver.opts;
        de_prob = de_prob,
        param_axis = de_wrapper.param_axis,
        t_domain = cd1_prob.t_domain,
        phase_space = cd1_prob.phase_space,
        kwargs...)
end

function LimitCycleProblem(point::SpecialPoint{Continuous},
                           opts::ContinuationOptions = ContinuationOptions();
                           num_mesh::Int = error("num_mesh required"),
                           degree::Int = error("degree required"),
                           amp = opts.h0,
                           kwargs...)
    @assert point.point_type == Codim1.PointTypes.hopf

    center = point.u[1:end - 1]
    t0 = point.u[end]

    vals, vecs = eigen(Codim1.ds_jacobian(point.J))
    _, idx0 = findmin(abs.(real.(vals)))  # real closest to zero
    v0 = normalize(vecs[:, idx0])
    w0 = imag(vals[idx0])
    if w0 < 0
        w0 = -w0
        v0 = conj(v0)
    end
    @assert w0 > opts.atol

    # TODO: use predictor-corrector rather than guessing amplitude
    vs = real(v0) .* amp
    vc = imag(v0) .* amp

    points = num_mesh * degree
    rads = range(0, step=2π / points, length=points)
    xs0 = center .+ (vs * sin.(rads)' .+
                     vc * cos.(rads)')

    l0 = 2π / w0  # period
    return LimitCycleProblem(;
        xs0 = xs0,
        l0 = l0,
        t0 = t0,
        num_mesh = num_mesh,
        degree = degree,
        kwargs...)
end

"""
    LimitCycleProblem(ode, param_axis, t_domain,
                      num_mesh, degree;
                      <keyword arguments>)

Solve `ode` and use it's solution as the initial limit cycle of the
continuation problem.

# Arguments
- `ode::ODEProblem`
- `param_axis::Lens`
- `t_domain::Tuple`
- `num_mesh::Integer`
- `degree::Integer`
"""
function LimitCycleProblem(de_prob::ODEProblem,
                           param_axis::Lens,
                           t_domain::Tuple,
                           num_mesh::Int,
                           degree::Int;
                           x0 = de_prob.u0,
                           l0 = de_prob.tspan[end] - de_prob.tspan[1],
                           time_offset = de_prob.tspan[1],
                           de_args = [],
                           de_opts = [],
                           kwargs...)
    ode = remake(
        de_prob;
        u0 = x0,
        tspan = (time_offset, time_offset + l0),
    )
    sol = solve(ode, de_args...; de_opts...)

    xs0 = similar(ode.u0, (length(ode.u0), num_mesh * degree))
    dt = l0 / (num_mesh * degree + 1)
    for i in 1:(num_mesh * degree)
        xs0[:, i] = sol(time_offset + (i - 1) * dt)
    end

    return LimitCycleProblem(;
        de_prob = ode,
        xs0 = xs0,
        l0 = l0,
        num_mesh = num_mesh,
        degree = degree,
        param_axis = param_axis,
        t_domain = t_domain,
        t0 = get(ode.p, param_axis),
        kwargs...)
end


text:
"""
Get and check item from the params, falling back on the definitions default.
"""
function get_item(key, def, key_path, param, classes)
    @debug "parsing $key_path ..."

    if haskey(param, key)
        return check_item(param[key], def, key_path; classes)
    else
        return def["default"]
    end
end

"""
Get the bin key as an integer, or error if it is not convertable.
"""
get_bin_int(i::Int) = i
function get_bin_int(s)
    try
        return parse(Int, s)
    catch
        error("Bins must be integers")
    end
end

"""
Get and validate a type == bin definition
"""
function get_bins(key, def, key_path, param, classes)
    @debug "parsing $key_path ..."

    if !haskey(param, key)
        return def["default"]
    end

    bins = merge_rec(def["default"], param[key])

    parsed_bins = Dict()
    for (bin, val) in bins
        bin_int = get_bin_int(bin)

        for (field, defn) in def["fields"]
            @assert haskey(val, field) "$field must be in $val [$key_path.$bin]"
            val[field] = check_item(val[field], defn, "$key_path.$bin.$field"; classes)
        end

        parsed_bins[bin_int] = val
    end

    return parsed_bins
end

"""
Get and validate a type == definition definition
"""
function get_defn(key, def, key_path, param, classes)
    @debug "parsing $key_path ..."

    if !haskey(param, key) || param[key] == Dict()
        parsed = def["default"]
    else
        parsed = param[key]
    end

    for (k, val) in parsed
        for (field, defn) in def["fields"]
            if !haskey(val, field)
                if haskey(defn, "default")
                    val[field] = defn["default"]
                else # no key in param and no default available, error
                    @assert haskey(val, field) "$field must be in $val [$key_path]"
                end
            end
        end
    end

    return parsed
end

"""
Recursively parse the passed params, using the definitions to validate
and provide defaults.
"""
function parse_params(
    defs::AbstractDict,
    params::AbstractDict;
    key_path::String = "",
    classes::AbstractDict = Dict(),
)
    @debug "parsing $key_path ..."

    # handles case of bin or def as direct default item
    if haskey(defs, "default")
        if defs["type"] == "bin"
            return get_bins("dummy", defs, key_path, Dict("dummy" => params), classes)
        else
            defs["type"] == "definition"
            return get_defn("dummy", defs, key_path, Dict("dummy" => params), classes)
        end
    end

    parsed = Dict()
    # all v are dicts at this point
    for (k, def) in defs
        kp = "$key_path.$k"
        if haskey(def, "default")
            dtype = def["type"]
            if dtype == "sub-dict"
                parsed[k] = parse_subdict(
                    def["default"],
                    def["keys"],
                    get(params, k, Dict()),
                    kp,
                    classes,
                )
            elseif dtype == "bin"
                parsed[k] = get_bins(k, def, kp, params, classes)
            elseif dtype == "definition"
                parsed[k] = get_defn(k, def, kp, params, classes)
            else
                parsed[k] = get_item(k, def, kp, params, classes)
            end
        else
            parsed[k] = parse_params(def, get(params, k, Dict()); key_path = kp, classes)
        end
    end

    return parsed
end

# params is a scalar, return it
parse_params(
    defs::AbstractDict,
    params;
    key_path::String = "",
    classes::AbstractDict = Dict(),
) = params

"""
Parse a type == sub-dict definition
"""
function parse_subdict(default, keys, params, key_path, classes)
    @debug "parsing $key_path ..."

    parsed = Dict()
    key, rem_keys = headtail(keys...)
    for val in keys_or_vals(classes[key])
        kp = "$key_path.$val"
        val_params = get(params, val, Dict())
        parsed[val] = parse_params(default, val_params; key_path = kp, classes)
        if length(rem_keys) > 0
            merge!(parsed[val], parse_subdict(default, rem_keys, val_params, kp, classes))
        end
    end

    return parsed
end

"""
Parse the classes definition first as it is needed in parsing the full params.
"""
function parse_classes(defs::AbstractDict, params::AbstractDict)
    if "classes" in keys(defs)
        return parse_params(defs["classes"], get(params, "classes", Dict()))
    else
        return Dict()
    end
end


text:
module PtrModif
  function divrem(a,b)
    r = Ref(CxxPtr{MyData}(C_NULL))
    q = PtrModif.divrem(CxxPtr(a),CxxPtr(b),r)
    return (q,r[])
  end

  include(joinpath(@__DIR__, "testcommon.jl"))
  @wrapmodule CxxWrap.CxxWrapCore.libpointer_modification

  function __init__()
    @initcxx
  end

end

using BenchmarkTools
using CxxWrap
using Test

@testset "$(basename(@__FILE__)[1:end-3])" begin

let d = PtrModif.MyData()
  PtrModif.setvalue!(d, 10)
  @test PtrModif.readpointerptr(Ref(CxxPtr(d))) == 10
  PtrModif.setvalue!(d, 20)
  @test PtrModif.readpointerref(Ref(CxxPtr(d))) == 20
end

let nd = Ref(CxxPtr{PtrModif.MyData}(C_NULL))
  @test isnull(nd[])
  PtrModif.writepointerref!(nd)
  @test !isnull(nd[])
  @test PtrModif.value(nd[][]) == 30
  PtrModif.setvalue!(nd[][], 40)
  @test PtrModif.value(nd[][]) == 40
  PtrModif.delete(nd[]) # Some manual management needed here
end

# Simulation of the Issue #133 use case
let a = PtrModif.MyData(9), b = PtrModif.MyData(2)
  (q,r) = PtrModif.divrem(a,b)
  @test PtrModif.value.((q,r[])) == (4,1)
  PtrModif.delete(r)
  (q,r) = PtrModif.divrem(q[],b)
  @test PtrModif.value(q) == 2
  @test isnull(r)
end

let a = PtrModif.MyData(9), b = PtrModif.MyData(2)
  (q,r) = PtrModif.prettydivrem(CxxPtr(a), CxxPtr(b))
  @test PtrModif.value.((q,r)) == (4,1)
end

end

GC.gc()
@test PtrModif.alive_count() == 0

# Must be after the GC test, because BenchmarkTools keeps a reference to the value

let a = PtrModif.MyData(11)
  println("value timing:")
  @btime PtrModif.value($a)
end


text:
"""
    activate_inner!(com, constraint::ActivatorConstraint)

Activate the inner constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_inner!(com, constraint::ActivatorConstraint)
    inner_constraint = constraint.inner_constraint
    if !constraint.inner_activated && inner_constraint.impl.activate
        !activate_constraint!(com, inner_constraint, inner_constraint.fct, inner_constraint.set) && return false
        constraint.inner_activated = true
        constraint.inner_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    activate_anti_inner!(com, constraint::ActivatorConstraint)

Activate the anti constraint of `constraint` when not activated yet.
Saves at which stage it was activated.
"""
function activate_anti_inner!(com, constraint::ActivatorConstraint)
    anti_constraint = constraint.anti_constraint
    if !constraint.anti_inner_activated && anti_constraint.impl.activate
        !activate_constraint!(com, anti_constraint, anti_constraint.fct, anti_constraint.set) && return false
        constraint.anti_inner_activated = true
        constraint.anti_inner_activated_in_backtrack_idx = com.c_backtrack_idx
    end
    return true
end

"""
    update_best_bound_constraint!(com::CS.CoM,
        constraint::ActivatorConstraint,
        fct::Union{MOI.VectorOfVariables, VAF{T}},
        set::IS,
        vidx::Int,
        lb::Int,
        ub::Int
    ) where {T<:Real}

Update the bound constraint associated with this constraint. This means that the `bound_rhs` bounds will be changed according to
the possible values the table constraint allows. `vidx`, `lb` and `ub` don't are not considered atm.
Additionally only a rough estimated bound is used which can be computed relatively fast.
This method calls the inner_constraint method if it exists and the indicator is activated.
"""
function update_best_bound_constraint!(
    com::CS.CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    vidx::Int,
    lb::Int,
    ub::Int,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    activator_vidx = constraint.indices[1]
    search_space = com.search_space
    activator_var = search_space[activator_vidx]
    if inner_constraint.impl.update_best_bound
        if CS.issetto(activator_var, Int(constraint.activate_on))
            return update_best_bound_constraint!(
                com,
                inner_constraint,
                inner_constraint.fct,
                inner_constraint.set,
                vidx,
                lb,
                ub,
            )
        else
            # if not activated (for example in a different subtree) we reset the bounds
            for rhs in constraint.bound_rhs
                rhs.lb = typemin(Int64)
                rhs.ub = typemax(Int64)
            end
        end
    end
    return true
end

function single_reverse_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    var::Variable,
    backtrack_idx::Int,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    # the variable must be part of the inner constraint
    if inner_constraint.impl.single_reverse_pruning &&
       (var.idx != constraint.indices[1] || constraint.activator_in_inner)
        single_reverse_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            var,
            backtrack_idx,
        )
    end
end

function reverse_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    backtrack_id::Int,
) where {
    T<:Real,
}
    # check if inner constraint should be deactived again
    if constraint.inner_activated && backtrack_id == constraint.inner_activated_in_backtrack_idx
        constraint.inner_activated = false
        constraint.inner_activated_in_backtrack_idx = 0
    end
    if constraint isa ReifiedConstraint
        if constraint.anti_inner_activated && backtrack_id == constraint.anti_inner_activated_in_backtrack_idx
            constraint.anti_inner_activated = false
            constraint.anti_inner_activated_in_backtrack_idx = 0
        end
    end

    inner_constraint = constraint.inner_constraint
    if inner_constraint.impl.reverse_pruning
        reverse_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            backtrack_id,
        )
    end
end

function restore_pruning_constraint!(
    com::CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
    prune_steps::Union{Int,Vector{Int}},
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    if inner_constraint.impl.restore_pruning
        restore_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
            prune_steps,
        )
    end
end

function finished_pruning_constraint!(
    com::CS.CoM,
    constraint::ActivatorConstraint,
    fct::Union{MOI.VectorOfVariables,VAF{T}},
    set,
) where {
    T<:Real,
}
    inner_constraint = constraint.inner_constraint
    if inner_constraint.impl.finished_pruning
        finished_pruning_constraint!(
            com,
            inner_constraint,
            inner_constraint.fct,
            inner_constraint.set,
        )
    end
end


text:
include("const.jl")

struct Model
    state::Array{Float64,3}
    state_tmp::Array{Float64,3}
    flux::Array{Float64,3}
    tend::Array{Float64,3}
    hy_dens_cell::Vector{Float64}
    hy_dens_theta_cell::Vector{Float64}
    hy_dens_int::Vector{Float64}
    hy_dens_theta_int::Vector{Float64}
    hy_pressure_int::Vector{Float64}
    weather_type::String
end

struct Grid
    nx::Int64
    nz::Int64
    nt::Int64
    dx::Float64
    dz::Float64
    dt::Float64
end

function init(nx_glob, nz_glob, sim_time, weather_type)
    r, hr, ht, u, w, t = 0.0, 0.0, 0.0, 0.0, 0.0, 0.0
    #  ierr = MPI_Init(argc,argv)

    #Set the cell grid size
    dx = xlen / nx_glob
    dz = zlen / nz_glob

    ##############################/
    # BEGIN MPI DUMMY SECTION
    # TODO: (1) GET NUMBER OF MPI RANKS
    #       (2) GET MY MPI RANK ID (RANKS ARE ZERO-BASED INDEX)
    #       (3) COMPUTE MY BEGINNING "I" INDEX (1-based index)
    #       (4) COMPUTE HOW MANY X-DIRECTION CELLS MY RANK HAS
    #       (5) FIND MY LEFT AND RIGHT NEIGHBORING RANK IDs
    ##############################/
    nranks = 1
    myrank = 0
    i_beg = 0
    nx = nx_glob
    left_rank = 0
    right_rank = 0
    #######################
    # END MPI DUMMY SECTION
    #######################
    println("num_vars = $NUM_VARS")
    #Vertical direction isn't MPI-ized, so the rank's local values = the global values
    k_beg = 0
    nz = nz_glob
    #masterproc = (myrank == 0)

    #Allocate the model data
    model = Model(
        zeros(nx + 2 * hs, nz + 2 * hs, NUM_VARS),  # state
        zeros(nx + 2 * hs, nz + 2 * hs, NUM_VARS),  # state_tmp
        zeros(nx + 1, nz + 1, NUM_VARS),            # flux
        zeros(nx, nz, NUM_VARS),                    # tend
        zeros(nz + 2 * hs),                         # hy_dens_cell
        zeros(nz + 2 * hs),                         # hy_dens_theta_cell
        zeros(nz + 1),                              # hy_dens_int
        zeros(nz + 1),                              # hy_dens_theta_int
        zeros(nz + 1),                              # hy_pressure_int
        weather_type,   # thermal, injection, mountain_waves, collision, turbulence, density_current
    )
    #Define the maximum stable time step based on an assumed maximum wind speed
    dt = min(dx, dz) / max_speed * cfl
    nt = round(sim_time / dt)

    grid = Grid(nx, nz, nt, dx, dz, dt)

    #Set initial elapsed model time and output_counter to zero
    etime = 0.0
    output_counter = 0.0

    #If I'm the master process in MPI, display some grid information
    #  if (masterproc) {
    println("nx_glob, nz_glob: $nx $nz")
    println("dx,dz: $dx $dz")
    println("dt: $dt")
    #  }
    #Want to make sure this info is displayed before further output
    #ierr = MPI_Barrier(MPI_COMM_WORLD)

    # Define quadrature weights and points
    nqpoints = 3
    qpoints = [
        0.112701665379258311482073460022,
        0.500000000000000000000000000000,
        0.887298334620741688517926539980,
    ]

    qweights = [
        0.277777777777777777777777777779,
        0.444444444444444444444444444444,
        0.277777777777777777777777777779,
    ]

    #####################################
    # Initialize the cell-averaged fluid state via Gauss-Legendre quadrature
    #####################################
    ########################/
    # TODO: MAKE THESE 2 LOOPS A PARALLEL_FOR
    ########################/
    for k = 1:nz+2*hs
        for i = 1:nx+2*hs
            #Use Gauss-Legendre quadrature to initialize a hydrostatic balance + temperature perturbation
            for kk = 1:nqpoints
                for ii = 1:nqpoints
                    #Compute the x,z location within the global domain based on cell and quadrature index
                    x = (i_beg + i - 1 - hs + 0.5) * dx + (qpoints[ii] - 0.5) * dx
                    z = (k_beg + k - 1 - hs + 0.5) * dz + (qpoints[kk] - 0.5) * dz

                    #Set the fluid state based on the user's specification
                    r, t, u, w, hr, ht = initialCondition(weather_type, x, z)
                    #Store into the fluid state array
                    model.state[i, k, ID_DENS] += r * qweights[ii] * qweights[kk]
                    model.state[i, k, ID_UMOM] += (r + hr) * u * qweights[ii] * qweights[kk]
                    model.state[i, k, ID_WMOM] += (r + hr) * w * qweights[ii] * qweights[kk]
                    model.state[i, k, ID_RHOT] +=
                        ((r + hr) * (t + ht) - hr * ht) * qweights[ii] * qweights[kk]
                end
            end
            for ll = 1:NUM_VARS
                model.state_tmp[i, k, ll] = model.state[i, k, ll]
            end
        end
    end
    #Compute the hydrostatic background state over vertical cell averages
    ########################/
    # TODO: MAKE THIS LOOP A PARALLEL_FOR
    ########################/
    for k = 1:nz+2*hs
        for kk = 1:nqpoints
            z = (k_beg + k - 1 - hs + 0.5) * dz
            #Set the fluid state based on the user's specification
            r, t, u, w, hr, ht = initialCondition(weather_type, 0.0, z)
            model.hy_dens_cell[k] = model.hy_dens_cell[k] + hr * qweights[kk]
            model.hy_dens_theta_cell[k] =
                model.hy_dens_theta_cell[k] + hr * ht * qweights[kk]
        end
    end
    #Compute the hydrostatic background state at vertical cell interfaces
    ########################/
    # TODO: MAKE THIS LOOP A PARALLEL_FOR
    ########################/
    for k = 1:nz+1
        z = (k_beg + k - 1) * dz
        #Set the fluid state based on the user's specification
        r, t, u, w, hr, ht = initialCondition(weather_type, 0.0, z)
        model.hy_dens_int[k] = hr
        model.hy_dens_theta_int[k] = hr * ht
        model.hy_pressure_int[k] = C0 * (hr * ht)^gamma
    end

    return model, grid
end

function initialCondition(funcname, x, z)
    func = Symbol(funcname)
    r, t, u, w, hr, ht = eval(:($func($x, $z)))
    return r, t, u, w, hr, ht
end

#This test case is initially balanced but injects fast, cold air from the left boundary near the model top
#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function injection(x, z)
    hr, ht = hydro_const_theta(z)
    r = 0.0
    t = 0.0
    u = 0.0
    w = 0.0
    return r, t, u, w, hr, ht
end

#Initialize a density current (falling cold thermal that propagates along the model bottom)
#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function density_current(x, z)
    hr, ht = hydro_const_theta(z)
    r = 0.0
    t = 0.0
    u = 0.0
    w = 0.0
    t = t + sample_ellipse_cosine(x, z, -20.0, xlen / 2, 5000.0, 4000.0, 2000.0)
    return r, t, u, w, hr, ht
end

#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function turbulence(x, z)
    hr, ht = hydro_const_theta(z)
    r = 0.0
    t = 0.0
    u = 0.0
    w = 0.0
    # call random_number(u);
    # call random_number(w);
    # u = (u-0.5)*20;
    # w = (w-0.5)*20;
    return r, t, u, w, hr, ht
end

#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function mountain_waves(x, z)
    hr, ht = hydro_const_bvfreq(z, 0.02)
    r = 0.0
    t = 0.0
    u = 15.0
    w = 0.0
    return r, t, u, w, hr, ht
end

#Rising thermal
#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function thermal(x, z)
    hr, ht = hydro_const_theta(z)
    r = 0.0
    t = 0.0
    u = 0.0
    w = 0.0
    t = t + sample_ellipse_cosine(x, z, 3.0, xlen / 2, 2000.0, 2000.0, 2000.0)
    #println(t)
    return r, t, u, w, hr, ht
end

#Colliding thermals
#x and z are input coordinates at which to sample
#r,u,w,t are output density, u-wind, w-wind, and potential temperature at that location
#hr and ht are output background hydrostatic density and potential temperature at that location
function collision(x, z, r, u, w, t, hr, ht)
    hr, ht = hydro_const_theta!(z)
    r = 0.0
    t = 0.0
    u = 0.0
    w = 0.0
    t = t + sample_ellipse_cosine(x, z, 20.0, xlen / 2, 2000.0, 2000.0, 2000.0)
    t = t + sample_ellipse_cosine(x, z, -20.0, xlen / 2, 8000.0, 2000.0, 2000.0)
end

#Establish hydrostatic balance using constant potential temperature (thermally neutral atmosphere)
#z is the input coordinate
#r and t are the output background hydrostatic density and potential temperature
function hydro_const_theta(z)
    theta0 = 300.0  #Background potential temperature
    exner0 = 1.0    #Surface-level Exner pressure
    #Establish hydrostatic balance first using Exner pressure
    ht = theta0                                  #Potential Temperature at z
    exner = exner0 - grav * z / (cp * theta0)   #Exner pressure at z
    p = p0 * exner^(cp / rd)                 #Pressure at z
    rt = (p / C0)^(1.0 / gamma)              #rho*theta at z
    hr = rt / ht                                  #Density at z
    return hr, ht
end

#Establish hydrstatic balance using constant Brunt-Vaisala frequency
#z is the input coordinate
#bv_freq0 is the constant Brunt-Vaisala frequency
#r and t are the output background hydrostatic density and potential temperature
function hydro_const_bvfreq(z, bv_freq0)
    theta0 = 300.0  #Background potential temperature
    exner0 = 1.0    #Surface-level Exner pressure
    ht = theta0 * exp(bv_freq0 * bv_freq0 / grav * z)                                    #Pot temp at z
    exner =
        exner0 - grav * grav / (cp * bv_freq0 * bv_freq0) * (ht - theta0) / (ht * theta0) #Exner pressure at z
    p = p0 * exner^(cp / rd)                                                         #Pressure at z
    rt = (p / C0)^(1.0 / gamma)                                                  #rho*theta at z
    hr = rt / ht
    return hr, ht                                                                    #Density at z
end

#Sample from an ellipse of a specified center, radius, and amplitude at a specified location
#x and z are input coordinates
#amp,x0,z0,xrad,zrad are input amplitude, center, and radius of the ellipse
function sample_ellipse_cosine(x, z, amp, x0, z0, xrad, zrad)
    #Compute distance from bubble center
    dist =
        sqrt(
            ((x - x0) / xrad) * ((x - x0) / xrad) + ((z - z0) / zrad) * ((z - z0) / zrad),
        ) * pi / 2.0
    #If the distance from bubble center is less than the radius, create a cos**2 profile
    if dist <= pi / 2.0
        return amp * cos(dist)^2.0
    else
        return 0.0
    end
end


text:
# Basic types for inferring a network

"""
Node with metadata

Fields:
* `label`: unique identifying label
* `binned_values`: data values discretized into bins
* `number_of_bins`: no. bins the data were discretized into
* `probabilities`: probability distribution across the bins
"""
struct Node
    label::String
    binned_values::Array{Int64}
    number_of_bins::Int64
    probabilities::Array{Float64}
end

# Constructs a Node from a line of a data file. line should be an array with
# the label as the first element, then the raw data values.
function Node(line::AbstractArray, discretizer, estimator, number_of_bins)

    label = string(line[1])
    raw_values = Array{Float64}(line[2:end])

    # Raw values are mapped to their bin IDs
    binned_values = zeros(Int, length(raw_values))

    # If the discretizer is Bayesian blocks, number_of_bins will be
    # overwritten by the ideal number of bins. Otherwise, it will remain
    # the same as the value passed in.
    number_of_bins = get_bin_ids!(raw_values, discretizer, number_of_bins, binned_values)

    probabilities = get_probabilities(estimator, get_frequencies_from_bin_ids(binned_values, number_of_bins))

    return Node(label, binned_values, number_of_bins, probabilities)

end

# Type for caching information between pairs of nodes:
# - mi: mutual information
# - si: specific information
struct NodePair
    mi::Float64
    si::Array{Float64}
end

"""
Undirected edge

Fields:
* `nodes`: the two nodes, in an arbitrary order
* `weight`: weight indicating confidence of edge existing in the true network
Weights are used to rank the edges, and different algorithms may have a
different scale. The relative weights within one inferred network are
therefore more meaningful than the absolute weight out of context.
"""
struct Edge
    nodes::Array{Node}
    weight::Float64
end


text:
# This computes the IPM
#
# sup_{g\in G} |E_P g(X) - E_Q g(Z)|
#
# by using a RKHS that contains only mean-zero functions under P.
# Given some base kernel k(x,x'), we construct a new kernel
# via the mapping
#
# k_0(x, x') = <\grad_x, \grad_x' k(x,x')> +
#   <u(x), \grad_x' k(x,x')> + <u(x'), \grad_x k(x,x')> +
#   <u(x), u(x')> k(x,x')
#
# where u(x) = \grad_x log p(x) is the score function.
#
# This is useful, because if H_0 is the RKHS with kernel k_0,
# we can solve the problem
#
# sup_{h\in H_0} \sum q_i h(x_i) - 1/2 ||h||_{H0}^2
#
# explicity with solution h(x) = \sum q_i k_0(x, x_i). This leaves
# the Stein discrepancy at \sum q_i q_j k_0(x_i, x_j).
#
# Args:
#  points  - the support points comprising the sample distribution
#  weights - the weights of the discrete sample
#  gradlogdensity - the gradlogdensity of the target distribution
#  kernel - the base kernel to make RKHS mean zero under target
#  checkpoints - if not empty, expects an array of ints that
#                mark intermediate n to compute the kernel. This is
#                more optimized than calling this function len(checkpoints) times
#  ncores - the number of cores to break the chunks into
function langevin_kernel_discrepancy(; points=[],
                                     weights=fill(1/size(points,1), size(points,1)),
                                     gradlogdensity=nothing,
                                     kernel=nothing,
                                     checkpoints=[],
                                     ncores=nprocs())
    # make sure kernel is defined
    if !isa(kernel, SteinKernel)
        error("Must supply a subtype of SteinKernel")
    end
    # make sure points is a 2-d array
    if length(size(points)) == 1
        points = reshape(points, (length(points), 1))
    end
    # start the timer
    tic()
    ## Extract inputs
    weights = vec(weights)
    n = length(weights)
    d = size(points,2)
    nocheckpoints = (length(checkpoints) == 0)
    # initialize checkpoints and sort
    if nocheckpoints
        sortedcheckpoints = [n]
        checkpointsindex = [1]
    else
        checkpointsindex = sortperm(checkpoints)
        sortedcheckpoints = checkpoints[checkpointsindex]
    end
    # Compute score functions for each observation
    gradlogdensities = compute_gradlogdensities(gradlogdensity, points)
    # chunk up the data in big batches
    datachunks = Any[
        Any[kernel, points, gradlogdensities, weights, sortedcheckpoints, i, ncores]
        for i in 1:ncores
    ]
    # parallelize the work
    worker_discrepancies = pmap((args) -> compute_kernel_sums(args...), datachunks)
    # merge the worker results and reweight discrepancies
    checkpoint_discrepancies = merge_worker_discrepancies(
        worker_discrepancies
    )
    # renormalize the discrepancies
    checkpoint_discrepancies = renormalize_discrepancies(
        checkpoint_discrepancies,
        weights,
        sortedcheckpoints
    )
    # put the discrepancies back in original order
    discrepancies = Array{Float64}(length(checkpoint_discrepancies))
    discrepancies[checkpointsindex] = checkpoint_discrepancies
    # stop the timer
    solvetime = toc()

    if nocheckpoints
        return LangevinKernelResult(
            points,
            weights,
            discrepancies[1],
            solvetime
        );
    else
        return LangevinKernelCheckpointsResult(
            points,
            weights,
            discrepancies,
            checkpoints,
            solvetime
        );
    end
end

# This is a helper function to compute the score functions
# used in the final kernel evaluations.
function compute_gradlogdensities(gradlogdensity::Function,
                                  points::Array{Float64,2})
    n, d = size(points)
    gradlogdensities = @parallel (vcat) for i=1:n
        gradlogdensity(points[i,:])'
    end
    gradlogdensities
end

# Computes the sum of discrepancies for indices (i,j) s.t.
# i % gap == istart and j >= i.
#
# This method is a bit clever b/c of checkpoints. It returns
# and array the size of checkpoints of the discrepancies
# that fall into each class. That is, suppose
#
# checkpoints = [n1, n2, n]
#
# Then an array of size 3 will be returned, where the first value
# is all (i, j) pairs above but i <= j <= n1. The second value will
# be all (i, j) pairs above but n1 < i <= j <= n2, and so on.
#
# @params
# kernel: the Stein kernel
# points: the xi
# gradlogdensities: the gradient of log densities
# weights: the weights
# checkpoints: SORTED array of n's to compute discrepancy for first n points
# istart: the lower bound
# gap: the gap to jump
# @returns: array of partial kernel sums
function compute_kernel_sums(kernel::SteinKernel,
                             points::Array{Float64,2},
                             gradlogdensities::Array{Float64,2},
                             weights::Array{Float64,1},
                             checkpoints::Array{Int},
                             istart::Int,
                             gap::Int)
    n = length(weights)
    discrepancies = [0.0 for _ in checkpoints]
    for i in istart:gap:n
        dindex = findfirst(checkpoints .>= i)
        if dindex == 0
            # biggest checkpoint smaller than n, no op
            break
        end
        for j in i:n
            if checkpoints[dindex] < j
                # crossed into a new cross section
                dindex += 1
            end
            m = (i == j) ? 1 : 2
            discrepancies[dindex] += m * compute_kernel_pair(
                 kernel,
                 points[i,:],
                 gradlogdensities[i,:],
                 weights[i],
                 points[j,:],
                 gradlogdensities[j,:],
                 weights[j])
        end
    end
    discrepancies
end

# This is a helper function farmed out to each worker
# that computes k(xi, xj) for the Stein kernel.
# @params
# kernel: the Stein kernel
# xi: the point xi
# ui: the score function at xi
# wi: the weight for xi
# xj: the point xj
# uj: the score function at xj
# wj: the weight for xj
function compute_kernel_pair(
        kernel::SteinKernel,
        xi::Array{Float64,1},
        ui::Array{Float64,1},
        wi::Float64,
        xj::Array{Float64,1},
        uj::Array{Float64,1},
        wj::Float64)
    kij = k0(kernel, xi, xj, ui, uj)
    wi * wj * kij
end

# Private helper function to merge discrepancies from workers
function merge_worker_discrepancies(
    worker_discrepancies::Array{Array{Float64,1},1})
    # merge all the workers
    checkpoint_discrepancies = reduce(+, worker_discrepancies)
    checkpoint_discrepancies
end

function renormalize_discrepancies(
    checkpoint_discrepancies::Array{Float64,1},
    weights::Array{Float64,1},
    checkpoints::Array{Int})
    # sum all the pairs
    cumulative_discrepancies = cumsum(checkpoint_discrepancies)
    # renormalize checkpoint weights to sum to 1
    checkpointweights = map(i -> sum(weights[1:i]), checkpoints)
    cumulative_discrepancies = cumulative_discrepancies ./ checkpointweights.^2
    cumulative_discrepancies
end