text:
function kernel_invocations_gemm(n)

  reps = parse(Int64, ENV["LAMP_REPS"])

  A = randn(n, n)
  B = randn(n, n)
  C = randn(n, n)

  Benchmarker.add_data(csv, "gemm_implicit", Benchmarker.measure(reps, gemm_implicit, A, B, C))
  Benchmarker.add_data(csv, "gemm_implicit_noup", Benchmarker.measure(reps, gemm_implicit_noup, A, B))
  Benchmarker.add_data(csv, "gemm_implicit_compact", Benchmarker.measure(reps, gemm_implicit_compact, A, B, C))
  Benchmarker.add_data(csv, "gemm_implicit_coeff", Benchmarker.measure(reps, gemm_implicit_coeff, A, B, C))
  Benchmarker.add_data(csv, "gemm_implicit_double_coeff", Benchmarker.measure(reps, gemm_implicit_double_coeff, A, B, C))
  Benchmarker.add_data(csv, "gemm_explicit", Benchmarker.measure(reps, gemm_explicit, A, B, C))
  Benchmarker.add_data(csv, "gemm_explicit_noup", Benchmarker.measure(reps, gemm_explicit_noup, A, B))
  Benchmarker.add_data(csv, "gemm_explicit_semi", Benchmarker.measure(reps, gemm_explicit_semi, A, B, C))

  @test(isapprox(gemm_implicit(A,B,C)[2], gemm_explicit(A,B,C)[2]))
  @test(isapprox(gemm_implicit(A,B,C)[2], gemm_implicit_compact(A,B,C)[2]))
  @test(isapprox(gemm_implicit(A,B,C)[2], gemm_explicit_semi(A,B,C)[2]))
  @test(isapprox(gemm_implicit_noup(A,B)[2], gemm_explicit_noup(A,B)[2]))

  A = LowerTriangular(randn(n, n))
  Benchmarker.add_data(csv, "trmm_implicit", Benchmarker.measure(reps, trmm_implicit, A, B))
  Benchmarker.add_data(csv, "trmm_implicit_compact", Benchmarker.measure(reps, trmm_implicit_compact, A, B))
  Benchmarker.add_data(csv, "trmm_explicit", Benchmarker.measure(reps, trmm_explicit, A, B))

  @test istril(A)
  @test(isapprox(trmm_implicit(A,B)[2], trmm_explicit(A,B)[2]))

  A = Diagonal(randn(n, n))
  B = randn(n, n)

  Benchmarker.add_data(csv, "diagmm", Benchmarker.measure(reps, diagmm, A, B))

  @test isdiag(A)
end

