text:
# ------------ Constrain in call hierarchy ------------ #

struct ConstrainedHierarchicalSelection{T <: ConstrainedSelectQuery} <: ConstrainedSelection
    tree::Dict{Address, ConstrainedSelection}
    query::T
    ConstrainedHierarchicalSelection() = new{ConstrainedByAddress}(Dict{Address, ConstrainedHierarchicalSelection}(), ConstrainedByAddress())
    ConstrainedHierarchicalSelection(csa::T) where T <: ConstrainedSelectQuery = new{T}(Dict{Address, ConstrainedHierarchicalSelection}(), csa)
end

has_sub(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Address = haskey(chs.tree, addr)
has_sub(chs::ConstrainedHierarchicalSelection, addr::Tuple{}) = false
has_sub(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = has_sub(chs, addr[1])
function has_sub(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && has_sub(get_sub(chs, addr[1]), addr[2 : end])
end

function get_sub(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Address
    haskey(chs.tree, addr) && return chs.tree[addr]
    return ConstrainedEmptySelection()
end
get_sub(chs::ConstrainedHierarchicalSelection, addr::Tuple{}) = ConstrainedEmptySelection()
get_sub(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = get_sub(chs, addr[1])
function get_sub(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Tuple
    haskey(chs.tree, addr[1]) && return get_sub(chs.tree[addr[1]], addr[2 : end])
    return ConstrainedEmptySelection()
end

set_sub!(chs::ConstrainedHierarchicalSelection, addr::T, sub::K) where {T <: Address, K <: ConstrainedSelection} = chs.tree[addr] = sub
set_sub!(chs::ConstrainedHierarchicalSelection, addr::Tuple{}, sub::K) where K <: ConstrainedSelection = nothing
set_sub!(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}, sub::K) where {T <: Address, K <: ConstrainedSelection} = set_sub!(chs, addr[1], sub)
function set_sub!(chs::ConstrainedHierarchicalSelection, addr::T, sub::K) where {T <: Tuple, K <: ConstrainedSelection}
    has_sub(chs, addr[1]) && set_sub!(chs.tree[addr[1]], addr[2 : end], sub)
end

has_top(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Address = has_top(chs.query, addr)
function has_top(chs::ConstrainedHierarchicalSelection, addr::Tuple{}) where T <: Address
    return false
end
function has_top(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address
    return has_top(chs.query, addr[1])
end
function has_top(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && return has_top(get_sub(chs, addr[1]), addr[2 : end])
end

get_top(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Address = get_top(chs.query, addr)
get_top(chs::ConstrainedHierarchicalSelection, addr::Tuple{}) = nothing
get_top(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = get_top(chs, addr[1])
function get_top(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && return get_top(chs.tree[addr[1]], addr[2 : end])
    error("ConstrainedHierarchicalSelection (get_top): invalid index at $addr.")
end

haskey(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Address = has_top(chs, addr) || has_sub(chs, addr)
haskey(chs::ConstrainedHierarchicalSelection, addr::Tuple{}) = false
haskey(chs::ConstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = haskey(chs, addr[1])
function haskey(chs::ConstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && haskey(get_sub(chs, addr[1]), addr[2 : end])
end

function getindex(chs::ConstrainedHierarchicalSelection, addr...)
    has_top(chs, addr) && return get_top(chs, addr)
    has_sub(chs, addr) && return get_sub(chs, addr)
    error("ConstrainedHierarchicalSelection (getindex): no query or subselection at $addr.")
end

push_query!(toplevel::Set, k::K, q::T) where {K, T <: Address} = push!(toplevel, (k, q))
push_query!(toplevel::Set, k::K, q::T) where {K, T <: Tuple} = push!(toplevel, (k, q...))

function dump_queries(chs::ConstrainedHierarchicalSelection)
    toplevel = dump_queries(chs.query)
    for (k, v) in chs.tree
        q = dump_queries(v)
        for (t, l) in q
            push!(toplevel, (k, t...) => l)
        end
    end
    return toplevel
end

isempty(chs::ConstrainedHierarchicalSelection) = begin
    !isempty(chs.query) && return false
    for (k, v) in chs.tree
        !isempty(chs.tree[k]) && return false
    end
    return true
end

function ==(sel1::ConstrainedHierarchicalSelection, sel2::ConstrainedHierarchicalSelection)
    sel1.query == sel2.query || return false
    keys(sel1.tree) == keys(sel2.tree) || return false
    for (k, v) in sel1.tree
        v == get_sub(sel2, k) || return false
    end
    return true
end

# Used to merge observations.
function merge!(sel1::ConstrainedHierarchicalSelection, sel2::ConstrainedHierarchicalSelection)
    merge!(sel1.query, sel2.query)
    for k in keys(sel2.tree)
        if haskey(sel1.tree, k)
            merge!(sel1.tree[k], sel2.tree[k])
        else
            sel1.tree[k] = sel2.tree[k]
        end
    end
end
function merge(cl::T, sel::ConstrainedHierarchicalSelection) where T <: CallSite
    cl_selection = get_selection(cl)
    merge!(cl_selection, sel)
    return cl_selection
end

# Used to build.
push!(sel::ConstrainedHierarchicalSelection, addr::T, val) where T <: Address = push!(sel.query, addr, val)
push!(sel::ConstrainedHierarchicalSelection, addr::Tuple{}, val) = nothing
push!(sel::ConstrainedHierarchicalSelection, addr::Tuple{T}, val) where T <: Address = push!(sel, addr[1], val)
function push!(sel::ConstrainedHierarchicalSelection, addr::T, val) where T <: Tuple
    if has_sub(sel, addr[1])
        push!(get_sub(sel, addr[1]), addr[2 : end], val)
    else
        new = ConstrainedHierarchicalSelection()
        push!(new, addr[2 : end], val)
        sel.tree[addr[1]] = new
    end
end

# Used for functional filter querying.
function filter(k_fn::Function, v_fn::Function, chs::ConstrainedHierarchicalSelection) where T <: Address
    top = ConstrainedHierarchicalSelection(filter(k_fn, v_fn, chs.query))
    for (k, v) in chs.tree
        top.tree[k] = filter(k_fn, v_fn, v)
    end
    isempty(top) && return ConstrainedEmptySelection()
    return top
end

# To and from arrays.
function fill_array!(val::T, arr::Vector{K}, f_ind::Int) where {K, T <: ConstrainedHierarchicalSelection}
    sorted_toplevel_keys = sort(collect(addresses(val.query)))
    sorted_tree_keys  = sort(collect(keys(val.tree)))
    idx = f_ind
    for k in sorted_toplevel_keys
        v = get_top(val, k)
        n = fill_array!(v, arr, idx)
        idx += n
    end
    for k in sorted_tree_keys
        n = fill_array!(get_sub(val, k), arr, idx)
        idx += n
    end
    idx - f_ind
end

function from_array(schema::T, arr::Vector{K}, f_ind::Int) where {K, T <: ConstrainedHierarchicalSelection}
    sel = ConstrainedHierarchicalSelection()
    sorted_toplevel_keys = sort(collect(addresses(schema.query)))
    sorted_tree_keys  = sort(collect(keys(schema.tree)))
    idx = f_ind
    for k in sorted_toplevel_keys
        (n, v) = from_array(get_top(schema, k), arr, idx)
        idx += n
        push!(sel, k, v)
    end
    for k in sorted_tree_keys
        (n, v) = from_array(get_sub(schema, k), arr, idx)
        idx += n
        sel.tree[k] = v
    end
    (idx - f_ind, sel)
end

# Used for pretty printing.
function collect!(args...) end
function collect!(par::T, addrs::Vector{Any}, chd::Dict{Any, Any}, chs::ConstrainedHierarchicalSelection) where T <: Tuple
    collect!(par, addrs, chd, chs.query)
    for (k, v) in chs.tree
        collect!((par..., k), addrs, chd, v)
    end
end
function collect!(addrs::Vector{Any}, chd::Dict{Any, Any}, chs::ConstrainedHierarchicalSelection)
    collect!(addrs, chd, chs.query)
    for (k, v) in chs.tree
        collect!((k, ), addrs, chd, v)
    end
end

function collect(chs::ConstrainedHierarchicalSelection)
    addrs = Any[]
    chd = Dict{Any, Any}()
    collect!(addrs, chd, chs)
    return addrs, chd
end

# Pretty printing.
function Base.display(chs::ConstrainedHierarchicalSelection; show_values = true)
    println("  __________________________________\n")
    println("             Constrained\n")
    addrs, chd = collect(chs)
    if show_values
        for a in addrs
            println(" $(a) = $(chd[a])")
        end
    else
        for a in addrs
            println(" $(a)")
        end
    end
    println("  __________________________________\n")
end

# ------------ CHS from vectors ----------- #

function ConstrainedHierarchicalSelection(a::Vector{Pair{K, J}}) where {K <: Tuple, J}
    top = ConstrainedHierarchicalSelection()
    for (addr, val) in a
        push!(top, addr, val)
    end
    return top
end

# ------------ CHS from traces ------------ #

site_push!(chs::ConstrainedHierarchicalSelection, addr::T, cs::ChoiceSite) where T <: Address = push!(chs, addr, cs.val)
site_push!(chs::ConstrainedHierarchicalSelection, addr::T, cs::CallSite) where T <: Address = set_sub!(chs, addr, get_selection(cs))
function push!(chs::ConstrainedHierarchicalSelection, tr::VectorizedTrace)
    for (k, cs) in enumerate(tr.subrecords)
        site_push!(chs, k, cs)
    end
end
function push!(chs::ConstrainedHierarchicalSelection, tr::HierarchicalTrace)
    for k in keys(tr.calls)
        site_push!(chs, k, tr.calls[k])
    end
    for k in keys(tr.choices)
        site_push!(chs, k, tr.choices[k])
    end
end
function get_selection(tr::VectorizedTrace)
    top = ConstrainedHierarchicalSelection()
    push!(top, tr)
    return top
end
function get_selection(tr::HierarchicalTrace)
    top = ConstrainedHierarchicalSelection()
    push!(top, tr)
    return top
end
function get_selection(cl::VectorizedCallSite)
    top = ConstrainedHierarchicalSelection()
    push!(top, cl.trace)
    return top
end
function get_selection(cl::HierarchicalCallSite)
    top = ConstrainedHierarchicalSelection()
    push!(top, cl.trace)
    return top
end

# ------------ Unconstrained selection in call hierarchy ------------ #

struct UnconstrainedHierarchicalSelection{T <: UnconstrainedSelectQuery} <: UnconstrainedSelection
    tree::Dict{Address, UnconstrainedHierarchicalSelection}
    query::T
    UnconstrainedHierarchicalSelection() = new{UnconstrainedByAddress}(Dict{Address, UnconstrainedHierarchicalSelection}(), UnconstrainedByAddress())
end

function get_sub(uhs::UnconstrainedHierarchicalSelection, addr::T) where T <: Address
    haskey(uhs.tree, addr) && return uhs.tree[addr]
    return UnconstrainedEmptySelection()
end
get_sub(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{}) = UnconstrainedEmptySelection()
get_sub(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = get_sub(uhs, addr[1])
function get_sub(uhs::UnconstrainedHierarchicalSelection, addr::T) where T <: Tuple
    haskey(uhs.tree, addr[1]) && return get_sub(uhs.tree[addr[1]], addr[2 : end])
    return UnconstrainedEmptySelection()
end

has_top(uhs::UnconstrainedHierarchicalSelection, addr::T) where T <: Address = has_top(uhs.query, addr)
has_top(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{}) where T <: Address = false
has_top(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = has_top(uhs, addr[1])
function has_top(uhs::UnconstrainedHierarchicalSelection, addr::T) where T <: Tuple
    haskey(uhs.tree, addr[1]) && has_top(uhs.tree[addr[1]], addr[2 : end])
end

has_sub(chs::UnconstrainedHierarchicalSelection, addr::T) where T <: Address = haskey(chs.tree, addr)
has_sub(chs::UnconstrainedHierarchicalSelection, addr::Tuple{}) = false
has_sub(chs::UnconstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = has_sub(chs, addr[1])
function has_sub(chs::UnconstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && has_sub(get_sub(chs, addr[1]), addr[2 : end])
end

function set_sub!(uhs::UnconstrainedHierarchicalSelection, addr::T, sub::K) where {T <: Address, K <: UnconstrainedSelection}
    uhs.tree[addr] = sub
end
set_sub!(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{}, sub::K) where K <: UnconstrainedSelection = nothing
function set_sub!(uhs::UnconstrainedHierarchicalSelection, addr::Tuple{T}, sub::K) where {T <: Address, K <: UnconstrainedSelection}
    set_sub!(uhs, addr[1], sub)
end
function set_sub!(uhs::UnconstrainedHierarchicalSelection, addr::T, sub::K) where {T <: Tuple, K <: UnconstrainedSelection}
    haskey(uhs.tree, addr[1]) && set_sub!(uhs.tree[1], addr[2 : end], sub)
end

function dump_queries(uhs::UnconstrainedHierarchicalSelection)
    toplevel = dump_queries(uhs.query)
    for (k, v) in uhs.tree
        toplevel = union(toplevel, dump_queries(v))
    end
    return toplevel
end

haskey(chs::UnconstrainedHierarchicalSelection, addr::T) where T <: Address = has_top(chs, addr) || has_sub(chs, addr)
haskey(chs::UnconstrainedHierarchicalSelection, addr::Tuple{}) = false
haskey(chs::UnconstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = haskey(chs, addr[1])
function haskey(chs::UnconstrainedHierarchicalSelection, addr::T) where T <: Tuple
    has_sub(chs, addr[1]) && haskey(get_sub(chs, addr[1]), addr[2 : end])
end

isempty(uhs::UnconstrainedHierarchicalSelection) = begin
    isempty(uhs.query) || return false
    for (k, v) in uhs.tree
        isempty(uhs.tree[k]) || return false
    end
    return true
end

function ==(sel1::UnconstrainedHierarchicalSelection, sel2::UnconstrainedHierarchicalSelection)
    sel1.query == sel2.query || return false
    keys(sel1.tree) == keys(sel2.tree) || return false
    for (k, v) in sel1.tree
        v == get_sub(sel2, k) || return false
    end
    return true
end

# Used to merge observations.
function merge!(sel1::UnconstrainedHierarchicalSelection, sel2::UnconstrainedHierarchicalSelection)
    merge!(sel1.query, sel2.query)
    for k in keys(sel2.tree)
        if haskey(sel1.tree, k)
            merge!(sel1.tree[k], sel2.tree[k])
        else
            sel1.tree[k] = sel2.tree[k]
        end
    end
end
function merge(cl::T, sel::UnconstrainedHierarchicalSelection) where T <: CallSite
    cl_selection = get_selection(cl)
    un_sel = selection(map(collect(dump_queries(cl_selection))) do k
                           (k, )
                       end)
    merge!(un_sel, sel)
    return un_sel
end

# Used to build.
push!(sel::UnconstrainedHierarchicalSelection, addr::T) where T <: Address = push!(sel.query, addr)
push!(sel::UnconstrainedHierarchicalSelection, addr::Tuple{}) = nothing
push!(sel::UnconstrainedHierarchicalSelection, addr::Tuple{T}) where T <: Address = push!(sel, addr[1])
function push!(sel::UnconstrainedHierarchicalSelection, addr::Tuple)
    if haskey(sel.tree, addr[1])
        push!(get_sub(sel, addr[1]), addr[2 : end])
    else
        new = UnconstrainedHierarchicalSelection()
        push!(new, addr[2 : end])
        sel.tree[addr[1]] = new
    end
end

# ------------ UHS from vectors ------------ #

function UnconstrainedHierarchicalSelection(a::Vector{K}) where K <: Tuple
    top = UnconstrainedHierarchicalSelection()
    for addr in a
        push!(top, addr)
    end
    return top
end

# Used in functional filter querying.
function filter(k_fn::Function, chs::UnconstrainedHierarchicalSelection) where T <: Address
    top = UnconstrainedHierarchicalSelection(filter(k_fn, chs.query))
    for (k, v) in chs.tree
        top.tree[k] = filter(k_fn, v)
    end
    isempty(top) && return UnconstrainedEmptySelection()
    return top
end

# Used in pretty printing.
function collect!(par::T, addrs::Vector{Any}, chs::UnconstrainedHierarchicalSelection) where T <: Any
    collect!(par, addrs, chs.query)
    for (k, v) in chs.tree
        collect!((par..., k), addrs, v)
    end
end

function collect!(addrs::Vector{Any}, chs::UnconstrainedHierarchicalSelection)
    collect!(addrs, chs.query)
    for (k, v) in chs.tree
        collect!((k, ), addrs, v)
    end
end

function collect(chs::UnconstrainedHierarchicalSelection)
    addrs = Any[]
    collect!(addrs, chs)
    return addrs
end

# Pretty printing.
function Base.display(chs::UnconstrainedHierarchicalSelection)
    println("  __________________________________\n")
    println("              Selection\n")
    addrs = collect(chs)
    for a in addrs
        println(" $(a)")
    end
    println("  __________________________________\n")
end


text:
# ---
# title: 349. Intersection of Two Arrays
# id: problem349
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Two Pointers, Binary Search, Sort
# link: <https://leetcode.com/problems/intersection-of-two-arrays/description/>
# hidden: true
# ---
#
# Given two arrays, write a function to compute their intersection.
#
# **Example 1:**
#
#
#
#     Input: nums1 = [1,2,2,1], nums2 = [2,2]
#     Output: [2]
#
#
# **Example 2:**
#
#
#
#     Input: nums1 = [4,9,5], nums2 = [9,4,9,8,4]
#     Output: [9,4]
#
# **Note:**
#
#   * Each element in the result must be unique.
#   * The result can be in any order.
#
#
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
module ColumnwiseLUSolver

export ManyColumnLU, SingleColumnLU

using ..Mesh.Grids
using ..Mesh.Grids: VerticalDirection
using ..Mesh.Topologies
using ..DGmethods
using ..DGmethods: BalanceLaw, DGModel, num_state, num_diffusive
using ..LinearSolvers
const LS = LinearSolvers
using ..MPIStateArrays
using LinearAlgebra
using GPUifyLoops

include("ColumnwiseLUSolver_kernels.jl")

abstract type AbstractColumnLUSolver <: AbstractLinearSolver  end

"""
    ManyColumnLU()

This solver is used for systems that are block diagonal where each block is
associated with a column of the mesh.  The systems are solved using a
non-pivoted LU factorization.
"""
struct ManyColumnLU <: AbstractColumnLUSolver  end

"""
    SingleColumnLU()

This solver is used for systems that are block diagonal where each block is
associated with a column of the mesh.  Moreover, each block is assumed to be
the same.  The systems are solved using a non-pivoted LU factorization.
"""
struct SingleColumnLU <: AbstractColumnLUSolver  end

struct ColumnwiseLU{F, AT}
  f::F
  A::AT
end

function LS.prefactorize(op, solver::AbstractColumnLUSolver, Q, args...)
  dg = op.f!

  # TODO: can we get away with just passing the grid?
  A = banded_matrix(op, dg, similar(Q), similar(Q), args...;
                    single_column = typeof(solver) <: SingleColumnLU)

  band_lu!(A, dg)

  ColumnwiseLU(dg, A)
end

function LS.linearsolve!(clu::ColumnwiseLU{F}, ::AbstractColumnLUSolver,
                         Q, Qrhs, args...) where {F <: DGModel}
  dg = clu.f
  A = clu.A
  Q .= Qrhs

  band_forward!(Q, A, dg)
  band_back!(Q, A, dg)
end

"""
    band_lu!(A, dg::DGModel)

"""
function band_lu!(A, dg::DGModel)
  bl = dg.balancelaw
  grid = dg.grid
  topology = grid.topology
  @assert isstacked(topology)
  @assert typeof(dg.direction) <: VerticalDirection

  FT = eltype(A)
  device = typeof(A) <: Array ? CPU() : CUDA()

  nstate = num_state(bl, FT)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqj = dimensionality(grid) == 2 ? 1 : Nq

  eband = num_diffusive(bl, FT) == 0 ? 1 : 2

  nrealelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nrealelem, nvertelem)

  threads = (Nq, Nqj)
  blocks = nhorzelem

  if ndims(A) == 2
    # single column case
    #
    # TODO Would it be faster to copy the matrix to the host and factorize it
    # there?
    threads = (1, 1)
    blocks = 1
    A = reshape(A, 1, 1, size(A)..., 1)
  end

  @launch(device, threads=threads, blocks=blocks,
          band_lu_knl!(A, Val(Nq), Val(threads[1]), Val(threads[2]),
                       Val(nstate), Val(nvertelem), Val(blocks),
                       Val(eband)))
end

function band_forward!(Q, A, dg::DGModel)
  bl = dg.balancelaw
  grid = dg.grid
  topology = grid.topology
  @assert isstacked(topology)
  @assert typeof(dg.direction) <: VerticalDirection

  FT = eltype(A)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  nstate = num_state(bl, FT)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqj = dimensionality(grid) == 2 ? 1 : Nq

  eband = num_diffusive(bl, FT) == 0 ? 1 : 2

  nrealelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nrealelem, nvertelem)

  @launch(device, threads=(Nq, Nqj), blocks=nhorzelem,
          band_forward_knl!(Q.data, A, Val(Nq), Val(Nqj), Val(nstate),
                            Val(nvertelem), Val(nhorzelem), Val(eband)))
end

function band_back!(Q, A, dg::DGModel)
  bl = dg.balancelaw
  grid = dg.grid
  topology = grid.topology
  @assert isstacked(topology)
  @assert typeof(dg.direction) <: VerticalDirection

  FT = eltype(A)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  nstate = num_state(bl, FT)
  N = polynomialorder(grid)
  Nq = N + 1
  Nqj = dimensionality(grid) == 2 ? 1 : Nq

  eband = num_diffusive(bl, FT) == 0 ? 1 : 2

  nrealelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nrealelem, nvertelem)

  @launch(device, threads=(Nq, Nqj), blocks=nhorzelem,
          band_back_knl!(Q.data, A, Val(Nq), Val(Nqj), Val(nstate),
                         Val(nvertelem), Val(nhorzelem), Val(eband)))
end

"""
    banded_matrix(dg::DGModel, [Q::MPIStateArray, dQ::MPIStateArray,
                  single_column=false])

Forms the banded matrices for each the column operator defined by the `DGModel`
dg.  If `single_column=false` then a banded matrix is stored for each column and
if `single_column=true` only the banded matrix associated with the first column
of the first element is stored. The bandwidth of the DG column banded matrix is
`p = q = (polynomialorder + 1) * nstate * nvertelem - 1` with `p` and `q` being
the upper and lower bandwidths.

The banded matrices are stored in the LAPACK band storage format
<https://www.netlib.org/lapack/lug/node124.html>.

The banded matrices are returned as an arrays where the array type matches that
of `Q`. If `single_column=false` then the returned array has 5 dimensions, which
are:
- first horizontal column index
- second horizontal column index
- band index (-q:p)
- vertical DOF index with state `s`, vertical DOF index `k`, and vertical
  element `ev` mapping to `s + nstate * (k - 1) + nstate * nvertelem * (ev - 1)`
- horizontal element index

If the `single_column=true` then the returned array has 2 dimensions which are
the band index and the vertical DOF index.
"""
function banded_matrix(dg::DGModel, Q::MPIStateArray = MPIStateArray(dg),
                       dQ::MPIStateArray = MPIStateArray(dg);
                       single_column = false)
  banded_matrix((dQ, Q) -> dg(dQ, Q, nothing, 0; increment=false),
                dg, Q, dQ; single_column = single_column)
end

"""
    banded_matrix(f!::Function, dg::DGModel,
                  Q::MPIStateArray = MPIStateArray(dg),
                  dQ::MPIStateArray = MPIStateArray(dg), args...;
                  single_column = false, args...)

Forms the banded matrices for each the column operator defined by the linear
operator `f!` which is assumed to have the same banded structure as the
`DGModel` dg.  If `single_column=false` then a banded matrix is stored for each
column and if `single_column=true` only the banded matrix associated with the
first column of the first element is stored. The bandwidth of the DG column
banded matrix is `p = q = (polynomialorder + 1) * nstate * nvertelem - 1` with
`p` and `q` being the upper and lower bandwidths.

The banded matrices are stored in the LAPACK band storage format
<https://www.netlib.org/lapack/lug/node124.html>.

The banded matrices are returned as an arrays where the array type matches that
of `Q`. If `single_column=false` then the returned array has 5 dimensions, which
are:
- first horizontal column index
- second horizontal column index
- band index (-q:p)
- vertical DOF index with state `s`, vertical DOF index `k`, and vertical
  element `ev` mapping to `s + nstate * (k - 1) + nstate * nvertelem * (ev - 1)`
- horizontal element index

If the `single_column=true` then the returned array has 2 dimensions which are
the band index and the vertical DOF index.

Here `args` are passed to `f!`.
"""
function banded_matrix(f!, dg::DGModel,
                       Q::MPIStateArray = MPIStateArray(dg),
                       dQ::MPIStateArray = MPIStateArray(dg),
                       args...; single_column = false)
  bl = dg.balancelaw
  grid = dg.grid
  topology = grid.topology
  @assert isstacked(topology)
  @assert typeof(dg.direction) <: VerticalDirection

  FT = eltype(Q.data)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  nstate = num_state(bl, FT)
  N = polynomialorder(grid)
  Nq = N + 1

  # p is lower bandwidth
  # q is upper bandwidth
  eband = num_diffusive(bl, FT) == 0 ? 1 : 2
  p = q = nstate * Nq * eband - 1

  nrealelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nrealelem, nvertelem)

  dim = dimensionality(grid)

  Nqj = dim == 2 ? 1 : Nq

  # first horizontal DOF index
  # second horizontal DOF index
  # band index -q:p
  # vertical DOF index
  # horizontal element index
  A = if single_column
    similar(Q.data, p + q + 1, Nq * nstate * nvertelem)
  else
    similar(Q.data, Nq, Nqj, p + q + 1, Nq * nstate * nvertelem,
            nhorzelem)
  end
  fill!(A, zero(FT))

  # loop through all DOFs in a column and compute the matrix column
  for ev = 1:nvertelem
    for s = 1:nstate
      for k = 1:Nq
        # Set a single 1 per column and rest 0
        @launch(device, threads=(Nq, Nqj, Nq), blocks=(nvertelem, nhorzelem),
                knl_set_banded_data!(bl, Val(dim), Val(N), Val(nvertelem),
                                     Q.data, k, s, ev, 1:nhorzelem,
                                     1:nvertelem))

        # Get the matrix column
        f!(dQ, Q, args...)

        # Store the banded matrix
        @launch(device, threads=(Nq, Nqj, Nq),
                blocks=(2 * eband + 1, nhorzelem),
                knl_set_banded_matrix!(bl, Val(dim), Val(N), Val(nvertelem),
                                       Val(p), Val(q), Val(eband+1),
                                       A, dQ.data, k, s, ev, 1:nhorzelem,
                                       -eband:eband))
      end
    end
  end
  A
end

"""
    banded_matrix_vector_product!(dg::DGModel, A, dQ::MPIStateArray,
                                  Q::MPIStateArray)

Compute a matrix vector product `dQ = A * Q` where `A` is assumed to be a matrix
created using the `banded_matrix` function.

This function is primarily for testing purposes.
"""
function banded_matrix_vector_product!(dg::DGModel, A, dQ::MPIStateArray,
                                       Q::MPIStateArray)
  bl = dg.balancelaw
  grid = dg.grid
  topology = grid.topology
  @assert isstacked(topology)
  @assert typeof(dg.direction) <: VerticalDirection

  FT = eltype(Q.data)
  device = typeof(Q.data) <: Array ? CPU() : CUDA()

  eband = num_diffusive(bl, FT) == 0 ? 1 : 2
  nstate = num_state(bl, FT)
  N = polynomialorder(grid)
  Nq = N + 1
  p = q = nstate * Nq * eband - 1

  nrealelem = length(topology.elems)
  nvertelem = topology.stacksize
  nhorzelem = div(nrealelem, nvertelem)

  dim = dimensionality(grid)

  Nqj = dim == 2 ? 1 : Nq

  @launch(device, threads=(Nq, Nqj, Nq),
          blocks=(nvertelem, nhorzelem),
          knl_banded_matrix_vector_product!(bl, Val(dim), Val(N),
                                            Val(nvertelem), Val(p), Val(q),
                                            dQ.data, A, Q.data, 1:nhorzelem,
                                            1:nvertelem))
end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
function _row_echelon(A; L = [])
    A, l = SymbolicWedderburn.row_echelon_form(A)
    @test all([x ∈ l for x in L]) && length(L) <= length(l)
    for (i, j) in enumerate(l)
        @test all([iszero(A[k, j]) for k = 1:size(A, 1) if k != i])
    end
    @test all([isone(A[i, j]) for (i, j) in enumerate(l)])
end

function _nullspace(A)
    V = SymbolicWedderburn.left_nullspace(A)
    W = SymbolicWedderburn.right_nullspace(A)
    @test all(iszero.(V*A))
    @test all(iszero.(A*W))
end

function _left_eigen(A)
    E = SymbolicWedderburn.left_eigen(A)
    for (val, space) in E
        @test space * A == space .* val
    end
end

@testset "RowEchelon" begin
    @testset "examples" begin
        A = GF{11}.([ 0 1 2; 1 2 1; 2 7 8])
        L = [1, 2]
        _row_echelon(A; L = L)

        A = GF{7}.([0 1; 1 2; 0 5])
        L = [1, 2]
        _row_echelon(A; L = L)

        A = GF{2}.([ 0 1 0 0; 0 0 1 1; 0 1 0 1; 0 0 0 0])
        L = [2, 3, 4]
        _row_echelon(A, L = L)

        #Handbook of Computational Group Theory p.252
        A = GF{3}.([1 0 2 0 0 0 0 0;
                    2 1 0 0 0 1 2 0;
                    0 2 1 0 0 0 1 2;
                    0 0 0 0 0 2 0 1])

        N = GF{3}.([1 0 2 0 0 0 0 0;
                    0 1 2 0 0 0 2 1;
                    0 0 0 0 0 1 0 2;
                    0 0 0 0 0 0 0 0])

     @test SymbolicWedderburn.row_echelon_form(A)[1] == N

    end
    @testset "random" begin
        for p in [2, 5, 7]
            for s in [3, 5, 6]
                for t in [3, 5, 9]
                    _row_echelon(GF{p}.(rand(Int8, s, t)))
                end
            end
        end
    end
end

@testset "NullSpace" begin
    @testset "examples" begin
        A = [2 1 0 2; 0 0 1 3; 0 1 0 3; 0 0 0 0]
        for p in [2, 3, 5]
            _nullspace(GF{p}.(A))
        end

        A  = [1 2 0 1 4; 0 1 2 1 3; 1 2 1 3 5; 0 1 2 2 7]
        for p in [2, 3, 5]
            _nullspace(GF{p}.(A))
        end

        A = GF{2}.([ 0 1 0 0; 0 0 1 1; 0 1 0 1; 0 0 0 0])
        _nullspace(A)

        #Handbook of Computational Group Theory p.252
        A = GF{3}.([1 1 1 0 0 0 0 0 2 0 2 0;
                    1 1 1 0 0 0 0 0 0 2 0 2;
                    1 1 1 0 0 0 0 0 2 0 2 0;
                    0 0 0 0 0 0 0 0 0 2 0 2;
                    0 0 0 0 0 0 0 0 1 0 1 0;
                    0 0 0 0 0 1 1 1 0 1 0 1;
                    0 0 0 0 0 1 1 1 1 0 1 0;
                    0 0 0 0 0 1 1 1 0 1 0 1])
        N = GF{3}.([1 0 2 0 0 0 0 0;
                    0 1 2 0 0 0 2 1;
                    0 0 0 1 1 0 2 1;
                    0 0 0 0 0 1 0 2])

        K = SymbolicWedderburn.left_nullspace(A)
        @test SymbolicWedderburn.row_echelon_form(K)[1] == N

   end
    @testset "random" begin
        for p in [2, 5, 7]
            for s in [3, 5, 6]
                for t in [3, 5, 9]
                    _nullspace(GF{p}.(rand(Int8, s, t)))
                end
            end
        end
    end
end

@testset "LeftEigen" begin
    @testset "examples" begin
        #Handbook of Computational Group Theory p.259
        M2 = GF{7}.([0 1 0 0;
                     3 2 0 0;
                     0 0 3 0;
                     0 0 0 3])
        M3 = GF{7}.([0 0 1 0;
                     0 0 3 0;
                     0 0 0 4;
                     4 4 0 0])
        M4 = GF{7}.([0 0 0 1;
                     0 0 0 3;
                     4 4 0 0;
                     0 0 4 0])

        E2 = SymbolicWedderburn.left_eigen(M2)
        @test size(E2[GF{7}(3)], 1) == 3

        E3 = SymbolicWedderburn.left_eigen(M3)
        E4 = SymbolicWedderburn.left_eigen(M4)
        for val in GF{7}.([0, -3, -5, -6])
            @test haskey(E3, val)
            @test haskey(E4, val)
        end

        @test E3[GF{7}(0)] == GF{7}.([1 2 0 0])
        @test E3[GF{7}(-3)] == GF{7}.([1 1 1 1])
        @test E3[GF{7}(-5)] == GF{7}.([1 1 2 4])
        @test E3[GF{7}(-6)] == GF{7}.([1 1 4 2])

    end

    @testset "random" begin
        for p in [2, 5, 7]
            for s in [3, 5, 6, 9]
                _left_eigen(GF{p}.(rand(Int8, s, s)))
            end
        end
    end

end

@testset "EigSpaceDec" begin
      @testset "examples" begin
        #Handbook of Computational Group Theory p.259
        M3 = GF{7}.([0 0 1 0;
                     0 0 3 0;
                     0 0 0 4;
                     4 4 0 0])
        esd = SymbolicWedderburn.EigenSpaceDecomposition(M3)
        @test SymbolicWedderburn.isdiag(esd)

        #Handbook of Computational Group Theory p.262
        M6 = GF{13}.([0 0 0 0 0 1;
                      0 0 1 0 0 0;
                      0 2 0 0 0 1;
                      0 0 0 0 2 0;
                      0 0 0 2 0 0;
                      2 0 1 0 0 0])
        esd = SymbolicWedderburn.EigenSpaceDecomposition(M6)

        @test length(esd) == 4
        @test esd[4] == GF{13}.([1 1 6 0 0 6])
        @test esd[3] == GF{13}.([1 12 1 0 0 12; 0 0 0 1 12 0])
        @test esd[2] == GF{13}.([1 1 1 0 0 1; 0 0 0 1 1 0])
        @test esd[1] == GF{13}.([1 12 6 0 0 7])

        M4 = GF{13}.([0 0 0 1 0 0;
                      0 0 0 0 1 0;
                      0 0 0 2 0 0;
                      0 3 0 0 0 3;
                      3 0 3 0 0 0;
                      0 0 0 0 2 0])

        esd = SymbolicWedderburn.refine(esd, M4)
        @test SymbolicWedderburn.isdiag(esd)
        @test esd[1] == GF{13}.([1 12 6 0 0 7])
        @test esd[2] == GF{13}.([1 1 1 1 1 1])
        @test esd[3] == GF{13}.([1 1 1 12 12 1])
        @test esd[4] == GF{13}.([1 12 1 8 5 12])
        @test esd[5] == GF{13}.([1 12 1 5 8 12])
        @test esd[6] == GF{13}.([1 1 6 0 0 6])

    end
end


text:
#!/usr/bin/env julia

import GR
const gr3 = GR.gr3

function draw(mesh::Cint)
    gr3.clear()
    gr3.drawmesh(mesh, 1, (0,0,0), (0,0,1), (0,1,0), (1,1,1), (1,1,1))
    GR.clearws()
    gr3.drawimage(0, 1, 0, 1, 500, 500, gr3.DRAWABLE_GKS)
    GR.updatews()
end

data = open(stream -> read!(stream, Array{UInt16}(undef, 93, 64, 64)), "mri.raw")
data = min.(data, 2000) / 2000.0 * typemax(UInt16)
data = convert(Array{UInt16, 3}, floor.(data))
data = permutedims(data, [3, 2, 1])

GR.setviewport(0, 1, 0, 1)
gr3.cameralookat(-3, 2, -2, 0, 0, 0, 0, 0, -1)

mesh = gr3.createisosurfacemesh(data, (2.0/63, 2.0/63, 2.0/92), (-1.0, -1.0, -1.0), 40000)

GR.setcolormap(1)

gr3.setlightsources(2, [0, -1, 1, 1, -1, 1], [0.45, 0.7, 0.9, 0.93, 0.91, 0.2])
directions, colors = gr3.getlightsources(2)
@show directions, colors

draw(mesh)


text:
using SparseArrays

include("GloVe.jl")

cooc_file = ARGS[1]
pert_file = ARGS[2]
out_file = ARGS[3]
vocab_file = ARGS[4]

vocab, ivocab = GloVe.load_vocab(vocab_file)
V = length(vocab)

δX = GloVe.load_cooc(pert_file, V)
@show nnz(δX)
println("here")

GloVe.perturb_coocs(cooc_file, out_file, -δX)


text:
module DOMTests

using Test
using Compat

import Documenter.Utilities.DOM: DOM, @tags, HTMLDocument

@tags div ul li p

@testset "DOM" begin
    for tag in (:div, :ul, :li, :p)
        TAG = @eval $tag
        @test isa(TAG, DOM.Tag)
        @test TAG.name === tag
    end

    @test div().name === :div
    @test div().text == ""
    @test isempty(div().nodes)
    @test isempty(div().attributes)

    @test div("...").name === :div
    @test div("...").text == ""
    @test length(div("...").nodes) === 1
    @test div("...").nodes[1].text == "..."
    @test div("...").nodes[1].name === Symbol("")
    @test isempty(div("...").attributes)

    @test div[".class"]("...").name === :div
    @test div[".class"]("...").text == ""
    @test length(div[".class"]("...").nodes) === 1
    @test div[".class"]("...").nodes[1].text == "..."
    @test div[".class"]("...").nodes[1].name === Symbol("")
    @test length(div[".class"]("...").attributes) === 1
    @test div[".class"]("...").attributes[1] == (:class => "class")
    @test div[:attribute].attributes[1] == (:attribute => "")
    @test div[:attribute => "value"].attributes[1] == (:attribute => "value")

    let d = div(ul(map(li, [string(n) for n = 1:10])))
        @test d.name === :div
        @test d.text == ""
        @test isempty(d.attributes)
        @test length(d.nodes) === 1
        let u = d.nodes[1]
            @test u.name === :ul
            @test u.text == ""
            @test isempty(u.attributes)
            @test length(u.nodes) === 10
            for n = 1:10
                let v = u.nodes[n]
                    @test v.name === :li
                    @test v.text == ""
                    @test isempty(v.attributes)
                    @test length(v.nodes) === 1
                    @test v.nodes[1].name === Symbol("")
                    @test v.nodes[1].text == string(n)
                    @test !isdefined(v.nodes[1], :attributes)
                    @test !isdefined(v.nodes[1], :nodes)
                end
            end
        end
    end

    @tags script style img

    @test string(div(p("one"), p("two"))) == "

one

two

"
    @test string(div[:key => "value"])    == "

"
    @test string(p(" < > & ' \" "))       == "
 < > & ' "
"
    @test string(img[:src => "source"])   == ""
    @test string(img[:none])              == ""
    @test string(script(" < > & ' \" "))  == ""
    @test string(style(" < > & ' \" "))   == ""
    @test string(script)                  == "


text:
#
# Spline
#

struct Spline
  x :: Vector{Float64}
  y :: Vector{Float64}
end

function fitspline(X :: AbstractVector, Y :: AbstractVector, options :: Options)
  X, Y = checkdata(X,Y)
  t = 1:length(X)
  A = hcat(X,Y)
  itp = Interpolations.scale(interpolate(A,
            (BSpline(Interpolations.Cubic(Natural(OnGrid()))), NoInterp())), t, 1:2)
  tfine = 1:length(X)/options.fine:length(X)
  x, y = [itp(t,1) for t in tfine], [itp(t,2) for t in tfine]
  return Spline(x,y)
end
fitspline(X :: AbstractVector, Y :: AbstractVector) = fitspline(X,Y,Options())

function Base.show( io :: IO, fit :: Spline )
  println("")
  println(" -------- Spline fit --------------------------- ")
  println("")
  println(" x spline: x = [",fit.x[1],", ",fit.x[2],"...")
  println(" y spline: y = [",fit.y[1],", ",fit.y[2],"...")
  println("")
  println(" ----------------------------------------------- ")
end

export fitspline


text:
@defcomp scenconverter begin
    regions = Index()

    population = Parameter(index=[time,regions])
    income = Parameter(index=[time,regions])
    energuse = Parameter(index=[time,regions])
    emission = Parameter(index=[time,regions])

    pop0 = Variable(index=[regions])
    scenpgrowth = Variable(index=[time,regions])
    gdp0 = Variable(index=[regions])
    scenypcgrowth = Variable(index=[time,regions])
    scenaeei = Variable(index=[time,regions])
    scenacei = Variable(index=[time,regions])
    energint0 = Variable(index=[regions])
    emissint0 = Variable(index=[regions])
    energint = Variable(index=[time,regions])
    emissint = Variable(index=[time,regions])

    function run_timestep(p, v, d, t)

        if is_first(t)
            for r in d.regions
                v.pop0[r] = p.population[t, r]
                v.gdp0[r] = p.income[t, r]
                v.energint0[r] = p.energuse[t,r] / p.income[t, r]
                v.emissint0[r]  = p.emission[t, r] / p.energuse[t, r]

                v.energint[t, r] = v.energint0[r]
                v.emissint[t, r] = v.emissint0[r]
            end
        end

        for r in d.regions
            if t.t < 1050
                v.scenpgrowth[t, r] = (p.population[t+1, r] / p.population[t, r] - 1.) * 100.
                v.scenypcgrowth[t, r] = (p.income[t+1, r] / p.income[t, r] / (1 + 0.01 * v.scenpgrowth[t, r]) - 1.) * 100.
            end

            v.energint[t, r] = p.energuse[t, r] / p.income[t,r]
            v.emissint[t, r] = p.emission[t, r] / p.energuse[t, r]

            if !is_first(t)
                v.scenaeei[t, r] = -(v.energint[t, r] / v.energint[t - 1, r] - 1.) * 100.
                v.scenacei[t, r] = -(v.emissint[t, r] / v.emissint[t - 1, r] - 1.) * 100.
            end
        end
    end
end


text:
# precise values were computed with Maple 2018.2.1
@testset "Type stability" begin
  @testset for F = [Float32, Float64, BigFloat]
    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(1)], 1, 0, F(1))
    @test isa(y, Symmetric{F, Matrix{F}})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(1)], 1, 0, F(0))
    @test isa(y, Symmetric{F, Matrix{F}})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(0), F(1)], 1, 0, F(1))
    @test isa(y, Symmetric{F, Matrix{F}})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(1), F(0)], 1, 0, F(1))
    @test isa(y, Symmetric{F, Matrix{F}})

    x, y = SimpleBirthDeathProcess.gradient_hessian([F(0), F(0)], 1, 0, F(1))
    @test isa(y, Symmetric{F, Matrix{F}})
  end
end

@testset "Parameters are zero" begin
  x, y = SimpleBirthDeathProcess.gradient_hessian([1.0, 1.0], 1, 0, 0.0)
  @test y[1, 1] === 0.0
  @test y[1, 2] === 0.0
  @test y[2, 1] === 0.0
  @test y[2, 2] === 0.0

  x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, 0.0], 1, 0, 1.0)
  @test y[1, 1] === NaN
  @test y[1, 2] === NaN
  @test y[2, 1] === NaN
  @test y[2, 2] === NaN

  x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, 0.0], 2, 2, 1.0)
  @test y[1, 1] === 0.0
  @test y[1, 2] === 4.0
  @test y[2, 1] === 4.0
  @test y[2, 2] === 0.0
end

@testset "Death rate μ = 0" begin
  u = ((1, 0), (1, 1), (1, 100), (1000, 0), (1000, 500), (1000, 998), (1000, 999), (1000, 1000), (1000, 2000))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λλ1 = [999999.99999991667, 99.999999916666667, 0.99999991666667083, 0.99999166670833317e-2, 0.99916708316804727e-4, 999999.99916666667, 99.999166670833317, 0.99916708316804727, 0.92067359420779232e-2, 0.45404052350475412e-6, 999999.91666667083, 99.916708316804727, 0.92067359420779232, 0.45404052350475412e-4, 0.37200759760208360e-43, 999991.66670833317, 92.067359420779232, 0.45404052350475412e-2, 0.37200759760208360e-41, 999167.08316804727, 0.45404052350475412, 0.37200759760208360e-39]
  λμ1 = fill(Float64(NaN), 22)
  μμ1 = fill(Float64(-Inf), 22)

  # case i = 1, j = 1
  λλ2 = zeros(Float64, 22)
  λμ2 = [0.99999933333358333e-6, 0.99993333583326667e-6, 0.99933358326668055e-6, 0.99335826680531781e-6, 0.93576803208889390e-6, 0.99993333583326667e-2, 0.99335826680531781e-2, 0.93576803208889390e-2, 0.52848223531423071e-2, 0.19990012015452253e-3, 0.99933358326668055, 0.93576803208889390, 0.52848223531423071, 0.19990012015452253e-1, 0.20000000000000000e-3, 99.335826680531781, 52.848223531423071, 1.9990012015452253, 0.20000000000000000e-1, 9357.6803208889390, 199.90012015452253, 2.0000000000000000]
  μμ2 = [-0.66666600000036667e-12, -0.66660000366652223e-10, -0.66600036652226745e-9, -0.66003652267341542e-8, -0.60352662965245339e-7, -0.66660000366652223e-6, -0.66003652267341542e-4, -0.60352662965245339e-3, -0.25781166884100533e-2, -0.19981839986871934e-3, -0.66600036652226745e-3, -0.60352662965245339e-1, -0.25781166884100533, -0.19981839986871934e-1, -0.20000000000000000e-3, -0.66003652267341542, -25.781166884100533, -1.9981839986871934, -0.20000000000000000e-1, -603.52662965245339, -199.81839986871934, -2.0000000000000000]

  # case i = 1, j = 100
  λλ3 = [-98999999.999991750, -9899.9999917500000, -98.999991750000412, -0.98999175004124984, -0.98917541233636680e-2, -98999999.917500000, -9899.9175004124984, -98.917541233636680, -0.91146685826571440, -0.44950011826970658e-4, -98999991.750000412, -9891.7541233636680, -91.146685826571440, -0.44950011826970658e-2, -0.36828752162606276e-41, -98999175.004124984, -9114.6685826571440, -0.44950011826970658, -0.36828752162606276e-39, -98917541.233636680, -44.950011826970658, -0.36828752162606276e-37]
  λμ3 = [0.42249966333346371e-4, 0.42246633463704967e-4, 0.42216346367467352e-4, 0.41914633723504859e-4, 0.39010408984157111e-4, 0.42246633463704967, 0.41914633723504859, 0.39010408984157111, 0.18835038709940091, 0.23990607963035859e-3, 42.216346367467352, 39.010408984157111, 18.835038709940091, 0.23990607963035859e-1, 0.20000000000000000e-3, 4191.4633723504859, 1883.5038709940091, 2.3990607963035859, 0.20000000000000000e-1, 390104.08984157111, 239.90607963035859, 2.0000000000000000]
  μμ3 = [0.82499663333665875e-5, 0.82466336658564841e-5, 0.82163665689906419e-5, 0.79166403062136151e-5, 0.51980671565871416e-5, 0.82466336658564841e-1, 0.79166403062136151e-1, 0.51980671565871416e-1, -0.51661751030422087e-1, -0.23577920519729713e-3, 8.2163665689906419, 5.1980671565871416, -5.1661751030422087, -0.23577920519729713e-1, -0.20000000000000000e-3, 791.66403062136151, -516.61751030422087, -2.3577920519729713, -0.20000000000000000e-1, 51980.671565871416, -235.77920519729713, -2.0000000000000000]

  # case i = 1000, j = 0
  λλ4 = [999999999.99991667, 99999.999916666667, 999.99991666667083, 9.9999166670833317, 0.99916708316804727e-1, 999999999.16666667, 99999.166670833317, 999.16708316804727, 9.2067359420779232, 0.45404052350475412e-3, 999999916.66667083, 99916.708316804727, 920.67359420779232, 0.45404052350475412e-1, 0.37200759760208360e-40, 999991666.70833317, 92067.359420779232, 4.5404052350475412, 0.37200759760208360e-38, 999167083.16804727, 454.04052350475412, 0.37200759760208360e-36]
  λμ4 = [-0.41666700000012917e-3, -0.41670000129170000e-3, -0.41700012920000678e-3, -0.42001295006791023e-3, -0.45132568999644072e-3, -4.1670000129170000, -4.2001295006791023, -4.5132568999644072, -9.2067359420779232, -19823.819669366568, -417.00012920000678, -451.32568999644072, -920.67359420779232, -1982381.9669366568, -0.26612359703979741e45, -42001.295006791023, -92067.359420779232, -198238196.69366568, -0.26612359703979741e47, -4513256.8999644072, -19823819669.366568, -0.26612359703979741e49]
  μμ4 = [-0.83333666666995834e-4, -0.83366669958516674e-4, -0.83666996016739077e-4, -0.86699767392810193e-4, -0.12014913769930394e-3, -0.83366669958516674, -0.86699767392810193, -1.2014913769930394, -10.318188478047674, -48472466.608935374, -83.666996016739077, -120.14913769930394, -1031.8188478047674, -4847246660.8935374, -0.72259737681257493e86, -8669.9767392810193, -103181.88478047674, -484724666089.35374, -0.72259737681257493e88, -1201491.3769930394, -48472466608935.374, -0.72259737681257493e90]

  # case i = 1000, j = 500
  λλ5 = [499999999.99995833, 49999.999958333333, 499.99995833333542, 4.9999583335416658, 0.49958354158402363e-1, 499999999.58333333, 49999.583335416658, 499.58354158402363, 4.6033679710389616, 0.22702026175237706e-3, 499999958.33333542, 49958.354158402363, 460.33679710389616, 0.22702026175237706e-1, 0.18600379880104180e-40, 499995833.35416658, 46033.679710389616, 2.2702026175237706, 0.18600379880104180e-38, 499583541.58402363, 227.02026175237706, 0.18600379880104180e-36]
  λμ5 = [-0.12103378413496805e-2, -0.12103873445002588e-2, -0.12108376564837180e-2, -0.12153688999561409e-2, -0.12635399336281270e-2, -12.103873445002588, -12.153688999561409, -12.635399336281270, -20.949909052664360, -39687.566126019797, -1210.8376564837180, -1263.5399336281270, -2094.9909052664360, -3968756.6126019797, -0.53278050790131786e45, -121536.88999561409, -209499.09052664360, -396875661.26019797, -0.53278050790131786e47, -12635399.336281270, -39687566126.019797, -0.53278050790131786e49]
  μμ5 = [-0.41667166665481152e-4, -0.41716654811796964e-4, -0.42165481427042550e-4, -0.46548388668952514e-4, -0.80071308902245171e-4, -0.41716654811796964, -0.46548388668952514, -0.80071308902245171, 4.9900093079764706, 851359.41856749361, -42.165481427042550, -80.071308902245171, 499.00093079764706, 85135941.856749361, 0.11631264148506534e85, -4654.8388668952514, 49900.093079764706, 8513594185.6749361, 0.11631264148506534e87, -800713.08902245171, 851359418567.49361, 0.11631264148506534e89]

  # case i = 1000, j = 998
  λλ6 = [1999999.9999998333, 199.99999983333333, 1.9999998333333417, 0.19999833334166663e-1, 0.19983341663360945e-3, 1999999.9983333333, 199.99833334166663, 1.9983341663360945, 0.18413471884155846e-1, 0.90808104700950824e-6, 1999999.8333333417, 199.83341663360945, 1.8413471884155846, 0.90808104700950824e-4, 0.74401519520416719e-43, 1999983.3334166663, 184.13471884155846, 0.90808104700950824e-2, 0.74401519520416719e-41, 1998334.1663360945, 0.90808104700950824, 0.74401519520416719e-39]
  λμ6 = [-0.99800083399958283, -0.99800090242833598, -0.99800174883367213, -0.99803244356437850, -1.0005638889995014, -9980.0090242833598, -9980.3244356437850, -10005.638889995014, -12624.487641335677, -19804194.946321424, -998001.74883367213, -1000563.8889995014, -1262448.7641335677, -1980419494.6321424, -0.26585747344275761e48, -99803244.356437850, -126244876.41335677, -198041949463.21424, -0.26585747344275761e50, -10005638889.995014, -19804194946321.424, -0.26585747344275761e52]
  μμ6 = fill(Float64(NaN), 22)

  # case i = 1000, j = 999
  λλ7 = [999999.99999991667, 99.999999916666667, 0.99999991666667083, 0.99999166670833317e-2, 0.99916708316804727e-4, 999999.99916666667, 99.999166670833317, 0.99916708316804727, 0.92067359420779232e-2, 0.45404052350475412e-6, 999999.91666667083, 99.916708316804727, 0.92067359420779232, 0.45404052350475412e-4, 0.37200759760208360e-43, 999991.66670833317, 92.067359420779232, 0.45404052350475412e-2, 0.37200759760208360e-41, 999167.08316804727, 0.45404052350475412, 0.37200759760208360e-39]
  λμ7 = fill(Float64(NaN), 22)
  μμ7 = fill(Float64(-Inf), 22)

  # case i = 1000, j = 1000
  λλ8 = zeros(Float64, 22)
  λμ8 = [0.99999999933358333, 0.99999993583333327, 0.99999958333328056, 1.0000183334055557, 1.0024346558790135, 9999.9993583333327, 10000.183334055557, 10024.346558790135, 12635.053587747725, 19804195.346120856, 999999.58333328056, 1002434.6558790135, 1263505.3587747725, 1980419534.6120856, 0.26585747344275761e48, 100001833.34055557, 126350535.87747725, 198041953461.20856, 0.26585747344275761e50, 10024346558.790135, 19804195346120.856, 0.26585747344275761e52]
  μμ8 = [-0.50066550000008370e-6, -0.50000550083336772e-2, -0.49999958333347847, -50.000723337672237, -5008.3279637162128, -50.000550083336772, -500007.23337672237, -50083279.637162128, -5898717394.4429521, -24253782099472752., -499999.58333347847, -5008327963.7162128, -589871739444.29521, -0.24253782099472752e19, -0.36129760523281962e95, -5000072333.7672237, -58987173944429.521, -0.24253782099472752e21, -0.36129760523281962e97, -50083279637162.128, -0.24253782099472752e23, -0.36129760523281962e99]

  # case i = 1000, j = 2000
  λλ9 = [-999999999.99991667, -99999.999916666667, -999.99991666667083, -9.9999166670833317, -0.99916708316804727e-1, -999999999.16666667, -99999.166670833317, -999.16708316804727, -9.2067359420779232, -0.45404052350475412e-3, -999999916.66667083, -99916.708316804727, -920.67359420779232, -0.45404052350475412e-1, -0.37200759760208360e-40, -999991666.70833317, -92067.359420779232, -4.5404052350475412, -0.37200759760208360e-38, -999167083.16804727, -454.04052350475412, -0.37200759760208360e-36]
  λμ9 = [0.24146676646692933e-2, 0.24145686709552363e-2, 0.24136692932358654e-2, 0.24047314317302164e-2, 0.23208587543955349e-2, 24.145686709552363, 24.047314317302164, 23.208587543955349, 19.751121085572186, 19784.611039706006, 2413.6692932358654, 2320.8587543955349, 1975.1121085572186, 1978461.1039706006, 0.26559188156119642e45, 240473.14317302164, 197511.21085572186, 197846110.39706006, 0.26559188156119642e47, 23208587.543955349, 19784611039.706006, 0.26559188156119642e49]
  μμ9 = [0.83332333331012451e-4, 0.83233310123964467e-4, 0.82331011901146989e-4, 0.73100692111297332e-4, -0.40456453246960823e-4, 0.83233310123964467, 0.73100692111297332, -0.40456453246960823, -41.832586470892298, -144784500.99608062, 82.331011901146989, -40.456453246960823, -4183.2586470892298, -14478450099.608062, -0.21562651734461104e87, 7310.0692111297332, -418325.86470892298, -1447845009960.8062, -0.21562651734461104e89, -404564.53246960823, -144784500996080.62, -0.21562651734461104e91]

  λλ = [λλ1 λλ2 λλ3 λλ4 λλ5 λλ6 λλ7 λλ8 λλ9]
  λμ = [λμ1 λμ2 λμ3 λμ4 λμ5 λμ6 λμ7 λμ8 λμ9]
  μμ = [μμ1 μμ2 μμ3 μμ4 μμ5 μμ6 μμ7 μμ8 μμ9]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, 0.0], i, j, t)
      @test y[1, 1] ≈ λλ[a, b]

      if (b == 1) || (b == 7)
        @test isnan(y[1, 2])
        @test isnan(y[2, 1])
        @test isinf(y[2, 2])
      elseif b == 6
        @test y[1, 2] ≈ λμ[a, b]
        @test y[2, 1] ≈ λμ[a, b]
        @test isnan(y[2, 2])
      else
        @test y[1, 2] ≈ λμ[a, b]
        @test y[2, 1] ≈ λμ[a, b]
        @test y[2, 2] ≈ μμ[a, b]
      end
    end
  end
end

@testset "Birth rate λ = 0" begin
  u = ((1, 0), (1, 1), (1, 2), (1, 3), (1, 100),
       (1000, 0), (1000, 500), (1000, 1000), (1000, 1001), (1000, 1002), (1000, 2000))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λλ1 = [.83333000000329166e-7, .83300003291483341e-7, .83000328983405698e-7, .80032734054953023e-7, .53115351712650373e-7, .83300003291483341e-3, .80032734054953023e-3, .53115351712650373e-3, -.49579428628294984e-3, -.36324045786442217e-6, .83000328983405698e-1, .53115351712650373e-1, -.49579428628294984e-1, -.36324045786442217e-4, -.36456744565004192e-43, 8.0032734054953023, -4.9579428628294984, -.36324045786442217e-2, -.36456744565004192e-41, 531.15351712650373, -.36324045786442217, -.36456744565004192e-39]
  λμ1 = [.41666633333346250e-6, .41663333462496667e-6, .41633346246667344e-6, .41334621673433880e-6, .38459233284917391e-6, .41663333462496667e-2, .41334621673433880e-2, .38459233284917391e-2, .18491471186490768e-2, .40410060076602078e-6, .41633346246667344, .38459233284917391, .18491471186490768, .40410060076602078e-4, .36825032086630255e-43, 41.334621673433880, 18.491471186490768, .40410060076602078e-2, .36825032086630255e-41, 3845.9233284917391, .40410060076602078, .36825032086630255e-39]
  μμ1 = [-999999.99999991667, -99.999999916666667, -.99999991666667083, -.99999166670833317e-2, -.99916708316804727e-4, -999999.99916666667, -99.999166670833317, -.99916708316804727, -.92067359420779232e-2, -.45404052350475412e-6, -999999.91666667083, -99.916708316804727, -.92067359420779232, -.45404052350475412e-4, -.37200759760208360e-43, -999991.66670833317, -92.067359420779232, -.45404052350475412e-2, -.37200759760208360e-41, -999167.08316804727, -.45404052350475412, -.37200759760208360e-39]

  # case i = 1, j = 1
  λλ2 = [-.66666600000036667e-12, -.66660000366652223e-10, -.66600036652226745e-9, -.66003652267341542e-8, -.60352662965245339e-7, -.66660000366652223e-6, -.66003652267341542e-4, -.60352662965245339e-3, -.25781166884100533e-2, -.19981839986871934e-3, -.66600036652226745e-3, -.60352662965245339e-1, -.25781166884100533, -.19981839986871934e-1, -.20000000000000000e-3, -.66003652267341542, -25.781166884100533, -1.9981839986871934, -.20000000000000000e-1, -603.52662965245339, -199.81839986871934, -2.0000000000000000]
  λμ2 = [.99999933333358333e-6, .99993333583326667e-6, .99933358326668055e-6, .99335826680531781e-6, .93576803208889390e-6, .99993333583326667e-2, .99335826680531781e-2, .93576803208889390e-2, .52848223531423071e-2, .19990012015452253e-3, .99933358326668055, .93576803208889390, .52848223531423071, .19990012015452253e-1, .20000000000000000e-3, 99.335826680531781, 52.848223531423071, 1.9990012015452253, .20000000000000000e-1, 9357.6803208889390, 199.90012015452253, 2.0000000000000000]
  μμ2 = zeros(Float64, 22)

  # case i = 1, j = 2
  λλ3 = fill(Float64(-Inf), 22)
  λμ3 = fill(Float64(NaN), 22)
  μμ3 = [999999.99999991667, 99.999999916666667, .99999991666667083, .99999166670833317e-2, .99916708316804727e-4, 999999.99916666667, 99.999166670833317, .99916708316804727, .92067359420779232e-2, .45404052350475412e-6, 999999.91666667083, 99.916708316804727, .92067359420779232, .45404052350475412e-4, .37200759760208360e-43, 999991.66670833317, 92.067359420779232, .45404052350475412e-2, .37200759760208360e-41, 999167.08316804727, .45404052350475412, .37200759760208360e-39]

  # case i = 1, j = 3
  λλ4 = fill(Float64(NaN), 22)
  λμ4 = [-.38333346666676750e-5, -.38334666767501333e-5, -.38346676751333886e-5, -.38467676338860596e-5, -.39768889048486304e-5, -.38334666767501333e-1, -.38467676338860596e-1, -.39768889048486304e-1, -.63698294237298154e-1, -118.94351610015747, -3.8346676751333886, -3.9768889048486304, -6.3698294237298154, -11894.351610015747, -.15967415822387845e43, -384.67676338860596, -636.98294237298154, -1189435.1610015747, -.15967415822387845e45, -39768.889048486304, -118943516.10015747, -.15967415822387845e47]
  μμ4 = [1999999.9999998333, 199.99999983333333, 1.9999998333333417, .19999833334166663e-1, .19983341663360945e-3, 1999999.9983333333, 199.99833334166663, 1.9983341663360945, .18413471884155846e-1, .90808104700950824e-6, 1999999.8333333417, 199.83341663360945, 1.8413471884155846, .90808104700950824e-4, .74401519520416719e-43, 1999983.3334166663, 184.13471884155846, .90808104700950824e-2, .74401519520416719e-41, 1998334.1663360945, .90808104700950824, .74401519520416719e-39]

  # case i = 1, j = 100
  λλ5 = [-.82500336666999132e-5, -.82533669991503648e-5, -.82836999317087239e-5, -.85900099084723869e-5, -.11968402995881409e-4, -.82533669991503648e-1, -.85900099084723869e-1, -.11968402995881409, -1.0404678434243096, -4895677.3769611642, -8.2836999317087239, -11.968402995881409, -104.04678434243096, -489567737.69611642, -.72981723836443202e85, -859.00099084723869, -10404.678434243096, -48956773769.611642, -.72981723836443202e87, -119684.02995881409, -4895677376961.1642, -.72981723836443202e89]
  λμ5 = [-.42270441829945015e-4, -.42273774960361360e-4, -.42304087875901366e-4, -.42608382463667796e-4, -.45770937268056287e-4, -.42273774960361360, -.42608382463667796, -.45770937268056287, -.93172486665707299, -2002.6105575176441, -42.304087875901366, -45.770937268056287, -93.172486665707299, -200261.05575176441, -.26883914394836677e44, -4260.8382463667796, -9317.2486665707299, -20026105.575176441, -.26883914394836677e46, -457709.37268056287, -2002610557.5176441, -.26883914394836677e48]
  μμ5 = [98999999.999991750, 9899.9999917500000, 98.999991750000412, .98999175004124984, .98917541233636680e-2, 98999999.917500000, 9899.9175004124984, 98.917541233636680, .91146685826571440, .44950011826970658e-4, 98999991.750000412, 9891.7541233636680, 91.146685826571440, .44950011826970658e-2, .36828752162606276e-41, 98999175.004124984, 9114.6685826571440, .44950011826970658, .36828752162606276e-39, 98917541.233636680, 44.950011826970658, .36828752162606276e-37]

  # case i = 1000, j = 0
  λλ6 = [.83333000000329166e-4, .83300003291483341e-4, .83000328983405698e-4, .80032734054953023e-4, .53115351712650373e-4, .83300003291483341, .80032734054953023, .53115351712650373, -.49579428628294984, -.36324045786442217e-3, 83.000328983405698, 53.115351712650373, -49.579428628294984, -.36324045786442217e-1, -.36456744565004192e-40, 8003.2734054953023, -4957.9428628294984, -3.6324045786442217, -.36456744565004192e-38, 531153.51712650373, -363.24045786442217, -.36456744565004192e-36]
  λμ6 = [.41666633333346250e-3, .41663333462496667e-3, .41633346246667344e-3, .41334621673433880e-3, .38459233284917391e-3, 4.1663333462496667, 4.1334621673433880, 3.8459233284917391, 1.8491471186490768, .40410060076602078e-3, 416.33346246667344, 384.59233284917391, 184.91471186490768, .40410060076602078e-1, .36825032086630255e-40, 41334.621673433880, 18491.471186490768, 4.0410060076602078, .36825032086630255e-38, 3845923.3284917391, 404.10060076602078, .36825032086630255e-36]
  μμ6 = [-999999999.99991667, -99999.999916666667, -999.99991666667083, -9.9999166670833317, -.99916708316804727e-1, -999999999.16666667, -99999.166670833317, -999.16708316804727, -9.2067359420779232, -.45404052350475412e-3, -999999916.66667083, -99916.708316804727, -920.67359420779232, -.45404052350475412e-1, -.37200759760208360e-40, -999991666.70833317, -92067.359420779232, -4.5404052350475412, -.37200759760208360e-38, -999167083.16804727, -454.04052350475412, -.37200759760208360e-36]

  # case i = 1000, j = 500
  λλ7 = [.41666166665514486e-4, .41616655144589631e-4, .41165514211314588e-4, .36551172143059331e-4, -.20145484770597893e-4, .41616655144589631, .36551172143059331, -.20145484770597893, -20.819207244228436, -72007626.637629603, 41.165514211314588, -20.145484770597893, -2081.9207244228436, -7200762663.7629603, -.10724035427265142e87, 3655.1172143059331, -208192.07244228436, -720076266376.29603, -.10724035427265142e89, -201454.84770597893, -72007626637629.603, -.10724035427265142e91]
  λμ7 = [.12063368253496155e-2, .12062873284900947e-2, .12058376393836500e-2, .12013686839547642e-2, .11594298763101958e-2, 12.062873284900947, 12.013686839547642, 11.594298763101958, 9.8629559705354181, 9872.5408342737259, 1205.8376393836500, 1159.4298763101958, 986.29559705354181, 987254.08342737259, .13253061369546797e45, 120136.86839547642, 98629.559705354181, 98725408.342737259, .13253061369546797e47, 11594298.763101958, 9872540834.2737259, .13253061369546797e49]
  μμ7 = [-499999999.99995833, -49999.999958333333, -499.99995833333542, -4.9999583335416658, -.49958354158402363e-1, -499999999.58333333, -49999.583335416658, -499.58354158402363, -4.6033679710389616, -.22702026175237706e-3, -499999958.33333542, -49958.354158402363, -460.33679710389616, -.22702026175237706e-1, -.18600379880104180e-40, -499995833.35416658, -46033.679710389616, -2.2702026175237706, -.18600379880104180e-38, -499583541.58402363, -227.02026175237706, -.18600379880104180e-36]

  # case i = 1000, j = 1000
  λλ8 = [-.50066550000008370e-6, -.50000550083336772e-2, -.49999958333347847, -50.000723337672237, -5008.3279637162128, -50.000550083336772, -500007.23337672237, -50083279.637162128, -5898717394.4429521, -24253782099472752., -499999.58333347847, -5008327963.7162128, -589871739444.29521, -.24253782099472752e19, -.36129760523281962e95, -5000072333.7672237, -58987173944429.521, -.24253782099472752e21, -.36129760523281962e97, -50083279637162.128, -.24253782099472752e23, -.36129760523281962e99]
  λμ8 = [.99999999933358333, .99999993583333327, .99999958333328056, 1.0000183334055557, 1.0024346558790135, 9999.9993583333327, 10000.183334055557, 10024.346558790135, 12635.053587747725, 19804195.346120856, 999999.58333328056, 1002434.6558790135, 1263505.3587747725, 1980419534.6120856, .26585747344275761e48, 100001833.34055557, 126350535.87747725, 198041953461.20856, .26585747344275761e50, 10024346558.790135, 19804195346120.856, .26585747344275761e52]
  μμ8 = zeros(Float64, 22)

  # case i = 1000, j = 1001
  λλ9 = fill(Float64(-Inf), 22)
  λμ9 = fill(Float64(NaN), 22)
  μμ9 = [999999.99999991667, 99.999999916666667, .99999991666667083, .99999166670833317e-2, .99916708316804727e-4, 999999.99916666667, 99.999166670833317, .99916708316804727, .92067359420779232e-2, .45404052350475412e-6, 999999.91666667083, 99.916708316804727, .92067359420779232, .45404052350475412e-4, .37200759760208360e-43, 999991.66670833317, 92.067359420779232, .45404052350475412e-2, .37200759760208360e-41, 999167.08316804727, .45404052350475412, .37200759760208360e-39]

  # case i = 1000, j = 1002
  λλ10 = fill(Float64(NaN), 22)
  λμ10 = [-1.0020008340009172, -1.0020009025716693, -1.0020017511670057, -1.0020325568984007, -1.0045740280231603, -10020.009025716693, -10020.325568984007, -10045.740280231603, -12675.072001219609, -19883530.670621499, -1002001.7511670057, -1004574.0280231603, -1267507.2001219609, -1988353067.0621499, -.26692250007811088e48, -100203255.68984007, -126750720.01219609, -198835306706.21499, -.26692250007811088e50, -10045740280.231603, -19883530670621.499, -.26692250007811088e52]
  μμ10 = [1999999.9999998333, 199.99999983333333, 1.9999998333333417, .19999833334166663e-1, .19983341663360945e-3, 1999999.9983333333, 199.99833334166663, 1.9983341663360945, .18413471884155846e-1, .90808104700950824e-6, 1999999.8333333417, 199.83341663360945, 1.8413471884155846, .90808104700950824e-4, .74401519520416719e-43, 1999983.3334166663, 184.13471884155846, .90808104700950824e-2, .74401519520416719e-41, 1998334.1663360945, .90808104700950824, .74401519520416719e-39]

  # case i = 1000, j = 2000
  λλ11 = [-.83334333330979118e-4, -.83433309791728800e-4, -.84330979667575035e-4, -.93098458713979056e-4, -.16031102641766044e-3, -.83433309791728800, -.93098458713979056, -1.6031102641766044, 9.7824061341796240, 919500.10049597720, -84.330979667575035, -160.31102641766044, 978.24061341796240, 91950010.049597720, .11596326709759849e85, -9309.8458713979056, 97824.061341796240, 9195001004.9597720, .11596326709759849e87, -1603110.2641766044, 919500100495.97720, .11596326709759849e89]
  λμ11 = [-.24186696686692983e-2, -.24187686749654403e-2, -.24196692984358704e-2, -.24287317357315902e-2, -.25250688354042896e-2, -24.187686749654403, -24.287317357315902, -25.250688354042896, -41.874457251134705, -79335.364991461609, -2419.6692984358704, -2525.0688354042896, -4187.4457251134705, -7933536.4991461609, -.10650271681332433e46, -242873.17357315902, -418744.57251134705, -793353649.91461609, -.10650271681332433e48, -25250688.354042896, -79335364991.461609, -.10650271681332433e50]
  μμ11 = [999999999.99991667, 99999.999916666667, 999.99991666667083, 9.9999166670833317, .99916708316804727e-1, 999999999.16666667, 99999.166670833317, 999.16708316804727, 9.2067359420779232, .45404052350475412e-3, 999999916.66667083, 99916.708316804727, 920.67359420779232, .45404052350475412e-1, .37200759760208360e-40, 999991666.70833317, 92067.359420779232, 4.5404052350475412, .37200759760208360e-38, 999167083.16804727, 454.04052350475412, .37200759760208360e-36]

  λλ = [λλ1 λλ2 λλ3 λλ4 λλ5 λλ6 λλ7 λλ8 λλ9 λλ10 λλ11]
  λμ = [λμ1 λμ2 λμ3 λμ4 λμ5 λμ6 λμ7 λμ8 λμ9 λμ10 λμ11]
  μμ = [μμ1 μμ2 μμ3 μμ4 μμ5 μμ6 μμ7 μμ8 μμ9 μμ10 μμ11]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, μ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([0.0, μ], i, j, t)
      @test y[2, 2] ≈ μμ[a, b]

      if (b == 3) || (b == 9)
        @test isinf(y[1, 1])
        @test isnan(y[1, 2])
        @test isnan(y[2, 1])
      elseif (b == 4) || (b == 10)
        @test isnan(y[1, 1])
        @test y[1, 2] ≈ λμ[a, b]
        @test y[2, 1] ≈ λμ[a, b]
      else
        @test y[1, 1] ≈ λλ[a, b]
        @test y[1, 2] ≈ λμ[a, b]
        @test y[2, 1] ≈ λμ[a, b]
      end
    end
  end
end

@testset "Equal rates" begin
  # formulas are subject to catastrophic cancellation for particular
  # combinations of parameters. We are willing to accept a lower relative
  # tolerance than the default of ϵ = sqrt(eps()) ≈ 1.49e-8
  ϵ = 1.0e-3

  u = ((1, 0), (1, 1), (100, 0), (100, 1), (100, 50), (100, 100), (100, 200))
  v = ((0.001, 0.001), (0.001, 0.1), (0.001, 1.0), (0.001, 10.0), (0.001, 100.0),
       (  0.1, 0.001), (  0.1, 0.1), (  0.1, 1.0), (  0.1, 10.0), (  0.1, 100.0),
       (  1.0, 0.001), (  1.0, 0.1), (  1.0, 1.0), (  1.0, 10.0), (  1.0, 100.0),
       ( 10.0, 0.001), ( 10.0, 0.1), ( 10.0, 1.0), ( 10.0, 10.0),
       (100.0, 0.001), (100.0, 0.1), (100.0, 1.0))

  # case i = 1, j = 0
  λλ1 = [.83333000000583333e-7, .83300005832500108e-7, .83000582501082002e-7, .80057510701565206e-7, .55096418732782369e-7, .83300005832500108e-3, .80057510701565206e-3, .55096418732782369e-3, -.20833333333333333e-3, -.13085399449035813e-3, .83000582501082002e-1, .55096418732782369e-1, -.20833333333333333e-1, -.13085399449035813e-1, -.16256576152664771e-2, 8.0057510701565206, -2.0833333333333333, -1.3085399449035813, -.16256576152664771, 550.96418732782369, -130.85399449035813, -16.256576152664771]
  λμ1 = [.41666600000091667e-6, .41660000916550014e-6, .41600091550141500e-6, .41009051400189524e-6, .35812672176308540e-6, .41660000916550014e-2, .41009051400189524e-2, .35812672176308540e-2, .14583333333333333e-2, .17217630853994490e-3, .41600091550141500, .35812672176308540, .14583333333333333, .17217630853994490e-1, .16746724177368232e-2, 41.009051400189524, 14.583333333333333, 1.7217630853994490, .16746724177368232, 3581.2672176308540, 172.17630853994490, 16.746724177368232]
  μμ1 = [-999999.99999991667, -99.999999916699994, -.99999991699941750, -.99999199424892984e-2, -.99944903581267218e-4, -999999.99916699994, -99.999199424892984, -.99944903581267218, -.10208333333333333e-1, -.23085399449035813e-3, -999999.91699941750, -99.944903581267218, -1.0208333333333333, -.23085399449035813e-1, -.17256576152664771e-2, -999991.99424892984, -102.08333333333333, -2.3085399449035813, -.17256576152664771, -999449.03581267218, -230.85399449035813, -17.256576152664771]

  # case i = 1, j = 1
  λλ2 = [-.66666550000166666e-12, -.66655001666450027e-10, -.66550166450266350e-9, -.65516452635362546e-8, -.56473829201101928e-7, -.66655001666450027e-6, -.65516452635362546e-4, -.56473829201101928e-3, -.20833333333333333e-2, -.19283746556473829e-2, -.66550166450266350e-3, -.56473829201101928e-1, -.20833333333333333, -.19283746556473829, -.16991798189719962, -.65516452635362546, -20.833333333333333, -19.283746556473829, -16.991798189719962, -564.73829201101928, -1928.3746556473829, -1699.1798189719962]
  λμ2 = [.99999866666850000e-6, .99986668499766695e-6, .99866849766949667e-6, .98684769467045714e-6, .88292011019283747e-6, .99986668499766695e-2, .98684769467045714e-2, .88292011019283747e-2, .45833333333333333e-2, .20110192837465565e-2, .99866849766949667, .88292011019283747, .45833333333333333, .20110192837465565, .17001601150214031, 98.684769467045714, 45.833333333333333, 20.110192837465565, 17.001601150214031, 8829.2011019283747, 2011.0192837465565, 1700.1601150214031]
  μμ2 = [-.66666550000166666e-12, -.66655001666450027e-10, -.66550166450266350e-9, -.65516452635362546e-8, -.56473829201101928e-7, -.66655001666450027e-6, -.65516452635362546e-4, -.56473829201101928e-3, -.20833333333333333e-2, -.19283746556473829e-2, -.66550166450266350e-3, -.56473829201101928e-1, -.20833333333333333, -.19283746556473829, -.16991798189719962, -.65516452635362546, -20.833333333333333, -19.283746556473829, -16.991798189719962, -564.73829201101928, -1928.3746556473829, -1699.1798189719962]

  # case i = 100, j = 0
  λλ3 = [.83333000000583333e-5, .83300005832500108e-5, .83000582501082002e-5, .80057510701565206e-5, .55096418732782369e-5, .83300005832500108e-1, .80057510701565206e-1, .55096418732782369e-1, -.20833333333333333e-1, -.13085399449035813e-1, 8.3000582501082002, 5.5096418732782369, -2.0833333333333333, -1.3085399449035813, -.16256576152664771, 800.57510701565206, -208.33333333333333, -130.85399449035813, -16.256576152664771, 55096.418732782369, -13085.399449035813, -1625.6576152664771]
  λμ3 = [.41666600000091667e-4, .41660000916550014e-4, .41600091550141500e-4, .41009051400189524e-4, .35812672176308540e-4, .41660000916550014, .41009051400189524, .35812672176308540, .14583333333333333, .17217630853994490e-1, 41.600091550141500, 35.812672176308540, 14.583333333333333, 1.7217630853994490, .16746724177368232, 4100.9051400189524, 1458.3333333333333, 172.17630853994490, 16.746724177368232, 358126.72176308540, 17217.630853994490, 1674.6724177368232]
  μμ3 = [-99999999.999991667, -9999.9999916699994, -99.999991699941750, -.99999199424892984, -.99944903581267218e-2, -99999999.916699994, -9999.9199424892984, -99.944903581267218, -1.0208333333333333, -.23085399449035813e-1, -99999991.699941750, -9994.4903581267218, -102.08333333333333, -2.3085399449035813, -.17256576152664771, -99999199.424892984, -10208.333333333333, -230.85399449035813, -17.256576152664771, -99944903.581267218, -23085.399449035813, -1725.6576152664771]

  # case i = 100, j = 1
  λλ4 = [.82499663333922499e-5, .82466339224158443e-5, .82163921659426155e-5, .79191419141914191e-5, .53980716253443526e-5, .82466339224158443e-1, .79191419141914191e-1, .53980716253443526e-1, -.22708333333333333e-1, -.14882920110192837e-1, 8.2163921659426155, 5.3980716253443526, -2.2708333333333333, -1.4882920110192837, -.33085808580858086, 791.91419141914191, -227.08333333333333, -148.82920110192837, -33.085808580858086, 53980.716253443526, -14882.920110192837, -3308.5808580858086]
  λμ4 = [.42249932666759250e-4, .42243267592382181e-4, .42182759132309582e-4, .41585808580858086e-4, .36337465564738292e-4, .42243267592382181, .41585808580858086, .36337465564738292, .14895833333333333, .19056473829201102e-1, 42.182759132309582, 36.337465564738292, 14.895833333333333, 1.9056473829201102, .33580858085808581, 4158.5808580858086, 1489.5833333333333, 190.56473829201102, 33.580858085808581, 363374.65564738292, 19056.473829201102, 3358.0858085808581]
  μμ4 = [-98999999.999991750, -9899.9999917533661, -98.999991783607834, -.98999208085808581, -.98946019283746556e-2, -98999999.917533661, -9899.9208085808581, -98.946019283746556, -1.0127083333333333, -.24782920110192837e-1, -98999991.783607834, -9894.6019283746556, -101.27083333333333, -2.4782920110192837, -.34075808580858086, -98999208.085808581, -10127.083333333333, -247.82920110192837, -34.075808580858086, -98946019.283746556, -24782.920110192837, -3407.5808580858086]

  # case i = 100, j = 50
  λλ5 = [.41666166666272575e-5, .41616662724514970e-5, .41166271338474794e-5, .36626128375626596e-5, -.12400880426595393e-5, .41616662724514970e-1, .36626128375626596e-1, -.12400880426595393e-1, -.17778913406221861, -.34764084409279366e-1, 4.1166271338474794, -1.2400880426595393, -17.778913406221861, -3.4764084409279366, -.45271312799763001, 366.26128375626596, -1777.8913406221861, -347.64084409279366, -45.271312799763001, -12400.880426595393, -34764.084409279366, -4527.1312799763001]
  λμ5 = [.11887244901946751e-3, .11886254761605637e-3, .11877240850878170e-3, .11785835861931029e-3, .10746652249990810e-3, 1.1886254761605637, 1.1785835861931029, 1.0746652249990810, .34849721940534291, .37203150470509560e-1, 118.77240850878170, 107.46652249990810, 34.849721940534291, 3.7203150470509560, .45522926644795513, 11785.835861931029, 3484.9721940534291, 372.03150470509560, 45.522926644795513, 1074665.2249990810, 37203.150470509560, 4552.2926644795513]
  μμ5 = [-49999999.999995833, -4999.9999958383337, -49.999995883372866, -.49999633738716244, -.50012400880426595e-2, -49999999.958383337, -4999.9633738716244, -50.012400880426595, -.67778913406221861, -.39764084409279366e-1, -49999995.883372866, -5001.2400880426595, -67.778913406221861, -3.9764084409279366, -.45771312799763001, -49999633.738716244, -6777.8913406221861, -397.64084409279366, -45.771312799763001, -50012400.880426595, -39764.084409279366, -4577.1312799763001]

  # case i = 100, j = 100
  λλ6 = [-.11664999933374163e-9, -.50643336453393712e-6, -.49390384192982280e-4, -.18466221277045134e-2, -.47453404521680332e-3, -.50643336453393712e-2, -18.466221277045134, -4.7453404521680332, -.49790935414088870, -.50407709465625055e-1, -49.390384192982280, -474.53404521680332, -49.790935414088870, -5.0407709465625055, -.57215735744432853, -184662.21277045134, -4979.0935414088870, -504.07709465625055, -57.215735744432853, -4745340.4521680332, -50407.709465625055, -5721.5735744432853]
  λμ6 = [.99999997666950010e-2, .99989870499564435e-2, .99008855740672263e-2, .51311135979071411e-2, .50021154820563956e-3, 99.989870499564435, 51.311135979071411, 5.0021154820563956, .50041560398306440, .50432661763749555e-1, 9900.8855740672263, 500.21154820563956, 50.041560398306440, 5.0432661763749555, .57217652502882766, 513111.35979071411, 5004.1560398306440, 504.32661763749555, 57.217652502882766, 5002115.4820563956, 50432.661763749555, 5721.7652502882766]
  μμ6 = [-.11664999933374163e-9, -.50643336453393712e-6, -.49390384192982280e-4, -.18466221277045134e-2, -.47453404521680332e-3, -.50643336453393712e-2, -18.466221277045134, -4.7453404521680332, -.49790935414088870, -.50407709465625055e-1, -49.390384192982280, -474.53404521680332, -49.790935414088870, -5.0407709465625055, -.57215735744432853, -184662.21277045134, -4979.0935414088870, -504.07709465625055, -57.215735744432853, -4745340.4521680332, -50407.709465625055, -5721.5735744432853]

  # case i = 100, j = 200
  λλ7 = [-99999999.999991667, -9999.9999916766676, -99.999991766760046, -.99999267620357329, -.10002603980125785e-1, -99999999.916766676, -9999.9267620357329, -100.02603980125785, -1.3571893101381411, -.79101645877138889e-1, -99999991.766760046, -10002.603980125785, -135.71893101381411, -7.9101645877138889, -.80092426825143571, -99999267.620357329, -13571.893101381411, -791.01645877138889, -80.092426825143571, -100026039.80125785, -79101.645877138889, -8009.2426825143571]
  λμ7 = [.23968626864655441e-3, .23966646554377574e-3, .23948615802693786e-3, .23765513126851455e-3, .21660834710927332e-3, 2.3966646554377574, 2.3765513126851455, 2.1660834710927332, .69687980914836946, .73955062212881197e-1, 239.48615802693786, 216.60834710927332, 69.687980914836946, 7.3955062212881197, .79592638194259420, 23765.513126851455, 6968.7980914836946, 739.55062212881197, 79.592638194259420, 2166083.4710927332, 73955.062212881197, 7959.2638194259420]
  μμ7 = [.83332333332402011e-5, .83233324017634264e-5, .82332399539510388e-5, .73237964267106197e-5, -.26039801257846584e-5, .83233324017634264e-1, .73237964267106197e-1, -.26039801257846584e-1, -.35718931013814113, -.69101645877138889e-1, 8.2332399539510388, -2.6039801257846584, -35.718931013814113, -6.9101645877138889, -.79092426825143571, 732.37964267106197, -3571.8931013814113, -691.01645877138889, -79.092426825143571, -26039.801257846584, -69101.645877138889, -7909.2426825143571]

  λλ = [λλ1 λλ2 λλ3 λλ4 λλ5 λλ6 λλ7]
  λμ = [λμ1 λμ2 λμ3 λμ4 λμ5 λμ6 λμ7]
  μμ = [μμ1 μμ2 μμ3 μμ4 μμ5 μμ6 μμ7]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, λ], i, j, t)
      @test isapprox(y[1, 1], λλ[a, b], rtol=ϵ)
      @test isapprox(y[2, 1], λμ[a, b], rtol=ϵ)
      @test isapprox(y[1, 2], λμ[a, b], rtol=ϵ)
      @test isapprox(y[2, 2], μμ[a, b], rtol=ϵ)
    end
  end
