text:
# To be included in MERAKit.jl.

"""
    ModifiedBinaryOp{T}

The operator type for `ModifiedBinaryLayer`.

Due to the structure of a `ModifiedBinaryMERA`, there's an alternating pattern of two-site
translation invariance at each layer. Because of this, every operator at a layer is defined
by two tensors, one for each possible position. `ModifiedBinaryOp` encapsulates this. The
two positions are called `mid` and `gap`, where `mid` refers to the position that has a
disentangler directly below it, and `gap` to the position where the disentangler below is
missing. Both `op.mid` and `op.gap` are of type `T`.
"""
struct ModifiedBinaryOp{T}
    mid::T
    gap::T
end

ModifiedBinaryOp(op::T) where T = ModifiedBinaryOp{T}(op, op)
ModifiedBinaryOp(op::ModifiedBinaryOp) = op  # Makes some method signatures simpler to write
Base.convert(::Type{<: ModifiedBinaryOp}, op::AbstractTensorMap) = ModifiedBinaryOp(op)

function Base.convert(::Type{ModifiedBinaryOp{T}}, op::ModifiedBinaryOp) where {T}
    return ModifiedBinaryOp{T}(convert(T, op.mid), convert(T, op.gap))
end

Base.iterate(op::ModifiedBinaryOp) = (op.mid, Val(1))
Base.iterate(op::ModifiedBinaryOp, state::Val{1}) = (op.gap, Val(2))
Base.iterate(op::ModifiedBinaryOp, state::Val{2}) = nothing
Base.length(op::ModifiedBinaryOp) = 2

Base.eltype(op::ModifiedBinaryOp) = promote_type((eltype(x) for x in op)...)
Base.copy(op::ModifiedBinaryOp) = ModifiedBinaryOp((deepcopy(x) for x in op)...)
Base.adjoint(op::ModifiedBinaryOp) = ModifiedBinaryOp((x' for x in op)...)
Base.imag(op::ModifiedBinaryOp) = ModifiedBinaryOp((imag(x) for x in op)...)
Base.real(op::ModifiedBinaryOp) = ModifiedBinaryOp((real(x) for x in op)...)
# The below could just as well use op.gap.
Base.one(op::ModifiedBinaryOp) = ModifiedBinaryOp(one(op.mid))

TensorKit.storagetype(::Type{ModifiedBinaryOp{T}}) where {T} = storagetype(T)

function pad_with_zeros_to(op::ModifiedBinaryOp, args...)
    return ModifiedBinaryOp((pad_with_zeros_to(x, args...) for x in op)...)
end

function gershgorin_bounds(op::ModifiedBinaryOp)
    lb_mid, ub_mid = gershgorin_bounds(op.mid)
    lb_gap, ub_gap = gershgorin_bounds(op.gap)
    return min(lb_mid, lb_gap), max(ub_mid, ub_gap)
end

support(op::ModifiedBinaryOp) = support(op.mid)  # Could equally well be op.gap.

function expand_support(op::ModifiedBinaryOp, n::Integer)
    mid = expand_support(op.mid, n)
    gap = expand_support(op.gap, n)
    return ModifiedBinaryOp(mid, gap)
end

function Base.similar(op::ModifiedBinaryOp, ::Type{elementT} = eltype(op)) where {elementT}
    mid = similar(op.mid, elementT)
    gap = similar(op.gap, elementT)
    return ModifiedBinaryOp(mid, gap)
end

TensorKit.space(op::ModifiedBinaryOp) = space(op.gap)
TensorKit.space(op::ModifiedBinaryOp, n) = space(op.gap, n)
TensorKit.domain(op::ModifiedBinaryOp) = domain(op.gap)
TensorKit.codomain(op::ModifiedBinaryOp) = codomain(op.gap)

function Base.:(==)(op1::ModifiedBinaryOp, op2::ModifiedBinaryOp)
    return (op1.mid == op2.mid) && (op1.gap == op2.gap)
end

function Base.isapprox(op1::ModifiedBinaryOp, op2::ModifiedBinaryOp)
    return (op1.mid ≈ op2.mid) && (op1.gap ≈ op2.gap)
end

# Pass element-wise arithmetic down onto the AbstractTensorMaps. Promote AbstractTensorMaps
# to ModifiedBinaryOps if necessary.
for op in (:+, :-, :/, :*)
    eval(:(
        Base.$(op)(x::ModifiedBinaryOp, y::ModifiedBinaryOp)
        = ModifiedBinaryOp(($(op)(xi, yi) for (xi, yi) in zip(x, y))...)
    ))
    eval( :(
        Base.$(op)(x::AbstractTensorMap, y::ModifiedBinaryOp)
        = $(op)(ModifiedBinaryOp(x), y)
    ))
    eval(:(
        Base.$(op)(x::ModifiedBinaryOp, y::AbstractTensorMap)
        = $(op)(x, ModifiedBinaryOp(y))
    ))
end

Base.:*(op::ModifiedBinaryOp, a::Number) = ModifiedBinaryOp(op.mid * a, op.gap * a)
Base.:*(a::Number, op::ModifiedBinaryOp) = op*a
Base.:/(op::ModifiedBinaryOp, a::Number) = ModifiedBinaryOp(op.mid / a, op.gap / a)

function Base.copy!(op1::ModifiedBinaryOp, op2::ModifiedBinaryOp)
    copy!(op1.mid, op2.mid)
    copy!(op1.gap, op2.gap)
    return op1
end

function Base.fill!(op::ModifiedBinaryOp, a::Number)
    fill!(op.mid, a)
    fill!(op.gap, a)
    return op
end

function LinearAlgebra.dot(op1::ModifiedBinaryOp, op2::ModifiedBinaryOp)
    dotmid = dot(op1.mid, op2.mid)
    dotgap = dot(op1.gap, op2.gap)
    return (dotmid + dotgap) / 2.0
end

function LinearAlgebra.dot(op1::ModifiedBinaryOp, t2::AbstractTensorMap)
    return dot(op1, ModifiedBinaryOp(t2))
end

function LinearAlgebra.dot(t1::AbstractTensorMap, op2::ModifiedBinaryOp)
    return dot(ModifiedBinaryOp(t1), op2)
end

function LinearAlgebra.rmul!(op::ModifiedBinaryOp, a::Number)
    rmul!(op.mid, a)
    rmul!(op.gap, a)
    return op
end

function LinearAlgebra.lmul!(a::Number, op::ModifiedBinaryOp)
    lmul!(a, op.mid)
    lmul!(a, op.gap)
    return op
end

function BLAS.axpby!(a::Number, X::ModifiedBinaryOp, b::Number, Y::ModifiedBinaryOp)
    axpby!(a, X.mid, b, Y.mid)
    axpby!(a, X.gap, b, Y.gap)
    return Y
end

function BLAS.axpby!(a::Number, X::AbstractTensorMap, b::Number, Y::ModifiedBinaryOp)
    return axpby!(a, ModifiedBinaryOp(X), b, Y)
end

function BLAS.axpby!(a::Number, X::ModifiedBinaryOp, b::Number, Y::AbstractTensorMap)
    return axpby!(a, X, b, ModifiedBinaryOp(Y))
end

function BLAS.axpy!(a::Number, X::ModifiedBinaryOp, Y::ModifiedBinaryOp)
    axpy!(a, X.mid, Y.mid)
    axpy!(a, X.gap, Y.gap)
    return Y
end

function BLAS.axpy!(a::Number, X::AbstractTensorMap, Y::ModifiedBinaryOp)
    return axpy!(a, ModifiedBinaryOp(X), Y)
end

function BLAS.axpy!(a::Number, X::ModifiedBinaryOp, Y::AbstractTensorMap)
    return axpy!(a, X, ModifiedBinaryOp(Y))
end

function LinearAlgebra.mul!(C::ModifiedBinaryOp, A::ModifiedBinaryOp, B::Number)
    mul!(C.mid, A.mid, B)
    mul!(C.gap, A.gap, B)
    return C
end

function LinearAlgebra.mul!(C::ModifiedBinaryOp, A::Number, B::ModifiedBinaryOp)
    mul!(C.mid, A, B.mid)
    mul!(C.gap, A, B.gap)
    return C
end

function LinearAlgebra.mul!(C::AbstractTensorMap, A::ModifiedBinaryOp, B::Number)
    return mul!(ModifiedBinaryOp(C), A, B)
end

function LinearAlgebra.mul!(C::ModifiedBinaryOp, A::AbstractTensorMap, B::Number)
    return mul!(C, ModifiedBinaryOp(A), B)
end

function LinearAlgebra.mul!(C::AbstractTensorMap, A::Number, B::ModifiedBinaryOp)
    return mul!(ModifiedBinaryOp(C), A, B)
end

function LinearAlgebra.mul!(C::ModifiedBinaryOp, A::Number, B::AbstractTensorMap)
    return mul!(C, A, ModifiedBinaryOp(B))
end

LinearAlgebra.tr(op::ModifiedBinaryOp) = (tr(op.mid) + tr(op.gap)) / 2.0

# The entropy of the density matrix is the average over the two different density matrices.
function densitymatrix_entropy(rho::ModifiedBinaryOp)
    rho_mid, rho_gap = rho
    S_mid, S_gap = densitymatrix_entropy(rho_mid), densitymatrix_entropy(rho_gap)
    S = (S_mid + S_gap) / 2.0
    return S
end


text:
const Shrinkage = Union{Symbol, Real}
abstract type LinearShrinkageTarget end

# Taxonomy from http://strimmerlab.org/publications/journals/shrinkcov2005.pdf
# page 13

"""
    DiagonalUnitVariance

Target for linear shrinkage: unit matrix.
A subtype of `LinearShrinkageTarget` where

* ``F_{ij}=1`` if ``i=j`` and
* ``F_{ij}=0`` otherwise
"""
struct DiagonalUnitVariance <: LinearShrinkageTarget end

"""
    DiagonalCommonVariance

Target for linear shrinkage: unit matrix multiplied by average variance of
variables.
A subtype of `LinearShrinkageTarget` where

* ``F_{ij}=v`` if ``i=j`` with ``v=\\mathrm{tr}(S)/p`` and
* ``F_{ij}=0`` otherwise
"""
struct DiagonalCommonVariance <: LinearShrinkageTarget end

"""
    DiagonalUnequalVariance

Target for linear shrinkage: diagonal of covariance matrix.
A subtype of `LinearShrinkageTarget` where

* ``F_{ij}=s_{ij}`` if ``i=j`` and
* ``F_{ij}=0`` otherwise
"""
struct DiagonalUnequalVariance <: LinearShrinkageTarget end

"""
    CommonCovariance

Target for linear shrinkage: see `target_C`.
A subtype `LinearShrinkageTarget` where

* ``F_{ij}=v`` if ``i=j`` with ``v=\\mathrm{tr}(S)/p`` and
* ``F_{ij}=c`` with ``c=\\sum_{i\\neq j} S_{ij}/(p(p-1))`` otherwise
"""
struct CommonCovariance <: LinearShrinkageTarget end

"""
    PerfectPositiveCorrelation

Target for linear shrinkage: see `target_E`.
A subtype of `LinearShrinkageTarget` where

* ``F_{ij}=S_{ij}`` if ``i=j`` and
* ``F_{ij}=\\sqrt{S_{ii}S_{jj}}`` otherwise
"""
struct PerfectPositiveCorrelation <: LinearShrinkageTarget end

"""
    ConstantCorrelation

Target for linear shrinkage: see `target_F`.
A subtype of `LinearShrinkageTarget` where

* ``F_{ij}=S_{ij}`` if ``i=j`` and
* ``F_{ij}=\\overline{r}\\sqrt{S_{ii}S_{jj}}`` otherwise where
  ``\\overline{r}`` is the average sample correlation
"""
struct ConstantCorrelation <: LinearShrinkageTarget end

"""
    LinearShrinkage(target, shrinkage; corrected=false)

Linear shrinkage estimator described by equation
``(1 - \\lambda) S + \\lambda F`` where ``S`` is standard covariance matrix,
``F`` is shrinkage target described by argument `target` and ``\\lambda`` is a
shrinkage parameter, either given explicitly in `shrinkage` or automatically
determined according to one of the supported methods.

The corrected estimator is used if `corrected` is true.
"""
struct LinearShrinkage{T<:LinearShrinkageTarget, S<:Shrinkage} <: CovarianceEstimator
    target::T
    shrinkage::S
    corrected::Bool
    function LinearShrinkage(t::TT, s::SS; corrected=false) where {TT<:LinearShrinkageTarget, SS<:Real}
        0 ≤ s ≤ 1 || throw(ArgumentError("Shrinkage value should be between 0 and 1. Got $s."))
        new{TT, SS}(t, s, corrected)
    end
    function LinearShrinkage(t::TT, s::Symbol=:auto; corrected=false) where TT <: LinearShrinkageTarget
        s ∈ (:auto, :lw, :ss, :rblw, :oas) || throw(ArgumentError("Shrinkage method $s not supported."))
        new{TT, Symbol}(t, s, corrected)
    end
end

LinearShrinkage(;
    target::LinearShrinkageTarget=DiagonalUnitVariance(),
    shrinkage::Shrinkage,
    corrected::Bool=false) = LinearShrinkage(target, shrinkage, corrected=corrected)

"""
    cov(lse::LinearShrinkage, X; dims=1)

Linear shrinkage covariance estimator for matrix `X` along dimension `dims`.
Computed using the method described by `lse`.
"""
function cov(lse::LinearShrinkage, X::AbstractMatrix{<:Real};
             dims::Int=1, mean=nothing)

    dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))

    Xc   = (dims == 1) ? copy(X) : copy(transpose(X))
    n, p = size(Xc)
    # sample covariance of size (p x p)
    S = cov(SimpleCovariance(corrected=lse.corrected), X; dims=dims, mean=mean)

    # NOTE: don't need to check if mean is proper as this is already done above
    if mean === nothing
        Xc .-= Statistics.mean(Xc, dims=1)
    elseif mean isa AbstractArray
        if dims == 1
            Xc .-= mean
        else
            Xc .-= mean'
        end
    end

    return linear_shrinkage(lse.target, Xc, S, lse.shrinkage, n, p, lse.corrected)