function diagmm(A::Diagonal{Float64}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,2))
  A_copy = Matrix(copy(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A_copy * B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_implicit(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A * B + result

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_implicit_coeff(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = 3.0 * A * B + result

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_implicit_double_coeff(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = 3.0 * A * B + 3.0 * result

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_implicit_compact(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result += A * B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_implicit_noup(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A * B

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function trmm_implicit_compact(A::LowerTriangular{Float64}, B::DenseArray{Float64, 2})
  result = copy(B)
  A_copy = Matrix(copy(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result *= A_copy

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function trmm_implicit(A::LowerTriangular{Float64}, B::DenseArray{Float64, 2})
  result = copy(B)
  A_copy = Matrix(copy(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = A_copy * result

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_explicit(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('N', 'N', 1.0, A, B, 1.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function trmm_explicit(A::LowerTriangular{Float64}, B::DenseArray{Float64, 2})
  result = copy(B)
  A_copy = Matrix(copy(A))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.trmm!('L', 'L', 'N', 'N', 1.0, A_copy, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_explicit_noup(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2})
  result = zeros(size(A,1), size(B,2))
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  LinearAlgebra.BLAS.gemm!('N', 'N', 1.0, A, B, 0.0, result)

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end

function gemm_explicit_semi(A::DenseArray{Float64, 2}, B::DenseArray{Float64, 2}, C::DenseArray{Float64, 2})
  result = copy(C)
  start::Float64 = 0.0
  finish::Float64 = 0.0
  Benchmarker.cachescrub()
  Benchmarker.gcscrub()
  start = time_ns()

  result = LinearAlgebra.BLAS.gemm('N', 'N', 1.0, A, B) + result

  finish = time_ns()
  return ((finish-start)*1e-9, result)
end


text:
"""
    DiscreteDynamics

A dynamics model of the form

```math
x_{k+1} = f(x_k, u_k)
```

where ``x`` is the `n`-dimensional state vector and ``u`` is the `m`-dimensional
control vector. The following methods should be defined on any model:

    state_dim(model)
    control_dim(model)

The output dimension `output_dim` is automatically defined to be equal to the
state dimension.

## Defining the dynamics
To define the dynamics using out-of-place (returning an SVector) evaluation,
the user must define one of the following:

    discrete_dynamics(model, z::AbstractKnotPoint)
    discrete_dynamics(model, x, u, t, dt)

where `t` is the current time and `dt` is the time step between ``x_k`` and ``x_{k+1}``.

The user shouldn't assume that the inputs are static arrays, although for best
computational performance these methods should be called by passing in static arrays.

To define the dynamics using in-place evaluation, the user must define one of the following:

    discrete_dynamics!(model, xn, z::AbstractKnotPoint)
    discrete_dynamics!(model, xn, x, u, t, dt)

Either of these function can be called by dispatching on the [`FunctionSignature`](@ref)
using:

    discrete_dynamics!(sig, model, xn, z::AbstractKnotPoint)

## Discretizing continuous dynamics
This type is the abstract type representing all discrete dynamics systems, whether
or not they are approximations of continuous systems. Users can choose to implement
their dynamics directly in discrete time using this interface, or can use
[`DiscretizedDynamics`](@ref) to discretize a continuous time system.

## Implicit dynamics functions
Sometimes continuous time systems are approximated using implicit integrators of the
form

```math
d(x_{k+1}, u_{k+1}, x_k, u_k) = 0
```

To evaluate these types of systems, use the [`dynamics_error`](@ref) and
[`dynamics_error!`](@ref) methods, which will return the output of `f` for
two consecutive knot points. The default functions such as `discrete_dynamics` or
`jacobian!` are not defined for these types of systems.

## Simulating the dynamics
Use [`propagate_dynamics!`](@ref) to conveniently save the next state directly into
the state vector of another [`KnotPoint`](@ref).
"""
abstract type DiscreteDynamics <: AbstractModel end
@inline evaluate(model::DiscreteDynamics, z::AbstractKnotPoint) =
    discrete_dynamics(model, z::AbstractKnotPoint)
@inline evaluate!(model::DiscreteDynamics, xn, z::AbstractKnotPoint) =
    discrete_dynamics!(model, xn, z::AbstractKnotPoint)

"""
    discrete_dynamics(model, z::AbstractKnotPoint)
    discrete_dynamics(model, x, u, t, dt)

Evaluate the discrete time dynamics, returning the output ``x_{k+1}``. For best
performance, the output should be a `StaticArrays.SVector`. This method is called
when using the `StaticReturn` [`FunctionSignature`](@ref).

Calling `evaluate` on a [`DiscreteDynamics`](@ref) model will call this function.
"""
discrete_dynamics(model::DiscreteDynamics, z::AbstractKnotPoint) =
    discrete_dynamics(model, state(z), control(z), time(z), timestep(z))
discrete_dynamics(model::DiscreteDynamics, x, u, t, dt) =
    error("Discrete dynamics not defined yet.")

"""
    discrete_dynamics!(model, xn, z::AbstractKnotPoint)
    discrete_dynamics!(model, xn, x, u, t, dt)

Evaluate the discrete time dynamics, storing the output in `xn`.
This method is called when using the `InPlace` [`FunctionSignature`](@ref).

Calling `evaluate!` on a [`DiscreteDynamics`](@ref) model will call this function.
"""
discrete_dynamics!(model::DiscreteDynamics, xn, z::AbstractKnotPoint) =
    discrete_dynamics!(model, xn, state(z), control(z), time(z), timestep(z))
discrete_dynamics!(model::DiscreteDynamics, xn, x, u, t, dt) =
    error("In-place discrete dynamics not defined yet.")

"""
    discrete_dynamics!(sig, xn, z::AbstractKnotPoint)

Evaluate the discrete dynamics function, storing the output in `xn`, using the
[`FunctionSignature`](@ref) `sig` to determine which method to call.
"""
discrete_dynamics!(::InPlace, model::DiscreteDynamics, xn, z::AbstractKnotPoint) =
    discrete_dynamics!(model, xn, z)
discrete_dynamics!(::StaticReturn, model::DiscreteDynamics, xn, z::AbstractKnotPoint) =
    xn .= discrete_dynamics(model, z)

jacobian!(model::DiscreteDynamics, J, y, x, u, p) = jacobian!(model, J, y, x, u, p.t, p.dt)
jacobian!(model::DiscreteDynamics, J, y, x, u, t, dt) =
    error("User-defined discrete dynamics Jacobian not defined.")

"""
    dynamics_error(model, z2, z1)

Evaluate the dynamics error between two knot points for a [`DiscreteDynamics`](@ref)
model. In general, this function takes the form:

```math
d(x_{k+1}, u_{k+1}, x_k, u_k) = 0
```

For explicit integration methods of the form ``x_{k+1} = f(x_k, u_k)``, this is just
```math
d(x_k, u_k) - x_{k+1} = 0
```

This is the method that should be used with implicit integration methods.

This form is the out-of-place form that should, in general, return a
`StaticVectors.SVector`. See [`dynamics_error!`](@ref) for the in-place method.

The Jacobian of this function is evaluated by calling [`dynamics_error_jacobian!`](@ref).
"""
dynamics_error(model::DiscreteDynamics, z2::AbstractKnotPoint, z1::AbstractKnotPoint) =
    discrete_dynamics(model, z1) - state(z2)

"""
    dynamics_error!(model, y2, y1, z2, z1)

Evaluate the dynamics error between two knot points for a [`DiscreteDynamics`](@ref) model.
The output is stored in `y2`, and `y1` is provided as a extra input that can be used
to store temporary results. Any specific usage of this variable is left to the user.

See [`dynamics_error`](@ref) for more details on this function.
"""
function dynamics_error!(
    model::DiscreteDynamics,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    discrete_dynamics!(model, y2, z1)
    y2 .-= state(z2)
    return nothing
end

"""
    dynamics_error_jacobian!(sig, diff, model, J2, J1, y2, y1, z2, z1)

Evaluate the Jacobian of [`dynamics_error`](@ref). The derivative with respect to
the first knotpoint `z1` should be stored in `J1`, and the derivative with respect to
the second knotpoint `z2` should be stored in `J2`. The variables `y2` and `y1` should be
vectors of the size of the output dimension (usually the state dimension) and are provided
as extra cache variables whose usage is left to be determined by the user. For example,
a user could choose to evaluate the dynamics error and store the result in one of these
variables to evaluate the error and it's Jacobian at the same time.

Both `J2` and `J1` must have dimensions `(p, n + m)` where `n`, `m`, and `p` are the output
of `state_dim`, `control_dim`, and `output_dim`. Usually `p = n`.

Note that for explicit integration methods of the form `x_{k+1} = f(x_k, u_k)` `J2` should
be equal to `[-I(n) zeros(n,m)]`.

## Implementing on a custom type
To implement this function on a new [`DiscreteDynamics`](@ref) model, define methods
for the [`FunctionSignature`](@ref) and [`DiffMethod`](@ref) of your choice. These are
automatically defined when using [`@autodiff`](@ref) on a [`DiscreteDynamics`](@ref) model.
To implement the `UserDefined` [`DiffMethod`](@ref), implement the following method:

    dynamics_error_jacobian!(model, J2, J1, y2, y1,
                             z2::AbstractKnotPoint, z1::AbstractKnotPoint)

"""
dynamics_error_jacobian!(
    ::FunctionSignature,
    ::UserDefined,
    model::DiscreteDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = dynamics_error_jacobian!(model, J2, J1, y2, y1, z2, z1)

dynamics_error_jacobian!(
    model::DiscreteDynamics,
    J2,
    J1,
    y2,
    y1,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
) = error("User-defined dynamics error Jacobian not defined for $(typeof(model)).")

"""
    propagate_dynamics!(sig, model, z2, z1)

Save the output of either [`discrete_dynamics`](@ref) or [`discrete_dynamics!`](@ref)
evaluated using `z1` into the state vector of `z2`. Useful for simulating discrete systems
forward in time.
"""
function propagate_dynamics!(
    ::InPlace,
    model::DiscreteDynamics,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    discrete_dynamics!(model, state(z2), z1)
    return nothing
end

function propagate_dynamics!(
    ::StaticReturn,
    model::DiscreteDynamics,
    z2::AbstractKnotPoint,
    z1::AbstractKnotPoint,
)
    setstate!(z2, discrete_dynamics(model, z1))
    return nothing
end

ext: jl


text:
module NeuralALU

using Flux
include("nalu.jl")

export NeuralAccumulator

end # module

max_stars_repo_name: dillondaudert/neural-alu

max_issues_repo_name: dillondaudert/neural-alu

max_forks_repo_name: dillondaudert/neural-alu


text:
"""
Hequation.jl

This file contains the function/Jacobian evaluations for
the Chandrasekhar H-equation examples and everything you should
need to run them.

If you only want to run the examples, you should not have to look
at the code.
"""
### And now for the functions ...
"""
function heqJ!(FP,F,x,pdata)

The is the Jacobian evaluation playing by nsol rules. The
precomputed data is a big deal for this one.
"""
function heqJ!(FP::Array{T,2}, F, x, pdata) where {T<:Real}
    pseed = pdata.pseed
    mu = pdata.mu
    n = length(x)
    #
    # Look at the formula in the notebook and you'll see what I did here.
    #
    pmu = pdata.pmu
    Gfix = pdata.gtmp
    @views Gfix .= x - F
    @views Gfix .= -(Gfix .* Gfix .* pmu)
    @views @inbounds for jfp = 1:n
        FP[:, jfp] .= Gfix[:, 1] .* pseed[jfp:jfp+n-1]
        FP[jfp, jfp] = 1.0 + FP[jfp, jfp]
    end
end

"""
heqf!(F,x,pdata)

The function evaluation as per nsold rules.

The precomputed data is a big deal for this example. In particular,
the output pdata.FFB from plan_fft! goes to the fixed point map
computation. Things get very slow if you do not use plan_fft or plan_fft!
"""
function heqf!(F, x, pdata)
    HeqFix!(F, x, pdata)
    #
    # naked BLAS call to fix the allocation blues
    #
    # Using any variation of F.=x-F really hurts
    #
    axpby!(1.0, x, -1.0, F)
end

"""
function HeqFix!(Gfix,x,pdata)
The fixed point map. Gfix goes directly into the function and
Jacobian evaluations for the nonlinear equations formulation.

The precomputed data is a big deal for this example. In particular,
the output pdata.FFA from plan_fft goes to the fixed point map
computation. Things get very slow if you do not use plan_fft.
"""
function HeqFix!(Gfix, x, pdata)
    n = length(x)
    Gfix .= x
    heq_hankel!(Gfix, pdata)
    Gfix .*= pdata.pmu
    Gfix .= 1.0 ./ (1.0 .- Gfix)
end

"""
heqinit(x0::Array{T,1}, c) where T :< Real

Initialize H-equation precomputed data.
"""
function heqinit(x0::Array{T,1}, c) where {T<:Real}
    n = length(x0)
    cval = ones(1)
    cval[1] = c
    vsize = (n,)
    bsize = (2 * n,)
    ssize = (2 * n - 1,)
    FFA = plan_fft(ones(bsize))
    mu = collect(0.5:1:n-0.5)
    pmu = mu * c
    mu = mu / n
    hseed = zeros(ssize)
    for is = 1:2*n-1
        hseed[is] = 1.0 / is
    end
    hseed = (0.5 / n) * hseed
    pseed = hseed
    gtmp = zeros(vsize)
    rstore = zeros(bsize)
    zstore = zeros(bsize) * (1.0 + im)
    hankel = zeros(bsize) * (1.0 + im)
    FFB = plan_fft!(zstore)
    bigseed = zeros(bsize)
    @views bigseed .= [hseed[n:2*n-1]; 0; hseed[1:n-1]]
    @views hankel .= conj(FFA * bigseed)
    return (
        cval = cval,
        mu = mu,
        hseed = hseed,
        pseed = pseed,
        gtmp = gtmp,
        pmu = pmu,
        rstore = rstore,
        zstore = zstore,
        hankel = hankel,
        FFB = FFB,
    )
end

"""
setc!(pdata, cin)

If you are varying c in a computation, this function
lets you set it.
"""
function setc!(pdata, cin)
    c = pdata.cval[1]
    cfix = cin / c
    pdata.pmu .*= cfix
    pdata.cval[1] = cin
end

"""
heq_hankel!(b,pdata)
Multiply an nxn Hankel matrix with seed in R^(2N-1) by a vector b
FFA is what you get with plan_fft before you start computing
"""
function heq_hankel!(b, pdata)
    reverse!(b)
    heq_toeplitz!(b, pdata)
end

"""
heq_toeplitz!(b,pdata)
Multiply an nxn Toeplitz matrix with seed in R^(2n-1) by a vector b
"""
function heq_toeplitz!(b, pdata)
    n = length(b)
    y = pdata.rstore
    y .*= 0.0
    @views y[1:n] = b
    heq_cprod!(y, pdata)
    @views b .= y[1:n]
end

"""
heq_cprod!(b,pdata)
Circulant matrix-vector product with FFT
compute u = C b

Using in-place FFT
"""

function heq_cprod!(b, pdata)
    xb = pdata.zstore
    xb .*= 0.0
    xb .+= b
    pdata.FFB \ xb
    hankel = pdata.hankel
    xb .*= hankel
    pdata.FFB * xb
    b .= real.(xb)
end

"""
Alternative formulation for CI. Tuned to match the paper

author="P. B. Bosma and W. A. DeRooij",
title="Efficient Methods to Calculate Chandrasekhar's H-functions",
journal="Astron. Astrophys.",
volume=126,
year=1983,
pages=283

"""
function heqbos!(F, x, pdata)
    c = pdata
    n = length(x)
    mu = 0.5:1:n-0.5
    mu = mu / n
    h = 1.0 / n
    cval = sqrt(1.0 - c)
    A = zeros(n, n)
    for j = 1:n
        for i = 1:n
            A[i, j] = mu[j] / (mu[i] + mu[j])
        end
    end
    A = (c / 2) * h * A
    F .= (A * x)
    for ig = 1:n
        F[ig] = 1.0 / (cval + F[ig])
    end
    F .= x - F
end


text:
using Test
using Random
using ForwardDiff2
using ForwardDiff2: Dual, Tag, value, partials, dualrun
using ChainRulesCore: Zero

Random.seed!(1234)

using StaticArrays

using Cassette

using NaNMath, SpecialFunctions
using DiffRules
using Calculus

using ForwardDiff2

macro dtest(expr)
    :(@test dualrun(()->$(esc(expr))))
end

macro dtest2(expr)
    :(@test dualrun(()->dualrun(()->$(esc(expr)))))
end

macro dtest3(expr)
    :(@test dualrun(()->dualrun(()->dualrun(()->$(esc(expr))))))
end

macro drun1(expr)
    :(dualrun(()->$(esc(expr))))
end

macro drun2(expr)
    :(dualrun(()->dualrun(()->$(esc(expr)))))
end

macro drun3(expr)
    :(dualrun(()->dualrun(()->dualrun(()->$(esc(expr))))))
end

macro dtest_broken(expr)
    :(@test_broken dualrun(()->$(esc(expr))))
end

const Partials{N,V} = SArray{Tuple{N},V,1,N}

const Tag1 = Tag{Nothing}
const Tag2 = Tag{Tag1}
const Tag3 = Tag{Tag2}

samerng() = MersenneTwister(1)

# By lower-bounding the Int range at 2, we avoid cases where differentiating an
# exponentiation of an Int value would cause a DomainError due to reducing the
# exponent by one
intrand(V) = V == Int ? rand(2:10) : rand(V)

dual_isapprox(::Partials{0,T}, ::Partials{0,S}) where {T,S} = true
dual_isapprox(a, b) = isapprox(a, b)
dual_isapprox(a::Dual{T,T1,T2}, b::Dual{T,T3,T4}) where {T,T1,T2,T3,T4} = dual_isapprox(value(a), value(b)) && all(dual_isapprox.(partials(a), partials(b)))
dual_isapprox(a::Dual{T,T1,T2}, b::Dual{T3,T4,T5}) where {T,T1,T2,T3,T4,T5} = error("Tags don't match")

dual1(primal, partial) = dualrun(()->Dual(primal, partial))
dual2(primal, partial) = dualrun(()->dual1(primal, partial))
dual3(primal, partial) = dualrun(()->dual2(primal, partial))

_mul_partials(p1, p2, v1, v2) = p1 * v1 + p2 * v2
_div_partials(a, b, aval, bval) = _mul_partials(a, b, inv(bval), -(aval / (bval*bval)))

function approx_dtest(x, y)
    x = partials(x)
    if x isa Zero
        @test iszero(y)
    else
        @dtest x[1] ≈ y
    end
end

const Partials{N,V} = SVector{N,V}

for N in (3), M in (4), V in (Int, Float32, Float64)
    println("  ...testing Dual{..,$V,$N} and Dual{..,Dual{..,$V,$M},$N}")

    global PARTIALS = Partials{N,V}(ntuple(n -> intrand(V), N))
    global PRIMAL = intrand(V)
    global FDNUM = dual1(PRIMAL, PARTIALS)

    global PARTIALS2 = Partials{N,V}(ntuple(n -> intrand(V), N))
    global PRIMAL2 = intrand(V)
    global FDNUM2 = dual1(PRIMAL2, PARTIALS2)

    global PARTIALS3 = Partials{N,V}(ntuple(n -> intrand(V), N))
    global PRIMAL3 = intrand(V)
    global FDNUM3 = dual1(PRIMAL3, PARTIALS3)

    global M_PARTIALS = Partials{M,V}(ntuple(m -> intrand(V), M))
    global NESTED_FDNUM = dual2(dual1(PRIMAL, M_PARTIALS), PARTIALS)

    global M_PARTIALS2 = Partials{M,V}(ntuple(m -> intrand(V), M))
    global NESTED_FDNUM2 = dual2(dual1(PRIMAL2, M_PARTIALS2), PARTIALS2)

    ################
    # Constructors #
    ################

    @test Dual(PRIMAL, PARTIALS) === Dual{Nothing}(PRIMAL, PARTIALS)

    @test typeof(NESTED_FDNUM) == Dual{Tag2,Dual{Tag1,V,Partials{M,V}},Partials{N,V}}

    #############
    # Accessors #
    #############

    @test value(PRIMAL) == PRIMAL
    @test value(FDNUM) == PRIMAL
    @test value(NESTED_FDNUM) === dual1(PRIMAL, M_PARTIALS)

    @test partials(FDNUM) == PARTIALS
    @test partials(NESTED_FDNUM) === PARTIALS

    for i in 1:N
        @test partials(FDNUM)[i] == PARTIALS[i]
        for j in 1:M
            @test partials(NESTED_FDNUM) == PARTIALS
        end
    end

    @test ForwardDiff2.npartials(FDNUM) == N
    @test ForwardDiff2.npartials(NESTED_FDNUM) == N

    #####################
    # Generic Functions #
    #####################

    @test FDNUM === copy(FDNUM)
    @test NESTED_FDNUM === copy(NESTED_FDNUM)

    if V != Int
        @test eps(FDNUM) === eps(PRIMAL)
        @test eps(typeof(FDNUM)) === eps(V)
        @test eps(NESTED_FDNUM) === eps(PRIMAL)
        @test eps(typeof(NESTED_FDNUM)) === eps(V)

        @test floor(Int, FDNUM) === floor(Int, PRIMAL)
        @test floor(Int, FDNUM2) === floor(Int, PRIMAL2)
        @test floor(Int, NESTED_FDNUM) === floor(Int, PRIMAL)

        @test floor(FDNUM) === floor(PRIMAL)
        @test floor(FDNUM2) === floor(PRIMAL2)
        @test floor(NESTED_FDNUM) === floor(PRIMAL)

        @test ceil(Int, FDNUM) === ceil(Int, PRIMAL)
        @test ceil(Int, FDNUM2) === ceil(Int, PRIMAL2)
        @test ceil(Int, NESTED_FDNUM) === ceil(Int, PRIMAL)

        @test ceil(FDNUM) === ceil(PRIMAL)
        @test ceil(FDNUM2) === ceil(PRIMAL2)
        @test ceil(NESTED_FDNUM) === ceil(PRIMAL)

        @test trunc(Int, FDNUM) === trunc(Int, PRIMAL)
        @test trunc(Int, FDNUM2) === trunc(Int, PRIMAL2)
        @test trunc(Int, NESTED_FDNUM) === trunc(Int, PRIMAL)

        @test trunc(FDNUM) === trunc(PRIMAL)
        @test trunc(FDNUM2) === trunc(PRIMAL2)
        @test trunc(NESTED_FDNUM) === trunc(PRIMAL)

        @test round(Int, FDNUM) === round(Int, PRIMAL)
        @test round(Int, FDNUM2) === round(Int, PRIMAL2)
        @test round(Int, NESTED_FDNUM) === round(Int, PRIMAL)

        @test round(FDNUM) === round(PRIMAL)
        @test round(FDNUM2) === round(PRIMAL2)
        @test round(NESTED_FDNUM) === round(PRIMAL)

        @test Base.rtoldefault(typeof(FDNUM)) ≡ Base.rtoldefault(typeof(PRIMAL))
        @test Base.rtoldefault(typeof(NESTED_FDNUM)) ≡ Base.rtoldefault(typeof(PRIMAL))
        @dtest Dual(PRIMAL-eps(V), PARTIALS) ≈ FDNUM
    end

    @test hash(FDNUM) === hash(PRIMAL)
    @test hash(FDNUM, hash(PRIMAL)) === hash(PRIMAL, hash(PRIMAL))
    @test hash(NESTED_FDNUM) === hash(PRIMAL)
    @test hash(NESTED_FDNUM, hash(PRIMAL)) === hash(PRIMAL, hash(PRIMAL))

    TMPIO = IOBuffer()
    write(TMPIO, FDNUM)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(FDNUM)) === FDNUM
    seekstart(TMPIO)
    write(TMPIO, FDNUM2)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(FDNUM2)) === FDNUM2
    seekstart(TMPIO)
    write(TMPIO, NESTED_FDNUM)
    seekstart(TMPIO)
    @test read(TMPIO, typeof(NESTED_FDNUM)) === NESTED_FDNUM
    close(TMPIO)

    @test zero(NESTED_FDNUM) === dual2(dual1(zero(PRIMAL), zero(M_PARTIALS)), zero(PARTIALS))
    @test zero(typeof(NESTED_FDNUM)) === zero(NESTED_FDNUM)
    @test zero(FDNUM) === dual1(zero(PRIMAL), zero(PARTIALS))
    @test zero(typeof(FDNUM)) === zero(FDNUM)

    @test one(NESTED_FDNUM) === dual2(dual1(one(PRIMAL), zero(M_PARTIALS)), zero(PARTIALS))
    @test one(typeof(NESTED_FDNUM)) === one(NESTED_FDNUM)
    @test one(FDNUM) === dual1(one(PRIMAL), zero(PARTIALS))
    @test one(typeof(FDNUM)) === one(FDNUM)

    if V <: Integer
        @test rand(samerng(), FDNUM) == rand(samerng(), value(FDNUM))
        @test rand(samerng(), NESTED_FDNUM) == rand(samerng(), value(NESTED_FDNUM))
    elseif V <: AbstractFloat
        for f in (rand, randn, randexp)
            @test f(samerng(), typeof(NESTED_FDNUM)) === dual2(dual1(f(samerng(), V), zero(Partials{M,V})), zero(Partials{N,V}))
            @test f(samerng(), typeof(FDNUM)) === dual1(f(samerng(), V), zero(Partials{N,V}))
        end
    end

    # Predicates #
    #------------#

    @test ForwardDiff2.isconstant(zero(FDNUM))
    @test ForwardDiff2.isconstant(one(FDNUM))
    @test ForwardDiff2.isconstant(FDNUM) == (N == 0)

    @test ForwardDiff2.isconstant(zero(NESTED_FDNUM))
    @test ForwardDiff2.isconstant(one(NESTED_FDNUM))
    @test ForwardDiff2.isconstant(NESTED_FDNUM) == (N == 0)

    @dtest isequal(FDNUM, dual1(PRIMAL, PARTIALS2))
    @dtest isequal(PRIMAL, PRIMAL2) == isequal(FDNUM, FDNUM2)

    @dtest isequal(NESTED_FDNUM, dual2(dual1(PRIMAL, M_PARTIALS2), PARTIALS2))
    @dtest isequal(PRIMAL, PRIMAL2) == isequal(NESTED_FDNUM, NESTED_FDNUM2)

    @dtest FDNUM == dual1(PRIMAL, PARTIALS2)
    @dtest (PRIMAL == PRIMAL2) == (FDNUM == FDNUM2)
    @dtest (PRIMAL == PRIMAL2) == (NESTED_FDNUM == NESTED_FDNUM2)

    @dtest isless(dual1(1, PARTIALS), dual1(2, PARTIALS2))
    @dtest !(isless(dual1(1, PARTIALS), dual1(1, PARTIALS2)))
    @dtest !(isless(dual1(2, PARTIALS), dual1(1, PARTIALS2)))

    @dtest isless(dual2(dual1(1, M_PARTIALS), PARTIALS), dual2(dual1(2, M_PARTIALS2), PARTIALS2))
    @dtest !(isless(dual2(dual1(1, M_PARTIALS), PARTIALS), dual2(dual1(1, M_PARTIALS2), PARTIALS2)))
    @dtest !(isless(dual2(dual1(2, M_PARTIALS), PARTIALS), dual2(dual1(1, M_PARTIALS2), PARTIALS2)))

    @dtest dual1(1, PARTIALS) < dual1(2, PARTIALS2)
    @dtest !(dual1(1, PARTIALS) < dual1(1, PARTIALS2))
    @dtest !(dual1(2, PARTIALS) < dual1(1, PARTIALS2))

    @dtest dual2(dual1(1, M_PARTIALS), PARTIALS) < dual2(dual1(2, M_PARTIALS2), PARTIALS2)
    @dtest !(dual2(dual1(1, M_PARTIALS), PARTIALS) < dual2(dual1(1, M_PARTIALS2), PARTIALS2))
    @dtest !(dual1(dual1(2, M_PARTIALS), PARTIALS) < dual2(dual1(1, M_PARTIALS2), PARTIALS2))

    @dtest dual1(1, PARTIALS) <= dual1(2, PARTIALS2)
    @dtest dual1(1, PARTIALS) <= dual1(1, PARTIALS2)
    @dtest !(dual1(2, PARTIALS) <= dual1(1, PARTIALS2))

    @dtest dual2(dual1(1, M_PARTIALS), PARTIALS) <= dual2(dual1(2, M_PARTIALS2), PARTIALS2)
    @dtest dual2(dual1(1, M_PARTIALS), PARTIALS) <= dual2(dual1(1, M_PARTIALS2), PARTIALS2)
    @dtest !(dual2(dual1(2, M_PARTIALS), PARTIALS) <= dual2(dual1(1, M_PARTIALS2), PARTIALS2))

    @dtest dual1(2, PARTIALS) > dual1(1, PARTIALS2)
    @dtest !(dual1(1, PARTIALS) > dual1(1, PARTIALS2))
    @dtest !(dual1(1, PARTIALS) > dual1(2, PARTIALS2))

    @dtest dual2(dual1(2, M_PARTIALS), PARTIALS) > dual2(dual1(1, M_PARTIALS2), PARTIALS2)
    @dtest !(dual2(dual1(1, M_PARTIALS), PARTIALS) > dual2(dual1(1, M_PARTIALS2), PARTIALS2))
    @dtest !(dual2(dual1(1, M_PARTIALS), PARTIALS) > dual2(dual1(2, M_PARTIALS2), PARTIALS2))

    @dtest dual1(2, PARTIALS) >= dual1(1, PARTIALS2)
    @dtest dual1(1, PARTIALS) >= dual1(1, PARTIALS2)
    @dtest !(dual1(1, PARTIALS) >= dual1(2, PARTIALS2))

    @dtest dual2(dual1(2, M_PARTIALS), PARTIALS) >= dual2(dual1(1, M_PARTIALS2), PARTIALS2)
    @dtest dual2(dual1(1, M_PARTIALS), PARTIALS) >= dual2(dual1(1, M_PARTIALS2), PARTIALS2)
    @dtest !(dual2(dual1(1, M_PARTIALS), PARTIALS) >= dual2(dual1(2, M_PARTIALS2), PARTIALS2))

    @test isnan(dual1(NaN, PARTIALS))
    @test !(isnan(FDNUM))

    @test isnan(dual2(dual1(NaN, M_PARTIALS), PARTIALS))
    @test !(isnan(NESTED_FDNUM))

    @test isfinite(FDNUM)
    @test !(isfinite(dual1(Inf, PARTIALS)))

    @test isfinite(NESTED_FDNUM)
    @test !(isfinite(dual2(dual1(NaN, M_PARTIALS), PARTIALS)))

    @test isinf(dual1(Inf, PARTIALS))
    @test !(isinf(FDNUM))

    @test isinf(dual2(dual1(Inf, M_PARTIALS), PARTIALS))

    @test isreal(FDNUM)

    @test isinteger(dual1(1.0, PARTIALS))
    @test isinteger(FDNUM) == (V == Int)

    @test isinteger(dual2(dual1(1.0, M_PARTIALS), PARTIALS))
    @test isinteger(NESTED_FDNUM) == (V == Int)

    @test iseven(dual1(2, 1))
    @test !(iseven(dual1(1, 1)))

    @test iseven(dual2(dual1(2, 1), 1))
    @test !(iseven(dual2(dual1(1, 1), 1)))

    @test isodd(dual1(1, 1))
    @test !(isodd(dual1(2, 1)))

    @test isodd(dual2(dual1(1, 1), 1))
    @test !(isodd(dual2(dual1(2, 1), 1)))

    ##############
    # Arithmetic #
    ##############

    # Addition/Subtraction #
    #----------------------#

    @dtest FDNUM + FDNUM2 === Dual(value(FDNUM) + value(FDNUM2), partials(FDNUM) + partials(FDNUM2))
    @dtest FDNUM + PRIMAL === Dual(value(FDNUM) + PRIMAL, partials(FDNUM))
    @dtest PRIMAL + FDNUM === Dual(value(FDNUM) + PRIMAL, partials(FDNUM))
    @test @drun2(NESTED_FDNUM + NESTED_FDNUM2) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) + value(NESTED_FDNUM2), partials(NESTED_FDNUM) + partials(NESTED_FDNUM2))
    @test @drun2(NESTED_FDNUM + PRIMAL) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) + PRIMAL, partials(NESTED_FDNUM))
    @test @drun2(PRIMAL + NESTED_FDNUM) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) + PRIMAL, partials(NESTED_FDNUM))

    @dtest FDNUM - FDNUM2 === Dual(value(FDNUM) - value(FDNUM2), partials(FDNUM) - partials(FDNUM2))
    @dtest FDNUM - PRIMAL === Dual(value(FDNUM) - PRIMAL, partials(FDNUM))
    @dtest PRIMAL - FDNUM === Dual(PRIMAL - value(FDNUM), -(partials(FDNUM)))
    @dtest -(FDNUM) === Dual(-(value(FDNUM)), -(partials(FDNUM)))

    @test @drun2(NESTED_FDNUM - NESTED_FDNUM2) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) - value(NESTED_FDNUM2), partials(NESTED_FDNUM) - partials(NESTED_FDNUM2))
    @test @drun2(NESTED_FDNUM - PRIMAL) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) - PRIMAL, partials(NESTED_FDNUM))
    @test @drun2(PRIMAL - NESTED_FDNUM) === @drun1 Dual{Tag2}(PRIMAL - value(NESTED_FDNUM), -(partials(NESTED_FDNUM)))
    @test @drun2(-(NESTED_FDNUM)) === @drun1 Dual{Tag2}(-(value(NESTED_FDNUM)), -(partials(NESTED_FDNUM)))

    # Multiplication #
    #----------------#

    @test dual_isapprox(@drun1(FDNUM * FDNUM2), Dual{Tag1}(value(FDNUM) * value(FDNUM2), _mul_partials(partials(FDNUM), partials(FDNUM2), value(FDNUM2), value(FDNUM))))
    @test @drun1(FDNUM * PRIMAL) === Dual{Tag1}(value(FDNUM) * PRIMAL, partials(FDNUM) * PRIMAL)
    @test @drun1(PRIMAL * FDNUM) === Dual{Tag1}(value(FDNUM) * PRIMAL, partials(FDNUM) * PRIMAL)

    @test @drun2(NESTED_FDNUM * NESTED_FDNUM2) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) * value(NESTED_FDNUM2), _mul_partials(partials(NESTED_FDNUM), partials(NESTED_FDNUM2), value(NESTED_FDNUM2), value(NESTED_FDNUM)))
    @test @drun2(NESTED_FDNUM * PRIMAL) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) * PRIMAL, partials(NESTED_FDNUM) * PRIMAL)
    @test @drun2(PRIMAL * NESTED_FDNUM) === @drun1 Dual{Tag2}(value(NESTED_FDNUM) * PRIMAL, partials(NESTED_FDNUM) * PRIMAL)

    # Division #
    #----------#

    #=
    #TODO: what are these tests???
    if M > 0 && N > 0
        @dtest Dual{1}(FDNUM) / Dual{1}(PRIMAL) === Dual{1}(FDNUM / PRIMAL)
        @dtest Dual{1}(PRIMAL) / Dual{1}(FDNUM) === Dual{1}(PRIMAL / FDNUM)
        @dtest_broken Dual{1}(FDNUM) / FDNUM2 === Dual{1}(FDNUM / FDNUM2)
        @dtest_broken FDNUM / Dual{1}(FDNUM2) === Dual{1}(FDNUM / FDNUM2)
        # following may not be exact, see #264
        @dtest Dual{1}(FDNUM / PRIMAL, FDNUM2 / PRIMAL) ≈ Dual{1}(FDNUM, FDNUM2) / PRIMAL
    end
    =#

    @test dual_isapprox(@drun1(FDNUM / FDNUM2), dual1(value(FDNUM) / value(FDNUM2), _div_partials(partials(FDNUM), partials(FDNUM2), value(FDNUM), value(FDNUM2))))
    @test dual_isapprox(@drun1(FDNUM / PRIMAL), dual1(value(FDNUM) / PRIMAL, partials(FDNUM) / PRIMAL))
    @test dual_isapprox(@drun1(PRIMAL / FDNUM), dual1(PRIMAL / value(FDNUM), (-(PRIMAL) / value(FDNUM)^2) * partials(FDNUM)))

    @test dual_isapprox(@drun2(NESTED_FDNUM / NESTED_FDNUM2), @drun1 Dual{Tag2}(value(NESTED_FDNUM) / value(NESTED_FDNUM2), _div_partials(partials(NESTED_FDNUM), partials(NESTED_FDNUM2), value(NESTED_FDNUM), value(NESTED_FDNUM2))))
    @test dual_isapprox(@drun2(NESTED_FDNUM / PRIMAL), @drun1 Dual{Tag2}(value(NESTED_FDNUM) / PRIMAL, partials(NESTED_FDNUM) / PRIMAL))
    @test dual_isapprox(@drun2(PRIMAL / NESTED_FDNUM), @drun1 Dual{Tag2}(PRIMAL / value(NESTED_FDNUM), (-(PRIMAL) / value(NESTED_FDNUM)^2) * partials(NESTED_FDNUM)))

    # Exponentiation #
    #----------------#

    @dtest dual_isapprox(FDNUM^FDNUM2, exp(FDNUM2 * log(FDNUM)))
    @dtest dual_isapprox(FDNUM^PRIMAL, exp(PRIMAL * log(FDNUM)))
    @dtest dual_isapprox(PRIMAL^FDNUM, exp(FDNUM * log(PRIMAL)))

    @dtest2 dual_isapprox(NESTED_FDNUM^NESTED_FDNUM2, exp(NESTED_FDNUM2 * log(NESTED_FDNUM)))
    @dtest2 dual_isapprox(NESTED_FDNUM^PRIMAL, exp(PRIMAL * log(NESTED_FDNUM)))
    @dtest2 dual_isapprox(PRIMAL^NESTED_FDNUM, exp(NESTED_FDNUM * log(PRIMAL)))

    @dtest_broken partials(NaNMath.pow(Dual(-2.0, 1.0), Dual(2.0, 0.0)))[1] == -4.0

    ###################################
    # General Mathematical Operations #
    ###################################

    @dtest conj(FDNUM) == FDNUM
    @dtest conj(NESTED_FDNUM) === NESTED_FDNUM

    @dtest transpose(FDNUM) === FDNUM
    @dtest transpose(NESTED_FDNUM) === NESTED_FDNUM

    @dtest abs(-FDNUM) === FDNUM
    @dtest abs(FDNUM) === FDNUM
    @test @drun2 abs(-NESTED_FDNUM) === NESTED_FDNUM
    @test @drun2 abs(NESTED_FDNUM) === NESTED_FDNUM

    if V != Int
        for (M, f, arity) in DiffRules.diffrules()
            in(f, (:hankelh1, :hankelh1x, :hankelh2, :hankelh2x, :/, :rem2pi)) && continue
            println("       ...auto-testing $(M).$(f) with $arity arguments")
            if arity == 1
                deriv = DiffRules.diffrule(M, f, :x)
                modifier = in(f, (:asec, :acsc, :asecd, :acscd, :acosh, :acoth)) ? one(V) : zero(V)
                @eval begin
                    x = rand() + $modifier
                    dx = dualrun(()->$M.$f(Dual(x, one(x))))
                    @dtest value(dx) == $M.$f(x)
                    approx_dtest(dx, $deriv)
                end
            elseif arity == 2
                derivs = DiffRules.diffrule(M, f, :x, :y)
                @eval begin
                    x, y = rand(1:10), rand()
                    dx = dualrun(()->$M.$f(Dual(x, one(x)), y))
                    dy = dualrun(()->$M.$f(x, Dual(y, one(y))))
                    actualdx = $(derivs[1])
                    actualdy = $(derivs[2])
                    @dtest value(dx) == $M.$f(x, y)
                    @dtest value(dy) == value(dx)
                    if isnan(actualdx)
                        @dtest isnan(partials(dx)[1])
                    else
                        approx_dtest(dx, actualdx)
                    end
                    if isnan(actualdy)
                        @dtest isnan(partials(dy)[1])
                    else
                        approx_dtest(dy, actualdy)
                    end
                end
            end
        end
    end

    # Special Cases #
    #---------------#

    @test_broken dual_isapprox(hypot(FDNUM, FDNUM2, FDNUM), sqrt(2*(FDNUM^2) + FDNUM2^2))
    @test_broken dual_isapprox(hypot(FDNUM, FDNUM2, FDNUM3), sqrt(FDNUM^2 + FDNUM2^2 + FDNUM3^2))

    @dtest all(map(dual_isapprox, value(sincos(FDNUM)), (sin(FDNUM), cos(FDNUM))))

    if V === Float32
        @dtest typeof(sqrt(FDNUM)) === typeof(FDNUM)
    end

    for f in (fma, muladd)
        @dtest dual_isapprox(f(FDNUM, FDNUM2, FDNUM3),   Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PRIMAL2*PARTIALS + PARTIALS3))
        @dtest dual_isapprox(f(FDNUM, FDNUM2, PRIMAL3),  Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PRIMAL2*PARTIALS))
        @dtest dual_isapprox(f(PRIMAL, FDNUM2, FDNUM3),  Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2 + PARTIALS3))
        @dtest dual_isapprox(f(PRIMAL, FDNUM2, PRIMAL3), Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL*PARTIALS2))
        @dtest dual_isapprox(f(FDNUM, PRIMAL2, FDNUM3),  Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL2*PARTIALS + PARTIALS3))
        @dtest dual_isapprox(f(FDNUM, PRIMAL2, PRIMAL3), Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PRIMAL2*PARTIALS))
        @dtest dual_isapprox(f(PRIMAL, PRIMAL2, FDNUM3), Dual(f(PRIMAL, PRIMAL2, PRIMAL3), PARTIALS3))
    end