end

@testset "Unequal rates" begin
  # TODO: Improve numerical accuracy of Hessian formulas

  # Hessian formulas are subject to catastrophic cancellation for particular
  # combinations of parameters. We are willing to accept a slightly lower
  # relative tolerance than the default of ϵ = sqrt(eps()) ≈ 1.49e-8
  ϵ = 1.0e-5

  u = ((1, 0), (1, 1), (100, 0), (100, 1), (100, 2), (100, 50), (100, 100), (100, 200))

  # θ = (λ - μ) * t
  # ω = log(λ / μ)
  # compute the gradient at (approximate) points (θ, ω) at different 't':
  # (-20, -10), (-10,  -7), (-5, -2), (-1, -5), (-0.5, -0.1), (-0.0001, -1)
  # (0.001, 1), (0.5, 0.1), ( 1,  5), ( 5,  2), (  10,    7), (     20, 10)
  v = ((0.001,      0.5, 20_001.0), (0.001,    10.0, 10_010.0), (0.001,   800.0, 5_800.0), (0.001,     7.0, 1_007.0), (0.001,  4_700.0, 5_200.0), (0.001,     0.15,     0.25),
       (0.001,     0.25,     0.15), (0.001, 5_200.0,  4_700.0), (0.001, 1_007.0,     7.0), (0.001, 5_800.0,   800.0), (0.001, 10_010.0,    10.0), (0.001, 20_001.0,      0.5),
       (  1.0,   0.0005,   20.005), (  1.0,   0.005,   10.005), (  1.0,     0.8,     5.8), (  1.0,   0.007,   1.007), (  1.0,      4.7,     5.2), (  1.0, 0.000154, 0.000254),
       (  1.0, 0.000254, 0.000154), (  1.0,     5.2,      4.7), (  1.0,   1.007,   0.007), (  1.0,     5.8,     0.8), (  1.0,   10.005,   0.005), (  1.0,   20.005,   0.0005),
       (100.0,   9.0e-6, 0.200009), (100.0,  9.0e-5,  0.10009), (100.0,  0.0078,  0.0578), (100.0,  7.0e-5, 0.01007), (100.0,   0.0475,  0.0525), (100.0, 1.542e-6, 2.542e-6),
       (100.0, 2.542e-6, 1.542e-6), (100.0,  0.0525,   0.0475), (100.0, 0.01007,  0.7e-4), (100.0,  0.0578,  0.0078), (100.0,  0.10009,  0.9e-4), (100.0, 0.200009,   9.0e-6))

  # case i = 1, j = 0
  λλ1 = [-.18540697809772868e-14, -.36287761340522706e-10, -.35694642215581045e-8, -.49562055156282736e-7, -.22340285563093801e-7, .83250030158085537e-7, .83283350162017842e-7, -.17703794411927834e-7, .78210530081870600e-7, .24204345766267671e-7, .99346812294128723e-8, .24997479586823461e-8, -.18467093505545231e-8, -.36305894497732242e-4, -.35694642215581045e-2, -.49562055156282736e-1, -.22340285563093801e-1, .83248697900264523e-1, .83282017637691697e-1, -.17703794411927834e-1, .78210530081870600e-1, .24204345766267671e-1, .99446311200498346e-2, .24987484167511895e-2, -.18549547957153271e-4, -.36291386522748753, -35.816673053497619, -495.62055156282736, -221.94997463167606, 832.48631287863005, 832.81951011965006, -176.35195117846158, 782.10530081870600, 243.83556320691192, 99.366700140161320, 24.997729541354437]
  λμ1 = [.19519184279360775e-14, .40365616213720209e-10, .43734903460821636e-8, .18370756603107506e-6, .34047356153559954e-7, .41653336899111878e-6, .41653336899111878e-6, .34047356153559954e-7, .18370756603107506e-6, .43734903460821636e-8, .40365616213720209e-10, .19519184279360775e-14, .19441489238350697e-8, .40387826021297718e-4, .43734903460821636e-2, .18370756603107506, .34047356153559954e-1, .41653070380522335, .41653070380522335, .34047356153559954e-1, .18370756603107506, .43734903460821636e-2, .40387826021297718e-4, .19441489238350697e-8, .19528507819547967e-4, .40370056237308282, 43.912346777538440, 1837.0756603107506, 337.08429969520399, 4165.3057054669798, 4165.3057054669798, 337.08429969520399, 1837.0756603107506, 43.912346777538440, .40370056237308282, .19528507819547967e-4]
  μμ1 = [-.20600664039392114e-14, -.45348730637067794e-10, -.55221702860509977e-8, -.90793510982310441e-6, -.54686042932637893e-7, -15.999999916716650, -44.444444361194414, -.67609638211350931e-7, -.20408212827361279e-1, -.15660694642215581e-5, -.10000000036287761e-1, -4.0000000000000019, -.20518426093300468e-8, -.45376374953288565e-4, -.55221702860509977e-2, -.90793510982310441, -.54686042932637893e-1, -15500030.916779982, -42165626.497962299, -.67609638211350931e-1, -20408.212827361279, -1.5660694642215581, -40000.000036305894, -4000000.0000000018, -.20610511451423516e-4, -.45354256856637602, -55.490354825678691, -9079.3510982310441, -539.16374256168153, -154756502190.19402, -420563183739.23737, -665.16327103056803, -204082128.27361279, -16472.371571146857, -123456790.48637066, -12345679012.345698]

  # case i = 1, j = 1
  λλ2 = [-.49997495972676260e-8, -.19981767417980095e-7, -.73479582682732096e-7, -.25777692189698083e-6, -.20928439306066007e-6, -.16660626716225975e-9, -.99966259297650355e-10, -.20001141075832814e-6, -.22317514206741620e-8, -.17931962707080546e-7, -.39829436473304624e-10, -.24997176337179830e-12, -.49977503487673977e-2, -.19981803684294514e-1, -.73479582682732096e-1, -.25777692189698083, -.20928439306066007, -.16927078280428639e-3, -.10263130794993936e-3, -.20001141075832814, -.22317514206741620e-2, -.17931962707080546e-1, -.19929655199380054e-4, -.24982184631764809e-6, -49.999995876023296, -199.81774668344547, -735.03988850315411, -2577.7692189698083, -2089.9381686080768, -1.6940400760732269, -1.0276455940332091, -1998.7421217016479, -22.317514206741620, -175.73541598233504, -.35851867266785165, -.44996942185081166e-2]
  λμ2 = [.49997497929649199e-8, .19989923127726490e-7, .75087634931780214e-7, .52606794364656549e-6, .23269853424159238e-6, .99973340456557089e-6, .99973340456557089e-6, .23269853424159238e-6, .52606794364656549e-6, .75087634931780214e-7, .19989923127726490e-7, .49997497929649199e-8, .49977505436465443e-2, .19989967547341645e-1, .75087634931780214e-1, .52606794364656549, .23269853424159238, .99972807419378003, .99972807419378003, .23269853424159238, .52606794364656549, .75087634931780214e-1, .19989967547341645e-1, .49977505436465443e-2, 49.999997833943020, 199.89932007773666, 751.23123595123575, 5260.6794364656549, 2320.2068187351327, 9997.2780767672930, 9997.2780767672930, 2320.2068187351327, 5260.6794364656549, 751.23123595123575, 199.89932007773666, 49.999997833943020]
  μμ2 = [-.24997176337179830e-12, -.39829436473304624e-10, -.17931962707080546e-7, -.22317514206741620e-8, -.20001141075832814e-6, -.99966259297650355e-10, -.16660626716225975e-9, -.20928439306066007e-6, -.25777692189698083e-6, -.73479582682732096e-7, -.19981767417980095e-7, -.49997495972676260e-8, -.24982184631764809e-6, -.19929655199380054e-4, -.17931962707080546e-1, -.22317514206741620e-2, -.20001141075832814, -.10263130794993936e-3, -.16927078280428639e-3, -.20928439306066007, -.25777692189698083, -.73479582682732096e-1, -.19981803684294514e-1, -.49977503487673977e-2, -.44996942185081166e-2, -.35851867266785165, -175.73541598233504, -22.317514206741620, -1998.7421217016479, -1.0276455940332091, -1.6940400760732269, -2089.9381686080768, -2577.7692189698083, -735.03988850315411, -199.81774668344547, -49.999995876023296]

  # case i = 100, j = 0
  λλ3 = [-.18540697809772868e-12, -.36287761340522706e-8, -.35694642215581045e-6, -.49562055156282736e-5, -.22340285563093801e-5, .83250030158085537e-5, .83283350162017842e-5, -.17703794411927834e-5, .78210530081870600e-5, .24204345766267671e-5, .99346812294128723e-6, .24997479586823461e-6, -.18467093505545231e-6, -.36305894497732242e-2, -.35694642215581045, -4.9562055156282736, -2.2340285563093801, 8.3248697900264523, 8.3282017637691697, -1.7703794411927834, 7.8210530081870600, 2.4204345766267671, .99446311200498346, .24987484167511895, -.18549547957153271e-2, -36.291386522748753, -3581.6673053497619, -49562.055156282736, -22194.997463167606, 83248.631287863005, 83281.951011965006, -17635.195117846158, 78210.530081870600, 24383.556320691192, 9936.6700140161320, 2499.7729541354437]
  λμ3 = [.19519184279360775e-12, .40365616213720209e-8, .43734903460821636e-6, .18370756603107506e-4, .34047356153559954e-5, .41653336899111878e-4, .41653336899111878e-4, .34047356153559954e-5, .18370756603107506e-4, .43734903460821636e-6, .40365616213720209e-8, .19519184279360775e-12, .19441489238350697e-6, .40387826021297718e-2, .43734903460821636, 18.370756603107506, 3.4047356153559954, 41.653070380522335, 41.653070380522335, 3.4047356153559954, 18.370756603107506, .43734903460821636, .40387826021297718e-2, .19441489238350697e-6, .19528507819547967e-2, 40.370056237308282, 4391.2346777538440, 183707.56603107506, 33708.429969520399, 416530.57054669798, 416530.57054669798, 33708.429969520399, 183707.56603107506, 4391.2346777538440, 40.370056237308282, .19528507819547967e-2]
  μμ3 = [-.20600664039392114e-12, -.45348730637067794e-8, -.55221702860509977e-6, -.90793510982310441e-4, -.54686042932637893e-5, -1599.9999916716650, -4444.4444361194414, -.67609638211350931e-5, -2.0408212827361279, -.15660694642215581e-3, -1.0000000036287761, -400.00000000000019, -.20518426093300468e-6, -.45376374953288565e-2, -.55221702860509977, -90.793510982310441, -5.4686042932637893, -1550003091.6779982, -4216562649.7962299, -6.7609638211350931, -2040821.2827361279, -156.60694642215581, -4000000.0036305894, -400000000.00000018, -.20610511451423516e-2, -45.354256856637602, -5549.0354825678691, -907935.10982310441, -53916.374256168153, -15475650219019.402, -42056318373923.737, -66516.327103056803, -20408212827.361279, -1647237.1571146857, -12345679048.637066, -1234567901234.5698]

  # case i = 100, j = 1
  λλ4 = [-.49999331501759428e-8, -.23574255790691843e-7, -.42685654061698444e-6, -.51644203823689717e-5, -.24209726638069463e-5, .82415863793833059e-5, .82449516997804687e-5, -.19526870575391837e-5, .77406107266845152e-5, .23782982681534188e-5, .98349361227540105e-6, .24747479793778890e-6, -.49979331729931026e-2, -.23576087239570006e-1, -.42685654061698444, -5.1644203823689717, -2.4209726638069463, 8.2414518213433835, 8.2448171148235281, -1.9526870575391837, 7.7406107266845152, 2.3782982681534188, .98449855122973424, .24737584343652144, -50.001832281271054, -235.74621934096673, -4280.8905207994184, -51644.203823689717, -24062.985657144007, 82414.450934908302, 82448.103856251323, -19457.585288369344, 77406.107266845152, 23963.985341501945, 9836.9447952033028, 2474.7707248998707]
  λμ4 = [.49999430328892856e-8, .23986119132884790e-7, .50806317919391441e-6, .18713116980722997e-4, .36033867934440278e-5, .42236536934686330e-4, .42236536934686330e-4, .36033867934440278e-5, .18713116980722997e-4, .50806317919391441e-6, .23986119132884790e-7, .49999430328892856e-8, .49979430143900040e-2, .23988362323450119e-1, .50806317919391441, 18.713116980722997, 3.6033867934440278, 42.236267750910892, 42.236267750910892, 3.6033867934440278, 18.713116980722997, .50806317919391441, .23988362323450119e-1, .49979430143900040e-2, 50.001931156217156, 239.86567575267186, 5098.5535669275413, 187131.16980722997, 35691.552488560328, 422362.54291799830, 422362.54291799830, 35691.552488560328, 187131.16980722997, 5098.5535669275413, 239.86567575267186, 50.001931156217156]
  μμ4 = [-.45391833736178023e-12, -.45293537695430162e-8, -.56462682102612932e-6, -.89887807623908010e-4, -.56139296610894795e-5, -1583.9999917550483, -4399.9999917584136, -.69026385759844022e-5, -2.0204133276856885, -.15511435654061698e-3, -.99000002357425579, -396.00000000499993, -.45295426464132272e-6, -.45121907755749480e-2, -.56462682102612932, -89.887807623908010, -5.6139296610894795, -1534503060.7613209, -4174397023.2984369, -6.9026385759844022, -2020413.3276856885, -155.11435654061698, -3960000.0235760872, -396000000.00499793, -.65401348521990447e-2, -45.259232960739078, -5669.2805437245254, -898878.07623908010, -55375.952635308120, -15320893716830.236, -41635755190186.194, -67941.102000634312, -20204133276.856885, -1631499.8254320420, -12222222457.968442, -1222222222272.2241]

  # case i = 100, j = 2
  λλ5 = [-3.9675807968071457, -.98353120353130817e-3, -.93436160596652992e-6, -.53803379208370003e-5, -.25641683453558282e-5, .81581697428904629e-5, .81615683835772432e-5, -.20753341824275296e-5, .76602474641244655e-5, .24904882015469049e-5, .11492595666259066e-5, .25111612718525662e-6, -3967703.0111576486, -1373.4068402261176, -.93436160596652992, -5.3803379208370003, -2.5641683453558282, 8.1580338525946443, 8.1614324661036148, -2.0753341824275296, 7.6602474641244655, 2.4904882015469049, 1.0969702424230145, .25100355437830446, -12289934549.448664, -10469328.673405890, -9453.3223546715762, -53803.379208370003, -25485.939233486575, 81580.270581297882, 81614.256702798768, -20676.383767596238, 76602.474641244655, 25074.631735591372, 11420.112145144349, 2495.0827706767988]
  λμ5 = [.76889029587052205e-5, .19071245160493535e-4, .66971956611129955e-6, .19080980991878433e-4, .37548311718115354e-5, .42839938990356436e-4, .42839938990356436e-4, .37548311718115354e-5, .19080980991878433e-4, .66971956611129955e-6, .19071245160493535e-4, .76889029587052205e-5, 7.6600737564213785, 26.738615177135929, .66971956611129955, 19.080980991878433, 3.7548311718115354, 42.839667141388342, 42.839667141388342, 3.7548311718115354, 19.080980991878433, .66971956611129955, 26.738615177135929, 7.6600737564213785, 23853.824847955337, 203169.20945850300, 6767.9349043015786, 190809.80991878433, 37193.808024389854, 428396.53549018413, 428396.53549018413, 37193.808024389854, 190809.80991878433, 6767.9349043015786, 203169.20945850300, 23853.824847955337]
  μμ5 = [.61406253478791166e-8, .17121663054101246e-6, -.42271037158032461e-6, -.88982025246563085e-4, -.56995945374571154e-5, -1567.9999918384316, -4355.5555473973858, -.70005649048850270e-5, -2.0000053803379208, -.15405936160596653e-3, -.98098353120353131, -395.96758079680715, .61260084561121748e-2, .11794950791270846, -.42271037158032461, -88.982025246563085, -5.6995945374571154, -1519003029.8446435, -4132231396.8006438, -7.0005649048850270, -2000005.3803379208, -154.05936160596653, -3921373.4068402261, -395967703.01115765, 45.303255793941769, 1637.7283596890351, -4259.3082316025070, -889820.25246563085, -56231.939323151794, -15166137214641.069, -41215192006448.651, -68920.842280577988, -20000053803.379208, -1620235.7023678208, -12109234760.772171, -1222166477759.3252]

  # case i = 100, j = 50
  λλ6 = [-194.55066906840326, -.32599119052765900, -.59562897793560737e-4, -.18391819420165313e-3, -.69550320976721772e-5, .41541626622441978e-5, .41591666999777895e-5, -.56621815351973071e-5, .39730457963362265e-5, .45952933846232730e-5, .48914913939068294e-5, .40533252781548184e-6, -194555671.02751812, -1107799.7010150164, -59.562897793560737, -183.91819420165313, -6.9550320976721772, 4.1539625581147943, 4.1589666766154148, -5.6621815351973071, 3.9730457963362265, 4.5952933846232730, 5.7117999855054703, .40471034145320188, -602347475865.08869, -3938116283.8473973, -622628.12490853171, -1839181.9420165313, -68888.336164667527, 41539.525528751507, 41589.566754141351, -56199.674371113878, 39730.457963362265, 46407.016944665743, 50179.519986744363, 3402.1876864777841]
  λμ6 = [.33317274458778183e-3, .44385064744836505e-3, .37450539963530036e-5, .93657072446996322e-4, .72379155404313330e-5, .11885254406529681e-3, .11885254406529681e-3, .72379155404313330e-5, .93657072446996322e-4, .37450539963530036e-5, .44385064744836505e-3, .33317274458778183e-3, 332.06430172421837, 1080.7328668039804, 3.7450539963530036, 93.657072446996322, 7.2379155404313330, 118.85214383842703, 118.85214383842703, 7.2379155404313330, 93.657072446996322, 3.7450539963530036, 1080.7328668039804, 332.06430172421837, 1087004.0275664352, 5096622.6033692642, 38985.506197141825, 936570.72446996322, 71661.415336059599, 1188521.2382696557, 1188521.2382696557, 71661.415336059599, 936570.72446996322, 38985.506197141825, 5096622.6033692642, 1087004.0275664352]
  μμ6 = [.28034502687804434e-6, .43924898959043323e-5, .31089675820073395e-5, -.45334236198912524e-4, -.75112939612328100e-5, -799.99999584083330, -2222.2222180680596, -.92184997300850337e-5, -1.0205920814595078, -.13768789779356074e-3, -.82599119052765900, -394.55066906840326, .27977281802351194, 5.2122996107553141, 3.1089675820073395, -45.334236198912524, -7.5112939612328100, -775001545.84413333, -2108281324.9065873, -9.2184997300850337, -1020592.0814595078, -137.68789779356074, -3107799.7010150164, -394555671.02751812, 2152.3001788844897, 45188.507851307978, 31440.721043036213, -453342.36198912524, -74340.263940274876, -7737825109561.1097, -21028159187046.659, -91049.000984612125, -10205920814.595078, -1444455.8698131997, -10110955790.020237, -1219631426482.3726]

  # case i = 100, j = 100
  λλ7 = [-393.30301526995897, -.75067304123836232, -.12976722601415863e-3, -.83684474263298963e-1, -.10898868297846672e-4, -.31391886379568040e-5, -.11339879474519501e-5, -.88737425000684134e-5, -.17163597983169786e-5, .54393570309882585e-5, .66725812303610013e-5, .53330165133917883e-6, -393311664.04246537, -2710860.9816446026, -129.76722601415863, -83684.474263298963, -10.898868297846672, -3.2401039151135311, -1.1949764514339901, -8.8737425000684134, -1.7163597983169786, 5.4393570309882585, 7.9628330322091234, .53234782806720267, -1217245446832.7660, -9140030225.4016853, -1359386.6056766866, -836844742.63298963, -107790.43672296675, -32451.914930594876, -11980.677587265358, -87943.539671499445, -17163.597983169786, 54980.474902023189, 68669.189150430979, 4226.4315192583257]
  λμ7 = [.60335441556889053e-3, .64632743105857367e-3, .53357963207583444e-5, .66244180694482306e-3, .98451257920285936e-5, .99962258261056460e-2, .99962258261056460e-2, .98451257920285936e-5, .66244180694482306e-3, .53357963207583444e-5, .64632743105857367e-3, .60335441556889053e-3, 601.54206498387442, 1595.1226482633815, 5.3357963207583444, 662.44180694482306, 9.8451257920285936, 9996.0638622623216, 9996.0638622623216, 9.8451257920285936, 662.44180694482306, 5.3357963207583444, 1595.1226482633815, 601.54206498387442, 2053234.1133155038, 7434229.0224092325, 55597.472644074633, 6624418.0694482306, 97468.451955022586, 99960556.802966576, 99960556.802966576, 97468.451955022586, 6624418.0694482306, 55597.472644074633, 7434229.0224092325, 2053234.1133155038]
  μμ7 = [.53330165133917883e-6, .66725812303610013e-5, .54393570309882585e-5, -.17163597983169786e-5, -.88737425000684134e-5, -.11339879474519501e-5, -.31391886379568040e-5, -.10898868297846672e-4, -.83684474263298963e-1, -.12976722601415863e-3, -.75067304123836232, -393.30301526995897, .53234782806720267, 7.9628330322091234, 5.4393570309882585, -1.7163597983169786, -8.8737425000684134, -1.1949764514339901, -3.2401039151135311, -10.898868297846672, -83684.474263298963, -129.76722601415863, -2710860.9816446026, -393311664.04246537, 4226.4315192583257, 68669.189150430979, 54980.474902023189, -17163.597983169786, -87943.539671499445, -11980.677587265358, -32451.914930594876, -107790.43672296675, -836844742.63298963, -1359386.6056766866, -9140030225.4016853, -1217245446832.7660]

  # case i = 100, j = 200
  λλ8 = [-791.26350380522894, -1.6545039909253089, -.27576180173812562e-3, -2.0412004584074806, -.18405950401239727e-4, -4444.4444361361200, -1599.9999916816669, -.14989002171344369e-4, -.90662317822383813e-4, .62316153778040884e-5, .87927240289525069e-5, .72800789309328876e-6, -791277158.13444722, -6225521.3405212886, -275.76180173812562, -2041200.4584074806, -18.405950401239727, -4216562649.8131755, -1550003091.6882670, -14.989002171344369, -90.662317822383813, 6.2316153778040884, 10.432127529655573, .72667274027502166, -2447641185301.2709, -20252974646.459562, -2892977.0438997491, -20412004584.074806, -181784.92565875428, -42056318374093.327, -15475650219122.223, -148342.72153978802, -906623.17822383813, 63017.354597360837, 90453.344152094253, 5600.9319099144138]
  λμ8 = [.10146238694520783e-2, .88826688763116462e-3, .74762452185149416e-5, .18950502993135776e-3, .14433391907803636e-4, .23964645763314889e-3, .23964645763314889e-3, .14433391907803636e-4, .18950502993135776e-3, .74762452185149416e-5, .88826688763116462e-3, .10146238694520783e-2, 1011.9998102448781, 2162.5688488720194, 7.4762452185149416, 189.50502993135776, 14.433391907803636, 239.64565713157803, 239.64565713157803, 14.433391907803636, 189.50502993135776, 7.4762452185149416, 2162.5688488720194, 1011.9998102448781, 3664966.7648956677, 10199412.682756421, 77833.944005514866, 1895050.2993135776, 142898.17635367199, 2396456.1710623868, 2396456.1710623868, 142898.17635367199, 1895050.2993135776, 77833.944005514866, 10199412.682756421, 3664966.7648956677]
  μμ8 = [.97798289496816377e-6, .97907270249575009e-5, .92042669830359552e-5, .79522461681136880e-5, -.11290777319273363e-4, .83183330924519572e-5, .83083244055993453e-5, -.13879015136414014e-4, -.38413187686838619e-3, -.11951180173812562e-3, -.65450399092530894, -391.26350380522894, .97654778713440154, 11.431128279155885, 9.2042669830359552, 7.9522461681136880, -11.290777319273363, 8.3179330287097254, 8.3079241680953625, -13.879015136414014, -384.13187686838619, -119.51180173812562, -2225521.3405212886, -391277158.13444722, 8100.7069251010026, 100435.36842296702, 92949.946400619898, 79522.461681136880, -112061.54240146603, 83179.130254442682, 83079.041561417928, -137463.59601886508, -3841318.7686838619, -1249321.5540904131, -7907295634.1138830, -1213073284066.7030]

  λλ = [λλ1 λλ2 λλ3 λλ4 λλ5 λλ6 λλ7 λλ8]
  λμ = [λμ1 λμ2 λμ3 λμ4 λμ5 λμ6 λμ7 λμ8]
  μμ = [μμ1 μμ2 μμ3 μμ4 μμ5 μμ6 μμ7 μμ8]

  @testset for (b, (i, j)) = enumerate(u)
    @testset for (a, (t, λ, μ)) = enumerate(v)
      x, y = SimpleBirthDeathProcess.gradient_hessian([λ, μ], i, j, t)
      @test isapprox(y[1, 1], λλ[a, b], rtol=ϵ)
      @test isapprox(y[2, 1], λμ[a, b], rtol=ϵ)
      @test isapprox(y[1, 2], λμ[a, b], rtol=ϵ)
      @test isapprox(y[2, 2], μμ[a, b], rtol=ϵ)
    end
  end
end


text:
@testset "Outer product" begin
    @test A₃ ∧ 5 == 5 ∧ A₃ == 5A₃
    @test ∧(v1, v3, v2) == -1(v1 ∧ v2 ∧ v3) == -1v123

    for r ∈ 0:4
        Aᵣ = As[r+1]
        @test a ∧ Aᵣ ≈ 0.5 * (a * Aᵣ + (-1) ^ r * Aᵣ * a)
    end

    @test A ∧ (B + C) ≈ A ∧ B + A ∧ C
    @test (B + C) ∧ A ≈ B ∧ A + C ∧ A
end