end

##############################################################################
# Helper functions for Linear Shrinkage estimators

"""
    rescale(M, d)

Internal function to scale the rows and the columns of a square matrix `M`
according to the elements in `d`.
This is useful when dealing with the standardised data matrix which can be
written `Xs=Xc*D` where `Xc` is the centered data matrix and `D=Diagonal(d)`
so that its simple covariance is `D*S*D` where `S` is the simple covariance
of `Xc`. Such `D*M*D` terms appear often in the computations of optimal
shrinkage λ.

* Space complexity: ``O(p^2)``
* Time complexity: ``O(p^2)``
"""
rescale(M::AbstractMatrix, d::AbstractVector) = M .* d .* d'

"""
    rescale!(M, d)

Same as `rescale` but in place (no allocation).
"""
rescale!(M::AbstractMatrix, d::AbstractVector) = (M .*= d .* d')

"""
    uccov(X)

Internal function to compute `X*X'/n` where `n` is the number of rows of `X`.
This corresponds to the uncorrected covariance of `X` if `X` is centered.
This operation appears often in the computations of optimal shrinkage λ.

* Space complexity: ``O(p^2)``
* Time complexity: ``O(2np^2)``
"""
uccov(X::AbstractMatrix) = (X' * X) / size(X, 1)

"""
    sumij(S)

Internal function to compute the sum of elements of a square matrix `S`.
A keyword `with_diag` can be passed to indicate whether to include or not the
diagonal of `S` in the sum. Both cases happen often in the computations of
optimal shrinkage λ.

* Space complexity: ``O(1)``
* Time complexity: ``O(p^2)``
"""
function sumij(S::AbstractMatrix; with_diag=false)
    acc = sum(S)
    with_diag || return acc - tr(S)
    return acc
end

"""
    sumij2(S)

Internal function identical to `sumij` except that it passes the function
`abs2` to the sum so that it is the sum of the elements of `S` squared
which is computed. This is significantly more efficient than using
`sumij(S.^2)` for large matrices as it allocates very little.

* Space complexity: ``O(1)``
* Time complexity: ``O(2p^2)``
"""
function sumij2(S::AbstractMatrix; with_diag=false)
    acc = sum(abs2, S)
    with_diag || return acc - sum(abs2, diag(S))
    return acc
end

"""
    sum_fij(Xc, S, n, κ)

Internal function corresponding to ``∑_{i≂̸j}f_{ij}`` that appears in
http://strimmerlab.org/publications/journals/shrinkcov2005.pdf p.11.

* Space complexity: ``O(np + 2p^2)``
* Time complexity: ``O(2np^2)``
"""
function sum_fij(Xc, S, n, κ)
    sd  = sqrt.(diag(S))
    M   = ((Xc.^3)' * Xc) ./ sd
    M .-= κ .* S .* sd
    M .*= sd'
    return sumij(M) / (n * κ)
end
##############################################################################

"""
    linear_shrinkage(target, Xc, S, λ, n, p, corrected)

Performs linear shrinkage with target of type `target` for data matrix `Xc`
of size `n` by `p` with covariance matrix `S` and shrinkage parameter `λ`.
Calculates corrected covariance if `corrected` is true.
"""
linear_shrinkage

## TARGET A

function linear_shrinkage(::DiagonalUnitVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    return linshrink(I, S, λ)
end

"""
    linear_shrinkage(::DiagonalUnitVariance, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a `DiagonalUnitVariance`.
"""
function linear_shrinkage(::DiagonalUnitVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F   = I
    κ   = n - Int(corrected)
    γ   = κ/n
    Xc² = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        ΣS² = sumij2(S, with_diag=true)
        λ   = sumij(uccov(Xc²), with_diag=true) / γ^2 - ΣS²
        λ  /= κ * (ΣS² - 2tr(S) + p)
    elseif λ == :ss
        # use the standardised data matrix
        d   = 1.0 ./ vec(sum(Xc², dims=1))
        S̄   = rescale(S, sqrt.(d)) # this has diagonal 1/κ
        ΣS̄² = sumij2(S̄, with_diag=true)
        λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
        λ  /= κ * ΣS̄² - p / κ
    else
        throw(ArgumentError("Unsupported shrinkage method for target DiagonalUnitVariance: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink(F, S, λ)
end

## TARGET B

target_B(S::AbstractMatrix, p::Int) = tr(S)/p * I

function linear_shrinkage(::DiagonalCommonVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    return linshrink(target_B(S, p), S, λ)
end

"""
    linear_shrinkage(::DiagonalCommonVariance, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a `DiagonalCommonVariance`.
"""
function linear_shrinkage(::DiagonalCommonVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F   = target_B(S, p)
    κ   = n - Int(corrected)
    γ   = κ/n
    Xc² = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        v   = F.λ # tr(S)/p
        ΣS² = sumij2(S, with_diag=true)
        λ   = sumij(uccov(Xc²), with_diag=true) / γ^2 - ΣS²
        λ  /= κ * (ΣS² - p*v^2)
    elseif λ == :ss
        # use the standardised data matrix
        d   = 1.0 ./ vec(sum(Xc², dims=1))
        S̄   = rescale(S, sqrt.(d)) # this has diagonal 1/κ
        v̄   = κ # tr(S̄)/p
        ΣS̄² = sumij2(S̄, with_diag=true)
        λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
        λ  /= κ * ΣS̄² - p/κ
    elseif λ == :rblw
        # https://arxiv.org/pdf/0907.4698.pdf equations 17, 19
        trS² = sum(abs2, S)
        tr²S = tr(S)^2
        # note: using corrected or uncorrected S does not change λ
        λ = ((n-2)/n * trS² + tr²S) / ((n+2) * (trS² - tr²S/p))
    elseif λ == :oas
        # https://arxiv.org/pdf/0907.4698.pdf equation 23
        trS² = sum(abs2, S)
        tr²S = tr(S)^2
        # note: using corrected or uncorrected S does not change λ
        λ = ((1.0-2.0/p) * trS² + tr²S) / ((n+1.0-2.0/p) * (trS² - tr²S/p))
    else
        throw(ArgumentError("Unsupported shrinkage method for target DiagonalCommonVariance: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink(F, S, λ)
end

## TARGET D

target_D(S::AbstractMatrix) = Diagonal(S)

function linear_shrinkage(::DiagonalUnequalVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    return linshrink(target_D(S), S, λ)
end

"""
    linear_shrinkage(::DiagonalUnequalVariance, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a `DiagonalUnequalVariance`.
"""
function linear_shrinkage(::DiagonalUnequalVariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F   = target_D(S)
    κ   = n - Int(corrected)
    γ   = κ/n
    Xc² = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        ΣS² = sumij2(S)
        λ   = sumij(uccov(Xc²)) / γ^2 - ΣS²
        λ  /= κ * ΣS²
    elseif λ == :ss
        # use the standardised data matrix
        d   = 1.0 ./ vec(sum(Xc², dims=1))
        ΣS̄² = sumij2(rescale(S, sqrt.(d)))
        λ   = sumij(rescale!(uccov(Xc²), d)) / γ^2 - ΣS̄²
        λ  /= κ * ΣS̄²
    else
        throw(ArgumentError("Unsupported shrinkage method for target DiagonalUnequalVariance: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink(F, S, λ)
end

## TARGET C

function target_C(S::AbstractMatrix, p::Int)
    v  = tr(S)/p
    c  = sumij(S; with_diag=false) / (p * (p - 1))
    F  = fill(c, (p, p))
    F -= Diagonal(F)
    F += v * I
    return F, v, c
end

function linear_shrinkage(::CommonCovariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    F, _, _ = target_C(S, p)
    return linshrink!(F, S, λ)
end

"""
    linear_shrinkage(::CommonCovariance, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a `CommonCovariance`.
"""
function linear_shrinkage(::CommonCovariance, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F, v, c = target_C(S, p)
    κ   = n - Int(corrected)
    γ   = κ/n
    Xc² = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        ΣS² = sumij2(S, with_diag=true)
        λ   = sumij(uccov(Xc²), with_diag=true) / γ^2 - ΣS²
        λ  /= κ * (ΣS² - p*(p-1)*c^2 - p*v^2)
    elseif λ == :ss
        d   = 1.0 ./ vec(sum(Xc², dims=1))
        S̄   = rescale(S, sqrt.(d))
        ΣS̄² = sumij2(S̄, with_diag=true)
        λ   = sumij(rescale!(uccov(Xc²), d), with_diag=true) / γ^2 - ΣS̄²
        λ  /= κ * ΣS̄² - p/κ - κ * sumij(S̄; with_diag=false)^2 / (p * (p - 1))
    else
        throw(ArgumentError("Unsupported shrinkage method for target CommonCovariance: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink!(F, S, λ)
end

## TARGET E

function target_E(S::AbstractMatrix)
    d = sqrt.(diag(S))
    return d*d'
end

function linear_shrinkage(::PerfectPositiveCorrelation, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    return linshrink!(target_E(S), S, λ)
end

"""
    linear_shrinkage(::PerfectPositiveCorrelation, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a
`PerfectPositiveCorrelation`.
"""
function linear_shrinkage(::PerfectPositiveCorrelation, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F   = target_E(S)
    κ   = n - Int(corrected)
    γ   = κ/n
    Xc² = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        ΣS² = sumij2(S)
        λ   = (sumij(uccov(Xc²)) / γ^2 - ΣS²) / κ
        λ  -= sum_fij(Xc, S, n, κ)
        λ  /= sumij2(S - F)
    elseif λ == :ss
        d   = 1.0 ./ vec(sum(Xc², dims=1))
        s   = sqrt.(d)
        S̄   = rescale(S, s)
        ΣS̄² = sumij2(S̄)
        λ   = (sumij(rescale!(uccov(Xc²), d)) / γ^2 - ΣS̄²) / κ
        λ  -= sum_fij(Xc .* s', S̄, n, κ)
        F̄  = target_E(S̄)
        λ  /= sumij2(S̄ - F̄)
    else
        throw(ArgumentError("Unsupported shrinkage method for target PerfectPositiveCorrelation: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink!(F, S, λ)
end

## TARGET F

function target_F(S::AbstractMatrix, p::Int)
    s  = sqrt.(diag(S))
    s_ = s*s'
    r̄  = (sum(S ./ s_) - p) / (p * (p - 1))
    F_ = r̄ * s_
    F  = F_ + (Diagonal(s_) - Diagonal(F_))
    return F, r̄
end

function linear_shrinkage(::ConstantCorrelation, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Real, n::Int, p::Int,
                          corrected::Bool)

    F, _ = target_F(S, p)
    return linshrink!(F, S, λ)
end

"""
    linear_shrinkage(::ConstantCorrelation, Xc, S, λ, n, p, corrected)

Compute the shrinkage estimator where the target is a `ConstantCorrelation`.
"""
function linear_shrinkage(::ConstantCorrelation, Xc::AbstractMatrix,
                          S::AbstractMatrix, λ::Symbol, n::Int, p::Int,
                          corrected::Bool)

    F, r̄ = target_F(S, p)
    κ    = n - Int(corrected)
    γ    = κ/n
    Xc²  = Xc.^2
    # computing the shrinkage
    if λ ∈ [:auto, :lw]
        ΣS² = sumij2(S)
        λ   = (sumij(uccov(Xc²)) / γ^2 - ΣS²) / κ
        λ  -= r̄ * sum_fij(Xc, S, n, κ)
        λ  /= sumij2(S - F)
    elseif λ == :ss
        d    = 1.0 ./ vec(sum(Xc², dims=1))
        s    = sqrt.(d)
        S̄    = rescale(S, s)
        F̄, r̄ = target_F(S̄, p)
        ΣS̄²  = sumij2(S̄)
        λ    = (sumij(rescale!(uccov(Xc²), d)) / γ^2 - ΣS̄²) / κ
        λ   -= r̄ * sum_fij(Xc .* s', S̄, n, κ)
        λ   /= sumij2(S̄ - F̄)
    else
        throw(ArgumentError("Unsupported shrinkage method for target ConstantCorrelation: $λ."))
    end
    λ = clamp(λ, 0.0, 1.0)
    return linshrink!(F, S, λ)
end

ext: jl


text:
using Enzyme
using Test

@testset "Threads $(Threads.nthreads())" begin
    function tasktest(M, x)
        xr = Ref(x)
        task = Threads.@spawn begin
            @inbounds M[1] = xr[]
        end
        @inbounds M[2] = x
        wait(task)
        nothing
    end

    R = Float64[0., 0.]
    dR = Float64[2., 3.]

    @test 5.0 ≈ Enzyme.autodiff(tasktest, Duplicated(R, dR), Active(2.0))[1]
    @test Float64[2.0, 2.0] ≈ R
    @test Float64[0.0, 0.0] ≈ dR

    Enzyme.fwddiff(tasktest, Duplicated(R, dR), Duplicated(2.0, 1.0))
    @test Float64[1.0, 1.0] ≈ dR

    function tasktest2(M, x)
        task = Threads.@spawn begin
           return
        end
        Base.wait(task)
        nothing
    end
    # The empty return previously resulted in an illegal instruction error
    @test 0.0 ≈ Enzyme.autodiff(tasktest2, Duplicated(R, dR), Active(2.0))[1]
    @test () === Enzyme.fwddiff(tasktest, Duplicated(R, dR), Duplicated(2.0, 1.0))
end

@testset "Advanced Threads $(Threads.nthreads())" begin
    function foo(y)
        Threads.@threads for i in 1:3
            y[i] *= 2
        end
        nothing
    end

    x = [1.0, 2.0, 3.0]
    dx = [1.0, 1.0, 1.0]
    Enzyme.autodiff(foo, Duplicated(x, dx))
    @test 2.0 ≈ x[1]
    @test 4.0 ≈ x[2]
    @test 6.0 ≈ x[3]
    @test 2.0 ≈ dx[1]
    @test 2.0 ≈ dx[2]
    @test 2.0 ≈ dx[3]
end


text:
# CLIMA solver configurations
#
# Contains helper functions to establish solver configurations to be
# used with the CLIMA driver.

"""
    CLIMA.SolverConfiguration

Parameters needed by `CLIMA.solve!()` to run a simulation.
"""
struct SolverConfiguration{FT}
    name::String
    mpicomm::MPI.Comm
    dg::DGModel
    Q::MPIStateArray
    t0::FT
    timeend::FT
    dt::FT
    init_on_cpu::Bool
    numberofsteps::Int
    init_args
    solver
end

"""
    DGmethods.courant(local_cfl, solver_config::SolverConfiguration;
                      Q=solver_config.Q, dt=solver_config.dt)

Returns the maximum of the evaluation of the function `local_courant`
pointwise throughout the domain with the model defined by `solver_config`. The
keyword arguments `Q` and `dt` can be used to call the courant method with a
different state `Q` or time step `dt` than are defined in `solver_config`.
"""
DGmethods.courant(
    f,
    sc::SolverConfiguration;
    Q = sc.Q,
    dt = sc.dt,
    simtime = gettime(sc.solver),
    direction = EveryDirection(),
) = DGmethods.courant(f, sc.dg, sc.dg.balancelaw, Q, dt, simtime, direction)

"""
    CLIMA.SolverConfiguration(
        t0::FT,
        timeend::FT,
        driver_config::DriverConfiguration,
        init_args...;
        init_on_cpu=false,
        ode_solver_type=driver_config.solver_type,
        ode_dt=nothing,
        modeldata=nothing,
        Courant_number=0.4,
        diffdir=EveryDirection(),
    )

Set up the DG model per the specified driver configuration, set up
the ODE solver, and return a `SolverConfiguration` to be used with
`CLIMA.invoke!()`.

# Arguments:
# - `t0::FT`: simulation start time.
# - `timeend::FT`: simulation end time.
# - `driver_config::DriverConfiguration`: from `AtmosLESConfiguration()`, etc.
# - `init_args...`: passed through to `init_state!()`.
# - `init_on_cpu=false`: run `init_state!()` on CPU?
# - `ode_solver_type=driver_config.solver_type`: override solver choice.
# - `ode_dt=nothing`: override timestep computation.
# - `modeldata=nothing`: passed through to `DGModel`.
# - `Courant_number=0.4`: for the simulation.
# - `diffdir=EveryDirection()`: restrict diffusivity direction.
# - `timeend_dt_adjust=true`: should `dt` be adjusted to hit `timeend` exactly
# - `CFL_direction=EveryDirection()`: direction for `calculate_dt`
"""
function SolverConfiguration(
    t0::FT,
    timeend::FT,
    driver_config::DriverConfiguration,
    init_args...;
    init_on_cpu = false,
    ode_solver_type = driver_config.solver_type,
    ode_dt = nothing,
    modeldata = nothing,
    Courant_number = nothing,
    diffdir = EveryDirection(),
    timeend_dt_adjust = true,
    CFL_direction = EveryDirection(),
) where {FT <: AbstractFloat}
    @tic SolverConfiguration

    bl = driver_config.bl
    grid = driver_config.grid
    numfluxnondiff = driver_config.numfluxnondiff
    numfluxdiff = driver_config.numfluxdiff
    gradnumflux = driver_config.gradnumflux

    # create DG model, initialize ODE state
    if Settings.restart_from_num > 0
        s_Q, s_aux, t0 = read_checkpoint(
            driver_config.name,
            driver_config.mpicomm,
            Settings.restart_from_num,
        )

        auxstate = restart_auxstate(bl, grid, s_aux)

        dg = DGModel(
            bl,
            grid,
            numfluxnondiff,
            numfluxdiff,
            gradnumflux,
            auxstate = auxstate,
            diffusion_direction = diffdir,
            modeldata = modeldata,
        )

        @info @sprintf("Restarting %s from time %8.2f", driver_config.name, t0)
        Q = restart_ode_state(dg, s_Q; init_on_cpu = init_on_cpu)
    else
        dg = DGModel(
            bl,
            grid,
            numfluxnondiff,
            numfluxdiff,
            gradnumflux,
            diffusion_direction = diffdir,
            modeldata = modeldata,
        )

        @info @sprintf("Initializing %s", driver_config.name,)
        Q = init_ode_state(dg, FT(0), init_args...; init_on_cpu = init_on_cpu)
    end
    update_aux!(dg, bl, Q, FT(0), dg.grid.topology.realelems)

    # create the linear model for IMEX solvers
    linmodel = nothing
    if isa(ode_solver_type, ExplicitSolverType)
        dtmodel = bl
    else # ode_solver_type === IMEXSolverType
        linmodel = ode_solver_type.linear_model(bl)
        dtmodel = linmodel
    end

    # default Courant number
    if Courant_number == nothing
        if ode_solver_type.solver_method == LSRK144NiegemannDiehlBusch
            Courant_number = FT(1.7)
        elseif ode_solver_type.solver_method == LSRK54CarpenterKennedy
            Courant_number = FT(0.3)
        else
            Courant_number = FT(0.5)
        end
    end

    # initial Δt specified or computed
    simtime = FT(0) # TODO: needs to be more general to account for restart:
    if ode_dt === nothing
        ode_dt = CLIMA.DGmethods.calculate_dt(
            dg,
            dtmodel,
            Q,
            Courant_number,
            simtime,
            CFL_direction,
        )
    end
    numberofsteps = convert(Int, cld(timeend, ode_dt))
    timeend_dt_adjust && (ode_dt = timeend / numberofsteps)

    # create the solver
    if isa(ode_solver_type, ExplicitSolverType)
        solver = ode_solver_type.solver_method(dg, Q; dt = ode_dt, t0 = t0)
    elseif isa(ode_solver_type, MultirateSolverType)
        fast_dg = DGModel(
            linmodel,
            grid,
            numfluxnondiff,
            numfluxdiff,
            gradnumflux,
            auxstate = dg.auxstate,
        )
        slow_model = RemainderModel(bl, (linmodel,))
        slow_dg = DGModel(
            slow_model,
            grid,
            numfluxnondiff,
            numfluxdiff,
            gradnumflux,
            auxstate = dg.auxstate,
        )
        slow_solver = ode_solver_type.slow_method(slow_dg, Q; dt = ode_dt)
        fast_dt = ode_dt / ode_solver_type.timestep_ratio
        fast_solver = ode_solver_type.fast_method(fast_dg, Q; dt = fast_dt)
        solver =
            ode_solver_type.solver_method((slow_solver, fast_solver), t0 = t0)
    else # solver_type === IMEXSolverType
        vdg = DGModel(
            linmodel,
            grid,
            numfluxnondiff,
            numfluxdiff,
            gradnumflux,
            auxstate = dg.auxstate,
            direction = VerticalDirection(),
        )
        solver = ode_solver_type.solver_method(
            dg,
            vdg,
            ode_solver_type.linear_solver(),
            Q;
            dt = ode_dt,
            t0 = t0,
        )
    end

    @toc SolverConfiguration

    return SolverConfiguration(
        driver_config.name,
        driver_config.mpicomm,
        dg,
        Q,
        t0,
        timeend,
        ode_dt,
        init_on_cpu,
        numberofsteps,
        init_args,
        solver,
    )
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
"""
	dtw(sequence1::Matrix{Cdouble}, sequence2::Matrix{Cdouble}) -> (path1, path2)

Compute the set of indices `path1, path2` that give the alignment of two
sequences over their last dimension by dynamic time warping.

Returns two vectors representing the path.
"""
function dtw(
	sequence1::Matrix{Cdouble},
	sequence2::Matrix{Cdouble}
)

	length1 = size(sequence1)[end]
	length2 = size(sequence2)[end]
	cost = 1 .- sequence1'*sequence2

	# Fill the first column...
	for i in 2:length1
		cost[i, 1] = cost[i - 1, 1] + cost[i, 1]
	end

	# ...and the first row of the cost matrix.
	for j in 2:length2
		cost[1, j] = cost[1, j - 1] + cost[1, j]
	end

	# Compute the cumulated cost of the remaining elements
	for i in 2:length1
		for j in 2:length2
			cost[i, j] = cost[i, j] +
				min(cost[i - 1, j], cost[i - 1, j - 1], cost[i, j - 1])
		end
	end

	return trackback(cost)
end

"""
	trackback(cost::AbstractMatrix) -> (path1, path2)

Compute the lowest cost path from a cost array. Returns two vectors
representing the path.
"""
function trackback(cost::Matrix{Cdouble})

	i, j = size(cost)
	# Start to create the path from the last element of the cost matrix.
	path1, path2 = [i], [j]

	while i > 1 && j > 1
		least = argmin([cost[i-1,j-1], cost[i-1,j], cost[i,j-1]])

		# Everytime that the minimal cost is reached by going back along one
		# dimension, decrease the corresponding index.
		push!(path1, (least in [1 2] ? i-=1 : i))
		push!(path2, (least in [1 3] ? j-=1 : j))
	end

	# One of the indexes may still be greater than one if the path hits a border
	# of the cost matrix, in this case the path must be completed with a straight
	# line to the first element of the cost matrix.
	if i > 1
		append!(path1, i-1:-1:1)
		append!(path2, ones(length(path1) - length(path2)))
	end

	if j > 1
		append!(path2, j-1:-1:1)
		append!(path1, ones(length(path2) - length(path1)))
	end

	return reverse(path1), reverse(path2)
end

max_stars_repo_name: gabriel-ss/interpretation-switcher

max_issues_repo_name: gabriel-ss/interpretation-switcher

max_forks_repo_name: gabriel-ss/interpretation-switcher


text:
#=
Chr support

Copyright 2017-2018 Gandalf Software, Inc., Scott P. Jones
Licensed under MIT License, see LICENSE.md
In part based on code for Char in Julia
=#

"""Type of codeunits"""
codeunit(::Type{UniStr})    = UInt32 # Note, the real type could be UInt8, UInt16, or UInt32
codeunit(::Type{<:Str{C}}) where {C<:CSE} = codeunit(C)

codeunit(::MaybeSub{S}) where {S<:Str} = codeunit(S)
codeunit(::Type{<:MaybeSub{S}}) where {S<:Str} = codeunit(S)

eltype(::Type{<:Str{BinaryCSE}}) = UInt8

eltype(::Type{UniStr})                = UTF32Chr
eltype(::Type{<:Str{Text1CSE}})       = Text1Chr
eltype(::Type{<:Str{Text2CSE}})       = Text2Chr
eltype(::Type{<:Str{Text4CSE}})       = Text4Chr
eltype(::Type{<:Str{ASCIICSE}})       = ASCIIChr
eltype(::Type{<:Str{LatinCSE}})       = LatinChr
eltype(::Type{<:Str{_LatinCSE}})      = _LatinChr
eltype(::Type{<:Str{<:Union{UCS2CSE, _UCS2CSE}}}) = UCS2Chr
eltype(::Type{<:Str{<:Union{UTF8CSE, UTF16CSE, UTF32_CSEs}}}) = UTF32Chr

get_codeunit(dat, pos) = codeunit(dat, pos)
get_codeunit(pnt::Ptr{<:CodeUnitTypes}, pos) = unsafe_load(pnt, pos)
get_codeunit(dat::AbstractArray{<:CodeUnitTypes}, pos) = dat[pos]
get_codeunit(str::MaybeSub{<:Str}, pos::Integer) = @preserve str get_codeunit(pointer(str), pos)

get_codeunit(dat) = get_codeunit(dat, 1)
get_codeunit(pnt::Ptr{<:CodeUnitTypes}) = unsafe_load(pnt)

codeunit(str::Str, pos::Integer) = get_codeunit(str, pos)

set_codeunit!(pnt::Ptr{<:CodeUnitTypes}, pos, ch) = unsafe_store!(pnt, ch, pos)
set_codeunit!(dat::AbstractArray{<:CodeUnitTypes}, pos, ch) = (dat[pos] = ch)
set_codeunit!(dat::String, pos, ch) = unsafe_store!(pointer(dat), pos, ch)

set_codeunit!(pnt::Ptr{<:CodeUnitTypes}, ch) = unsafe_store!(pnt, ch)
set_codeunit!(dat::AbstractArray{<:CodeUnitTypes}, ch) = (dat[1] = ch)
set_codeunit!(dat::String, ch) = set_codeunit!(dat, 1, ch)


text:
## An image registration example (a difficult problem with many local minima)
#
# The goal is to discover the translation and rotation needed to align two images.
# You can visualize the images below with
#    using ImageView
#    imshow(img)
#    imshow(imgrot)
# You'll see that they differ by a big rotation.

# You could make this more efficient by using Fourier methods to compute the ideal
# translation given any rotation. But the point here is to illustrate the brute-force approach.

using Images, TestImages, CoordinateTransformations, QuadDIRECT, Test

# Takes a parameter vector `x` and returns a rigid transformation
function tfmx(x, img)
    dx, dy, θ = x
    return Translation(dx, dy) ∘ recenter(RotMatrix(θ), center(img))
end

# Compute the remaining mismatch after applying the transformation implied by `x` to `img`
function mismatch(x, img, target)
    tfm = tfmx(x, img)
    imgw = warp(img, tfm)
    mismatch(imgw, target)
end

function mismatch(img1, img2)
    # Compute the mismatch over the indices that the two images have in common
    inds1, inds2 = axes(img1), axes(img2)
    inds = (intersect(inds1[1], inds2[1]), intersect(inds1[2], inds2[2]))
    (isempty(inds[1]) || isempty(inds[2])) && return Inf
    val = 0.0
    counter = 0
    for x in inds[2], y in inds[1]
    v1, v2 = img1[y, x], img2[y, x]
    if !isnan(v1) && !isnan(v2)
        val += abs2(v1-v2)
        counter += 1
    end
    end
    counter == 0 && return Inf
    # This is like "nanmean"
    return val/counter
end

const img = float32.(testimage("lighthouse"))
x = [25, -17, 2*pi/3]
tfm = tfmx(x, img)
const imgrot = warp(img, tfm)

# Test that the ground truth gives a very low mismatch (actually, 0)
@test mismatch(x, img, imgrot) < 1e-3

f(x) = mismatch(x, img, imgrot)

upper = [30, 30, pi] # allow a search of 30 pixels in either direction, plus any rotation
lower = -upper
splits = ([-15, 0, 15], [-15, 0, 15], [-pi/2, 0, pi/2])

# Transformations that differ by less than a pixel or less than 0.02 radians are not
# worth comparing
minwidth = [1, 1, 0.02]

# Simple iterative analysis
root, x0 = analyze(f, splits, lower, upper; maxevals=10, minwidth=minwidth)
# Iterative refinement with control over convergence. `fvalue=0.02` forces it to keep looking until it finds something good.
@warn("This will take a long time")
root = analyze!(root, f, x0, splits, lower, upper; maxevals=10000, fvalue=0.01, rtol=0, minwidth=minwidth, print_interval=100)

box = minimum(root)
println("Ground truth: ", x)
println("Solution: ", position(box, x0))


text:
export skip_next

import Base: show

"""
    skip_next()

Creates a `skip_next` operator, which filters out all `next` messages by the source Observable by emitting only
`error` and `complete` messages.

# Producing

Stream of type `<: Subscribable{L}` where `L` refers to type of source stream

# Examples
```jldoctest
using Rocket

source = from([ 1, 2, 3 ])
subscribe!(source |> skip_next(), logger())
;

# output
[LogActor] Completed

```

See also: [`AbstractOperator`](@ref), [`InferableOperator`](@ref), [`ProxyObservable`](@ref), [`skip_error`](@ref), [`skip_complete`](@ref), [`logger`](@ref)
"""
skip_next() = SkipNextOperator()

struct SkipNextOperator <: InferableOperator end

function on_call!(::Type{L}, ::Type{L}, operator::SkipNextOperator, source) where L
    return proxy(L, source, SkipNextProxy())
end

operator_right(operator::SkipNextOperator, ::Type{L}) where L = L

struct SkipNextProxy <: ActorProxy end

actor_proxy!(::Type{L}, proxy::SkipNextProxy, actor::A) where { L, A } = SkipNextActor{L, A}(actor)

struct SkipNextActor{L, A} <: Actor{L}
    actor :: A
end

on_next!(actor::SkipNextActor{L}, data::L) where L = begin end
on_error!(actor::SkipNextActor, err)               = error!(actor.actor, err)
on_complete!(actor::SkipNextActor)                 = complete!(actor.actor)

Base.show(io::IO, ::SkipNextOperator)         = print(io, "SkipNextOperator()")
Base.show(io::IO, ::SkipNextProxy)            = print(io, "SkipNextProxy($L)")
Base.show(io::IO, ::SkipNextActor{L}) where L = print(io, "SkipNextActor($L)")


text:
module MyPkg2

include("func.jl")

export great,
       f,
       g,
       func

greet() = print("Hello World!")
f() = π

g() = 100

end # module


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

const Elasticity2DSurfaceElements = Union{Poi1,Seg2,Seg3}
const Elasticity2DVolumeElements = Union{Tri3,Tri6,Quad4,Quad8,Quad9}

function assemble!(assembly::Assembly, problem::Problem{Elasticity},
                   elements::Union{Vector{Element}, Vector{Element{T}}},
                   time, ::Type{Val{:plane_stress}}) where T
    for element in elements
        gdofs = get_gdofs(problem, element)
        Km, Kg, f = assemble(problem, element, time, Val{:plane})
        add!(assembly.K, gdofs, gdofs, Km)
        add!(assembly.Kg, gdofs, gdofs, Kg)
        add!(assembly.f, gdofs, f)
    end
end

function assemble!(assembly::Assembly, problem::Problem{Elasticity},
                   elements::Union{Vector{Element}, Vector{Element{T}}},
                   time, ::Type{Val{:plane_strain}}) where T
    for element in elements
        gdofs = get_gdofs(problem, element)
        Km, Kg, f = assemble(problem, element, time, Val{:plane})
        add!(assembly.K, gdofs, gdofs, Km)
        add!(assembly.Kg, gdofs, gdofs, Kg)
        add!(assembly.f, gdofs, f)
    end
end

""" Plane elasticity equations (plane stress, plane strain). """
function assemble(problem::Problem{Elasticity},
                  element::Element{El}, time,
                  ::Type{Val{:plane}}) where El<:Elasticity2DVolumeElements

    props = problem.properties
    dim = get_unknown_field_dimension(problem)
    nnodes = length(element)
    BL = zeros(3, dim*nnodes)
    BNL = zeros(4, dim*nnodes)
    Km = zeros(dim*nnodes, dim*nnodes)
    Kg = zeros(dim*nnodes, dim*nnodes)
    f = zeros(dim*nnodes)
    Dtan = zeros(3,3)

    for ip in get_integration_points(element)

        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)
        dN = element(ip, time, Val{:Grad})
        # kinematics

        gradu = element("displacement", ip, time, Val{:Grad})
        fill!(BL, 0.0)

        if props.finite_strain
            strain = 1/2*(gradu + gradu' + gradu'*gradu)
            F = I + gradu
            for i=1:size(dN, 2)
                BL[1, 2*(i-1)+1] += F[1,1]*dN[1,i]
                BL[1, 2*(i-1)+2] += F[2,1]*dN[1,i]
                BL[2, 2*(i-1)+1] += F[1,2]*dN[2,i]
                BL[2, 2*(i-1)+2] += F[2,2]*dN[2,i]
                BL[3, 2*(i-1)+1] += F[1,1]*dN[2,i] + F[1,2]*dN[1,i]
                BL[3, 2*(i-1)+2] += F[2,1]*dN[2,i] + F[2,2]*dN[1,i]
            end
        else # linearized strain
            strain = 1/2*(gradu + gradu')
            F = I
            for i=1:size(dN, 2)
                BL[1, 2*(i-1)+1] = dN[1,i]
                BL[2, 2*(i-1)+2] = dN[2,i]
                BL[3, 2*(i-1)+1] = dN[2,i]
                BL[3, 2*(i-1)+2] = dN[1,i]
            end
        end

        strain_vec = [strain[1,1]; strain[2,2]; strain[1,2]]

        # calculate stress
        E = element("youngs modulus", ip, time)
        nu = element("poissons ratio", ip, time)
        if props.formulation == :plane_stress
            D = E/(1.0 - nu^2) .* [
                1.0  nu 0.0
                nu  1.0 0.0
                0.0 0.0 (1.0-nu)/2.0]
        elseif props.formulation == :plane_strain
            D = E/((1.0+nu)*(1.0-2.0*nu)) .* [
                1.0-nu      nu               0.0
                    nu  1.0-nu               0.0
                   0.0     0.0  (1.0-2.0*nu)/2.0]
        else
            error("unknown plane formulation: $(props.formulation)")
        end

        # calculate stress

        if haskey(element, "plasticity")
            plastic_def = element("plasticity")[ip.id]

            calculate_stress! = plastic_def["type"]
            yield_surface_ = plastic_def["yield_surface"]
            params = plastic_def["params"]

            initialize_internal_params!(params, ip, Val{:type_2d})

            if time == 0.0
                error("Given step time = $(time). Please select time > 0.0")
            end

            t_last = ip("prev_time", time)
            update!(ip, "prev_time", time => t_last)
            dt = time - t_last

            stress_last = ip("stress", t_last)
            strain_last = ip("strain", t_last)

            dstrain_vec = strain_vec - strain_last
            stress_vec = [0.0, 0.0, 0.0]
            pstrain = zeros(3)
            calculate_stress!(stress_vec, stress_last, dstrain_vec, pstrain, D, params, Dtan, yield_surface_, time, dt, Val{:type_2d})
        else
            stress_vec = D * ([1.0, 1.0, 2.0] .* strain_vec)
            Dtan[:,:] = D[:,:]
        end

        :strain in props.store_fields && update!(ip, "strain", time => strain_vec)
        :stress in props.store_fields && update!(ip, "stress", time => stress_vec)
        :stress11 in props.store_fields && update!(ip, "stress11", time => stress_vec[1])
        :stress22 in props.store_fields && update!(ip, "stress22", time => stress_vec[2])
        :stress12 in props.store_fields && update!(ip, "stress12", time => stress_vec[3])

        Km += w*BL'*Dtan*BL

        # stress = [stress_vec[1] stress_vec[3]; stress_vec[3] stress_vec[2]]
        # cauchy_stress = F'*stress*F/det(F)
        # cauchy_stress = [cauchy_stress[1,1]; cauchy_stress[2,2]; cauchy_stress[1,2]]
        # update!(ip, "cauchy stress", time => cauchy_stress)

        # material stiffness end

        if props.geometric_stiffness
            # take geometric stiffness into account

            fill!(BNL, 0.0)
            for i=1:size(dN, 2)
                BNL[1, 2*(i-1)+1] = dN[1,i]
                BNL[2, 2*(i-1)+1] = dN[2,i]
                BNL[3, 2*(i-1)+2] = dN[1,i]
                BNL[4, 2*(i-1)+2] = dN[2,i]
            end

            S2 = zeros(2*dim, 2*dim)
            S2[1,1] = stress_vec[1]
            S2[2,2] = stress_vec[2]
            S2[1,2] = S2[2,1] = stress_vec[3]
            S2[3:4,3:4] = S2[1:2,1:2]

            Kg += w*BNL'*S2*BNL # geometric stiffness

        end

        # rhs, internal and external load

        f -= w*BL'*stress_vec

        if haskey(element, "displacement load")
            b = element("displacement load", ip, time)
            f += w*vec(b*N)
        end

        for i=1:dim
            if haskey(element, "displacement load $i")
                b = element("displacement load $i", ip, time)
                f[i:dim:end] += w*vec(b*N)
            end
        end

    end

    return Km, Kg, f
end

function assemble(problem::Problem{Elasticity},
                  element::Element{El},
                  time, ::Type{Val{:plane}}) where El<:Elasticity2DSurfaceElements

    props = problem.properties
    dim = get_unknown_field_dimension(problem)
    nnodes = length(element)
    Km = zeros(dim*nnodes, dim*nnodes)
    Kg = zeros(dim*nnodes, dim*nnodes)
    f = zeros(dim*nnodes)

    for ip in get_integration_points(element)

        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)

        if haskey(element, "displacement traction force")
            T = element("displacement traction force", ip, time)
            f += w*vec(T*N)
        end

        for i=1:dim
            # traction force for ith component
            if haskey(element, "displacement traction force $i")
                T = element("displacement traction force $i", ip, time)
                f[i:dim:end] += w*vec(T*N)
            end
        end

        if haskey(element, "nt displacement traction force")
            # traction force given in normal-tangential direction
            T = element("nt displacement traction force", ip, time)
            Q = element("normal-tangential coordinates", ip, time)
            f += w*vec(Q'*T*N)
        end

    end

    return Km, Kg, f
end


text:
using ProgressMeter
using LinearAlgebra

  """
      bootstrapPosteriorCorEstimation( predictions::Matrix{Float64}, t::Vector{Float64}, samplingBatchSize::Int64, nrRuns::Int64 )

      compute posterior p( h* = h | S ).
      #Arguments
      - `predictions::Matrix{Float64}`: each column is the prediction of one hypothesis.
      - `t::Vector{Float64}`:           label vector.
      - `samplingBatchSize::Int64`:     sample size per main iteration.
      - `nrRuns::Int64`:                number of main  iterations.
      #Return
      - `Vector{Float64}`:              posterior p( h* = h | S ).
  """
  @views function bootstrapPosteriorCorEstimation( predictions::Matrix{Float64}, t::Vector{Float64}, samplingBatchSize::Int64, nrRuns::Int64 )
    len      = size( predictions )[1];
    width    = size( predictions )[2];
    lenCache = samplingBatchSize;
    res      = zeros( Int64, width );
    @showprogress 1 "Computing..." for i in 1:1:nrRuns
      @inbounds samplingIndexCache = rand( 1:size( predictions, 1 ), lenCache );
      @inbounds label              = argmaxUniProb.argmaxProb( ( transpose( predictions[samplingIndexCache,:] ) * t[samplingIndexCache]  ) );
      @inbounds res[label]        += 1;
    end
    return  res ./ nrRuns;
  end

  """
      bootstrapPosteriorCorEstimation( predictions::Vector{Matrix{Float64}}, T::Matrix{Float64}, samplingFactor::Float64, nrRuns::Int64 )

      compute posterior p( h* = h | S ).
      #Arguments
      - `predictions::Matrix{Float64}`: each column is the prediction of one hypothesis.
      - `T::Matrix{Float64}`:           label matrix.
      - `samplingBatchSize::Int64`:     sample size per main iteration.
      - `nrRuns::Int64`:                number of main  iterations.
      #Return
      - `Vector{Float64}`:              posterior p( h* = h | S ).
  """
  @views function bootstrapPosteriorCorEstimation( predictions::Vector{Matrix{Float64}}, T::Matrix{Float64}, samplingFactor::Float64, nrRuns::Int64 )
    len      = size( T, 1 );
    width    = size( predictions, 1 );
    lenCache = round( Int64, samplingFactor * len );
    res      = zeros( Int64, width );
    cache    = zeros( Float64, width );
    @showprogress 1 "Computing..." for i in 1:1:nrRuns
      samplingIndexCache = rand( 1:len, lenCache );
      for index in 1:1:width
        @inbounds cache[index] = sum( predictions[index][samplingIndexCache,:] .* T[samplingIndexCache,:] )
      end
      label                 = argmaxUniProb.argmaxProb( cache );
      @inbounds res[label] += 1;
    end
    return  res ./ nrRuns;
  end


text:
using JLD2

using Oceananigans.Fields: FunctionField
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Coriolis: HydrostaticSphericalCoriolis, fᶠᶠᵃ
using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel, VectorInvariant
using Oceananigans.TurbulenceClosures: HorizontallyCurvilinearAnisotropicDiffusivity
using Oceananigans.AbstractOperations: KernelFunctionOperation, volume

function run_hydrostatic_free_turbulence_regression_test(grid, free_surface, arch; regenerate_data=false)

    #####
    ##### Constructing Grid and model
    #####

    model = HydrostaticFreeSurfaceModel(grid = grid,
                                architecture = arch,
                          momentum_advection = VectorInvariant(),
                                free_surface = free_surface,
                                    coriolis = HydrostaticSphericalCoriolis(),
                                     closure = HorizontallyCurvilinearAnisotropicDiffusivity(νh=1e+5, κh=1e+4))

    #####
    ##### Imposing initial conditions:
    #####    u = function of latitude
    #####    v = function of longitude
    #####    vertical shear for u-velocity
    #####

    step_function(x, d, c) = 1/2 * (1 + tanh((x - c) / d))
    polar_mask(y)          = step_function(y, -5, 40) * step_function(y, 5, -40)
    shear_func(x, y, z, p) = p.U * (0.5 + z / p.Lz) * polar_mask(y)

    set!(model, u = (λ, φ, z) -> polar_mask(φ) * exp(-φ^2 / 200),
                v = (λ, φ, z) -> polar_mask(φ) * sind(2λ))

    u, v, w = model.velocities
    U       = 0.1 * maximum(abs, u)
    shear   = FunctionField{Face, Center, Center}(shear_func, grid, parameters=(U=U, Lz=grid.Lz))
    u      .= u + shear

    # Time-scale for gravity wave propagation across the smallest grid cell
    # wave_speed is the hydrostatic (shallow water) gravity wave speed
    gravity    = model.free_surface.gravitational_acceleration
    wave_speed = sqrt(gravity * grid.Lz)

    minimum_Δx = grid.radius * cosd(maximum(abs, view(grid.φᵃᶜᵃ, 1:grid.Ny))) * deg2rad(grid.Δλ)
    minimum_Δy = grid.radius * deg2rad(grid.Δφ)

    wave_time_scale = min(minimum_Δx, minimum_Δy) / wave_speed

    # Δt based on wave propagation time scale
    Δt = 0.2 * wave_time_scale

    #####
    ##### Simulation setup
    #####

    stop_iteration = 20

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

    η = model.free_surface.η

    free_surface_str = string(typeof(model.free_surface).name.wrapper)
    x_topology_str = string(topology(grid, 1))
    output_prefix = "hydrostatic_free_turbulence_regression_$(x_topology_str)_$free_surface_str"

    if regenerate_data
        @warn "Generating new data for the Hydrostatic regression test."

        directory =  joinpath(dirname(@__FILE__), "data")
        outputs   = (; u, v, w, η)
        simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs,
                                                              dir = directory,
                                                         schedule = IterationInterval(stop_iteration),
                                                           prefix = output_prefix,
                                                            force = true)
    end

    # Let's gooooooo!
    run!(simulation)

    # Test results
    test_fields = (
        u = Array(interior(u)),
        v = Array(interior(v)),
        w = Array(interior(w)),
        η = Array(interior(η))
    )

    datadep_path = "regression_test_data/" * output_prefix * ".jld2"
    regression_data_path = @datadep_str datadep_path
    file = jldopen(regression_data_path)

    truth_fields = (
        u = file["timeseries/u/$stop_iteration"][:, :, :],
        v = file["timeseries/v/$stop_iteration"][:, :, :],
        w = file["timeseries/w/$stop_iteration"][:, :, :],
        η = file["timeseries/η/$stop_iteration"][:, :, :]
    )

    close(file)

    summarize_regression_test(test_fields, truth_fields)

    @test all(test_fields.u .≈ truth_fields.u)
    @test all(test_fields.v .≈ truth_fields.v)
    @test all(test_fields.w .≈ truth_fields.w)
    @test all(test_fields.η .≈ truth_fields.η)

    return nothing
end


text:
export conf_coord, cc_from_dict, calc_overlap!, calc_capt_coeff!

# define constants
ħ = 6.582119514E-16 # eV⋅s
kB  = 8.6173303E-5 # eV⋅K⁻¹

# occupation of highest eigenvalue for each potential
# (should be small to ensure convergence of partition function)
occ_cut_off = 1E-5

"""
Stores two `potential`s with e-ph coupling constant `W` to calculate the capture coefficient `capt_coeff`(`temperature`).

## Fields

- `name` -- the name of a configuration coordinate.
- `V1` and `V2` -- the initial and fianal `potential`s.
- `W` -- the e-ph coupling matrix element.
- `g` -- the degeneracy.
- `temperature` -- the temperature range where `capt_coeff` is calculated.
- `capt_coeff` -- the capture coefficient.

## Constructor
    conf_coord(pot_i::potential, pot_f::potential)

"""
mutable struct conf_coord
    # Configuration coordinate
    name::String
    # potentials
    V1::potential; V2::potential
    # e-ph coupling matrix element; degeneracy
    W::Float64; g::Int
    # to calculate vibrational wave function overlap integral:
    # (initial state) phonon eigenvalue; phonon overlap; Gaussian function energy
    ϵ_matrix::Array{Float64,2}; overlap_matrix::Array{Float64,2}; δ_matrix::Array{Float64,2}
    temperature::Array{Float64,1}; capt_coeff::Array{Float64,1}
    partial_capt_coeff::Array{Float64, 3}
end
conf_coord(pot_i::potential, pot_f::potential) = conf_coord("", pot_i, pot_f, Inf, 1,
           Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0), Array{Float64}(undef, 0, 0),
           [], [], Array{Float64}(undef, 0, 0, 0))

"""
    calc_overlap!(cc::conf_coord; cut_off = 0.25, σ = 0.025)

Calculate phonon overlap between phonon wave functions 'potenrial.χ'.
If energy difference is larger then the cutoff (eV) `abs(cc.V1.ϵ[i] - cc.V2.ϵ[j]) > cut_off`,
the overlap will not be calculated.
Delta functions are replaced by a Gaussians function with widths `σ`.
"""
function calc_overlap!(cc::conf_coord; cut_off = 0.25, σ = 0.025)
    Q₀ = cc.V1.Q0
    ΔL = (maximum(cc.V1.Q) - minimum(cc.V1.Q))/length(cc.V1.Q)
    cc.overlap_matrix = zeros(length(cc.V1.ϵ), length(cc.V2.ϵ))
    cc.δ_matrix = zeros(length(cc.V1.ϵ), length(cc.V2.ϵ))

    for i in UnitRange(1, length(cc.V1.ϵ))
        for j in UnitRange(1, length(cc.V2.ϵ))
            Δϵ = abs(cc.V1.ϵ[i] - cc.V2.ϵ[j])
            if  Δϵ < cut_off
                integrand = (cc.V1.χ[i, :] .* (cc.V1.Q .- Q₀) .* cc.V2.χ[j, :])
                overlap = sum(integrand)*ΔL

                cc.overlap_matrix[i, j] = overlap
                cc.δ_matrix[i, j] = exp(-(Δϵ/σ)^2/2)/(σ*sqrt(2π))
            end
        end
    end
end

# calculating carrier capture rate in units cm³/s for a given:
# V: 	  volume of supercell [cm³]
# temperature: range from Tmin to Tmax for NT steps
"""
    calc_capt_coeff!(cc::conf_coord, V::Float64, temperature)

Calculte the capture coefficient `cc.capt_coeff` as a function of `temperature` which is a `UnitRange`.
`V` is a volume where the electron-phonon coupling matrix element `cc.W` is calculated.
The lowest thermal occupation number of the eigenstate must be lower than `occ_cut_off = 1E-5`.

    @assert occ_high < occ_cut_off "occ(ϵ_max, T_max): \$occ_high should be less than $occ_cut_off"

"""
function calc_capt_coeff!(cc::conf_coord, V::Float64, temperature)
    # TODO:       convergence over σ
    partial_capt_coeff = zeros(length(cc.V1.ϵ), length(cc.V2.ϵ), length(temperature))
    cc.capt_coeff = zeros(length(temperature))
    Z = zeros(length(temperature)) # partition function
    β = 1 ./ (kB * temperature)

    for ϵ in cc.V1.ϵ
        Z = Z + exp.(-β*ϵ)
    end

    for i in 1:length(cc.V1.ϵ)
        for j in 1:length(cc.V2.ϵ)
            ϵ = cc.V1.ϵ[i]
            overlap = cc.overlap_matrix[i, j]
            δ = cc.δ_matrix[i, j]
            occ = exp.(-β*ϵ) ./ Z
            partial_capt_coeff[i, j, :] = occ * overlap .* overlap * δ
        end
    end
    # term in capture rate expression prior to summation
    partial_capt_coeff = V*2π/ħ*cc.g*cc.W^2 * partial_capt_coeff

    replace!(partial_capt_coeff, NaN => 0)

    occ_high = exp(-β[end]*cc.V1.ϵ[end]) / Z[end]

    if isnan(occ_high)
        occ_high = 0
    end

    @assert occ_high < occ_cut_off "occ(ϵ_max, T_max): $occ_high should be less than $occ_cut_off"

    cc.capt_coeff = dropdims(sum(partial_capt_coeff, dims = (1, 2)), dims = (1, 2))
    replace!(cc.capt_coeff, 0=>1E-127)
    cc.partial_capt_coeff = partial_capt_coeff
    cc.temperature = temperature
end

# importing parameters
"""
Depreciated.
Construct `conf_coord` from two potentials `pot_i` (initila) and 'pot_f' (final) and configure dictionalry `cfg`.
"""
function cc_from_dict(pot_i, pot_f, cfg::Dict)::conf_coord
    cc = conf_coord(pot_i, pot_f)
    cc.name = "$(cfg["initial"]) => $(cfg["final"])"
    cc.W = cfg["W"]
    cc.g = cfg["g"]
    return cc
end


text:
const CHAIN_INDEX_DOC = """
`chain_index::Vector`: An optional vector of integers specifying which chain each
    step belongs to. For instance, `chain_index[step]` should return `2` if
    `log_likelihood[:, step]` belongs to the second chain.
"""

const LIKELIHOOD_FUNCTION_ARG = """
`ll_fun::Function`: A function taking a single data point and returning the log-likelihood
of that point. This function must take the form `f(θ[1], ..., θ[n], data)`, where `θ` is the
parameter vector. See also the `splat` keyword argument.
"""

const R_EFF_DOC = """
`r_eff::AbstractArray{T}`: An (optional) vector of relative effective sample sizes used
in ESS calculations. If left empty, calculated automatically using the FFTESS method
from InferenceDiagnostics.jl. See `relative_eff` to calculate these values.
"""

const ARGS = """
`args...`: Positional arguments to be passed to
"""

const KWARGS = """
`kwargs...`: Keyword arguments to be passed to
"""

###############
## FUNCTIONS ##
###############

"""
Assume that all objects belong to a single chain if chain index is missing. Inform user.
"""
function _assume_one_chain(matrix)
    @info "Chain information was not provided; " *
          "all samples are assumed to be drawn from a single chain."
    return ones(length(matrix))
end

"""
Convert a matrix+chain_index representation to a 3d array representation.
"""
function _convert_to_array(matrix::AbstractMatrix, chain_index::AbstractVector)
    indices = unique(chain_index)
    biggest_idx = maximum(indices)
    dims = size(matrix)
    if dims[2] ≠ length(chain_index)
        throw(ArgumentError("Some entries do not have a chain index."))
    elseif !issetequal(indices, 1:biggest_idx)
        throw(
            ArgumentError(
                "Indices must be numbered from 1 through the total number of chains."
            ),
        )
    else
        # Check how many elements are in each chain, assign to "counts"
        counts = count.(eachslice(chain_index .== indices'; dims=2))
        # check if all inputs are the same length
        if !all(==(counts[1]), counts)
            throw(ArgumentError("All chains must be of equal length."))
        end
    end
    new_ratios = similar(matrix, dims[1], dims[2] ÷ biggest_idx, biggest_idx)
    for i in 1:biggest_idx
        new_ratios[:, :, i] .= matrix[:, chain_index .== i]
    end
    return new_ratios
end


text:
##################################################
# Helper functions
##################################################

"""
    spline_catmullrom(u, y0, y1, y2, y3)

Smooth interpolation between `y1, y2`, using the Catmull-Rom spline.
The function uses `y0` and `y3` to ensure that the derivative at `y1`
and `y2` is continuous.  `u` is a parameter between `0` and `1`.

"""
spline_catmullrom(u, y0, y1, y2, y3) = y1 + 0.5u*(-y0 + y2 + u*(2y0 - 5y1 + 4y2 - y3 + u*(-y0 + 3y1 - 3y2 + y3)))

_orbitalrank(orbital_symbol) =
    startswith(string(orbital_symbol), "s") ? 1 :
    startswith(string(orbital_symbol), "p") ? 2 :
    startswith(string(orbital_symbol), "d") ? 3 :
    startswith(string(orbital_symbol), "f") ? 4 :
    error("invalid orbital symbol: $orbital_symbol")

##################################################


text:
"""
    $(TYPEDEF)

A type to represent memory-mapped TIFF data. Useful for opening and operating on
images too large to store in memory.

$(FIELDS)
"""
mutable struct DiskTaggedImage{T <: Colorant, O <: Unsigned, AA <: AbstractArray} <: AbstractDenseTIFF{T, 3}

    """
    Pointer to keep track of the backing file
    """
    file::TiffFile{O}
    """
    The associated tags for each slice in this array
    """
    ifds::Vector{IFD{O}}

    dims::NTuple{3, Int}

    """
    An internal cache to fill reading from disk
    """
    cache::AA

    """
    The index of the currently loaded slice
    """
    cache_index::Int

    """
    A flag tracking whether this file is editable
    """
    readonly::Bool

    function DiskTaggedImage(file::TiffFile{O}, ifds::Vector{IFD{O}}) where {O}
        ifd = ifds[1]
        dims = (nrows(ifd), ncols(ifd), length(ifds))
        cache = getcache(ifd)
        new{eltype(cache), O, typeof(cache)}(file, ifds, dims, cache, -1, false)
    end
end

Base.size(A::DiskTaggedImage) = A.dims

function Base.getindex(A::DiskTaggedImage{T, O, AA}, i1::Int, i2::Int, i::Int) where {T, O, AA}
    # check the loaded cache is already the correct slice
    if A.cache_index == i
        return A.cache[i2, i1]
    end

    ifd = A.ifds[i]

    # if the file isn't open, lets open a handle and update it
    if !isopen(A.file.io)
        path = A.file.filepath
        A.file.io = Stream(format"TIFF", open(path), path)
    end

    read!(A.cache, A.file, ifd)

    A.cache_index = i

    return A.cache[i2, i1]
end

function Base.setindex!(A::DiskTaggedImage, I...)
    error("This array is on disk and is read only. Convert to a mutable in-memory version by running "*
          "`copy(arr)`. \n\n𝗡𝗼𝘁𝗲: For large files this can be quite expensive. A future PR will add "*
          "support for reading and writing to/from disk.")
end


text:
include("hrb_utils.jl")

type RBMeta
  # Metadata attached to Rutherford-Boeing data.
  title :: AbstractString
  key   :: AbstractString

  mxtype :: AbstractString
  nrow :: Int
  ncol :: Int
  nnzero :: Int
  neltvl :: Int

  hermitian :: Bool
  assembled :: Bool
  pattern_only :: Bool

  ptrfmt :: AbstractString
  indfmt :: AbstractString
  valfmt :: AbstractString

  dattyp :: AbstractString
  positn :: Char
  orgniz :: Char
  caseid :: AbstractString
  numerf :: Char
  auxfm1 :: AbstractString
  auxfm2 :: AbstractString
  auxfm3 :: AbstractString
end

@compat RBDataType = Union{Array{Int,2}, Array{Float64,2}, Array{Complex128,2}, SparseMatrixCSC}

type RutherfordBoeingData

  meta :: RBMeta
  data :: RBDataType

  function RutherfordBoeingData(file_name :: AbstractString)
    rb = open(file_name)

    # Read header.
    seekstart(rb)
    line  = readline(rb)
    title = strip(line[1:72])      # A72
    key   = strip(line[73:end-1])  # A8
    buffer1 = readline(rb)         # A80
    buffer2 = readline(rb)         # A80

    if buffer2[3] in ['a', 'e']
      # Read a matrix.
      mxtype = buffer2[1:3]
      nrow, ncol, nnzero, neltvl = map(s -> parse(Int, s),
                                       split(chomp(buffer2[4:end])))

      pattern_only = (mxtype[1] in ['p', 'q'])
      if pattern_only
        ptrfmt, indfmt = split(strip(readline(rb)))
      else
        ptrfmt, indfmt, valfmt = split(strip(readline(rb)))
      end

      # Read integer data.
      np1 = (mxtype[2:3] == "re") ? (2*ncol + 1) : (ncol + 1)
      ip  = read_array(rb, np1, ptrfmt)
      ind = read_array(rb, nnzero, indfmt)

      nreal = neltvl > 0 ? neltvl : nnzero

      if pattern_only
        vals = ones(Float64, nreal)  # To be able to build a SparseMatrixCSC
      else
        vals = read_array(rb, nreal, valfmt)
      end
      # Ensure row indices are sorted in each column.
      sortsparse!(ip, ind, vals)
      data = SparseMatrixCSC(nrow, ncol, ip, ind, vals)

      meta = RBMeta(title, key, mxtype, nrow, ncol, nnzero, neltvl,
                    mxtype[2] == 's', mxtype[3] == 'a', pattern_only,
                    ptrfmt, indfmt, pattern_only ? "" : valfmt,
                    "", '\0', '\0', "", '\0', "", "", "")

    else

      # Read supplementary data.
      types_with_values = ["rhs", "sln", "est", "ipt", "icv"]

      dattyp = buffer1[1:3]
      positn = buffer1[4]
      orgniz = buffer1[5]
      caseid = buffer1[7:14]
      numerf = buffer1[16]
      intvals = map(s -> parse(Int, s), split(chomp(buffer1[17:end])))
      auxfmts = split(strip(buffer2))
      auxfm1 = auxfm2 = auxfm3 = ""
      if dattyp in types_with_values
        (nrow, ncol, nnzero) = intvals
      else
        (nrow, ncol) = intvals
      end
      if dattyp == "ord"
        nnzero = nrow * ncol
      end

      if (dattyp == "rhs" && orgniz == 's') || (dattyp in ["ord", "ipt", "icv"])
        # Read indices.
        auxfm = auxfmts[1]
        if (dattyp == "rhs" && orgniz == 's')
          auxfm1 = auxfmts[1]
          ip = read_array(rb, ncol+1, auxfm1)
          auxfm = auxfmts[2]
        end
        ind = read_array(rb, nnzero, auxfm)
      end

      if dattyp in types_with_values
        # Read values.
        if dattyp != "rhs"
          nnzero = nrow * ncol
        end
        if dattyp == "rhs" && orgniz == 's'
          auxfm = auxfmts[3]
        else
          auxfm = auxfmts[1]
        end

        vals = read_array(rb, nnzero, auxfm)

        if (dattyp == "rhs" && orgniz == 's') || (dattyp in ["ipt", "icv"])
          # Ensure row indices are sorted in each column.
          sortsparse!(ip, ind, vals)
          data = SparseMatrixCSC(nrow, ncol, ip, ind, vals)
        else
          data = reshape(vals, (nrow, ncol))
        end
      else
        data = reshape(ind, (nrow, ncol))
      end

      meta = RBMeta(title, key, "", nrow, ncol, nnzero, 0,
                    false, orgniz != 'e', false,
                    "", "", "",
                    dattyp, positn, orgniz, caseid, numerf,
                    auxfm1, auxfm2, auxfm3)
    end

    close(rb)

    new(meta, data)
  end
end

# Displaying Rutherford-Boeing data instances.

import Base.show, Base.print
function show(io :: IO, rb :: RutherfordBoeingData)
  typ = (rb.meta.mxtype == "") ? rb.meta.dattyp : rb.meta.mxtype
  s  = @sprintf("Rutherford-Boeing data %s of type %s\n", rb.meta.key, typ)
  s *= @sprintf("%d rows, %d cols, %d nonzeros\n", rb.meta.nrow, rb.meta.ncol, rb.meta.nnzero)
  print(io, s)
end

function print(io :: IO, rb :: RutherfordBoeingData)
  @printf("Rutherford-Boeing data %s\n", rb.meta.key)
  @printf("%s\n", rb.meta.title)
  @printf("%d rows, %d cols, %d nonzeros\n", rb.meta.nrow, rb.meta.ncol, rb.meta.nnzero)
  if rb.meta.mxtype != ""
    if rb.meta.mxtype[1] == 'P'
      dtype = "pattern only"
    elseif rb.meta.mxtype[1] == 'R'
      dtype = "real"
    else
      dtype = "complex"
    end
    herm = rb.meta.hermitian ? "hermitian" : "non-hermitian"
    assm = rb.meta.assembled ? "assembled" : "elemental"
    @printf("(%s, %s, %s)\n", dtype, herm, assm)
  end
  display(rb.data)
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    sssr, atom_sssrmem, bond_sssrmem,
    atom_sssrsizes, bond_sssrsizes,
    isringatom, isringbond,
    atom_sssrcount, bond_sssrcount,
    scaffoldmem, componentmem,
    atomsymbol, charge, multiplicity, bondorder,
    nodedegree, valence, lonepair, heavyatomcount,
    explicithcount, implicithcount, hcount, connectivity,
    ishdonor, hacceptorcount, ishacceptor, hdonorcount,
    isrotatable, rotatablecount,
    countatoms, molecularformula, empiricalformula,
    pielectron, hybridization,
    isaromaticring, isaromatic, isaromaticbond

# Molecular graph topology

sssr = mincyclenodes
atom_sssrmem = node_cyclemem
bond_sssrmem = edge_cyclemem
scaffoldmem = two_edge_membership
componentmem = connected_membership

function atom_sssrsizes(mol::UndirectedGraph)
    vec = [Set{Int}() for i in 1:nodecount(mol)]
    rings = sssr(mol)
    for (i, cs) in enumerate(atom_sssrmem(mol))
        for c in cs
            push!(vec[i], length(rings[c]))
        end
    end
    return vec
end

atom_sssrsizes(view::SubgraphView) = atom_sssrsizes(view.graph)

function bond_sssrsizes(mol::UndirectedGraph)
    vec = [Set{Int}() for i in 1:edgecount(mol)]
    rings = sssr(mol)
    for (i, cs) in enumerate(bond_sssrmem(mol))
        for c in cs
            push!(vec[i], length(rings[c]))
        end
    end
    return vec
end

bond_sssrsizes(view::SubgraphView) = bond_sssrsizes(view.graph)

# TODO: waiting for fix #28992
# @cache atom_isringmem(mol) = .!isempty.(node_cyclemem(mol))
isringatom(mol::UndirectedGraph) = [!i for i in isempty.(atom_sssrmem(mol))]
isringatom(view::SubgraphView) = isringatom(view.graph)

# @cache bond_isringmem(mol) = .!isempty.(node_cyclemem(mol))
isringbond(mol::UndirectedGraph) = [!i for i in isempty.(bond_sssrmem(mol))]
isringbond(view::SubgraphView) = isringbond(view.graph)

atom_sssrcount(mol::UndirectedGraph) = length.(atom_sssrmem(mol))
atom_sssrcount(view::SubgraphView) = atom_sssrcount(view.graph)

bond_sssrcount(mol::UndirectedGraph) = length.(bond_sssrmem(mol))
bond_sssrcount(view::SubgraphView) = bond_sssrcount(view.graph)

# Elemental properties

"""
    atomsymbol(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom symbols within the molecule that have ``n``
atoms.
"""
@cache atomsymbol(mol::GraphMol) = getproperty.(nodeattrs(mol), :symbol)
atomsymbol(view::SubgraphView) = atomsymbol(view.graph)

"""
    charge(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom charges within the molecule that have ``n``
atoms.
"""
@cache charge(mol::GraphMol) = getproperty.(nodeattrs(mol), :charge)
charge(view::SubgraphView) = charge(view.graph)

"""
    multiplicity(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of atom multiplicities within the molecule that
have ``n`` atoms (1: non-radical, 2: radical, 3: biradical).
"""
@cache multiplicity(mol::GraphMol) = getproperty.(nodeattrs(mol), :multiplicity)
multiplicity(view::SubgraphView) = multiplicity(view.graph)

"""
    bondorder(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of bond orders within the molecule that have ``n``
bonds.
"""
@cache bondorder(mol::GraphMol) = getproperty.(edgeattrs(mol), :order)
bondorder(view::SubgraphView) = bondorder(view.graph)

"""
    nodedegree(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of node degrees within the molecular graph that
have ``n`` atom nodes.

`nodedegree` values correspond to SMARTS `D` property.
"""
@cache nodedegree(mol::GraphMol) = [degree(mol, n) for n in 1:nodecount(mol)]
nodedegree(view::SubgraphView) = nodedegree(view.graph)

# Valence

"""
    lonepair(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size ``n`` vector of the number of lone pairs within the molecule that
have ``n`` atoms.

Note that implicit hydrogens are available for only organic atoms. The lonepair
value of inorganic atoms would be `nothing`. The result can take negative value if the atom has empty valence shells.
"""
function lonepair(mol::GraphMol)
    defs = Dict(
        :H => 0, :B => -1, :C => 0, :N => 1, :O => 2, :F => 3,
        :Si => 0, :P => 1, :S => 2, :Cl => 3,
        :As => 1, :Se => 2, :Br => 3, :I => 3
    )
    vec = Union{Int,Nothing}[]
    for (sym, chg) in zip(atomsymbol(mol), charge(mol))
        num = get(defs, sym, nothing)
        v = num === nothing ? nothing : num - chg
        push!(vec, v)
    end
    return vec
end

lonepair(view::SubgraphView) = lonepair(view.graph)

function apparentvalence(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    bondorder_ = bondorder(mol)
    for i in 1:nodecount(mol)
        for inc in incidences(mol, i)
            vec[i] += bondorder_[inc]
        end
    end
    return vec
end

apparentvalence(view::SubgraphView) = apparentvalence(view.graph)

"""
    valence(mol::GraphMol) -> Vector{Union{Int,Nothing}}

Return the size ``n`` vector of intrinsic valences (with considering implicit
hydrogens) within the molecule that have ``n`` atoms.

`valence` values correspond to SMARTS `v` property. `valence` values of inorganic atoms would be `nothing`.
"""
@cache function valence(mol::GraphMol)
    vec = Union{Int,Nothing}[]
    atomsymbol_ = atomsymbol(mol)
    lonepair_ = lonepair(mol)
    for i in 1:nodecount(mol)
        if lonepair_[i] === nothing
            push!(vec, nothing)
        elseif atomsymbol_[i] === :H
            push!(vec, 1)
        else
            push!(vec, 4 - abs(lonepair_[i]))
        end
    end
    return vec
end

valence(view::SubgraphView) = valence(view.graph)

"""
    explicithcount(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of adjacent explicit hydrogen nodes within the molecule that have ``n`` atoms.
"""
@cache function explicithcount(mol::GraphMol)
    vec = zeros(Int, nodecount(mol))
    atomsymbol_ = atomsymbol(mol)
    for i in 1:nodecount(mol)
        for adj in adjacencies(mol, i)
            if atomsymbol_[adj] == :H
                vec[i] += 1
            end
        end
    end
    return vec
end

explicithcount(view::SubgraphView) = explicithcount(view.graph)

"""
    implicithcount(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of implicit hydrogens within
the molecule that have ``n`` atoms.
"""
@cache function implicithcount(mol::GraphMol)
    hcnt = (v, av) -> v === nothing ? 0 : max(0, v - av)
    return hcnt.(valence(mol), apparentvalence(mol))
end

implicithcount(view::SubgraphView) = implicithcount(view.graph)

"""
    heavyatomcount(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of adjacent non-hydrogen atoms within
the molecule that have ``n`` atoms.
"""
heavyatomcount(mol::GraphMol) = nodedegree(mol) - explicithcount(mol)
heavyatomcount(view::SubgraphView) = heavyatomcount(view.graph)

"""
    hcount(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the total number of hydrogens attached to the atom within the molecule that have ``n`` atoms.

`hcount` values correspond to SMARTS `H` property.
"""
@cache hcount(mol::GraphMol) = explicithcount(mol) + implicithcount(mol)
hcount(view::SubgraphView) = hcount(view.graph)

"""
    connectivity(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the total number of adjacent atoms (including implicit hydrogens) within the molecule that have ``n`` atoms.

`connectivity` values correspond to SMARTS `X` property.
"""
@cache connectivity(mol::GraphMol) = nodedegree(mol) + implicithcount(mol)
connectivity(view::SubgraphView) = connectivity(view.graph)

# Hydrogen bond donor/acceptor

function ishacceptor(mol::GraphMol)
    ac = (sym, lp) -> lp === nothing ? false : sym in (:N, :O, :F) && lp > 0
    return ac.(atomsymbol(mol), lonepair(mol))
end

ishacceptor(view::SubgraphView) = ishacceptor(view.graph)

"""
    hacceptorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond acceptors (N, O and F).
"""
hacceptorcount(mol::GraphMol) = reduce(+, ishacceptor(mol); init=0)

function ishdonor(mol::GraphMol)
    dc = (sym, h) -> sym in (:N, :O) && h > 0
    return dc.(atomsymbol(mol), hcount(mol))
end

ishdonor(view::SubgraphView) = ishdonor(view.graph)

"""
    hdonorcount(mol::GraphMol) -> Int

Return the number of hydrogen bond donors (O and N attached to hydrogens).
"""
hdonorcount(mol::GraphMol) = reduce(+, ishdonor(mol); init=0)

# Rotatable bonds

"""
    isrotatable(mol::GraphMol)

Return whether the bonds are rotatable or not.
"""
@cache function isrotatable(mol::GraphMol)
    nodedegree_ = nodedegree(mol)
    isringbond_ = isringbond(mol)
    bondorder_ = bondorder(mol)
    vec = Bool[]
    for (i, (u, v)) in enumerate(edgesiter(mol))
        rot = (!isringbond_[i] && bondorder_[i] == 1
            && nodedegree_[u] != 1 && nodedegree_[v] != 1)
        push!(vec, rot)
    end
    return vec
end

isrotatable(view::SubgraphView) = isrotatable(view.graph)

"""
    rotatablecount(mol::GraphMol) -> Int

Return the number of rotatable bonds.
"""
rotatablecount(mol::GraphMol) = reduce(+, isrotatable(mol); init=0)

# Composition

"""
    countatoms(mol::GraphMol) -> Dict{Symbol,Int}

Count the number of atoms and return symbol => count dict.
"""
function countatoms(mol::GraphMol)
    counter = Dict{Symbol,Int}()
    for sym in atomsymbol(mol)
        if !haskey(counter, sym)
            counter[sym] = 1
        else
            counter[sym] += 1
        end
    end
    hcnt = reduce(+, implicithcount(mol); init=0)
    if hcnt > 0
        if !haskey(counter, :H)
            counter[:H] = hcnt
        else
            counter[:H] += hcnt
        end
    end
    return counter
end

function writeformula(counter::Dict{Symbol,Int})
    strs = []
    if haskey(counter, :C)
        push!(strs, "C")
        c = pop!(counter, :C)
        c > 1 && push!(strs, string(c))
        if haskey(counter, :H)
            push!(strs, "H")
            h = pop!(counter, :H)
            h > 1 && push!(strs, string(h))
        end
    end
    for sym in sort(collect(keys(counter)))
        push!(strs, string(sym))
        cnt = counter[sym]
        cnt > 1 && push!(strs, string(cnt))
    end
    return join(strs)
end

"""
    molecularformula(mol::GraphMol) -> String

Return the molecular formula in Hill system.
"""
function molecularformula(mol::GraphMol)
    counter = countatoms(mol)
    return writeformula(counter)
end

"""
    empiricalformula(mol::GraphMol) -> String

Return the empirical formula in Hill system.
"""
function empiricalformula(mol::GraphMol)
    counter = countatoms(mol)
    if length(counter) > 1
        dv = reduce(gcd, values(counter))
        for k in keys(counter)
            counter[k] = div(counter[k], dv)
        end
    end
    return writeformula(counter)
end

# Hybridization

"""
    pielectron(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of the number of ``\\pi`` electrons (including implicit hydrogens) within the molecule that have ``n`` atoms.

The ``\\pi`` electron enumelation is based on the following simple rules.

- Any atom incident to a double bond -> +1
- Any atom incident to two double bond -> +2
- Any atom incident to a triple bond -> +1
- Uncharged N and O not incident to any multiple bonds and adjacent to atoms described above -> +2

These rules are applied for only typical organic atoms. The values for inorganic atoms will be 0.
"""
@cache function pielectron(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    charge_ = charge(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    vec = zeros(Int, nodecount(mol))
    for i in 1:nodecount(mol)
        vec[i] = pie_[i]
        atomsymbol_[i] in (:N, :O) || continue
        for adj in adjacencies(mol, i)
            if pie_[i] == 0 && pie_[adj] > 0 && charge_[i] == 0
                vec[i] = 2
                break
            end
        end
    end
    return vec
end

pielectron(view::SubgraphView) = pielectron(view.graph)

"""
    hybridization(mol::GraphMol) -> Vector{Int}

Return the size ``n`` vector of orbital hybridization symbols (`:sp3`, `:sp2`, `:sp`, `:none`) within the molecule that have ``n`` atoms.

These hybridizations are for only typical organic atoms. Inorganic atoms and other orbitals like s, sp3d and sp3d2 will be `:none`.
"""
@cache function hybridization(mol::GraphMol)
    atomsymbol_ = atomsymbol(mol)
    connectivity_ = connectivity(mol)
    lonepair_ = lonepair(mol)
    pielectron_ = pielectron(mol)
    organic_symbols = [:B, :C, :N, :O, :F, :Si, :P, :S, :Cl, :As, :Se, :Br, :I]
    vec = Symbol[]
    for i in 1:nodecount(mol)
        if atomsymbol_[i] in organic_symbols
            val = connectivity_[i] + lonepair_[i]
            if val == 4
                if (atomsymbol_[i] in [:N, :O] && pielectron_[i] == 2)
                    push!(vec, :sp2)  # adjacent to conjugated bonds
                else
                    push!(vec, :sp3)
                end
            elseif val == 3
                push!(vec, :sp2)
            elseif val == 2
                push!(vec, :sp)
            else
                push!(vec, :none)
            end
        else
            push!(vec, :none)
        end
    end
    return vec
end

hybridization(view::SubgraphView) = hybridization(view.graph)

# Aromatic rings

function satisfyhuckel(mol::GraphMol, ring)
    # TODO: need improvement
    cnt = 0
    atomsymbol_ = atomsymbol(mol)
    nodedegree_ = nodedegree(mol)
    pie_ = apparentvalence(mol) - nodedegree(mol)
    lonepair_ = lonepair(mol)
    carbonylO = findall(
        (atomsymbol_ .== :O) .* (nodedegree_ .== 1) .* (pie_ .== 1))
    carbonylC = Int[]
    for o in carbonylO
        c = iterate(adjacencies(mol, o))[1]
        if atomsymbol_[c] == :C
            push!(carbonylC, c)
        end
    end
    for r in ring
        if r in carbonylC
            continue
        elseif pie_[r] == 1
            cnt += 1
        elseif lonepair_[r] === nothing
            return false
        elseif lonepair_[r] > 0
            cnt += 2
        elseif lonepair_[r] < 0
            continue
        else
            return false
        end
    end
    return cnt % 4 == 2
end

@cache function isaromaticring(mol::GraphMol)
    vec = falses(circuitrank(mol))
    for (i, ring) in enumerate(sssr(mol))
        if satisfyhuckel(mol, ring)
            vec[i] = true
        elseif nodeattrtype(mol) === SmilesAtom # SMILES aromatic atom
            sub = nodesubgraph(mol, Set(ring))
            if all(nodeattr(mol, n).isaromatic for n in nodeset(sub))
                vec[i] = true
            end
        end
    end
    return vec
end

"""
    isaromatic(mol::GraphMol)

Return the vector whether the atom belongs to an aromatic ring.

Note that aromaticity described here means simplified binary descriptor
(aromatic or not) based on classical Huckel's rule. This is intended for use in
some kind of pharmaceutical researches. Non-classical aromaticity such as
Moebius aromaticity is not considered.
"""
@cache function isaromatic(mol::GraphMol)
    aromatic = falses(nodecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for n in nodeset(sub)
            aromatic[n] = true
        end
    end
    return aromatic
end

isaromatic(view::SubgraphView) = isaromatic(view.graph)

@cache function isaromaticbond(mol::GraphMol)
    aromaticbond = falses(edgecount(mol))
    for ring in sssr(mol)[isaromaticring(mol)]
        sub = nodesubgraph(mol, Set(ring))
        for e in edgeset(sub)
            aromaticbond[e] = true
        end
    end
    return aromaticbond
end

isaromaticbond(view::SubgraphView) = isaromaticbond(view.graph)


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

Used by/in unclaimed high-level pnml labels.
Is an abstract syntax tree (ast) expressed in XML.
Note the structural similarity to [`PnmlLabel`](@ref) and [`AnyElement`](@ref)

# Extra
There are various defined structure ast variants:
  - Sort of a Place type [builtin, multi, product, user]
  - Term of Place HLMarking  [variable, operator]
  - Term of Transition Condition  [variable, operator]
  - Term of Arc Inscription [variable, operator]
  - Declarations of Declaration * [sort, variable, operator]
"""
struct Structure{T} #TODO
    tag::Symbol
    dict::T
    xml::XMLNode
end
Structure(p::Pair{Symbol,PnmlDict}, node::XMLNode) = Structure(p.first, p.second, node)

tag(s::Structure) = s.tag
dict(s::Structure) = s.dict
xmlnode(s::Structure) = s.xml


text:
@doc raw"""
`TableauVPRK`: Tableau of a Variational Partitioned Runge-Kutta method
```math
\begin{aligned}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , &
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, V_{n,j} , &
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, V_{n,i} , \\
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , &
P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, F_{n,j} - d_i \lambda , &
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, F_{n,i} , \\
&&
0 &= \sum \limits_{i=1}^{s} d_i V_i , &&
\end{aligned}
```
satisfying the symplecticity conditions
```math
\begin{aligned}
b_{i} \bar{a}_{ij} + b_{j} a_{ji} &= b_{i} b_{j} , &
\bar{b}_i &= b_i .
\end{aligned}
```
"""
struct TableauVPRK{T} <: AbstractTableau{T}
    @HeaderTableau

    q::Tableau{T}
    p::Tableau{T}

    R∞::Int

    d::Vector{T}

    function TableauVPRK{T}(name, o, q, p, R∞, d) where {T}
        @assert q.s == p.s == length(d)
        new(name, o, q.s, q, p, R∞, d)
    end

    function TableauVPRK{T}(name, o, q, p, R∞) where {T}
        @assert q.s == p.s
        new(name, o, q.s, q, p, R∞)
    end
end

function TableauVPRK(name::Symbol, order::Int, q::Tableau{T}, p::Tableau{T}, R∞::Int, d::Vector{T}) where {T}
    TableauVPRK{T}(name, order, q, p, R∞, d)
end

function TableauVPRK(name::Symbol, order::Int, q::Tableau{T}, p::Tableau{T}, R∞::Int) where {T}
    TableauVPRK{T}(name, order, q, p, R∞)
end

function TableauVPRK(name::Symbol, order::Int, q::Tableau{T}, R∞::Int, d::Vector{T}) where {T}
    TableauVPRK{T}(name, order, q, get_symplectic_conjugate_coefficients(q), R∞, d)
end

function TableauVPRK(name::Symbol, order::Int, q::Tableau{T}, R∞::Int) where {T}
    TableauVPRK{T}(name, order, q, get_symplectic_conjugate_coefficients(q), R∞)
end

Base.:(==)(tab1::TableauVPRK, tab2::TableauVPRK) = (tab1.o == tab2.o
                                                 && tab1.s == tab2.s
                                                 && tab1.q == tab2.q
                                                 && tab1.p == tab2.p
                                                 && tab1.R∞ == tab2.R∞)


text:
using Libdl

"""
    validate_libR(libR, raise=true)

Checks that the R library `libR` can be loaded and is satisfies version requirements.

If `raise` is set to `false`, then returns a boolean indicating success rather
than throwing exceptions.
"""
function validate_libR(libR, raise=true)
    if !isfile(libR)
        raise || return false
        error("Could not find library $libR. Make sure that R shared library exists.")
    end
    # Issue #143
    # On linux, sometimes libraries linked from libR (e.g. libRblas.so) won't open unless LD_LIBRARY_PATH is set correctly.
    libptr = try
        Libdl.dlopen(libR)
    catch er
        raise || return false
        Base.with_output_color(:red, stderr) do io
            print(io, "ERROR: ")
            showerror(io, er)
            println(io)
        end
        if Sys.iswindows()
            error("Try adding $(dirname(libR)) to the \"PATH\" environmental variable and restarting Julia.")
        else
            error("Try adding $(dirname(libR)) to the \"LD_LIBRARY_PATH\" environmental variable and restarting Julia.")
        end
    end
    # R_tryCatchError is only available on v3.4.0 or later.
    if Libdl.dlsym_e(libptr, "R_tryCatchError") == C_NULL
        msg = "R library $libR appears to be too old. RCall.jl requires R 3.4.0 or later."
        if raise
            error(msg)
        else
            @info msg
            return false
        end
    end
    Libdl.dlclose(libptr)
    return true
end

function locate_libR(Rhome, raise=true)
    @static if Sys.iswindows()
        libR = joinpath(Rhome, "bin", Sys.WORD_SIZE==64 ? "x64" : "i386", "R.dll")
    else
        libR = joinpath(Rhome, "lib", "libR.$(Libdl.dlext)")
    end
    return validate_libR(libR, raise) ? libR : ""
end


text:
################################################################################
#
#          EllCrv/EllCrvMap.jl : Abstract Variety Maps,
#               Elliptic Curve Maps, and
#               Elliptic Curve Blowups at subvarieties.
#
################################################################################

################################################################################
#
#  Submodules
#
################################################################################

include("EllCrv.jl")

################################################################################
#
#  Exports
#
################################################################################

export AbstractVarietyMap, AbstractCurveMap, EllCrvBlowup, EllCrvMap

################################################################################
#
#  Abstract Types/Structs
#
################################################################################

Map(::Type{T}) where T <: AbstractAlgebra.Map = AbstractAlgebra.Map(T)

mutable struct AbstractVarietyMap{T} <: AbstractAlgebra.Map{AbstractVariety, AbstractVariety}
    source::AbstractVariety
    target::AbstractVariety

    function AbstractVarietyMap(X::AbstractVariety, Y::AbstractVariety)
        z = new(X, Y)
        return z
    end
end

mutable struct AbstractCurveMap <: AbstractVarietyMap
    source::AbstractCurve
    target::AbstractCurve

    function AbstractCurveMap(X::AbstractCurve, Y::AbstractCurve)
        z = new(source, target)
        return z
    end
end

mutable struct EllCrvBlowup <: AbstractVarietyMap
    source::AbstractVariety
    target::EllCrv

    function EllCrvBlowup(Z::AbstractVariety, Y::EllCrv)
        z = new(source, target)
        return z
    end
end

mutable struct EllCrvMap <: AbstractCurveMap
    source::EllCrv
    target::EllCrv

    function EllCrvMap(E1::EllCrv, E2::EllCrv)
        z = new(source, target)
        return z
    end
end


text:
# df = CSV.File("top5000-lemmas.txt")|> DataFrame;
# top5000 = df[!, :lemma];
# top5000nv = df[(df[!, :PoS] .== "v") .| (df[!, :PoS] .== "n"), :lemma];

# top5000[endswith.(top5000, "s").&(.!endswith.(top5000, "ss"))] + "has", "is", "was"
const s_ending_words = Set(["numerous", "serious", "towards", "clothes", "focus", "bonus", "sometimes", "suspicious", "his",
        "species", "thanks", "hers", "as", "aids", "proceedings", "analysis", "bus", "us", "gas", "lots", "versus",
        "emphasis", "nervous", "upstairs", "dangerous", "curious", "delicious", "anonymous", "various", "yours",
        "lens", "physics", "dynamics", "indigenous", "besides", "anxious", "always", "works", "ours", "news", "bias",
        "gorgeous", "vs", "diabetes", "virus", "genius", "consensus", "hypothesis", "ambitious", "whereas", "jealous",
        "crisis", "tremendous", "economics", "canvas", "odds", "obvious", "jeans", "stimulus", "chaos", "corps",
        "mathematics", "statistics", "precious", "basis", "ridiculous", "perhaps", "headquarters", "lyrics", "yes",
        "tennis", "continuous", "diagnosis", "plus", "its", "earnings", "previous", "mysterious", "christmas",
        "politics", "ourselves", "themselves", "thus", "conscious", "sales", "religious", "generous", "enormous",
        "this", "status", "series", "olympics", "famous", "campus", "census", "ethics", "terms", "has", "is", "was"])

# top5000nv[endswith.(top5000nv, r"xe|che|she|oe|ie")]
const xe_ending_words = Set(["tie", "lie", "rookie", "pie", "cookie", "toe", "movie", "calorie", "headache", "die", "shoe"])

# top5000nv[endswith.(top5000nv, r"fe|f")]
const f_ending_words = Set(["beef", "stuff", "wolf", "gulf", "sheriff", "grief", "chief", "life", "shelf", "relief", "roof",
        "proof", "wildlife", "playoff", "wife", "belief", "leaf", "knife", "self", "staff", "thief", "chef", "half",
        "golf", "cliff", "cafe"])


text:
__precompile__()

module OrbitTomography

using LinearAlgebra
using Statistics
using SparseArrays
using Distributed
using Printf
import Base.Iterators: partition

using EFIT
using Equilibrium
using GuidingCenterOrbits
using HDF5
using JLD2, FileIO
using Clustering
using Images
using StatsBase
using FillArrays
using ProgressMeter
using NearestNeighbors
using SparseArrays
using StaticArrays
using NonNegLeastSquares
using HCubature
using Interpolations
using Distributed
using Optim
using Sobol
using ForwardDiff
using Roots
import IterTools: nth

const S3 = SVector{3}
const S4 = SVector{4}
const S33 = SMatrix{3,3}
const S44 = SMatrix{4,4}

const e0 = 1.60217733e-19 # Coulombs / Joules
const mu0 = 4*pi*1e-7 # N/A^2
const c0 = 2.99792458e8 # m/s

const mass_u = 1.6605402e-27 # kg
const e_amu = 5.48579909070e-4 # amu
const H1_amu = 1.007276466879 # amu
const H2_amu = 2.0141017778 # amu
const H3_amu = 3.01550071632 # amu
const He3_amu = 3.01602931914 # amu
const B5_amu = 10.81 # amu
const C6_amu = 12.011 # amu

include("polyharmonic.jl")
export PolyharmonicSpline

include("spectra.jl")
export InstrumentalResponse, kernel
export ExperimentalSpectra, TheoreticalSpectra

include("io.jl")
export FIDASIMSpectra
export AbstractDistribution, FIDASIMGuidingCenterFunction, FIDASIMGuidingCenterParticles, FIDASIMFullOrbitParticles
export read_fidasim_distribution, write_fidasim_distribution
export FIDASIMBeamGeometry, FIDASIMSpectraGeometry, FIDASIMNPAGeometry, write_fidasim_geometry
export FIDASIMPlasmaParameters

include("fidasim_utils.jl")
export split_spectra, merge_spectra, apply_instrumental!, merge_spectra_geometry
export split_particles, fbm2mc
export impurity_density, ion_density

include("orbits.jl")
export OrbitGrid, orbit_grid, segment_orbit_grid,combine_orbits, fbm2orbit, mc2orbit
export map_orbits, bin_orbits
export write_orbit_grid, read_orbit_grid
export orbit_index, orbit_matrix
export OrbitSpline

include("covariance.jl")
export epr_cov
export RepeatedBlockDiagonal, ep_cov, eprz_cov,transform_eprz_cov, eprz_kernel
export get_covariance, get_correlation, get_correlation_matrix, get_covariance_matrix
export get_global_covariance, get_global_covariance_matrix

include("weights.jl")
export AbstractWeight, FIDAOrbitWeight

include("tomography.jl")
export OrbitSystem, lcurve_point, lcurve, marginal_loglike, optimize_alpha!, estimate_rtol, optimize_parameters, inv_chol, solve

include("transforms.jl")
export EPDensity, local_distribution, RZDensity, rz_profile, EPRZDensity, eprz_distribution

include("analytic.jl")
export slowing_down, bimaxwellian, maxwellian

end # module


text:
""" SemanticModels

provides the AbstractModel type and constructors for building
metamodeling tooling for categorical model representations

"""
module SemanticModels

export model, AbstractModel

"""    AbstractModel

a placeholder struct to dispatch on how to parse the expression tree into a model.
"""
abstract type AbstractModel end

""" model(::AbstractModel, x)
dig into the expression that describes a model and break it down into components. This allows you to construct a structured representation of the modeling problem at the expression level. Just like how julia modeling frameworks build structured representations of the problems in data structures. This version builds them at the expression level.
The first argument is the type you want to construct, the second argument is the model structure that you want to analyze. For example
```
model(PetriModel, x::Petri.Model)::PetriModel
```
"""
function model(::Type{T}, x) where T<:AbstractModel
  error("NotImplemented: model(::$T, $(typeof(x)))")
end

include("CategoryTheory.jl")
include("PetriModels.jl")
include("RelOlogModels.jl")
include("GraphModels.jl")
include("OpenModels.jl")
#include("OpenPetris.jl")
include("PetriCospans.jl")
include("WiringDiagrams.jl")

include("exprmodels/ExprModels.jl")

end