end

@testset "Exponentiation of zero" begin
    x0 = 0.0
    x1 = Dual{Tag1}(x0, 1.0)
    x2 = Dual{Tag2}(x1, 1.0)
    x3 = Dual{Tag3}(x2, 1.0)
    pow = ^  # to call non-literal power
    @test_skip @dtest3 pow(x3, 2) === x3^2 == x3 * x3 # this triggers a Julia internal error
    @dtest2 pow(x2, 1) === x2^1 == x2
    @dtest pow(x1, 0) === x1^0 === Dual(1.0, 0.0)
end


text:
using GBDTs
using Random,Test

let
    X = [true, true, true, false, false]
    y = [x ? 1 : 2 for x in X]

    grammar = @grammar begin
        b = x | true | false
    end

    Random.seed!(0)
    p = MonteCarlo(10,5)
    members = collect(1:length(X))
    model = induce_tree(grammar, :b, p, X, y, 2)
    @test length(model) == 3
    @test id(model.tree) == 1
    @test label(model.tree) == 1
    @test !isleaf(model.tree)
    leafs = collect(Leaves(model.tree))
    @test length(leafs) == 2
    @test sort([n.label for n in leafs]) == [1, 2]
    y_bool = partition(X, members, :x, Main)
    @test X == y_bool
    members_true, members_false = members_by_bool(members, y_bool)
    @test members_true == findall(X)
    @test members_false == findall(!, X)
    y_pred = classify(model, X)
    @test y_pred == y
    mvec = node_members(model, X, members)
    @test mvec[1] == members
    @test mvec[2] == members_true
    @test mvec[3] == members_false
    @test model.tree == model[1]
    @test children_id(model.tree) == [2,3]
    @test get_expr(nothing) == nothing
    @test get_expr(model[1]) == :x
end


text:
plasma_frequency(n, q, m) = sqrt(4*pi*n*q^2 / m)
gyro_frequency(q, B, m) = q*B / (m * 3e10)
electron_collision_time(n, T, q, m, l) = 3*sqrt(m)*T^(3/2) / (4*sqrt(2*pi)*l*n*q^4)
ion_collision_time(n, T, q, m, l) = 3*sqrt(m)*T^(3/2) / (4*sqrt(pi)*l*n*q^4)
sound_speed(T, m) = sqrt(T / m)
plasma_beta(n, T, B) = 4*pi*n*T / B^2

#NOTE: input cgs units, except T0, input eV for T0
function dimensionless_parameters(a, R0, n0, T0, B0)

    # constants
    c = 3.00e10
    q = 4.8032e-10
    me = 9.1094e-28
    mp = 1.6726219e-24

    # convert T0 from eV to erg
    T0 = T0 * 1.60218e-12

    # assume deuterium ions
    mi = 2 * mp

    # important parameters
    wpe = plasma_frequency(n0, q, me)
    wci = gyro_frequency(q, B0, mi)
    tau_i = ion_collision_time(n0, T0, q, mi, 1)
    tau_e = electron_collision_time(n0, T0, q, me, 1)
    beta0 = plasma_beta(n0, T0, B0)
    cs0 = sound_speed(T0, mi)
    t0 = sqrt(a*R0/2) / cs0

    # parameters for gdb model
    am = R0*beta0/a
    ad = cs0^2 * t0 / (wci*a^2)
    ki = 2*3.9*t0*tau_i*T0 / (R0^2 * mi)
    ke = 2*3.2*t0*tau_e*T0 / (R0^2 * me)
    er = 2*a/R0
    eg = 0.08*tau_i / t0
    ev = cs0*t0/R0
    de2 = c / (a*wpe)
    eta = 0.51*t0*de2 / tau_e

    return Float64[am, ad, ki, ke, er, eg, ev, de2, eta]
end

function partial_s(f, f_x, f_y, psi_x, psi_y, Ds, am)

    return Ds*f + am*(psi_x .* f_y - psi_y .* f_x)
end

function partial_ss(f, f_s, f_x, f_y, f_xy, f_xx, f_yy, psi_x, psi_y, psi_xy, psi_xx, psi_yy, Dss, Ds, Dx, Dy, am)

    _a = Dss*f + am*Ds*(psi_x.*f_y - psi_y.*f_x) + am*(psi_x.*Dy*f_s - psi_y.*Dx*f_s)
    _b = (f_x.*psi_xy + f_xx.*psi_y - f_y.*psi_xx - psi_x.*f_xy).*psi_y
    _c = (f_x.*psi_yy - f_y.*psi_xy - f_yy.*psi_x + psi_y.*f_xy).*psi_x

    return _a + am^2 * (_b - _c)
end

function lnn_total_ion_derivative(phi_c, Pe_c, n, j_s, u_s, er, ev, ad)

    _a = er*phi_c - er*ad*Pe_c ./ n
    _b = er*ad*j_s ./ n - ev*u_s
    return _a + _b
end

function lnTe_total_electron_derivative_adiabatic(lnn_Dt, j, lnn_s, n, Te, lnTe_c, er, ad)

    _a = (2/3)*(lnn_Dt + er*ad*j.*lnn_s./n)
    _b = -(5/3)*er*ad*(Te .* lnTe_c)

    return _a + _b
end

function lnTi_total_ion_derivative_adiabatic(lnn_Dt, Ti, lnTi_c, er, ad)

    _a = (2/3)*lnn_Dt
    _b = (5/3)*er*ad*(Ti .* lnTi_c)

    return _a + _b
end

# NOTE: this one is written to be species-ambiguous since it applies to Te and Ti
function lnT_partial_thermal(lnT_s, lnT_ss, T, n, k)

    _a = (7/2)*lnT_s.^2 + lnT_ss
    return k*T.^(5/2) .* _a ./ n
end

function u_total_ion_derivative(Pe_s, Pi_s, G_s, n, Ti, u_c, ev, eg, er, ad)

    _a = -ev*(Pe_s + Pi_s) ./ n
    _b = -4*ev*eg*G_s ./ n
    _c = er*ad*(Ti .* u_c)

    return +(_a, _b, _c)
end

function general_partial_derivative(f_Dt, f_x, f_y, f_s, phi_x, phi_y, v)

    return f_Dt - phi_x.*f_y + phi_y.*f_x - v.*f_s
end

function w_partial_t(n, lnn_x, lnn_y, lnn_xy, lnn_xx, lnn_yy,
                     Pe_x, Pe_y, Pe_xy, Pe_xx, Pe_yy, Pe_xxx, Pe_yyy, Pe_xxy, Pe_xyy, Pe_c,
                     phi_x, phi_y, phi_xy, phi_xx, phi_yy, phi_xxx, phi_yyy, phi_xxy, phi_xyy,
                     G_c, Pi_c, j_s, ad, eg)

    _a = (lnn_x.*lnn_y + lnn_xy).*(Pe_y.*phi_y - Pe_x.*phi_x)
    _b = (lnn_x.^2 + lnn_xx).*Pe_x.*phi_y - (lnn_y.^2 + lnn_yy).*Pe_y.*phi_x
    _c = -Pe_x.*lnn_x.^2 .*phi_y + Pe_x.*lnn_x.*lnn_y.*phi_x
    _d = Pe_x.*lnn_x.*phi_xy - Pe_x.*lnn_y.*phi_xx + Pe_xx.*lnn_x.*phi_y
    _e = -Pe_xx.*lnn_y.*phi_x - Pe_xx.*phi_xy - Pe_xxx.*phi_y
    _f = -Pe_y.*lnn_x.*lnn_y.*phi_y + Pe_y.*lnn_x.*phi_yy
    _g = Pe_y.*lnn_y.^2 .*phi_x - Pe_y.*lnn_y.*phi_xy + Pe_yy.*lnn_x.*phi_y
    _h = -Pe_yy.*lnn_y.*phi_x + Pe_yy.*phi_xy + Pe_yyy.*phi_x + phi_x.*Pe_xxy
    _i = phi_xx.*Pe_xy - ad.*phi_y.*Pe_xyy - phi_yy.*Pe_xy
    _j = lnn_x.*phi_x.*phi_xy - lnn_x.*phi_xx.*phi_y + lnn_y.*phi_x.*phi_yy - lnn_y.*phi_y.*phi_xy
    _k = -lnn_x.*phi_xx.*phi_y + lnn_y.*phi_x.*phi_yy - lnn_y.*phi_y.*phi_xy
    _l = phi_x.*phi_yyy + phi_x.*phi_xxy - phi_xxx.*phi_y - phi_y.*phi_xyy

    _A = ad * +(_a, _b, _c, _d, _e, _f, _g, _h, _i)
    _B = n .* +(_j, _k, _l)

    return _A + _B - eg*G_c - Pe_c - Pi_c + j_s
end

function A_partial_t(phi_x, phi_y, phi_s, j, jn, jn_x, jn_y, jn_s, Pe_s, u, n, Te, de2, ev, ad, am, eta)

    _a = de2 * (phi_x .* jn_y - phi_y .* jn_x)
    _b = de2 * (u - jn / ev).*jn_s / ev
    _c = (phi_s - ad*Pe_s ./ n) / am
    _d = eta*(j ./ Te.^(3/2))

    return +(_a, _b, _c, _d)
end

function Abohm_def(Te, Ti, n, phi, er, ev, ad, de2)

    cs = sqrt.(Te + Ti)
    lambda = 2.695
    jbohm = ev * n.*cs.*(1 .- exp.(lambda .- phi ./(ad*Te))) / (ad*er)

    return -de2*jbohm
end

function G_def(Ti, phi_c, Pi_c, u_s, er, ev, ad)

    return Ti.^(5/2) .* (er*(phi_c + ad*Pi_c) - 4*ev*u_s)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

analyze stepper benchmark results
see stepper/README.md
=#

using CSV
using DataFrames
using Printf
import BenchmarkProfiles

keep_set = "various"

enhancements = [
    "basic",
    "prox",
    "TOA",
    "curve",
    "comb",
    ]

compare_pairs = [
    ["basic", "prox"],
    ["prox", "toa"],
    ["toa", "curve"],
    ["curve", "comb"],
    ["basic", "comb"],
    ]

# geomean shifts
time_shift = 1e-3
total_shift = 1e-4
piter_shift = 1e-5

# file locations
bench_file = joinpath(@__DIR__, "raw", "bench.csv")
output_dir = mkpath(joinpath(@__DIR__, "analysis"))
tex_dir = mkpath(joinpath(output_dir, "tex"))
stats_dir = mkpath(joinpath(output_dir, "stats"))
csv_dir = mkpath(joinpath(output_dir, "csvs"))

# get extra info about runs; uses hardcoded enhancement names
function extra_stats(all_df)
    all_df = transform(all_df, [:n, :p, :q] => ((x, y, z) -> x .+ y .+ z) => :npq)

    # basic or comb converged
    for enh in ("basic", "comb")
        enh_conv = filter(t -> ((t.enhancement == enh) && t.conv), all_df)
        enh_data = select(enh_conv,
            :npq, :iters, :solve_time,
            :iters => ByRow(log10) => :log_iters,
            :solve_time => ByRow(log10) => :log_solve_time,
            [:time_uprhs, :solve_time] => ((x, y) -> x ./ y) => :prop_rhs,
            )
        CSV.write(joinpath(csv_dir, enh * "conv.csv"), enh_data)
    end

    # basic/comb both converged
    two_solver = filter(t -> (t.enhancement in ("basic", "comb")), all_df)
    two_solver = combine(groupby(two_solver, :inst_key), names(all_df),
        :conv => all => :two_conv)
    two_solver_conv = filter(t -> t.two_conv, two_solver)
    rel_impr = (x -> (x[1] - x[2]) / x[1])
    two_solver_conv = combine(groupby(two_solver_conv, :inst_key),
        [:enhancement, :solve_time], :solve_time => rel_impr => :time_impr,
        [:enhancement, :iters], :iters => rel_impr => :iters_impr,
        )
    filter!(t -> (t.enhancement == "comb"), two_solver_conv)
    CSV.write(joinpath(csv_dir, "basiccombconv.csv"),
        select(two_solver_conv, :solve_time, :time_impr, :iters, :iters_impr))

    # logged stats for instances
    basic_df = filter(t -> (t.enhancement == "basic"), all_df)
    CSV.write(joinpath(csv_dir, "inst_stats.csv"), select(basic_df,
        :num_cones => ByRow(log10) => :log_numcones,
        :npq => ByRow(log10) => :log_npq,
        ),)

    # update examplestats.csv with unique cones and instance count for each example
    exstats = open(joinpath(stats_dir, "examplestats.csv"), "w")
    println(exstats, "example,num_instances,cones")
    for g in groupby(basic_df, :example, sort = true)
        ex = first(g).example
        n = nrow(g)
        cone_lists = [eval(Meta.parse.(x)) for x in g[!, :cone_types]]
        cones = unique(vcat(cone_lists...))
        println(exstats, "$ex,$n,$cones")
    end
    close(exstats)

    # count instances with non-default solver options
    nondef = filter!(!ismissing, basic_df[!, :nondefaults])
    count_nondefault = length(nondef)
    @info("$count_nondefault instances have non-default solver options")
    if !iszero(count_nondefault)
        @info("unique non-default solver options are:")
        for nd in unique(nondef)
            println(nd)
        end
    end

    # count number of instances in "every"
    count_every = nrow(unique(select(filter(t -> t.every_conv, all_df),
        :inst_key)))
    @info("$count_every instances are in the every set")

    return
end

# generic functions without hardcoded enhancement names

function post_process()
    all_df = preprocess_df()
    if isempty(all_df)
        error("input CSV has no $keep_set set instances with a full list of runs")
    end

    make_agg_tables(all_df)
    make_subtime_tables(all_df)

    for comp in compare_pairs, metric in [:solve_time, :iters]
        make_perf_profiles(all_df, comp, metric)
    end
    extra_stats(all_df)

    return
end

process_entry(x::Float64) = @sprintf("%.2f", x)

process_entry(x::Int) = string(x)

get_cap(x::AbstractVector, y::AbstractVector) = 2 * maximum(x[y])

function shifted_geomean(
    metric::AbstractVector{<:Real},
    conv::AbstractVector{Bool};
    shift::Real = 0,
    cap::AbstractVector = fill(Inf, length(conv)),
    use_cap::Bool = false,
    )
    if use_cap
        @assert length(cap) == length(metric)
        x = copy(metric)
        x[.!conv] .= cap[.!conv]
    else
        x = metric[conv]
    end
    @assert all(isfinite, x)
    return exp(sum(log, x .+ shift) / length(x)) - shift
end

# get enhancement name from solver options
function get_enhancement(x)
    sol_opt = eval(Meta.parse(x))
    @assert length(sol_opt) == 1
    return sol_opt[1]
end

function preprocess_df()
    all_df = CSV.read(bench_file, DataFrame)

    # only keep wanted instance set
    filter!(t -> (t.inst_set == keep_set), all_df)

    # get enhancement name from solver options
    transform!(all_df, :solver_options => ByRow(get_enhancement) => :enhancement)

    # check if any instances could be duplicates:
    possible_dupes = nonunique(all_df, [:enhancement, :example, :inst_data,
        :n, :p, :q, :nu, :cone_types, :num_cones, :max_q])
    if any(possible_dupes)
        df_dupes = unique!(all_df[possible_dupes, [:example, :model_type,
            :inst_set, :inst_num, :inst_data, :extender]])
        println()
        @warn("possible instance/option duplicates detected; inspect instance set " *
            "for duplicates of each of the below (unique) rows:")
        println("\n", df_dupes, "\n")
    end

    # get converged instances, identify by instance key
    ok_status = ["Optimal", "PrimalInfeasible", "DualInfeasible"]
    str_missing(s) = (ismissing(s) ? "" : string(s))
    transform!(all_df,
        :status => ByRow(x -> (!ismissing(x) && x in ok_status)) => :conv,
        # identify instances by example + model_type + inst_data + extender
        [:example, :model_type, :inst_data, :extender] =>
        ((s1, s2, s3, s4) -> s1 .* s2 .* s3 .* str_missing.(s4)) => :inst_key,
        )

    # assumes that nothing returned incorrect status, which is checked manually
    all_df = combine(groupby(all_df, :inst_key), names(all_df),
        :conv => all => :every_conv)

    # assert all instances are solved by at least one stepper (removes noise)
    none_df = combine(groupby(all_df, :inst_key), :inst_key,
        :conv => (x -> !any(x)) => :none_conv)
    @assert isempty(filter(t -> t.none_conv, none_df))

    # assert all instances took at least one iteration
    @assert minimum(all_df[!, :iters]) >= 1

    return all_df
end

function make_agg_tables(all_df)
    # calculate caps for replacing time/iters of unconverged instances
    # use double the largest value for the same instance across all steppers
    # assumes instances are in the same order for all steppers
    max_df = combine(groupby(all_df, :inst_key),
        [:solve_time, :conv] => get_cap => :max_time,
        [:iters, :conv] => get_cap => :max_iters,
        )
    cap(x::Symbol) = max_df[!, x]

    # collect aggregated summary statistics
    df_agg = combine(groupby(all_df, :enhancement),
        # geometric mean over instances where every stepper converged
        [:solve_time, :every_conv] => ((x, y) ->
            shifted_geomean(x, y, shift = time_shift)) => :time_geomean_everyconv,
        [:iters, :every_conv] => ((x, y) ->
            shifted_geomean(x, y, shift = 1)) => :iters_geomean_everyconv,
        # geometric mean over instances where this stepper converged
        [:solve_time, :conv] => ((x, y) ->
            shifted_geomean(x, y, shift = time_shift)) => :time_geomean_thisconv,
        [:iters, :conv] => ((x, y) ->
            shifted_geomean(x, y, shift = 1)) => :iters_geomean_thisconv,
        # geometric mean over all instances, use caps on unconverged isntances
        [:solve_time, :conv] => ((x, y) ->
            shifted_geomean(x, y, cap = cap(:max_time), use_cap = true,
            shift = time_shift)) => :time_geomean_all,
        [:iters, :conv] => ((x, y) ->
            shifted_geomean(x, y, shift = 1, cap = cap(:max_iters),
            use_cap = true)) => :iters_geomean_all,
        :status => (x -> count(isequal("Optimal"), x)) => :optimal,
        :status => (x -> count(isequal("PrimalInfeasible"), x)) => :priminfeas,
        :status => (x -> count(isequal("DualInfeasible"), x)) => :dualinfeas,
        :status => (x -> count(isequal("NumericalFailure"), x)) => :numerical,
        :status => (x -> count(isequal("SlowProgress"), x)) => :slowprogress,
        :status => (x -> count(isequal("TimeLimit"), x)) => :timelimit,
        :status => (x -> count(isequal("IterationLimit"), x)) => :iterationlimit,
        :status => length => :total,
        )

    sort!(df_agg, order(:enhancement, by = (x ->
        findfirst(isequal(x), lowercase.(enhancements)))))
    CSV.write(joinpath(stats_dir, "agg.csv"), df_agg)

    # prepare latex table

    # combine feasible and infeasible statuses
    transform!(df_agg, [:optimal, :priminfeas, :dualinfeas] =>
        ByRow((x...) -> sum(x)) => :converged)

    cols = [:converged, :iters_geomean_everyconv, :iters_geomean_thisconv,
        :iters_geomean_all, :time_geomean_everyconv, :time_geomean_thisconv,
        :time_geomean_all]
    sep = " & "
    tex = open(joinpath(tex_dir, "agg.tex"), "w")

    for i in 1:length(enhancements)
        row_str = enhancements[i]
        for c in cols
            subdf = df_agg[!, c]
            x = (startswith(string(c), "time") ? subdf[i] * 1000 : subdf[i])
            row_str *= sep * process_entry(x)
        end
        row_str *= " \\\\"
        println(tex, row_str)
    end
    close(tex)

    return
end

function make_subtime_tables(all_df)
    divfunc(x, y) = (x ./ y)
    init_cols = [:time_rescale, :time_initx, :time_inity, :time_unproc,
        :time_loadsys]
    transform!(all_df,
        [:time_upsys, :iters] => divfunc => :time_upsys_piter,
        [:time_uprhs, :iters] => divfunc => :time_uprhs_piter,
        [:time_getdir, :iters] => divfunc => :time_getdir_piter,
        [:time_search, :iters] => divfunc => :time_search_piter,
        init_cols => ((x...) -> sum(x)) => :time_init,
        )

    metrics = [:init, :lhs, :rhs, :direc, :search, :lhs_piter,
        :rhs_piter, :direc_piter, :search_piter]
    sets = [:_thisconv, :_everyconv, :_all]

    # calculate caps for replacing time/iters of unconverged instances
    # use double the largest value for the same instance across all steppers
    # assumes instances are in the same order for all steppers
    max_df = combine(groupby(all_df, :inst_key),
        [:time_init, :conv] => get_cap => :max_init,
        [:time_upsys, :conv] => get_cap => :max_upsys,
        [:time_uprhs, :conv] => get_cap => :max_uprhs,
        [:time_getdir, :conv] => get_cap => :max_getdir,
        [:time_search, :conv] => get_cap => :max_search,
        [:time_upsys_piter, :conv] => get_cap => :max_upsys_iter,
        [:time_uprhs_piter, :conv] => get_cap => :max_uprhs_iter,
        [:time_getdir_piter, :conv] => get_cap => :max_getdir_iter,
        [:time_search_piter, :conv] => get_cap => :max_search_iter,
        )
    cap(x::Symbol) = max_df[!, x]

    function get_subtime_df(set, convcol, use_cap)
        subtime_df = combine(groupby(all_df, :enhancement),
            [:time_init, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = total_shift, cap = cap(:max_init),
                use_cap = use_cap)) => Symbol(:init, set),
            [:time_upsys, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = total_shift, cap = cap(:max_upsys),
                use_cap = use_cap)) => Symbol(:lhs, set),
            [:time_uprhs, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = total_shift, cap = cap(:max_uprhs),
                use_cap = use_cap)) => Symbol(:rhs, set),
            [:time_getdir, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = total_shift, cap = cap(:max_getdir),
                use_cap = use_cap)) => Symbol(:direc, set),
            [:time_search, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = total_shift, cap = cap(:max_search),
                use_cap = use_cap)) => Symbol(:search, set),
            [:time_upsys_piter, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = piter_shift,
                cap = cap(:max_upsys_iter), use_cap = use_cap)) =>
                Symbol(:lhs_piter, set),
            [:time_uprhs_piter, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = piter_shift,
                cap = cap(:max_uprhs_iter), use_cap = use_cap)) =>
                Symbol(:rhs_piter, set),
            [:time_getdir_piter, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = piter_shift,
                cap = cap(:max_getdir_iter), use_cap = use_cap)) =>
                Symbol(:direc_piter, set),
            [:time_search_piter, convcol] => ((x, y) ->
                shifted_geomean(x, y, shift = piter_shift,
                cap = cap(:max_search_iter), use_cap = use_cap)) =>
                Symbol(:search_piter, set),
            )

        sort!(subtime_df, order(:enhancement,
            by = (x -> findfirst(isequal(x), lowercase.(enhancements)))))
        CSV.write(joinpath(stats_dir, "subtime" * string(set) * ".csv"),
            subtime_df)

        return subtime_df
    end

    # prepare latex table

    sep = " & "
    for s in sets
        if s == :_thisconv
            convcol = :conv
            use_cap = false
        elseif s == :_everyconv
            convcol = :every_conv
            use_cap = false
        elseif s == :_all
            convcol = :conv
            use_cap = true
        end
        subtime_df = get_subtime_df(s, convcol, use_cap)

        subtime_tex = open(joinpath(tex_dir, "subtime" * string(s) * ".tex"), "w")
        for i in 1:nrow(subtime_df)
            row_str = sep * enhancements[i]
            for m in metrics
                col = Symbol(m, s)
                row_str *= sep * process_entry(subtime_df[i, col] * 1000)
            end
            row_str *= " \\\\"
            println(subtime_tex, row_str)
        end
        close(subtime_tex)
    end

    return
end

function make_perf_profiles(all_df, comp, metric)
    pp = filter(t -> t.enhancement in comp, all_df)

    # BenchmarkProfiles.jl expects NaNs for failures
    select!(pp,
        :inst_key,
        :enhancement,
        [metric, :conv] => ByRow((x, y) -> (y ? x : NaN)) => metric,
        )
    wide_df = unstack(pp, :enhancement, metric)

    (x_plot, y_plot, max_ratio) = BenchmarkProfiles.performance_profile_data(
        Matrix{Float64}(wide_df[!, string.(comp)]), logscale = true)

    # make steps like :steppost in Plots
    for s in 1:2
        x = vcat(0, repeat(x_plot[s], inner = 2))
        y = vcat(0, 0, repeat(y_plot[s][1:(end - 1)], inner = 2), y_plot[s][end])
        CSV.write(joinpath(csv_dir, comp[s] * "_vs_" * comp[2 - s + 1] * "_" *
            string(metric) * ".csv"), DataFrame(x = x, y = y))
    end

    return
end

post_process()
;


text:
module RadioHead
export Coax, computeLength, Unit, toMeters, toFeet, toString, Foot, Meter, partial

abstract type Unit end

mutable struct Meter <: Unit
    value::Float32
end

mutable struct Foot <: Unit
    value::Float32
end

mutable struct Coax
    centerFreq::Float32
    velocityFactor::Float32
    # unit is meters
    wavelengthUnit
    function Coax(centerFreq::Float64, velocityFactor::Float64)
        # centerFreq is center of desired band in mHz
        wlu = (300 / centerFreq) * velocityFactor
        new(centerFreq, velocityFactor, wlu)
    end
end

function partial(divisor::Int64, c::Coax)
    return Meter(c.wavelengthUnit / divisor)
end

function computeLength(approxLength::Unit, c::Coax)
    # approxLength is in meters
    # waveLengths = ceil.(toMeters(approxLength).value ./ c.wavelengthUnit)
    waveLength = ceil(toMeters(approxLength).value / c.wavelengthUnit)
    coaxLength = Meter(waveLength * c.wavelengthUnit)
    tc = typeof(approxLength)
    if tc == Meter
        return coaxLength
    else
        return toFeet(coaxLength)
    end
end

function toMeters(u)
    # convert u to meters
    return Meter(conversion_factor(u, Meter) * u.value)
end

function toFeet(u)
    # convert u to feet
    return Foot(conversion_factor(u, Foot) * u.value)
end

# function convert_unit(u, new_u{T})
#     return T(conversion_factor(u, new_u) * u.value)
# end

function conversion_factor(from_u, to_u)
    x = Dict(
        (Meter, Foot) => 3.2808,
        (Foot, Meter) => 0.3048
    )
    if typeof(from_u) != to_u
        return x[(typeof(from_u), to_u)]
    else
        return 1.0
    end
end

function toString(t::Any)
    tc = typeof(t)
    if tc == Meter
        return string(t.value) * " m"
    elseif tc == Foot
        return string(t.value) * " ft"
    elseif tc == Coax
        return "cf " * string(t.centerFreq) * " vf " * string(t.velocityFactor)
    end
    return ""
end

end # module


text:
function read_binary_vec{T<:Real}(io::IO, ::Type{T})
    vec = T[]
    while !eof(io)
        push!(vec, read(io, T))
    end
    vec
end

function infer_shape{T}(vec::Vector{T}, shape::Tuple{Vararg{Int}})

    #=
    Setting a dimension to -1 will have us infer the other one
    =#

    l = length(vec)

    tot = 1
    dim_to_infer = 0
    for (i,v) in enumerate(shape)
        if v == -1
            @assert(dim_to_infer == 0)
            dim_to_infer = i
        else
            tot *= v
        end
    end

    if dim_to_infer != 0
        @assert(mod(l, tot) == 0)
        dims = collect(shape)
        dims[dim_to_infer] = div(l, tot)
        shape = tuple(dims...)
    end

    shape
end

function _convert_to_column_major_array!{T}(retval::Array{T, 1}, vec::Vector{T}, shape::Tuple{Int})
    copy!(retval, vec)
end
function _convert_to_column_major_array!{T}(retval::Array{T, 2}, vec::Vector{T}, shape::Tuple{Int, Int})
    n, m = shape

    count = 0
    for i in 1 : n
        for j in 1: m
            count += 1
            retval[i,j] = vec[count]
        end
    end

    retval
end
function _convert_to_column_major_array!{T}(retval::Array{T, 3}, vec::Vector{T}, shape::Tuple{Int, Int, Int})
    m, n, o = shape

    count = 0
    for i in 1 : m
        for j in 1 : n
            for k in 1 : o
                count += 1
                retval[i,j, k] = vec[count]
            end
        end
    end

    retval
end
function _convert_to_column_major_array!{T}(retval::Array{T, 4}, vec::Vector{T}, shape::Tuple{Int, Int, Int, Int})
    m, n, o, p = shape

    count = 0
    for i in 1 : m
        for j in 1 : n
            for k in 1 : o
                for l in 1 : p
                    count += 1
                    retval[i,j, k, l] = vec[count]
                end
            end
        end
    end

    retval
end

function convert_to_column_major_array!{T}(retval::Array{T, 1}, vec::Vector{T}, shape::Tuple{Int})
    @assert(shape[1] == length(vec))
    _convert_to_column_major_array(retval, vec, shape)
end
function convert_to_column_major_array!{T}(retval::Array{T, 2}, vec::Vector{T}, shape::Tuple{Int, Int})
    n, m = shape
    @assert(n*m == length(vec))
    @assert(size(retval) == shape)

    _convert_to_column_major_array(retval, vec, shape)
end
function convert_to_column_major_array!{T}(retval::Array{T, 3}, vec::Vector{T}, shape::Tuple{Int, Int, Int})
    m, n, o = shape
    @assert(n*m*o == length(vec))
    @assert(size(retval) == shape)

    _convert_to_column_major_array(retval, vec, shape)
end
function convert_to_column_major_array!{T}(retval::Array{T, 4}, vec::Vector{T}, shape::Tuple{Int, Int, Int, Int})
    m, n, o, p = shape
    @assert(n*m*o*p == length(vec))
    @assert(size(retval) == shape)

    _convert_to_column_major_array(retval, vec, shape)
end

function convert_to_column_major_array{T}(vec::Vector{T}, shape::Tuple{Int})
    @assert(shape[1] == length(vec))
    deepcopy(vec)
end
function convert_to_column_major_array{T}(vec::Vector{T}, shape::Tuple{Int, Int})
    n, m = shape
    @assert(n*m == length(vec))

    retval = Array{T}(n, m)

    _convert_to_column_major_array(retval, vec, shape)
end
function convert_to_column_major_array{T}(vec::Vector{T}, shape::Tuple{Int, Int, Int})
    m, n, o = shape
    @assert(n*m*o == length(vec))

    retval = Array{T}(m, n, o)

    _convert_to_column_major_array(retval, vec, shape)
end
function convert_to_column_major_array{T}(vec::Vector{T}, shape::Tuple{Int, Int, Int, Int})
    m, n, o, p = shape
    @assert(n*m*o*p == length(vec))

    retval = Array{T}(m, n, o, p)

    _convert_to_column_major_array(retval, vec, shape)
end


text:
"""
Two dimensional square lattice.
"""
struct SquareLattice <: AbstractCubicLattice
    L::Int
    sites::Int
    neighs::Matrix{Int} # row = up, right, down, left; col = siteidx
    neighs_cartesian::Array{Int, 3} # row (1) = up, right, down, left; cols (2,3) = cartesian siteidx
    lattice::Matrix{Int}
    bond_lattice::Array{Int, 3}
    bond_info::Matrix{Int}# num bond,src, trg,bond neighs(4)
    # for generic checkerboard decomposition
    # OPT: implement Assaad's two-group square lattice version
    n_bonds::Int
    bonds::Matrix{Int} # src, trg, type
    bond_checkerboard::Array{Int, 3}
end

# constructors
"""
    SquareLattice(L::Int)

Create a square lattice with linear dimension `L`.
"""
function SquareLattice(L::Int)
    sites = L^2
    lattice = convert(Array, reshape(1:L^2, (L, L)))
    bond_lattice = convert(Array, reshape(1:2*L^2, (L, L,2)))
    neighs, neighs_cartesian = build_neighbortable(SquareLattice, lattice, L)
     bond_info , bond_checkerboard = build_b_neighbortable(SquareLattice, lattice,bond_lattice, L)
    # for generic checkerboard decomposition
    n_bonds = 2*sites
    bonds = zeros(n_bonds, 3)
    bondid = 1
    for src in lattice
        nup = neighs[1, src]
        bonds[bondid,:] .= [src,nup,0]
        bondid += 1

        nright = neighs[2, src]
        bonds[bondid,:] .= [src,nright,0]
        bondid += 1
    end

    return SquareLattice(L, sites, neighs, neighs_cartesian, lattice,
                bond_lattice, bond_info, n_bonds, bonds,bond_checkerboard)
end

function build_neighbortable(::Type{SquareLattice}, lattice, L)
    up = circshift(lattice,(-1,0))
    right = circshift(lattice,(0,-1))
    down = circshift(lattice,(1,0))
    left = circshift(lattice,(0,1))
    neighs = vcat(up[:]',right[:]',down[:]',left[:]')

    neighs_cartesian = Array{Int, 3}(undef, 4, L, L)
    neighs_cartesian[1,:,:] = up
    neighs_cartesian[2,:,:] = right
    neighs_cartesian[3,:,:] = down
    neighs_cartesian[4,:,:] = left
    return neighs, neighs_cartesian
end

function build_b_neighbortable(::Type{SquareLattice}, lattice, bond_lattice, L)
    num_check=div(L^2,2)
    bond_checkerboard=zeros(Int64,3,num_check,4)# 4 typse of Checkerboards
    bond_index = 1:L^2
    right = circshift(lattice,(0,-1))
    b_ul = bond_lattice[:,:,2]
    b_ur = circshift(bond_lattice[:,:,2],(0,-1))
    b_dl = circshift(bond_lattice[:,:,2],(1,0))
    b_dr = circshift(bond_lattice[:,:,2],(1,-1))
    b_n1 = vcat(bond_index[:]',lattice[:]',right[:]',b_ul[:]',b_ur[:]',b_dl[:]',b_dr[:]')
    for i=1:div(L,2)
        bond_checkerboard[:,(i-1)*L+1:i*L,1]=b_n1[1:3,2*(i-1)*L+1:(2*i-1)*L]#horizontal (right) odd column
        bond_checkerboard[:,(i-1)*L+1:i*L,2]=b_n1[1:3,(2*i-1)*L+1:2*i*L]#horizontal (right) even column
    end
    bond_index = L^2+1:2*L^2
    up = circshift(lattice,(-1,0))
    b_ul = circshift(bond_lattice[:,:,1],(-1,1))
    b_ur = circshift(bond_lattice[:,:,1],(-1,0))
    b_dl = circshift(bond_lattice[:,:,1],(0,1))
    b_dr = bond_lattice[:,:,1]
    b_n2 = vcat(bond_index[:]',lattice[:]',up[:]',b_ul[:]',b_ur[:]',b_dl[:]',b_dr[:]')
    for i=1:num_check
        bond_checkerboard[:,i,3]=b_n2[1:3,2*i-1] #vertical (up) odd
        bond_checkerboard[:,i,4]=b_n2[1:3,2*i] # vertical (up) even
    end
    bond_info=cat(b_n1,b_n2,dims=2)

    return  bond_info, bond_checkerboard
end

@inline nsites(s::SquareLattice) = s.sites
@inline n_bonds(s::SquareLattice) = s.n_bonds
@inline neighbors_lookup_table(s::SquareLattice) = s.neighs
@inline bond_lookup_table(s::SquareLattice) = s.bond_info
@inline bond_checkerboard_table(s::SquareLattice) = s.bond_checkerboard

max_stars_repo_name: yosefcap/MC

max_issues_repo_name: yosefcap/MC

max_forks_repo_name: yosefcap/MC


text:
import JLSO

import Base:
    iterate,
    length,
    eltype,
    getindex,
    show

using FilePathsBase: AbstractPath

using StringDistances:
    findnearest,
    Levenshtein

"""
    FailedTests

Container for failed tests from a [`@quickcheck`](@ref) run. Wrapper around
a `Dict{Symbol, Any}` used for dispatch.
"""
struct FailedTests <: AbstractDict{Symbol, Any}
    dict::Dict{Symbol, Any}
end

"""
    load(io)

Load a collection of failed test cases serialized by a
[`@quickcheck`](@ref) run.
 Argument `io` can be of type `IO`, `AbstractString` or `AbstractPath`.

# Examples
```jldoctest
julia> ft = JCheck.load("JCheck_test.jchk")
2 failing predicates:
Product commute
Is odd
```
"""
load(io::Union{IO, AbstractString, AbstractPath}) =
    FailedTests(JLSO.load(io))

## Iteration interface.

iterate(ft::FailedTests, state::Int = 1) = iterate(ft.dict, state)

eltype(::Type{FailedTests}) = Pair{Symbol, Any}

length(ft::FailedTests) = length(ft.dict)

## Indexing interface.
getindex(ft::FailedTests, key::Symbol) = getindex(ft.dict, key)

getindex(ft::FailedTests, key::AbstractString) =
    getindex(ft.dict, Symbol(key))

"""
    @getcases ft, desc...

Get the predicate with description `desc` and the valuations for which it
failed.

# Note
The predicate with description closest to the one given (in the sense of
the Levenshtein distance) will be returned; there is no need to pass the
exact description.

# Examples
``` jldoctest
julia> ft = JCheck.load("JCheck_test.jchk")
2 failing predicates:
Product commute
Is odd

julia> pred, valuations = @getcases ft iod
NamedTuple{(:predicate, :valuations), Tuple{Function, Vector{Tuple}}}((Serialization.__deserialized_types__.var"#3#4"(), Tuple[(0,), (-9223372036854775808,), (6444904272543528628,)]))

julia> map(x -> pred(x...), valuations)
3-element Vector{Bool}:
 0
 0
 0
```
"""
macro getcases(ft, desc...)
    _ft = esc(ft)
    _desc = join(String.(desc), ' ')

    :(getcases($_ft, $_desc))
end

function getcases(ft::FailedTests, desc::AbstractString)
    desc_matched = first(findnearest(desc,
                                     String.(keys(ft)),
                                     Levenshtein()))

    ft[desc_matched]
end

function show(io::IO, ::MIME"text/plain", ft::FailedTests)
    nbfailed = length(ft)

    header = "$nbfailed failing predicate" *
        (nbfailed > 1 ? "s" : "") *
        ":\n"

    print(io, header * join(String.(keys(ft)), "\n"))
end


text:
import HoTRG_lemon.SpinModule:
	initialize!

type Classical2dSquareLattice{T} <: Classical2dLattice{T}
	spinmodel::ClassicalSpinModel{T}

	tensorT::Array{T,4}
	tensorTtilde::Array{T,4}
 	function Classical2dSquareLattice{T}(spinmodel::ClassicalSpinModel{T})
		this = new{T}()
 		this.spinmodel = spinmodel
 		initialize!(this)
 		return this
 	end
end

Classical2dSquareLattice{T}(spinmodel::ClassicalSpinModel{T}) = Classical2dSquareLattice{T}(spinmodel)

#---
# initialize functions for constructor
function initialize!{T}(lattice::Classical2dSquareLattice{T})
	lattice.tensorT = initializeTensorT2D(lattice.spinmodel)
	lattice.tensorTtilde = initializeTensorT2D(lattice.spinmodel; tilde = true)
end

function initializeTensorT2D{T}(spinmodel::ClassicalSpinModel{T}; tilde::Bool = false)
	if tilde
		measureOperator = getMeasureOperator(spinmodel)
	else
		measureOperator = ones(T,getStates(spinmodel))
	end
	factorW = getFactorW(spinmodel)
	tensorT = makeTensorT(factorW, measureOperator)
	return tensorT
end

#---
# functions for Classical2dSquareLattice
function getCoarserate(lattice::Classical2dSquareLattice)
	return 2
end

#---
# functions for test
function testClassical2dSquareLattice()
	spinmodel = testIsingModel()
	return buildLattice("classical_2d_square", spinmodel)
end


text:
using Test, GasChromatographySimulator, PlutoUI

# general parameters
L = 10.0
d = 0.25e-3
df = 0.25e-6
sp = "SLB5ms"
gas = "He"
time_steps = [0.0, 60.0, 600.0, 300.0]
temp_steps = [40.0, 40.0, 300.0, 300.0]
pin_steps = [200.0, 200.0, 300.0, 300.0].*1000.0 .+ 101300
pout_steps = [0.0, 0.0, 0.0, 0.0].*1000.0
ΔT_steps = [20.0, 30.0, 50.0, 40.0]
x₀_steps = zeros(length(time_steps))
L₀_steps = L.*ones(length(time_steps))
α_steps = zeros(length(time_steps))
db_path = string(@__DIR__, "/data")
db = "Database_test.csv"

@testset "gradient function check" begin
    # Gradient function for d
    a_d = [d, 0.0, 2*L, 0.0]
    d_fun(x) = GasChromatographySimulator.gradient(x, a_d)
    @test d_fun(L) == d/2

    # Gradient function for T with all parameters equal zero
    a_gf = zeros(2,3)
    gf_fun1(x) = GasChromatographySimulator.gradient(x, a_gf)
    @test gf_fun1(0.0) == zeros(size(a_gf)[1])

    # Gradient function with changing α
    αs = [-5.0, 0.0, 0.0, 5.0]
    a_gf = [ΔT_steps x₀_steps L₀_steps αs]
    gf_fun2(x) = GasChromatographySimulator.gradient(x, a_gf)
    @test gf_fun2(L) == -ΔT_steps
end

@testset "parameters check" begin
    # default Options
    opt = GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "inlet", true)
    @test GasChromatographySimulator.Options() == opt

    # Column
    a_d = [d]
    a_df = [df]
    d_fun(x) = GasChromatographySimulator.gradient(x, a_d)
    df_fun(x) = GasChromatographySimulator.gradient(x, a_df)
    col = GasChromatographySimulator.Column(L, d_fun, a_d, df_fun, a_df, sp, gas)
    col_c = GasChromatographySimulator.constructor_System(L, a_d[1], a_df[1], sp, gas)
    @test col.d(col.L) == col_c.d(col_c.L)

    # Program
    a_gf = [ΔT_steps x₀_steps L₀_steps α_steps]
    gf(x) = GasChromatographySimulator.gradient(x, a_gf)
    T_itp = GasChromatographySimulator.temperature_interpolation(time_steps, temp_steps, gf, col.L)
    pin_itp = GasChromatographySimulator.steps_interpolation(time_steps, pin_steps)
    pout_itp = GasChromatographySimulator.steps_interpolation(time_steps, pout_steps)
    prog = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, gf, a_gf, T_itp, pin_itp, pout_itp)
    prog_c = GasChromatographySimulator.constructor_Program(time_steps, temp_steps, pin_steps, pout_steps, ΔT_steps, x₀_steps, L₀_steps, α_steps, opt.Tcontrol, col.L)
    prog_c_ng = GasChromatographySimulator.constructor_Program(time_steps, temp_steps, pin_steps, pout_steps, col.L) # without gradient
    @test prog.T_itp == prog_c.T_itp
    @test prog_c.T_itp != prog_c_ng.T_itp
    @test prog.Fpin_itp == prog_c.Fpin_itp
    @test prog.pout_itp == prog_c.pout_itp

    # conventional program
    TP0 = [40.0, 1.0, 5.0, 280.0, 2.0, 20.0, 320.0, 2.0]
    pP0 = [400000.0, 10.0, 5000.0, 500000.0, 20.0]
    FP0 = [1/(6e7), 60.0]
    ts1, Ts = GasChromatographySimulator.conventional_program(TP0)
    ts2, ps = GasChromatographySimulator.conventional_program(pP0)
    ts3, Fs = GasChromatographySimulator.conventional_program(FP0)
    @test ts1 == [0.0, 60.0, 2880.0, 120.0, 120.0, 120.0]
    @test Fs == 1/(6e7).*ones(2)

    TP = GasChromatographySimulator.temperature_program(ts1, Ts; time_unit="min")
    FP = GasChromatographySimulator.temperature_program(ts3, Fs; time_unit="s")
    @test TP == TP0
    @test FP./[1, 60] == FP0

    cts12 = GasChromatographySimulator.common_time_steps(ts1, ts2)
    cts21 = GasChromatographySimulator.common_time_steps(ts2, ts1)
    @test cts12 == cts21

    cts123 = GasChromatographySimulator.common_time_steps(cts12, ts3)
    cts13 = GasChromatographySimulator.common_time_steps(ts1, ts3)
    new_Fs = GasChromatographySimulator.new_value_steps(Fs, ts3, cts13)
    new_Ts = GasChromatographySimulator.new_value_steps(Ts, ts1, cts123)
    @test new_Fs == 1/(6e7).*ones(length(cts13))
    @test new_Ts == [40.0, 40.0, 85.0, 185.0, 280.0, 280.0, 280.0, 320.0, 320.0, 320.0]

    prog_conv = GasChromatographySimulator.Program([40.0, 1.0, 5.0, 280.0, 2.0, 20.0, 320.0, 2.0], [400000.0, 10.0, 5000.0, 500000.0, 20.0], L; pout="vacuum", time_unit="min")
    prog_conv_s_atm = GasChromatographySimulator.Program([40.0, 1.0*60.0, 5.0/60.0, 280.0, 2.0*60.0, 20.0/60.0, 320.0, 2.0*60.0], [400000.0, 10.0*60.0, 5000.0/60.0, 500000.0, 20.0*60.0], L; pout="atmosphere", time_unit="s")
    @test prog_conv.temp_steps == [40.0, 40.0, 85.0, 185.0, 280.0, 280.0, 280.0, 320.0, 320.0]
    @test prog_conv.time_steps == prog_conv_s_atm.time_steps
    @test prog_conv.pout_steps == prog_conv_s_atm.pout_steps .- 101300

    TP = [40.0, 1.0, 5.0, 200.0, 0.0, 10.0, 280.0, 2.0, 20.0, 320.0, 2.0]
    FpinP = [400000.0, 10.0, 5000.0, 500000.0, 20.0]
    poutP = [0.0, 100.0]
    ΔTP = [0.0, 5.0, 10.0, 60.0, 10.0, 5.0, 80.0, 0.0, -10.0, 0.0, 10.0]
    x₀P = [0.0, 100.0]
    L₀P = [L, 100.0]
    αP = [0.0, 10.0, -0.5, -5.0, 15.0, 1.0, 5.0, 10.0]

    prog_conv_grad = GasChromatographySimulator.Program(TP, FpinP, poutP, ΔTP, x₀P, L₀P, αP, "inlet", L; time_unit="min")
    @test prog_conv_grad.time_steps == [0.0, 60.0, 240.0, 300.0, 60.0, 540.0, 60.0, 240.0, 300.0, 180.0, 120.0, 360.0, 120.0, 120.0, 120.0, 180.0, 300.0, 2700.0]

    # Substance from the load-function-> 1st test
    db_path = string(@__DIR__, "/data")
    db = "Database_test.csv"
    solutes = ["C10", "C11"]
    init_t = zeros(length(solutes))
    init_τ = zeros(length(solutes))
    sub = GasChromatographySimulator.load_solute_database(db_path, db, col.sp, col.gas, solutes, init_t, init_τ)
    @test sub[1].name == "C10"
    @test sub[2].Tchar == 124.746 + 273.15

    # no stationary phase & other gas for diffusivity
    # H2
    sub_0 = GasChromatographySimulator.load_solute_database(db_path, db, "", "H2", solutes, init_t, init_τ)
    @test sub_0[1].Tchar == 1.0
    @test sub_0[2].ann == "no sp"
    @test round(sub_0[1].Dag; sigdigits=5) == 0.00011885
    # N2
    sub_0 = GasChromatographySimulator.load_solute_database(db_path, db, "", "N2", solutes, init_t, init_τ)
    @test round(sub_0[1].Dag; sigdigits=5) == 2.8398e-5
    # Ar
    sub_0 = GasChromatographySimulator.load_solute_database(db_path, db, "", "Ar", solutes, init_t, init_τ)
    @test round(sub_0[1].Dag; sigdigits=5) == 2.5268e-5
    # test for error-cases of GasChromatographySimulator.load_solute_database

    # parameters
    par = GasChromatographySimulator.Parameters(col, prog, sub, opt)
    @test par.col.L == L
    @test par.prog.gf(0.0) == par.prog.gf(L) + par.prog.a_gf[:,1]
    @test par.prog.T_itp(0.0, 0.0) == temp_steps[1] + 273.15
    @test par.prog.T_itp(L, sum(time_steps)) == temp_steps[end] - ΔT_steps[end] + 273.15
    @test par.sub[1].Dag == GasChromatographySimulator.diffusivity(142.28, 10, 22, 0, 0, 0, "He")
    @test par.sub[2].Dag == GasChromatographySimulator.diffusivity(156.31, 11, 24, 0, 0, 0, "He")
end

@testset "solving check" begin
    # test simulation with different settings:
    # - with/without gradient
    # - Tcontrol = "inlet"/"outlet" (with gradient)
    # - odesys = true/false
    # - alg = OwrenZen3/OwrenZen5
    # same settings for:
    # - the Column
    # - the Substances
    # - temperature and pressure program
    a_d = [d]
    a_df = [df]
    d_fun(x) = GasChromatographySimulator.gradient(x, a_d)
    df_fun(x) = GasChromatographySimulator.gradient(x, a_df)
    col = GasChromatographySimulator.Column(L, d_fun, a_d, df_fun, a_df, sp, gas)

    solutes = ["C10", "C11"]
    init_t = zeros(length(solutes))
    init_τ = zeros(length(solutes))
    sub = GasChromatographySimulator.load_solute_database(db_path, db, col.sp, col.gas, solutes, init_t, init_τ)

    # Tinlet = "inlet"/"outlet" and odesys = true/false
    opt = [ #GasChromatographySimulator.Options(alg=OwrenZen3()),
            GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "inlet", true),
            GasChromatographySimulator.Options(OwrenZen5(), 1e-6, 1e-3, "outlet", false)
            ]

    par_g = Array{GasChromatographySimulator.Parameters}(undef, length(opt))
    results_g = Array{Any}(undef, length(opt))
    for i=1:length(opt)
        prog_g = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, ΔT_steps, x₀_steps, L₀_steps, α_steps, opt[i].Tcontrol, col.L)
        par_g[i] = GasChromatographySimulator.Parameters(col, prog_g, sub, opt[i])
        results_g[i] = GasChromatographySimulator.simulate(par_g[i])
    end

    @test length(results_g[1]) == 2
    @test length(results_g[2]) == 3

    #@test isapprox(results_g[1][1].tR[1], 123.184, atol=1e-3)
    @test isapprox(results_g[1][1].tR[1], 123.186, atol=1e-3)
    @test isapprox(results_g[2][1].tR[1], 51.4564, atol=1e-4)

    #@test isapprox(results_g[1][1].τR[2], 0.548163, atol=1e-5)
    @test isapprox(results_g[1][1].τR[2], 0.547683, atol=1e-5)
    @test isapprox(results_g[2][1].τR[2], 0.524873, atol=1e-4)

    # sol_extraction()
    df_sol = GasChromatographySimulator.sol_extraction(results_g[1][2], par_g[1])
    @test df_sol.t[1][end] ≈ results_g[1][1].tR[1]
    @test df_sol.τ²[2][end] ≈ results_g[1][1].τR[2]^2

    df_sol = GasChromatographySimulator.sol_extraction(results_g[2][2], results_g[2][3], par_g[2])
    @test df_sol.t[1][end] ≈ results_g[2][1].tR[1]
    @test df_sol.τ²[2][end] ≈ results_g[2][1].τR[2]^2

    # ng = true
    opt_ng = GasChromatographySimulator.Options(ng=true)
    prog_o_ng = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, [zeros(length(time_steps)) x₀_steps L₀_steps α_steps], opt_ng.Tcontrol, col.L)
    par_o_ng = GasChromatographySimulator.Parameters(col, prog_o_ng, sub, opt_ng)
    results_o_ng = GasChromatographySimulator.simulate(par_o_ng)
    @test isapprox(results_o_ng[1].tR[1], 87.401, atol=1e-3)
    @test isapprox(results_o_ng[1].τR[2], 0.59028, atol=1e-5)

    # vis = "HP"
    #opt_vis = GasChromatographySimulator.Options(vis="HP")
    #prog_vis = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, [zeros(length(time_steps)) x₀_steps L₀_steps α_steps], opt_vis.Tcontrol, col.L)
    #par_vis = GasChromatographySimulator.Parameters(col, prog_vis, sub, opt_vis)
    #results_vis = GasChromatographySimulator.simulate(par_vis)
    #@test isapprox(results_vis[1].tR[1], 87.078, atol=1e-3)
    #@test isapprox(results_vis[1].τR[2], 0.59294, atol=1e-5)

    # vis = "HP" and control = "Flow"
    F_steps = [1.0, 1.0, 1.5, 1.2]./(60e6)
    opt_control = GasChromatographySimulator.Options(vis="HP", control="Flow")
    prog_control = GasChromatographySimulator.Program(time_steps, temp_steps, F_steps, pout_steps, [ΔT_steps x₀_steps L₀_steps α_steps], opt_control.Tcontrol, col.L)
    par_control = GasChromatographySimulator.Parameters(col, prog_control, sub, opt_control)
    results_control = GasChromatographySimulator.simulate(par_control)
    @test isapprox(results_control[1].tR[1], 184.280, atol=1e-3)
    @test isapprox(results_control[1].τR[2], 0.25771, atol=1e-5)

    # compare_peaklist from Misc.jl
    pl1 = results_g[1][1]
    pl2 = results_g[2][1]
    comp = GasChromatographySimulator.compare_peaklist(pl1, pl2)
    @test comp.ΔtR == pl1.tR .- pl2.tR

    # compare_measurement_simulation from Misc.jl
    meas = DataFrame(Name=["C9", "C10"], RT=[100.0, 123.0])
    pl1 = results_g[1][1]
    comp = GasChromatographySimulator.compare_measurement_simulation(meas, pl1)
    @test isnan(comp.simulated_tR[1])
end

@testset "plots check" begin

    a_d = [d]
    a_df = [df]
    d_fun(x) = GasChromatographySimulator.gradient(x, a_d)
    df_fun(x) = GasChromatographySimulator.gradient(x, a_df)
    col = GasChromatographySimulator.Column(L, d_fun, a_d, df_fun, a_df, sp, gas)

    solutes = ["C10", "C11"]
    init_t = zeros(length(solutes))
    init_τ = zeros(length(solutes))
    sub = GasChromatographySimulator.load_solute_database(db_path, db, col.sp, col.gas, solutes, init_t, init_τ)

    opt_ng = GasChromatographySimulator.Options(ng=true)
    prog_o_ng = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, [zeros(length(time_steps)) x₀_steps L₀_steps α_steps], opt_ng.Tcontrol, col.L)
    par_o_ng = GasChromatographySimulator.Parameters(col, prog_o_ng, sub, opt_ng)
    results_o_ng = GasChromatographySimulator.simulate(par_o_ng)

    # plot functions
    t_start = 0.0
    t_end = 200.0
    p_chrom, t, chrom = GasChromatographySimulator.plot_chromatogram(results_o_ng[1], (t_start, t_end))
    @test t[1] == t_start && t[end] == t_end
    tR = results_o_ng[1].tR[1]
    τR = results_o_ng[1].τR[1]
    t0 = round(tR; digits=1)
    @test chrom[findfirst(t.==t0)] == 1/sqrt(2*π*τR^2)*exp(-(t0-tR)^2/(2*τR^2))
    t2, chrom2 = GasChromatographySimulator.plot_chromatogram!(p_chrom, results_o_ng[1], (t_start, t_end); mirror=true)
    @test t == t2
    @test chrom == -chrom2
    @test p_chrom[1][1][:y] == - p_chrom[1][2][:y]
    @test p_chrom[1][1][:x] == p_chrom[1][2][:x]
    # mirrored the other way
    p_chrom, t, chrom = GasChromatographySimulator.plot_chromatogram(results_o_ng[1], (t_start, t_end); mirror=true)
    @test t[1] == t_start && t[end] == t_end
    tR = results_o_ng[1].tR[1]
    τR = results_o_ng[1].τR[1]
    t0 = round(tR; digits=1)
    @test chrom[findfirst(t.==t0)] == -1/sqrt(2*π*τR^2)*exp(-(t0-tR)^2/(2*τR^2))
    t2, chrom2 = GasChromatographySimulator.plot_chromatogram!(p_chrom, results_o_ng[1], (t_start, t_end); mirror=false)
    @test t == t2
    @test chrom == -chrom2
    @test p_chrom[1][1][:y] == - p_chrom[1][2][:y]
    @test p_chrom[1][1][:x] == p_chrom[1][2][:x]

    p_flow = GasChromatographySimulator.plot_flow(par_o_ng)
    #@test p_flow[1][:xaxis][:optimized_ticks][1][1] = 0.0
    #@test p_flow[1][:xaxis][:optimized_ticks][1][end] = 800.0
    @test isapprox(p_flow[1][1][:y][1], 11.9797, atol=1e-4)

    p_pres = GasChromatographySimulator.plot_pressure(par_o_ng)
    #@test p_pres[1][:xaxis][:optimized_ticks][1][1] = 0.0
    #@test p_pres[1][:xaxis][:optimized_ticks][1][end] = 800.0
    @test p_pres[1][1][:y][1] == par_o_ng.prog.Fpin_steps[1]
    @test p_pres[1][2][:y][end] == par_o_ng.prog.pout_steps[end]

    p_temp = GasChromatographySimulator.plot_temperature(par_o_ng)
    @test p_temp[1][1][:y][1] == par_o_ng.prog.temp_steps[1]
    @test p_temp[1][1][:y][end] == par_o_ng.prog.temp_steps[end]
    @test p_temp[1][2][:x][end] == sum(par_o_ng.prog.time_steps)

    p_temp = GasChromatographySimulator.plot_temperature(par_o_ng; selector="T(x)")
    @test p_temp[1][1][:y][1] == par_o_ng.prog.temp_steps[1]
    @test p_temp[1][end][:y][end] == par_o_ng.prog.temp_steps[end]
    @test p_temp[1][2][:x][end] == par_o_ng.col.L

    p_temp = GasChromatographySimulator.plot_temperature(par_o_ng; selector="T(x,t)")
    @test p_temp[1][1][:y][end] == sum(par_o_ng.prog.time_steps)
    @test p_temp[1][1][:x][end] == par_o_ng.col.L

    p_local = GasChromatographySimulator.local_plots("z", "z", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] == p_local[1][1][:y][end]
    p_local = GasChromatographySimulator.local_plots("t", "t", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] == p_local[1][1][:y][end]
    p_local = GasChromatographySimulator.local_plots("T", "T", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] == p_local[1][1][:y][end]
    p_local = GasChromatographySimulator.local_plots("τ", "τ", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] == p_local[1][1][:y][end]
    p_local = GasChromatographySimulator.local_plots("σ", "σ", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] === p_local[1][1][:y][end]
    p_local = GasChromatographySimulator.local_plots("u", "u", results_o_ng[2], par_o_ng)
    @test p_local[1][1][:x][end] === p_local[1][1][:y][end]
end

@testset "UI checks" begin
    col_val = GasChromatographySimulator.UI_Column(["HP5", "DB5"])
    @test typeof(col_val) == PlutoUI.CombineNotebook.CombinedBonds
    prog_val = GasChromatographySimulator.UI_Program()
    @test typeof(prog_val) == PlutoUI.CombineNotebook.CombinedBonds
    prog_val = GasChromatographySimulator.UI_Program(default=("0 60", "40 140", "100 100", "vacuum"))
    @test typeof(prog_val) == PlutoUI.CombineNotebook.CombinedBonds
    sub_val = GasChromatographySimulator.UI_Substance(["C6", "C7", "C8", "C9"])
    @test typeof(prog_val) == PlutoUI.CombineNotebook.CombinedBonds
    sub_val = GasChromatographySimulator.UI_Substance(["C6", "C7", "C8", "C9"]; default=(1:2, 0.1, 1.0))
    @test typeof(prog_val) == PlutoUI.CombineNotebook.CombinedBonds
    opt_val = GasChromatographySimulator.UI_Options()
    @test typeof(opt_val) == PlutoUI.CombineNotebook.CombinedBonds
    prog = GasChromatographySimulator.setting_prog(("0 60", "40 140", "100 100", "vacuum"), L)
    @test prog.pout_steps[1] == 0.0
    prog = GasChromatographySimulator.setting_prog(("0 60", "40 140", "100 100", "atmosphere"), L)
    @test prog.pout_steps[1] == 101300.0
end

@testset "misc checks" begin
    @test "C10" in GasChromatographySimulator.all_solutes(sp, DataFrame(CSV.File(string(db_path,"/",db), header=1, silencewarnings=true)))
    @test "C11" in GasChromatographySimulator.all_solutes(sp, DataFrame(CSV.File(string(db_path,"/",db), header=1, silencewarnings=true)))

    a_d = [d]
    a_df = [df]
    d_fun(x) = GasChromatographySimulator.gradient(x, a_d)
    df_fun(x) = GasChromatographySimulator.gradient(x, a_df)
    col = GasChromatographySimulator.Column(L, d_fun, a_d, df_fun, a_df, sp, gas)
    a_gf = [ΔT_steps x₀_steps L₀_steps α_steps]
    gf(x) = GasChromatographySimulator.gradient(x, a_gf)
    T_itp = GasChromatographySimulator.temperature_interpolation(time_steps, temp_steps, gf, col.L)
    pin_itp = GasChromatographySimulator.steps_interpolation(time_steps, pin_steps)
    pout_itp = GasChromatographySimulator.steps_interpolation(time_steps, pout_steps)
    # without gradient
    prog_ng = GasChromatographySimulator.Program(time_steps, temp_steps, pin_steps, pout_steps, col.L)

    # random time in the temperature program
    t = rand()*sum(time_steps)
    # temperature at random time at column inlet
    T_test = prog_ng.T_itp(0.0, t)

    # viscosity tests
    # vis = "Blumberg"
        # gas = "He"
        η_T = GasChromatographySimulator.viscosity(T_test, "He")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, col.gas)
        @test η_T ≈ η_t
        # gas = "H2"
        η_T = GasChromatographySimulator.viscosity(T_test, "H2")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "H2")
        @test η_T ≈ η_t
        # gas = "N2"
        η_T = GasChromatographySimulator.viscosity(T_test, "N2")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "N2")
        @test η_T ≈ η_t
        # gas = "Ar"
        η_T = GasChromatographySimulator.viscosity(T_test, "Ar")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "Ar")
        @test η_T ≈ η_t
    # vis = "HP
        # gas = "He"
        η_T = GasChromatographySimulator.viscosity(T_test, "He"; vis="HP")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "He"; vis="HP")
        @test η_T ≈ η_t
        # gas = "H2"
        η_T = GasChromatographySimulator.viscosity(T_test, "H2"; vis="HP")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "H2"; vis="HP")
        @test η_T ≈ η_t
        # gas = "N2"
        η_T = GasChromatographySimulator.viscosity(T_test, "N2"; vis="HP")
        η_t = GasChromatographySimulator.viscosity(0.0, t, prog_ng.T_itp, "N2"; vis="HP")
        @test η_T ≈ η_t

    tM_T = GasChromatographySimulator.holdup_time(T_test, prog_ng.Fpin_itp(t), prog_ng.pout_itp(t), col.L, col.a_d[1], col.gas) # only defined for non-gradient case
    tM_t = GasChromatographySimulator.holdup_time(t, prog_ng.T_itp, prog_ng.Fpin_itp, prog_ng.pout_itp, col.L, col.d, col.gas)
    tM_t_ng = GasChromatographySimulator.holdup_time(t, prog_ng.T_itp, prog_ng.Fpin_itp, prog_ng.pout_itp, col.L, col.d, col.gas; ng=true)
    @test tM_T ≈ tM_t
    @test tM_T ≈ tM_t_ng

    F_T = GasChromatographySimulator.flow(T_test, prog_ng.Fpin_itp(t), prog_ng.pout_itp(t), col.L, col.a_d[1], col.gas) # only defined for non-gradient case
    F_t = GasChromatographySimulator.flow(t, prog_ng.T_itp, prog_ng.Fpin_itp, prog_ng.pout_itp, col.L, col.d, col.gas)

    @test F_T  ≈  F_t

    # control = "Flow"
    F_steps = [1.0, 1.0, 1.5, 1.2]./(60e6)
    # without gradient
    prog_F = GasChromatographySimulator.Program(time_steps, temp_steps, F_steps, pout_steps, col.L)
    F_T = GasChromatographySimulator.flow(T_test, prog_F.Fpin_itp(t), prog_F.pout_itp(t), col.L, col.a_d[1], col.gas; control="Flow") # only defined for non-gradient case
    F_t = GasChromatographySimulator.flow(t, prog_F.T_itp, prog_F.Fpin_itp, prog_F.pout_itp, col.L, col.d, col.gas; control="Flow")
    pin(t) = GasChromatographySimulator.inlet_pressure(t, prog_F.T_itp, prog_F.Fpin_itp, prog_F.pout_itp, col.L, col.d, col.gas; ng=false, vis="Blumberg", control="Flow")

    @test F_T  ≈  F_t
    @test F_T ≈ GasChromatographySimulator.flow(T_test, pin(t), prog_F.pout_itp(t), col.L, col.a_d[1], col.gas; control="Pressure")

    tM_T = GasChromatographySimulator.holdup_time(T_test, prog_F.Fpin_itp(t), prog_F.pout_itp(t), col.L, col.a_d[1], col.gas; control="Flow") # only defined for non-gradient case
    tM_t = GasChromatographySimulator.holdup_time(t, prog_F.T_itp, prog_F.Fpin_itp, prog_F.pout_itp, col.L, col.d, col.gas; control="Flow") # prog_F is without gradient -> should be the same
    @test tM_T ≈ tM_t

    tM_T = GasChromatographySimulator.holdup_time(T_test, prog_F.Fpin_itp(t), prog_F.pout_itp(t), col.L, col.a_d[1], col.gas; control="Flow") # only defined for non-gradient case
    tM_t = GasChromatographySimulator.holdup_time(t, prog_F.T_itp, prog_F.Fpin_itp, prog_F.pout_itp, col.L, col.d, col.gas; control="Flow", ng=true)
    @test tM_T ≈ tM_t

    F_T = GasChromatographySimulator.flow(T_test, prog_ng.Fpin_itp(t), prog_ng.pout_itp(t), col.L, col.a_d[1], col.gas) # only defined for non-gradient case
    F_t = GasChromatographySimulator.flow(t, prog_ng.T_itp, prog_ng.Fpin_itp, prog_ng.pout_itp, col.L, col.d, col.gas)
    F_t_ng = GasChromatographySimulator.flow(t, prog_ng.T_itp, prog_ng.Fpin_itp, prog_ng.pout_itp, col.L, col.d, col.gas, ng=true)
    @test F_T  ≈  F_t
    @test F_T  ≈  F_t_ng

    common = GasChromatographySimulator.common(["aa", "bb", "cc"], ["dd", "aa", "ee"])
    @test common == ["aa"]
end

println("Test run successful.")


text:
# Functions relating to converting Nimrod data to an ASCII ESRI grid

"""
    toasciigrid(n::Nimrod, outfile::AbstractString)

Write an instance of `Nimrod` to an ESRI ASCII file.

Arguments
=========

- `n::Nimrod` -- The nimrod data structure to write out.
- `outfile::AbstractString` -- The name of the file to write to.

Returns
=======

- `nothing`

Notes
=====

- It is assumed that the spatial resolution is the same in the x and y
  directions.
- For instances of `Nimrod{AbstractFloat}`, the no data value is given by header
  field 38.
- For instances of `Nimrod{Integer}`, the no data value is given by header field
  25.
"""
toasciigrid(n::Nimrod, outfile::AbstractString) = open((f) -> toasciigrid(n, f), outfile, "w")
function toasciigrid{T<:Integer}(n::Nimrod{T}, outfile::IOStream)
    hdrstr = """nrows $(size(n.data)[1])
ncols $(size(n.data)[2])
xllcenter $(n.llcenter[1])
yllcenter $(n.llcenter[2])
cellsize $(n.xpix_size)
nodata_value $(n.hdr_dump[25])
"""

    write(outfile, hdrstr)
    writedlm(outfile, n.data, ' ')
    return
end

function toasciigrid{T<:AbstractFloat}(n::Nimrod{T}, outfile::IOStream)
    hdrstr = """nrows $(size(n.data)[1])
ncols $(size(n.data)[2])
xllcenter $(n.llcenter[1])
yllcenter $(n.llcenter[2])
cellsize $(n.xpix_size)
nodata_value $(n.hdr_dump[38])
"""

    write(outfile, hdrstr)
    writedlm(outfile, n.data, ' ')
    return
end