text:
"""
    MattssonSvärdNordström2004()

Coefficients of the SBP operators given in
- Mattsson, Svärd, Nordström (2004)
  Stable and Accurate Artificial Dissipation.
  Journal of Scientific Computing 21.1, pp. 57-79.
"""
struct MattssonSvärdNordström2004 <: SourceOfCoefficients end

function Base.show(io::IO, source::MattssonSvärdNordström2004)
    if get(io, :compact, false)
      summary(io, source)
    else
        print(io,
            "Mattsson, Svärd, Nordström (2004) \n",
            "  Stable and Accurate Artificial Dissipation. \n",
            "  Journal of Scientific Computing 21.1, pp. 57-79. \n",
            "See also (first and second derivatives) \n",
            "  Mattsson, Nordström (2004) \n",
            "  Summation by parts operators for finite difference approximations of second \n",
            "    derivatives. \n",
            "  Journal of Computational Physics 199, pp. 503-540.")
    end
end

function dissipation_coefficients(source::MattssonSvärdNordström2004, order::Int,
                                  grid, left_weights, right_weights, mode=FastMode())
    T = promote_type(eltype(grid), eltype(left_weights), eltype(right_weights))
    inv_left_weights = one(T) ./ left_weights
    inv_right_weights = one(T) ./ right_weights
    if order == 2
        coefficient_cache = MattssonSvärdNordström2004Cache2(inv_left_weights, inv_right_weights)
        b = fill(one(T), length(grid))
        b[1] = T(0)
    elseif order == 4
        coefficient_cache = MattssonSvärdNordström2004Cache4(inv_left_weights, inv_right_weights)
        b = fill(one(T), length(grid))
        b[1] = b[end] = T(0)
    elseif order == 6
        coefficient_cache = MattssonSvärdNordström2004Cache6(inv_left_weights, inv_right_weights)
        b = fill(one(T), length(grid))
        b[1] = b[2] = b[end] = T(0)
    elseif order == 8
        coefficient_cache = MattssonSvärdNordström2004Cache8(inv_left_weights, inv_right_weights)
        b = fill(one(T), length(grid))
        b[1] = b[2] = b[end] = b[end-1] = T(0)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end

    VarCoefDerivativeCoefficients(coefficient_cache, left_weights, right_weights, mode, order, 2, source), b
end

struct MattssonSvärdNordström2004Cache2{T,LeftWidth,RightWidth} <: AbstractCoefficientCache{T}
    inv_left_weights::SVector{LeftWidth,T}
    inv_right_weights::SVector{RightWidth,T}

    function MattssonSvärdNordström2004Cache2(inv_left_weights::SVector{LeftWidth,T}, inv_right_weights::SVector{RightWidth,T}) where {T,LeftWidth,RightWidth}
        boundary_length = 2

        if length(inv_left_weights) < boundary_length
            inv_left_weights = SVector(inv_left_weights..., ntuple(j->one(T), boundary_length-length(inv_left_weights))...)
        end
        if length(inv_right_weights) < boundary_length
            inv_right_weights = SVector(inv_right_weights..., ntuple(j->one(T), boundary_length-length(inv_right_weights))...)
        end

        new{T,length(inv_left_weights),length(inv_right_weights)}(inv_left_weights, inv_right_weights)
    end
end

lower_bandwidth(cache::MattssonSvärdNordström2004Cache2) = 1
upper_bandwidth(cache::MattssonSvärdNordström2004Cache2) = 1

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache2, u::AbstractVector, b::AbstractVector, α)
    @unpack inv_left_weights, inv_right_weights = cache
    α = -α # negative semidefinite operator

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    (b[1] + b[2]) * u[1]
                    + (b[1] - b[2]) * u[2]
                )
        dest[2] = α * inv_left_weights[2] * (
                    (b[1] - b[2]) * u[1]
                    + (b[1] + b[2] + b[3]) * u[2]
                    - b[3] * u[3]
                )
        for i in 3:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        -b[i] * u[i-1]
                        + (b[i] + b[i+1]) * u[i]
                        - b[i+1] * u[i+1]
                    )
        end

        dest[end] = α * inv_right_weights[1] * (
                        b[end] * u[end]
                        - b[end] * u[end-1]
                    )
        dest[end-1] = α * inv_right_weights[2] * (
                        -b[end-1] * u[end-2]
                        + (b[end-1] + b[end]) * u[end-1]
                        - b[end] * u[end]
                    )
        for i in 2:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        -b[end-i] * u[end-(i-1)]
                        + (b[end-i] + b[end-(i+1)]) * u[end-i]
                        - b[end-(i+1)] * u[end-(i+1)]
                    )
        end
    end
end

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache2, u::AbstractVector, b::AbstractVector, α, β)
    @unpack inv_left_weights, inv_right_weights = cache
    α = -α # negative semidefinite operator

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    (b[1] + b[2]) * u[1]
                    + (b[1] - b[2]) * u[2]
                ) + β*dest[1]
        dest[2] = α * inv_left_weights[2] * (
                    (b[1] - b[2]) * u[1]
                    + (b[1] + b[2] + b[3]) * u[2]
                    - b[3] * u[3]
                ) + β*dest[2]
        for i in 3:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        -b[i] * u[i-1]
                        + (b[i] + b[i+1]) * u[i]
                        - b[i+1] * u[i+1]
                    ) + β*dest[i]
        end

        dest[end] = α * inv_right_weights[1] * (
                        b[end] * u[end]
                        - b[end] * u[end-1]
                    ) + β*dest[end]
        dest[end-1] = α * inv_right_weights[2] * (
                        -b[end-1] * u[end-2]
                        + (b[end-1] + b[end]) * u[end-1]
                        - b[end] * u[end]
                    ) + β*dest[end-1]
        for i in 2:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        -b[end-i] * u[end-(i-1)]
                        + (b[end-i] + b[end-(i+1)]) * u[end-i]
                        - b[end-(i+1)] * u[end-(i+1)]
                    ) + β*dest[end-i]
        end
    end
end

@inline function convolve_interior_coefficients_loopbody(i, cache::MattssonSvärdNordström2004Cache2, u, b)
    @inbounds begin
        b_i = b[i]
        b_ip1 = b[i+1]

        retval = b_i * u[i-1] - (b_i + b_ip1) * u[i] + b_ip1 * u[i+1]
    end

    retval
end

struct MattssonSvärdNordström2004Cache4{T,LeftWidth,RightWidth} <: AbstractCoefficientCache{T}
    inv_left_weights::SVector{LeftWidth,T}
    inv_right_weights::SVector{RightWidth,T}

    function MattssonSvärdNordström2004Cache4(inv_left_weights::SVector{LeftWidth,T}, inv_right_weights::SVector{RightWidth,T}) where {T,LeftWidth,RightWidth}
        boundary_length = 3

        if length(inv_left_weights) < boundary_length
            inv_left_weights = SVector(inv_left_weights..., ntuple(j->one(T), boundary_length-length(inv_left_weights))...)
        end
        if length(inv_right_weights) < boundary_length
            inv_right_weights = SVector(inv_right_weights..., ntuple(j->one(T), boundary_length-length(inv_right_weights))...)
        end

        new{T,length(inv_left_weights),length(inv_right_weights)}(inv_left_weights, inv_right_weights)
    end
end

lower_bandwidth(cache::MattssonSvärdNordström2004Cache4) = 2
upper_bandwidth(cache::MattssonSvärdNordström2004Cache4) = 2

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache4, u::AbstractVector, b::AbstractVector, α)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2]) * u[1]
                    + 2*(b[1] + b[2]) * u[2]
                    - (b[1] + b[2]) * u[3]
                )
        dest[2] = α * inv_left_weights[2] * (
                    2*(b[1] + b[2]) * u[1]
                    - (4*b[1] + 4*b[2] + b[3]) * u[2]
                    + 2*(b[1] + b[2] + b[3]) * u[3]
                    - b[3] * u[4]
                )
        dest[3] = α * inv_left_weights[3] * (
                    - (b[1] + b[2]) * u[1]
                    + 2*(b[1] + b[2] + b[3]) * u[2]
                    - (b[1] + b[2] + 4*b[3] + b[4]) * u[3]
                    + 2*(b[3] + b[4]) * u[4]
                    - b[4] * u[5]
                )
        for i in 4:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        - b[i-1] * u[i-2]
                        + 2*(b[i-1] + b[i]) * u[i-1]
                        - (b[i-1] + 4*b[i] + b[i+1]) * u[i]
                        + 2*(b[i] + b[i+1]) * u[i+1]
                        - b[i+1] * u[i+2]
                    )
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end] + b[end-1]) * u[end]
                        + 2*(b[end] + b[end-1]) * u[end-1]
                        - (b[end] + b[end-1]) * u[end-2]
                    )
        dest[end-1] = α * inv_right_weights[2] * (
                        2*(b[end] + b[end-1]) * u[end]
                        - (b[end-2] + 4*b[end-1] + 4*b[end]) * u[end-1]
                        + 2*(b[end-2] + b[end-1] + b[end]) * u[end-2]
                        - b[end-2] * u[end-3]
                    )
        dest[end-2] = α * inv_right_weights[3] * (
                        - (b[end-1] + b[end]) * u[end]
                        + 2*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - (b[end-3] + 4*b[end-2] + b[end-1] + b[end]) * u[end-2]
                        + 2*(b[end-3] + b[end-2]) * u[end-3]
                        - b[end-3] * u[end-4]
                    )
        for i in 3:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        - b[end-(i-1)] * u[end-(i-2)]
                        + 2*(b[end-(i-1)] + b[i]) * u[end-(i-1)]
                        - (b[end-(i-1)] + 4*b[end-i] + b[end-(i+1)]) * u[end-i]
                        + 2*(b[end-i] + b[end-(i+1)]) * u[end-(i+1)]
                        - b[end-(i+1)] * u[end-(i+2)]
                    )
        end
    end
end

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache4, u::AbstractVector, b::AbstractVector, α, β)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2]) * u[1]
                    + 2*(b[1] + b[2]) * u[2]
                    - (b[1] + b[2]) * u[3]
                ) + β*dest[1]
        dest[2] = α * inv_left_weights[2] * (
                    2*(b[1] + b[2]) * u[1]
                    - (4*b[1] + 4*b[2] + b[3]) * u[2]
                    + 2*(b[1] + b[2] + b[3]) * u[3]
                    - b[3] * u[4]
                ) + β*dest[2]
        dest[3] = α * inv_left_weights[3] * (
                    - (b[1] + b[2]) * u[1]
                    + 2*(b[1] + b[2] + b[3]) * u[2]
                    - (b[1] + b[2] + 4*b[3] + b[4]) * u[3]
                    + 2*(b[3] + b[4]) * u[4]
                    - b[4] * u[5]
                ) + β*dest[3]
        for i in 4:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        - b[i-1] * u[i-2]
                        + 2*(b[i-1] + b[i]) * u[i-1]
                        - (b[i-1] + 4*b[i] + b[i+1]) * u[i]
                        + 2*(b[i] + b[i+1]) * u[i+1]
                        - b[i+1] * u[i+2]
                    ) + β*dest[i]
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end] + b[end-1]) * u[end]
                        + 2*(b[end] + b[end-1]) * u[end-1]
                        - (b[end] + b[end-1]) * u[end-2]
                    ) + β*dest[end]
        dest[end-1] = α * inv_right_weights[2] * (
                        + 2*(b[end] + b[end-1]) * u[end]
                        - (b[end-2] + 4*b[end-1] + 4*b[end]) * u[end-1]
                        + 2*(b[end-2] + b[end-1] + b[end]) * u[end-2]
                        - b[end-2] * u[end-3]
                    ) + β*dest[end-1]
        dest[end-2] = α * inv_right_weights[3] * (
                        - (b[end-1] + b[end]) * u[end]
                        + 2*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - (b[end-3] + 4*b[end-2] + b[end-1] + b[end]) * u[end-2]
                        + 2*(b[end-3] + b[end-2]) * u[end-3]
                        - b[end-3] * u[end-4]
                    ) + β*dest[end-2]
        for i in 3:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        - b[end-(i-1)] * u[end-(i-2)]
                        + 2*(b[end-(i-1)] + b[i]) * u[end-(i-1)]
                        - (b[end-(i-1)] + 4*b[end-i] + b[end-(i+1)]) * u[end-i]
                        + 2*(b[end-i] + b[end-(i+1)]) * u[end-(i+1)]
                        - b[end-(i+1)] * u[end-(i+2)]
                    ) + β*dest[end-i]
        end
    end
end

@inline function convolve_interior_coefficients_loopbody(i, cache::MattssonSvärdNordström2004Cache4, u, b)
    @inbounds begin
        b_im1 = b[i-1]
        b_i   = b[i]
        b_ip1 = b[i+1]

        retval = (
                    - b_im1 * u[i-2]
                    + 2*(b_im1 + b_i) * u[i-1]
                    - (b_im1 + 4b_i + b_ip1) * u[i]
                    + 2*(b_i + b_ip1) * u[i+1]
                    - b_ip1 * u[i+2]
                )
    end

    retval
end

struct MattssonSvärdNordström2004Cache6{T,LeftWidth,RightWidth} <: AbstractCoefficientCache{T}
    inv_left_weights::SVector{LeftWidth,T}
    inv_right_weights::SVector{RightWidth,T}

    function MattssonSvärdNordström2004Cache6(inv_left_weights::SVector{LeftWidth,T}, inv_right_weights::SVector{RightWidth,T}) where {T,LeftWidth,RightWidth}
        boundary_length = 4

        if length(inv_left_weights) < boundary_length
            inv_left_weights = SVector(inv_left_weights..., ntuple(j->one(T), boundary_length-length(inv_left_weights))...)
        end
        if length(inv_right_weights) < boundary_length
            inv_right_weights = SVector(inv_right_weights..., ntuple(j->one(T), boundary_length-length(inv_right_weights))...)
        end

        new{T,length(inv_left_weights),length(inv_right_weights)}(inv_left_weights, inv_right_weights)
    end
end

lower_bandwidth(cache::MattssonSvärdNordström2004Cache6) = 3
upper_bandwidth(cache::MattssonSvärdNordström2004Cache6) = 3

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache6, u::AbstractVector, b::AbstractVector, α)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 3*(b[1] + b[2] + b[3]) * u[2]
                    - 3*(b[1] + b[2] + b[3]) * u[3]
                    + (b[1] + b[2] + b[3]) * u[4]
                )
        dest[2] = α * inv_left_weights[2] * (
                    3*(b[1] + b[2] + b[3]) * u[1]
                    - (9*b[1] + 9*b[2] + 9*b[3] + b[4]) * u[2]
                    + 3*(3*b[1] + 3*b[2] + 3*b[3] + b[4]) * u[3]
                    - 3*(b[1] + b[2] + b[3] + b[4]) * u[4]
                    + b[4] * u[5]
                )
        dest[3] = α * inv_left_weights[3] * (
                    - 3*(b[1] + b[2] + b[3]) * u[1]
                    + 3*(3*b[1] + 3*b[2] + 3*b[3] + b[4]) * u[2]
                    - (9*b[1] + 9*b[2] + 9*b[3] + 9*b[4] + b[5]) * u[3]
                    + 3*(b[1] + b[2] + b[3] + 3*b[4] + b[5]) * u[4]
                    - 3*(b[4] + b[5]) * u[5]
                    + b[5] * u[6]
                )
        dest[4] = α * inv_left_weights[4] * (
                    (b[1] + b[2] + b[3]) * u[1]
                    - 3*(b[1] + b[2] + b[3] + b[4]) * u[2]
                    + 3*(b[1] + b[2] + b[3] + 3*b[4] + b[5]) * u[3]
                    - (b[1] + b[2] + b[3] + 9*b[4] + 9*b[5] + b[6]) * u[4]
                    + 3*(b[4] + 3*b[5] + b[6]) * u[5]
                    - 3*(b[5] + b[6]) * u[6]
                    + b[6] * u[7]
                )
        for i in 5:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        b[i-1] * u[i-3]
                        - 3*(b[i-1] + b[i]) * u[i-2]
                        + 3*(b[i-1] + 3*b[i] + b[i+1]) * u[i-1]
                        - (b[i-1] + 9*b[i] + 9*b[i+1] + b[i+2]) * u[i]
                        + 3*(b[i] + 3*b[i+1] + b[i+2]) * u[i+1]
                        - 3*(b[i+1] + b[i+2]) * u[i+2]
                        + b[i+2] * u[i+3]
                    )
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end] + b[end-1]) * u[end]
                        + 3*(b[end] + b[end-1]) * u[end-1]
                        - 3*(b[end] + b[end-1]) * u[end-2]
                        + (b[end] + b[end-1]) * u[end-3]
                    )
        dest[end-1] = α * inv_right_weights[2] * (
                        3*(b[end] + b[end-1]) * u[end]
                        - (b[end-2] + 9*b[end-1] + 9*b[end]) * u[end-1]
                        + 3*(b[end-2] + 3*b[end-1] + 3*b[end]) * u[end-2]
                        - 3*(b[end-2] + b[end-1] + b[end]) * u[end-3]
                        + b[end-2] * u[end-4]
                    )
        dest[end-2] = α * inv_right_weights[3] * (
                        - 3*(b[end] + b[end-1]) * u[end]
                        + 3*(b[end-2] + 3*b[end-1] + 3*b[end]) * u[end-1]
                        - (b[end-3] + 9*b[end-2] + 9*b[end-1] + 9*b[end]) * u[end-2]
                        + 3*(b[end-3] + 3*b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - 3*(b[end-3] + b[end-2]) * u[end-4]
                        + b[end-3] * u[end-5]
                    )
        dest[end-3] = α * inv_right_weights[4] * (
                        (b[end] + b[end-1]) * u[end]
                        - 3*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        + 3*(b[end-3] + 3*b[end-2] + b[end-1] + b[end]) * u[end-2]
                        - (b[end-4] + 9*b[end-3] + 9*b[end-2] + b[end-1] + b[end]) * u[end-3]
                        + 3*(b[end-4] + 3*b[end-3] + b[end-2]) * u[end-4]
                        - 3*(b[end-4] + b[end-3]) * u[end-5]
                        + b[end-4] * u[end-6]
                    )
        for i in 4:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        b[end-(i+1)] * u[end-(i+3)]
                        - 3*(b[end-(i+1)] + b[end-i]) * u[end-(i+2)]
                        + 3*(b[end-(i+1)] + 3*b[end-i] + b[end-(i-1)]) * u[end-(i+1)]
                        - (b[end-(i+1)] + 9*b[end-i] + 9*b[end-(i-1)] + b[end-(i-2)]) * u[end-i]
                        + 3*(b[end-i] + 3*b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-1)]
                        - 3*(b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-2)]
                        + b[end-(i-2)] * u[end-(i-3)]
                    )
        end
    end
end

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache6, u::AbstractVector, b::AbstractVector, α, β)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 3*(b[1] + b[2] + b[3]) * u[2]
                    - 3*(b[1] + b[2] + b[3]) * u[3]
                    + (b[1] + b[2] + b[3]) * u[4]
                ) + β*dest[1]
        dest[2] = α * inv_left_weights[2] * (
                    3*(b[1] + b[2] + b[3]) * u[1]
                    - (9*b[1] + 9*b[2] + 9*b[3] + b[4]) * u[2]
                    + 3*(3*b[1] + 3*b[2] + 3*b[3] + b[4]) * u[3]
                    - 3*(b[1] + b[2] + b[3] + b[4]) * u[4]
                    + b[4] * u[5]
                ) + β*dest[2]
        dest[3] = α * inv_left_weights[3] * (
                    - 3*(b[1] + b[2] + b[3]) * u[1]
                    + 3*(3*b[1] + 3*b[2] + 3*b[3] + b[4]) * u[2]
                    - (9*b[1] + 9*b[2] + 9*b[3] + 9*b[4] + b[5]) * u[3]
                    + 3*(b[1] + b[2] + b[3] + 3*b[4] + b[5]) * u[4]
                    - 3*(b[4] + b[5]) * u[5]
                    + b[5] * u[6]
                ) + β*dest[3]
        dest[4] = α * inv_left_weights[4] * (
                    (b[1] + b[2] + b[3]) * u[1]
                    - 3*(b[1] + b[2] + b[3] + b[4]) * u[2]
                    + 3*(b[1] + b[2] + b[3] + 3*b[4] + b[5]) * u[3]
                    - (b[1] + b[2] + b[3] + 9*b[4] + 9*b[5] + b[6]) * u[4]
                    + 3*(b[4] + 3*b[5] + b[6]) * u[5]
                    - 3*(b[5] + b[6]) * u[6]
                    + b[6] * u[7]
                ) + β*dest[4]
        for i in 5:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        b[i-1] * u[i-3]
                        - 3*(b[i-1] + b[i]) * u[i-2]
                        + 3*(b[i-1] + 3*b[i] + b[i+1]) * u[i-1]
                        - (b[i-1] + 9*b[i] + 9*b[i+1] + b[i+2]) * u[i]
                        + 3*(b[i] + 3*b[i+1] + b[i+2]) * u[i+1]
                        - 3*(b[i+1] + b[i+2]) * u[i+2]
                        + b[i+2] * u[i+3]
                    ) + β*dest[i]
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end] + b[end-1]) * u[end]
                        + 3*(b[end] + b[end-1]) * u[end-1]
                        - 3*(b[end] + b[end-1]) * u[end-2]
                        + (b[end] + b[end-1]) * u[end-3]
                    ) + β*dest[end]
        dest[end-1] = α * inv_right_weights[2] * (
                        3*(b[end] + b[end-1]) * u[end]
                        - (b[end-2] + 9*b[end-1] + 9*b[end]) * u[end-1]
                        + 3*(b[end-2] + 3*b[end-1] + 3*b[end]) * u[end-2]
                        - 3*(b[end-2] + b[end-1] + b[end]) * u[end-3]
                        + b[end-2] * u[end-4]
                    ) + β*dest[end-1]
        dest[end-2] = α * inv_right_weights[3] * (
                        - 3*(b[end] + b[end-1]) * u[end]
                        + 3*(b[end-2] + 3*b[end-1] + 3*b[end]) * u[end-1]
                        - (b[end-3] + 9*b[end-2] + 9*b[end-1] + 9*b[end]) * u[end-2]
                        + 3*(b[end-3] + 3*b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - 3*(b[end-3] + b[end-2]) * u[end-4]
                        + b[end-3] * u[end-5]
                    ) + β*dest[end-2]
        dest[end-3] = α * inv_right_weights[4] * (
                        (b[end] + b[end-1]) * u[end]
                        - 3*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        + 3*(b[end-3] + 3*b[end-2] + b[end-1] + b[end]) * u[end-2]
                        - (b[end-4] + 9*b[end-3] + 9*b[end-2] + b[end-1] + b[end]) * u[end-3]
                        + 3*(b[end-4] + 3*b[end-3] + b[end-2]) * u[end-4]
                        - 3*(b[end-4] + b[end-3]) * u[end-5]
                        + b[end-4] * u[end-6]
                    ) + β*dest[end-3]
        for i in 4:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        b[end-(i+1)] * u[end-(i+3)]
                        - 3*(b[end-(i+1)] + b[end-i]) * u[end-(i+2)]
                        + 3*(b[end-(i+1)] + 3*b[end-i] + b[end-(i-1)]) * u[end-(i+1)]
                        - (b[end-(i+1)] + 9*b[end-i] + 9*b[end-(i-1)] + b[end-(i-2)]) * u[end-i]
                        + 3*(b[end-i] + 3*b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-1)]
                        - 3*(b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-2)]
                        + b[end-(i-2)] * u[end-(i-3)]
                    ) + β*dest[end-i]
        end
    end
end

@inline function convolve_interior_coefficients_loopbody(i, cache::MattssonSvärdNordström2004Cache6, u, b)
    @inbounds begin
        b_im1 = b[i-1]
        b_i   = b[i]
        b_ip1 = b[i+1]
        b_ip2 = b[i+2]

        retval = (
                    b_im1 * u[i-3]
                    - 3*(b_im1 + b_i) * u[i-2]
                    + 3*(b_im1 + 3*b_i + b_ip1) * u[i-1]
                    - (b_im1 + 9*b_i + 9*b_ip1 + b_ip2) * u[i]
                    + 3*(b_i + 3*b_ip1 + b_ip2) * u[i+1]
                    - 3*(b_ip1 + b_ip2) * u[i+2]
                    + b_ip2 * u[i+3]
                )
    end

    retval
end

struct MattssonSvärdNordström2004Cache8{T,LeftWidth,RightWidth} <: AbstractCoefficientCache{T}
    inv_left_weights::SVector{LeftWidth,T}
    inv_right_weights::SVector{RightWidth,T}

    function MattssonSvärdNordström2004Cache8(inv_left_weights::SVector{LeftWidth,T}, inv_right_weights::SVector{RightWidth,T}) where {T,LeftWidth,RightWidth}
        boundary_length = 5

        if length(inv_left_weights) < boundary_length
            inv_left_weights = SVector(inv_left_weights..., ntuple(j->one(T), boundary_length-length(inv_left_weights))...)
        end
        if length(inv_right_weights) < boundary_length
            inv_right_weights = SVector(inv_right_weights..., ntuple(j->one(T), boundary_length-length(inv_right_weights))...)
        end

        new{T,length(inv_left_weights),length(inv_right_weights)}(inv_left_weights, inv_right_weights)
    end
end

lower_bandwidth(cache::MattssonSvärdNordström2004Cache8) = 4
upper_bandwidth(cache::MattssonSvärdNordström2004Cache8) = 4

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache8, u::AbstractVector, b::AbstractVector, α)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 4*(b[1] + b[2] + b[3]) * u[2]
                    - 6*(b[1] + b[2] + b[3]) * u[3]
                    + 4*(b[1] + b[2] + b[3]) * u[4]
                    - (b[1] + b[2] + b[3]) * u[5]
                )
        dest[2] = α * inv_left_weights[2] * (
                    4*(b[1] + b[2] + b[3]) * u[1]
                    - (16*b[1] + 16*b[2] + 16*b[3] + b[4]) * u[2]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + b[4]) * u[3]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 6*b[4]) * u[4]
                    + 4*(b[1] + b[2] + b[3] + b[4]) * u[5]
                    - b[4] * u[6]
                )
        dest[3] = α * inv_left_weights[3] * (
                    - 6*(b[1] + b[2] + b[3]) * u[1]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + b[4]) * u[2]
                    - (36*b[1] + 36*b[2] + 36*b[3] + 16*b[4] + b[5]) * u[3]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + 6*b[4] + b[5]) * u[4]
                    - (6*b[1] + 6*b[2] + 6*b[3] + 16*b[4] + 6*b[5]) * u[5]
                    + 4*(b[4] + b[5]) * u[6]
                    - b[5] * u[7]
                )
        dest[4] = α * inv_left_weights[4] * (
                    4*(b[1] + b[2] + b[3]) * u[1]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 6*b[4]) * u[2]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + 6*b[4] + b[5]) * u[3]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 36*b[4] + 16*b[5] + b[6]) * u[4]
                    + 4*(b[1] + b[2] + b[3] + 6*b[4] + 6*b[5] + b[6]) * u[5]
                    - (6*b[4] + 16*b[5] + 6*b[6]) * u[6]
                    + 4*(b[5] + b[6]) * u[7]
                    - b[6] * u[8]
                )
        dest[5] = α * inv_left_weights[5] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 4*(b[1] + b[2] + b[3] + b[4]) * u[2]
                    - (6*b[1] + 6*b[2] + 6*b[3] + 16*b[4] + 6*b[5]) * u[3]
                    + 4*(b[1] + b[2] + b[3] + 6*b[4] + 6*b[5] + b[6]) * u[4]
                    - (b[1] + b[2] + b[3] + 16*b[4] + 36*b[5] + 16*b[6] + b[7]) * u[5]
                    + 4*(b[4] + 6*b[5] + 6*b[6] + b[7]) * u[6]
                    - (6*b[5] + 16*b[6] + 6*b[7]) * u[7]
                    + 4*(b[6] + b[7]) * u[8]
                    - b[7] * u[9]
                )
        for i in 6:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        - b[i-2] * u[i-4]
                        + 4*(b[i-2] + b[i-1]) * u[i-3]
                        - (6*b[i-2] + 16*b[i-1] + 6*b[i]) * u[i-2]
                        + 4*(b[i-2] + 6*b[i-1] + 6*b[i] + b[i+1]) * u[i-1]
                        - (b[i-2] + 16*b[i-1] + 36*b[i] + 16*b[i+1] + b[i+2]) * u[i]
                        + 4*(b[i-1] + 6*b[i] + 6*b[i+1] + b[i+2]) * u[i+1]
                        - (6*b[i] + 16*b[i+1] + 6*b[i+2]) * u[i+2]
                        + 4*(b[i+1] + b[i+2]) * u[i+3]
                        - b[i+2] * u[i+4]
                    )
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - 6*(b[end-2] + b[end-1] + b[end]) * u[end-2]
                        + 4*(b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - (b[end-2] + b[end-1] + b[end]) * u[end-4]
                    )
        dest[end-1] = α * inv_right_weights[2] * (
                        4*(b[end-2] + b[end-1] + b[end]) * u[end]
                        - (b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-1]
                        + 4*(b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        - (6*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-3]
                        + 4*(b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        - b[end-3] * u[end-5]
                    )
        dest[end-2] = α * inv_right_weights[3] * (
                        - 6*(b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-1]
                        - (b[end-4] + 16*b[end-3] + 36*b[end-2] + 36*b[end-1] + 36*b[end]) * u[end-2]
                        + 4*(b[end-4] + 6*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-3]
                        - (6*b[end-4] + 16*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-4]
                        + 4*(b[end-4] + b[end-3]) * u[end-5]
                        - b[end-4] * u[end-6]
                    )
        dest[end-3] = α * inv_right_weights[4] * (
                        4*(b[end-2] + b[end-1] + b[end]) * u[end]
                        - (6*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-1]
                        + 4*(b[end-4] + 6*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        - (b[end-5] + 16*b[end-4] + 36*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-3]
                        + 4*(b[end-5] + 6*b[end-4] + 6*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        - (6*b[end-5] + 16*b[end-4] + 6*b[end-3]) * u[end-5]
                        + 4*(b[end-5] + b[end-4]) * u[end-6]
                        - b[end-5] * u[end-7]
                    )
        dest[end-4] = α * inv_right_weights[5] * (
                        - (b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - (6*b[end-4] + 16*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        + 4*(b[end-5] + 6*b[end-4] + 6*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - (b[end-6] + 16*b[end-5] + 36*b[end-4] + 16*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        + 4*(b[end-6] + 6*b[end-5] + 6*b[end-4] + b[end-3]) * u[end-5]
                        - (6*b[end-6] + 16*b[end-5] + 6*b[end-4]) * u[end-6]
                        + 4*(b[end-6] + b[end-5]) * u[end-7]
                        - b[end-6] * u[end-8]
                    )
        for i in 5:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        - b[end-(i+2)] * u[end-(i+4)]
                        + 4*(b[end-(i+2)] + b[end-(i+1)]) * u[end-(i+3)]
                        - (6*b[end-(i+2)] + 16*b[end-(i+1)] + 6*b[end-i]) * u[end-(i+2)]
                        + 4*(b[end-(i+2)] + 6*b[end-(i+1)] + 6*b[end-i] + b[end-(i-1)]) * u[end-(i+1)]
                        - (b[end-(i+2)] + 16*b[end-(i+1)] + 36*b[end-i] + 16*b[end-(i-1)] + b[end-(i-2)]) * u[end-i]
                        + 4*(b[end-(i+1)] + 6*b[end-i] + 6*b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-1)]
                        - (6*b[end-i] + 16*b[end-(i-1)] + 6*b[end-(i-2)]) * u[end-(i-2)]
                        + 4*(b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-3)]
                        - b[end-(i-2)] * u[end-(i-4)]
                    )
        end
    end
end

function convolve_boundary_coefficients!(dest::AbstractVector, cache::MattssonSvärdNordström2004Cache8, u::AbstractVector, b::AbstractVector, α, β)
    @unpack inv_left_weights, inv_right_weights = cache

    @inbounds begin
        dest[1] = α * inv_left_weights[1] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 4*(b[1] + b[2] + b[3]) * u[2]
                    - 6*(b[1] + b[2] + b[3]) * u[3]
                    + 4*(b[1] + b[2] + b[3]) * u[4]
                    - (b[1] + b[2] + b[3]) * u[5]
                ) + β*dest[1]
        dest[2] = α * inv_left_weights[2] * (
                    + 4*(b[1] + b[2] + b[3]) * u[1]
                    - (16*b[1] + 16*b[2] + 16*b[3] + b[4]) * u[2]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + b[4]) * u[3]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 6*b[4]) * u[4]
                    + 4*(b[1] + b[2] + b[3] + b[4]) * u[5]
                    - b[4] * u[6]
                ) + β*dest[2]
        dest[3] = α * inv_left_weights[3] * (
                    - 6*(b[1] + b[2] + b[3]) * u[1]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + b[4]) * u[2]
                    - (36*b[1] + 36*b[2] + 36*b[3] + 16*b[4] + b[5]) * u[3]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + 6*b[4] + b[5]) * u[4]
                    - (6*b[1] + 6*b[2] + 6*b[3] + 16*b[4] + 6*b[5]) * u[5]
                    + 4*(b[4] + b[5]) * u[6]
                    - b[5] * u[7]
                ) + β*dest[3]
        dest[4] = α * inv_left_weights[4] * (
                    4*(b[1] + b[2] + b[3]) * u[1]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 6*b[4]) * u[2]
                    + 4*(6*b[1] + 6*b[2] + 6*b[3] + 6*b[4] + b[5]) * u[3]
                    - (16*b[1] + 16*b[2] + 16*b[3] + 36*b[4] + 16*b[5] + b[6]) * u[4]
                    + 4*(b[1] + b[2] + b[3] + 6*b[4] + 6*b[5] + b[6]) * u[5]
                    - (6*b[4] + 16*b[5] + 6*b[6]) * u[6]
                    + 4*(b[5] + b[6]) * u[7]
                    - b[6] * u[8]
                ) + β*dest[4]
        dest[5] = α * inv_left_weights[5] * (
                    - (b[1] + b[2] + b[3]) * u[1]
                    + 4*(b[1] + b[2] + b[3] + b[4]) * u[2]
                    - (6*b[1] + 6*b[2] + 6*b[3] + 16*b[4] + 6*b[5]) * u[3]
                    + 4*(b[1] + b[2] + b[3] + 6*b[4] + 6*b[5] + b[6]) * u[4]
                    - (b[1] + b[2] + b[3] + 16*b[4] + 36*b[5] + 16*b[6] + b[7]) * u[5]
                    + 4*(b[4] + 6*b[5] + 6*b[6] + b[7]) * u[6]
                    - (6*b[5] + 16*b[6] + 6*b[7]) * u[7]
                    + 4*(b[6] + b[7]) * u[8]
                    - b[7] * u[9]
                ) + β*dest[5]
        for i in 6:length(inv_left_weights)
            dest[i] = α * inv_left_weights[i] * (
                        - b[i-2] * u[i-4]
                        + 4*(b[i-2] + b[i-1]) * u[i-3]
                        - (6*b[i-2] + 16*b[i-1] + 6*b[i]) * u[i-2]
                        + 4*(b[i-2] + 6*b[i-1] + 6*b[i] + b[i+1]) * u[i-1]
                        - (b[i-2] + 16*b[i-1] + 36*b[i] + 16*b[i+1] + b[i+2]) * u[i]
                        + 4*(b[i-1] + 6*b[i] + 6*b[i+1] + b[i+2]) * u[i+1]
                        - (6*b[i] + 16*b[i+1] + 6*b[i+2]) * u[i+2]
                        + 4*(b[i+1] + b[i+2]) * u[i+3]
                        - b[i+2] * u[i+4]
                    ) + β*dest[i]
        end

        dest[end] = α * inv_right_weights[1] * (
                        - (b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - 6*(b[end-2] + b[end-1] + b[end]) * u[end-2]
                        + 4*(b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - (b[end-2] + b[end-1] + b[end]) * u[end-4]
                    ) + β*dest[end]
        dest[end-1] = α * inv_right_weights[2] * (
                        4*(b[end-2] + b[end-1] + b[end]) * u[end]
                        - (b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-1]
                        + 4*(b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        - (6*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-3]
                        + 4*(b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        - b[end-3] * u[end-5]
                    ) + β*dest[end-1]
        dest[end-2] = α * inv_right_weights[3] * (
                        - 6*(b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-1]
                        - (b[end-4] + 16*b[end-3] + 36*b[end-2] + 36*b[end-1] + 36*b[end]) * u[end-2]
                        + 4*(b[end-4] + 6*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-3]
                        - (6*b[end-4] + 16*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-4]
                        + 4*(b[end-4] + b[end-3]) * u[end-5]
                        - b[end-4] * u[end-6]
                    ) + β*dest[end-2]
        dest[end-3] = α * inv_right_weights[4] * (
                        4*(b[end-2] + b[end-1] + b[end]) * u[end]
                        - (6*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-1]
                        + 4*(b[end-4] + 6*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        - (b[end-5] + 16*b[end-4] + 36*b[end-3] + 16*b[end-2] + 16*b[end-1] + 16*b[end]) * u[end-3]
                        + 4*(b[end-5] + 6*b[end-4] + 6*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        - (6*b[end-5] + 16*b[end-4] + 6*b[end-3]) * u[end-5]
                        + 4*(b[end-5] + b[end-4]) * u[end-6]
                        - b[end-5] * u[end-7]
                    ) + β*dest[end-3]
        dest[end-4] = α * inv_right_weights[5] * (
                        - (b[end-2] + b[end-1] + b[end]) * u[end]
                        + 4*(b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-1]
                        - (6*b[end-4] + 16*b[end-3] + 6*b[end-2] + 6*b[end-1] + 6*b[end]) * u[end-2]
                        + 4*(b[end-5] + 6*b[end-4] + 6*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-3]
                        - (b[end-6] + 16*b[end-5] + 36*b[end-4] + 16*b[end-3] + b[end-2] + b[end-1] + b[end]) * u[end-4]
                        + 4*(b[end-6] + 6*b[end-5] + 6*b[end-4] + b[end-3]) * u[end-5]
                        - (6*b[end-6] + 16*b[end-5] + 6*b[end-4]) * u[end-6]
                        + 4*(b[end-6] + b[end-5]) * u[end-7]
                        - b[end-6] * u[end-8]
                    ) + β*dest[end-4]
        for i in 5:length(inv_right_weights)-1
            dest[end-i] = α * inv_right_weights[i+1] * (
                        - b[end-(i+2)] * u[end-(i+4)]
                        + 4*(b[end-(i+2)] + b[end-(i+1)]) * u[end-(i+3)]
                        - (6*b[end-(i+2)] + 16*b[end-(i+1)] + 6*b[end-i]) * u[end-(i+2)]
                        + 4*(b[end-(i+2)] + 6*b[end-(i+1)] + 6*b[end-i] + b[end-(i-1)]) * u[end-(i+1)]
                        - (b[end-(i+2)] + 16*b[end-(i+1)] + 36*b[end-i] + 16*b[end-(i-1)] + b[end-(i-2)]) * u[end-i]
                        + 4*(b[end-(i+1)] + 6*b[end-i] + 6*b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-1)]
                        - (6*b[end-i] + 16*b[end-(i-1)] + 6*b[end-(i-2)]) * u[end-(i-2)]
                        + 4*(b[end-(i-1)] + b[end-(i-2)]) * u[end-(i-3)]
                        - b[end-(i-2)] * u[end-(i-4)]
                    ) + β*dest[end-i]
        end
    end
end

@inline function convolve_interior_coefficients_loopbody(i, cache::MattssonSvärdNordström2004Cache8, u, b)
    @inbounds begin
        b_im2 = b[i-2]
        b_im1 = b[i-1]
        b_i   = b[i]
        b_ip1 = b[i+1]
        b_ip2 = b[i+2]

        retval = (
                    - b_im2 * u[i-4]
                    + 4*(b_im2 + b_im1) * u[i-3]
                    - (6*b_im2 + 16*b_im1 + 6*b_i) * u[i-2]
                    + 4*(b_im2 + 6*b_im1 + 6*b_i + b_ip1) * u[i-1]
                    - (b_im2 + 16*b_im1 + 36*b_i + 16*b_ip1 + b_ip2) * u[i]
                    + 4*(b_im1 + 6*b_i + 6*b_ip1 + b_ip2) * u[i+1]
                    - (6*b_i + 16*b_ip1 + 6*b_ip2) * u[i+2]
                    + 4*(b_ip1 + b_ip2) * u[i+3]
                    - b_ip2 * u[i+4]
                )
    end

    retval
end

function first_derivative_coefficients(source::MattssonSvärdNordström2004, order::Int,
                                       T=Float64, mode=FastMode())
    if order == 2
        left_boundary = (
            DerivativeCoefficientRow{T,1,2}(SVector(-one(T), one(T))),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(1//2))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 4
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,4}(SVector(T(-24//17),
                                                    T(59//34),
                                                    T(-4//17),
                                                    T(-3//34) )),
            # q2
            DerivativeCoefficientRow{T,1,3}(SVector(T(-1//2),
                                                    T(0),
                                                    T(1//2))),
            # q3
            DerivativeCoefficientRow{T,1,5}(SVector(T(4//43),
                                                    T(-59//86),
                                                    T(0),
                                                    T(59//86),
                                                    T(-4//43))),
            # q4
            DerivativeCoefficientRow{T,1,6}(SVector(T(3//98),
                                                    T(0),
                                                    T(-59//98),
                                                    T(0),
                                                    T(32//49),
                                                    T(-4//49))),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(2//3), T(-1//12))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(17//48),
                                T(59//48),
                                T(43//48),
                                T(49//48) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 6
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,6}(SVector(T(-21600//13649),
                                                    T(104009//54596),
                                                    T(30443//81894),
                                                    T(-33311//27298),
                                                    T(16863//27298),
                                                    T(-15025//163788) )),
            # q2
            DerivativeCoefficientRow{T,1,6}(SVector(T(-104009//240260),
                                                    T(0),
                                                    T(-311//72078),
                                                    T(20229//24026),
                                                    T(-24337//48052),
                                                    T(36661//360390) )),
            # q3
            DerivativeCoefficientRow{T,1,6}(SVector(T(-30443//162660),
                                                    T(311//32532),
                                                    T(0),
                                                    T(-11155//16266),
                                                    T(41287//32532),
                                                    T(-21999//54220) )),
            # q4 #TODO
            DerivativeCoefficientRow{T,1,7}(SVector(T(33311//107180),
                                                    T(-20229//21436),
                                                    T(485//1398),
                                                    T(0),
                                                    T(4147//21436),
                                                    T(25427//321540),
                                                    T(72//5359) )),
            # q5
            DerivativeCoefficientRow{T,1,8}(SVector(T(-16863//78770),
                                                    T(24337//31508),
                                                    T(-41287//47262),
                                                    T(-4147//15754),
                                                    T(0),
                                                    T(342523//472620),
                                                    T(-1296//7877),
                                                    T(144//7877) )),
            # q6
            DerivativeCoefficientRow{T,1,9}(SVector(T(15025//525612),
                                                    T(-36661//262806),
                                                    T(21999//87602),
                                                    T(-25427//262806),
                                                    T(-342523//525612),
                                                    T(0),
                                                    T(32400//43801),
                                                    T(-6480//43801),
                                                    T(720//43801) ))
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(3//4), T(-3//20), T(1//60))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(13649//43200),
                                T(12013//8640),
                                T(2711//4320),
                                T(5359//4320),
                                T(7877//8640),
                                T(43801//43200) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    elseif order == 8
        left_boundary = (
            # q1
            DerivativeCoefficientRow{T,1,8}(SVector(T(-2540160//1498139),
                                                    T(5544277//5992556),
                                                    T(198794991//29962780),
                                                    T(-256916579//17977668),
                                                    T(20708767//1498139),
                                                    T(-41004357//5992556),
                                                    T(27390659//17977668),
                                                    T(-2323531//29962780) )),
            # q2
            DerivativeCoefficientRow{T,1,8}(SVector(T(-5544277//31004596),
                                                    T(0),
                                                    T(-85002381//22146140),
                                                    T(49607267//4429228),
                                                    T(-165990199//13287684),
                                                    T(7655859//1107307),
                                                    T(-7568311//4429228),
                                                    T(48319961//465068940) )),
            # q3
            DerivativeCoefficientRow{T,1,8}(SVector(T(-66264997//8719620),
                                                    T(9444709//415220),
                                                    T(0),
                                                    T(-20335981//249132),
                                                    T(32320879//249132),
                                                    T(-35518713//415220),
                                                    T(2502774//103805),
                                                    T(-3177073//1743924) )),
            # q4
            DerivativeCoefficientRow{T,1,8}(SVector(T(256916579//109619916),
                                                    T(-49607267//5219996),
                                                    T(61007943//5219996),
                                                    T(0),
                                                    T(-68748371//5219996),
                                                    T(65088123//5219996),
                                                    T(-66558305//15659988),
                                                    T(3870214//9134993) )),
            # q5
            DerivativeCoefficientRow{T,1,9}(SVector(T(-20708767//2096689),
                                                    T(165990199//3594324),
                                                    T(-96962637//1198108),
                                                    T(68748371//1198108),
                                                    T(0),
                                                    T(-27294549//1198108),
                                                    T(14054993//1198108),
                                                    T(-42678199//25160268),
                                                    T(-2592//299527) )),
            # q6
            DerivativeCoefficientRow{T,1,10}(SVector(T(13668119//8660148),
                                                     T(-850651//103097),
                                                     T(35518713//2061940),
                                                     T(-21696041//1237164),
                                                     T(9098183//1237164),
                                                     T(0),
                                                     T(-231661//412388),
                                                     T(7120007//43300740),
                                                     T(3072//103097),
                                                     T(-288//103097) )),
             # q7
             DerivativeCoefficientRow{T,1,11}(SVector(T(-27390659//56287644),
                                                      T(7568311//2680364),
                                                      T(-22524966//3350455),
                                                      T(66558305//8041092),
                                                      T(-14054993//2680364),
                                                      T(2084949//2680364),
                                                      T(0),
                                                      T(70710683//93812740),
                                                      T(-145152//670091),
                                                      T(27648//670091),
                                                      T(-2592//670091) )),
              # q8
              DerivativeCoefficientRow{T,1,12}(SVector(T(2323531//102554780),
                                                       T(-48319961//307664340),
                                                       T(9531219//20510956),
                                                       T(-3870214//5127739),
                                                       T(2246221//3238572),
                                                       T(-21360021//102554780),
                                                       T(-70710683//102554780),
                                                       T(0),
                                                       T(4064256//5127739),
                                                       T(-1016064//5127739),
                                                       T(193536//5127739),
                                                       T(-18144//5127739) )),
        )
        right_boundary = .- left_boundary
        upper_coef = SVector(T(4//5), T(-1//5), T(4//105), T(-1//280))
        central_coef = zero(T)
        lower_coef = -upper_coef
        left_weights = SVector( T(1498139//5080320),
                                T(1107307//725760),
                                T(20761//80640),
                                T(1304999//725760),
                                T(299527//725760),
                                T(103097//80640),
                                T(670091//725760),
                                T(5127739//5080320) )
        right_weights = left_weights
        left_boundary_derivatives = Tuple{}()
        right_boundary_derivatives = left_boundary_derivatives

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 1, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end

function second_derivative_coefficients(source::MattssonSvärdNordström2004, order::Int,
                                        T=Float64, mode=FastMode())
    if order == 2
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(1))
        central_coef = T(-2)
        lower_coef = upper_coef
        left_weights = SVector(T(1//2))
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,3}(SVector(T(-3//2),
                                                    T(2),
                                                    T(-1//2) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 4
        left_boundary = (
            #TODO:
            # d1
            DerivativeCoefficientRow{T,1,4}(SVector(T(2),
                                                    T(-5),
                                                    T(4),
                                                    T(-1) )),
            # d2
            DerivativeCoefficientRow{T,1,3}(SVector(T(1),
                                                    T(-2),
                                                    T(1) )),
            # d3
            DerivativeCoefficientRow{T,1,5}(SVector(T(-4//43),
                                                    T(59//43),
                                                    T(-110//43),
                                                    T(59//43),
                                                    T(-4//43) )),
            # d4
            DerivativeCoefficientRow{T,1,6}(SVector(T(-1//49),
                                                    T(0),
                                                    T(59//49),
                                                    T(-118//49),
                                                    T(64//49),
                                                    T(-4//49) ))
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(4//3), T(-1//12))
        central_coef = T(-5//2)
        lower_coef = upper_coef
        left_weights = SVector(T(17//48), T(59//48), T(43//48), T(49//48))
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,4}(SVector(T(-11/6),
                                                    T(3),
                                                    T(-3//2),
                                                    T(1//3) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 6
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,6}(SVector(T(114170//40947),
                                                    T(-438107//54596),
                                                    T(336409//40947),
                                                    T(-276997//81894),
                                                    T(3747//13649),
                                                    T(21035//163788) )),
            # d2
            DerivativeCoefficientRow{T,1,6}(SVector(T(6173//5860),
                                                    T(-2066//879),
                                                    T(3283//1758),
                                                    T(-303//293),
                                                    T(2111//3516),
                                                    T(-601//4395) )),
            # d3
            DerivativeCoefficientRow{T,1,6}(SVector(T(-52391//81330),
                                                    T(134603//32532),
                                                    T(-21982//2711),
                                                    T(112915//16266),
                                                    T(-46969//16266),
                                                    T(30409//54220) )),
            # d4
            DerivativeCoefficientRow{T,1,7}(SVector(T(68603//321540),
                                                    T(-12423//10718),
                                                    T(112915//32154),
                                                    T(-75934//16077),
                                                    T(53369//21436),
                                                    T(-54899//160770),
                                                    T(48//5359) )),
            # d5
            DerivativeCoefficientRow{T,1,8}(SVector(T(-7053//39385),
                                                    T(86551//94524),
                                                    T(-46969//23631),
                                                    T(53369//15754),
                                                    T(-87904//23631),
                                                    T(820271//472620),
                                                    T(-1296//7877),
                                                    T(96//7877) )),
            # d6
            DerivativeCoefficientRow{T,1,9}(SVector(T(21035//525612),
                                                    T(-24641//131403),
                                                    T(30409//87602),
                                                    T(-54899//131403),
                                                    T(820271//525612),
                                                    T(-117600//43801),
                                                    T(64800//43801),
                                                    T(-6480//43801),
                                                    T(480//43801) )),
        )
        right_boundary = left_boundary
        upper_coef = SVector(T(3//2), T(-3//20), T(1//90))
        central_coef = T(-49//18)
        lower_coef = upper_coef
        left_weights = SVector( T(13649//43200),
                                T(12013//8640),
                                T(2711//4320),
                                T(5359//4320),
                                T(7877//8640),
                                T(43801//43200) )
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,5}(SVector(T(-25//12),
                                                    T(4),
                                                    T(-3),
                                                    T(4//3),
                                                    T(-1//4) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    elseif order == 8
        left_boundary = (
            # d1
            DerivativeCoefficientRow{T,1,7}(SVector(T(4870382994799//1358976868290),
                                                    T(-893640087518//75498714905),
                                                    T(926594825119//60398971924),
                                                    T(-1315109406200//135897686829),
                                                    T(39126983272//15099742981),
                                                    T(12344491342//75498714905),
                                                    T(-451560522577//2717953736580) )),
            # d2
            DerivativeCoefficientRow{T,1,8}(SVector(T(333806012194//390619153855),
                                                    T(-154646272029//111605472530),
                                                    T(1168338040//33481641759),
                                                    T(82699112501//133926567036),
                                                    T(-171562838//11160547253),
                                                    T(-28244698346//167408208795),
                                                    T(11904122576//167408208795),
                                                    T(-2598164715//312495323084) )),
            # d3
            DerivativeCoefficientRow{T,1,8}(SVector(T(7838984095//52731029988),
                                                    T(1168338040//5649753213),
                                                    T(-88747895//144865467),
                                                    T(423587231//627750357),
                                                    T(-43205598281//22599012852),
                                                    T(4876378562//1883251071),
                                                    T(-5124426509//3766502142),
                                                    T(10496900965//39548272491) )),
            # d4
            DerivativeCoefficientRow{T,1,8}(SVector(T(-94978241528//828644350023),
                                                    T(82699112501//157837019052),
                                                    T(1270761693//13153084921),
                                                    T(-167389605005//118377764289),
                                                    T(48242560214//39459254763),
                                                    T(-31673996013//52612339684),
                                                    T(43556319241//118377764289),
                                                    T(-44430275135//552429566682) )),
            # d5
            DerivativeCoefficientRow{T,1,9}(SVector(T(1455067816//21132528431),
                                                    T(-171562838//3018932633),
                                                    T(-43205598281//36227191596),
                                                    T(48242560214//9056797899),
                                                    T(-52276055645//6037865266),
                                                    T(57521587238//9056797899),
                                                    T(-80321706377//36227191596),
                                                    T(8078087158//21132528431),
                                                    T(-1296//299527) )),
            # d6
            DerivativeCoefficientRow{T,1,10}(SVector(T(10881504334//327321118845),
                                                     T(-28244698346//140280479505),
                                                     T(4876378562//9352031967),
                                                     T(-10557998671//12469375956),
                                                     T(57521587238//28056095901),
                                                     T(-278531401019//93520319670),
                                                     T(73790130002//46760159835),
                                                     T(-137529995233//785570685228),
                                                     T(2048//103097),
                                                     T(-144//103097) )),
             # d7
             DerivativeCoefficientRow{T,1,11}(SVector(T(-135555328849//8509847458140),
                                                      T(11904122576//101307707835),
                                                      T(-5124426509//13507694378),
                                                      T(43556319241//60784624701),
                                                      T(-80321706377//81046166268),
                                                      T(73790130002//33769235945),
                                                      T(-950494905688//303923123505),
                                                      T(239073018673//141830790969),
                                                      T(-145152//670091),
                                                      T(18432//670091),
                                                      T(-1296//670091) )),
              # d8
              DerivativeCoefficientRow{T,1,12}(SVector(T(0),
                                                       T(-2598164715//206729925524),
                                                       T(10496900965//155047444143),
                                                       T(-44430275135//310094888286),
                                                       T(425162482//2720130599),
                                                       T(-137529995233//620189776572),
                                                       T(239073018673//155047444143),
                                                       T(-144648000000//51682481381),
                                                       T(8128512//5127739),
                                                       T(-1016064//5127739),
                                                       T(129024//5127739),
                                                       T(-9072//5127739) )),

        )
        right_boundary = left_boundary
        upper_coef = SVector(T(8//5), T(-1//5), T(8//315), T(-1//560))
        central_coef = T(-205//72)
        lower_coef = upper_coef
        left_weights = SVector( T(1498139//5080320),
                                T(1107307//725760),
                                T(20761//80640),
                                T(1304999//725760),
                                T(299527//725760),
                                T(103097//80640),
                                T(670091//725760),
                                T(5127739//5080320) )
        right_weights = left_weights
        left_boundary_derivatives = (
            DerivativeCoefficientRow{T,1,7}(SVector(T(-4723//2100),
                                                    T(839//175),
                                                    T(-157//35),
                                                    T(278//105),
                                                    T(-103//140),
                                                    T(-1//175),
                                                    T(6//175) )),
        )
        right_boundary_derivatives = (-left_boundary_derivatives[1], )

        DerivativeCoefficients(left_boundary, right_boundary,
                                left_boundary_derivatives, right_boundary_derivatives,
                                lower_coef, central_coef, upper_coef,
                                left_weights, right_weights, mode, 2, order, source)
    else
        throw(ArgumentError("Order $order not implemented/derived."))
    end
end


text:
"""
Delay coordinates embedding, `Dataset` structure and neighhborhoods.
Basic package used in the ecosystem of DynamicalSystems.jl.
"""
module DelayEmbeddings

include("dataset.jl")
include("reconstruction.jl")
include("various.jl")
include("neighborhoods.jl")

include("estimate_delay.jl")
include("estimate_dimension.jl")

end

ext: jl


text:
import Base.isnan
import Base.isempty

isnan(out::FBAOut)  = isnan(objval(out))
isempty(out::FBAOut) = isnan(out)


text:
using Test
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation

@testset "Integration Tests" begin
	@testset "M" begin
		@test Lar.M(0,0)==0.5
		@test Lar.M(1,0)==0.16666666666666666
		@test Lar.M(2,0)==0.08333333333333333
		@test Lar.M(3,0)==0.05
		@test Lar.M(1,1)==0.041666666666666685
		@test Lar.M(2,0)==0.08333333333333333
		@test Lar.M(2,1)==0.016666666666666663
		@test Lar.M(3,0)==0.05
		@test Lar.M(3,1)==0.008333333333333338
		@test Lar.M(3,2)==0.0023809523809523586
		@test Lar.M(3,3)==0.0008928571428571397
	end

	@testset "TT" begin
		tau=[0.0 1.0 0.0; 0.0 0.0 1.0; 0.0 0.0 0.0];
		@test Lar.TT(tau, 0,0,0)==0.5
		@test Lar.TT(tau, 1,0,0)==0.16666666666666666
		@test Lar.TT(tau, 1,1,0)==0.041666666666666685
		@test Lar.TT(tau, 1,1,1)==0.0
		@test Lar.TT(tau, 2,0,0)==0.08333333333333333
		@test Lar.TT(tau, 2,1,0)==0.016666666666666663
		@test Lar.TT(tau, 2,2,0)==0.005555555555555545
		@test Lar.TT(tau, 2,2,1)==0.0
		@test Lar.TT(tau, 2,2,2)==0.0
	end

	@testset "II" begin
		V = [0.0 1.0 0.0; 0.0 0.0 1.0; 0.0 0.0 0.0];
		FV = [[1,2,3]];
		P = V,FV;
		@test Lar.II(P, 0,0,0)==0.5
		@test Lar.II(P, 1,0,0)==0.16666666666666666
		@test Lar.II(P, 0,1,0)>=0.1666666666666666
		@test Lar.II(P, 0,0,1)==0.0
		@test Lar.II(P, 1,1,1)==0.0
	end

	@testset "III" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.III(P, 0,0,0)>0.166666666
		@test Lar.III(P, 0,0,0)<0.166666888
		@test Lar.III(P, 1,0,0)>0.041666666
		@test Lar.III(P, 1,0,0)<0.041666888
		@test Lar.III(P, 0,1,0)>0.041666666
		@test Lar.III(P, 0,1,0)<0.041666888
		@test Lar.III(P, 0,0,1)>0.041666666
		@test Lar.III(P, 0,0,1)<0.041666888
		@test Lar.III(P, 10,10,10)>1.3377e-11
		@test Lar.III(P, 10,10,10)<1.3388e-11
	end

	@testset "surface" begin
		V,FV = Lar.simplexGrid([1,1]);
		P = [V;[0 0 0 0]], FV
		@test Lar.surface(P)==1.0
		p = Lar.Struct([Lar.t(0.5,0.5,0), Lar.r(0,0,pi/4), P]);
		q = Lar.struct2lar(p);
		@test Lar.surface(q)>1.0000000
		@test Lar.surface(q)<1.0000222
	end

	@testset "volume" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.volume(P)>0.166666666
		@test Lar.volume(P)<0.166668888
	end

	@testset "firstMoment" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.firstMoment(P)[1]<0.0416667
		@test Lar.firstMoment(P)[1]>0.0416665

		@test Lar.firstMoment(P)[2]<0.0416667
		@test Lar.firstMoment(P)[2]>0.0416665

		@test Lar.firstMoment(P)[3]<0.0416667
		@test Lar.firstMoment(P)[3]>0.0416665

		@test abs(Lar.firstMoment(P)[1]-Lar.firstMoment(P)[2])<0.00001
		@test abs(Lar.firstMoment(P)[2]-Lar.firstMoment(P)[3])<0.00001
		@test abs(Lar.firstMoment(P)[3]-Lar.firstMoment(P)[1])<0.00001
	end

	@testset "secondMoment" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.secondMoment(P)[1]<0.0166666669
		@test Lar.secondMoment(P)[1]>0.0166666664

		@test Lar.secondMoment(P)[2]<0.0166666669
		@test Lar.secondMoment(P)[2]>0.0166666664

		@test Lar.secondMoment(P)[3]<0.0166666669
		@test Lar.secondMoment(P)[3]>0.0166666664

		@test abs(Lar.secondMoment(P)[1]-Lar.secondMoment(P)[2])<0.00001
		@test abs(Lar.secondMoment(P)[2]-Lar.secondMoment(P)[3])<0.00001
		@test abs(Lar.secondMoment(P)[3]-Lar.secondMoment(P)[1])<0.00001
	end

	@testset "inertiaProduct" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.inertiaProduct(P)[1]<0.00833666
		@test Lar.inertiaProduct(P)[1]>0.00833000

		@test Lar.inertiaProduct(P)[2]<0.00833666
		@test Lar.inertiaProduct(P)[2]>0.00833000

		@test Lar.inertiaProduct(P)[3]<0.00833666
		@test Lar.inertiaProduct(P)[3]>0.00833000
	end

	@testset "centroid" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.centroid(P)[1]<0.26
		@test Lar.centroid(P)[1]>0.24

		@test Lar.centroid(P)[2]<0.26
		@test Lar.centroid(P)[2]>0.24

		@test Lar.centroid(P)[3]<0.26
		@test Lar.centroid(P)[3]>0.24
	end

	@testset "inertiaMoment" begin
		V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];
		FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];
		P = V,FV;
		@test Lar.inertiaMoment(P)[1]<0.0333555
		@test Lar.inertiaMoment(P)[1]>0.0333111

		@test Lar.inertiaMoment(P)[2]<0.0333555
		@test Lar.inertiaMoment(P)[2]>0.0333111

		@test Lar.inertiaMoment(P)[3]<0.0333555
		@test Lar.inertiaMoment(P)[3]>0.0333111
	end

end


text:
export FspMatrix, matvec, matvec!

abstract type FspMatrix end

include("sparse/fspsparsematrix.jl")


text:
"""
    mutate_inputs(v::MutationVertex; kwargs...)

Apply the input size mutation from the given vertex.

In the context of a neural network, this means removing or adding columns or rows from a matrix.

Supplied `kwargs` will be passed on to computation.
"""
function mutate_inputs end

"""
    mutate_outputs(v::MutationVertex; kwargs...)

Apply the output size mutation from the given vertex.

In the context of a neural network, this means removing or adding columns or rows from a matrix.

Supplied `kwargs` will be passed on to computation.
"""
function mutate_outputs end

"""
    apply_mutation(v::AbstractVertex; kwargs...)

Apply the size mutation from the given vertex.

In the context of a neural network, this means removing or adding columns or rows from a matrix.

Supplied `kwargs` will be passed on to computation.
"""
function apply_mutation(v::AbstractVertex; kwargs...) end
function apply_mutation(v::MutationVertex; kwargs...)
    # Empty nin leads to stack overflow due to splatting below.
    !isempty(nin(v)) && mutate_inputs(v; kwargs...)
    mutate_outputs(v; kwargs...)
end
apply_mutation(g::CompGraph; kwargs...) = apply_mutation.(vertices(g); kwargs...)

## Input mutation. Basically just boilerplate-y traversal of the vertex composition hierachy
# until we hit the computation to mutate, then someone else will do the actual work
mutate_inputs(v::CompVertex, inputs...; kwargs...) = mutate_inputs(v.computation, inputs...; kwargs...)
mutate_inputs(v::AbstractVertex, inputs...; kwargs...) = mutate_inputs(base(v), inputs...; kwargs...)
mutate_inputs(v::AbstractVertex, s::IoChange; kwargs...) = mutate_inputs(v, in_inds(s)...; kwargs...)
mutate_inputs(v::AbstractVertex, s::IoIndices; kwargs...) = mutate_inputs(v, in_inds(s)...; kwargs...)
mutate_inputs(v::AbstractVertex, s::IoSize; kwargs...) = mutate_inputs(v, nin(s)...; kwargs...)

function mutate_inputs(v::MutationVertex; kwargs...)
    mutate_inputs(v, op(v); kwargs...)
    reset_in!(op(v))
end

# Noops
function mutate_inputs(f::Function, inputs...; kwargs...) end #Maybe bad idea as it hides errors
function mutate_inputs(v::AbstractVertex, s::NoOp; kwargs...) end

# Output mutations. Also just traversal of the vertex composition hierachy.
mutate_outputs(v::CompVertex, outputs; kwargs...) = mutate_outputs(v.computation, outputs; kwargs...)
mutate_outputs(v::AbstractVertex, outputs; kwargs...) = mutate_outputs(base(v), outputs; kwargs...)
mutate_outputs(v::AbstractVertex, s::IoChange; kwargs...) = mutate_outputs(v, out_inds(s); kwargs...)
mutate_outputs(v::AbstractVertex, s::IoIndices; kwargs...) = mutate_outputs(v, out_inds(s); kwargs...)
mutate_outputs(v::AbstractVertex, s::IoSize; kwargs...) = mutate_outputs(v, nout(s); kwargs...)

function mutate_outputs(v::MutationVertex; kwargs...)
    mutate_outputs(v, op(v); kwargs...)
    reset_out!(op(v))
end

# Noops
function mutate_outputs(f::Function, inputs...; kwargs...) end #Maybe bad idea as it hides errors
function mutate_outputs(v::AbstractVertex, s::NoOp; kwargs...) end


text:
using ControlSystems
using PredictiveControl
using Test

include( "../testUtils.jl" )

clqrUnstable  = createSampleLTISystem()
clqrStablized = createSampleLTISystem( usek = true )
clqrWithS     = createSampleLTISystem( uses = true, usek = true )

# These are all conditions that are not allowed
@test_throws DomainError conditionBound( clqrUnstable )
@test_throws DomainError conditionBound( clqrWithS )

# A test for a single-input pre-stabilized system
@test conditionBound( clqrStablized ) ≈ 9.986301 atol=1e-4

# Setup a test for a multi-input system (it means actual matrices are used in the computation)
A = [0.9 1.0;
     0.0 0.9]
B = [0.0 1.0;
     1.0 0.0]
C = [1.0 0.0;
     0.0 1.0]
D = [0.0 0.0;
     0.0 0.0]
Ts = 0.1
sys = StateSpace( A, B, C, D, Ts )

Q = [1.0 1.0;
     1.0 1.0]
R = [1.0 0.0
     0.0 1.0]

clqrTwoInputs = ConstrainedTimeInvariantLQR( sys, 5, Q, R, :Q )

@test conditionBound( clqrTwoInputs ) ≈ 12201.000000000648


text:
__precompile__(true)
module Bar
import Baz, Foo
this = "Bar"
end

max_stars_repo_name: nreigl/julia

max_issues_repo_name: BenChung/Julia-BTI

max_forks_repo_name: BenChung/Julia-BTI


text:
module SubspacePreservingSparsification

using LinearAlgebra
using SparseArrays
using SparseArrays: AbstractSparseMatrixCSC

include("math_functions.jl")
include("sparse_pattern.jl")
include("bin_pattern.jl")
include("compute.jl")

export pinv_qr
export sparsity_pattern
export bin_sparse_matrix!
export binned_minimization
export sparsify

end  # end of module


text:
# Autogenerated wrapper script for hidapi_jll for aarch64-linux-musl
export hidapi

using libusb_jll
JLLWrappers.@generate_wrapper_header("hidapi")
JLLWrappers.@declare_library_product(hidapi, "libhidapi-libusb.so.0")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        hidapi,
        "lib/libhidapi-libusb.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Dates
using ArgParse
using Printf
using Random
using GaussianProcesses

# before_all = Dates.now()

include("../btg.jl")
s = ArgParseSettings()
# The defaut setting: --test: multiple length scale, QMC
@add_arg_table! s begin
    "--test"
        help = "test or not"
        action = :store_true
    "--single"
        help = "use single length scale or not"
        action = :store_true
    "--sparse"
        help = "use SparseGrid or not"
        action = :store_true
    "--validate"
        help = "do cross validation or not"
        action = :store_true
    "--fast"
        help = "use fast validation or not"
        action = :store_true
    "--GP"
        help = "test GP model"
        action = :store_true
    "--logGP"
        help = "test log-GP model"
        action = :store_true
    "--posc"
        help = "another option with an argument"
        arg_type = Int
        default = 7
    "--singletest"
        help = "write log to single test"
        action = :store_true
    "--ntrain"
        help = "another option with an argument"
        arg_type = Int
        default = 200
    "--ntest"
        help = "another option with an argument"
        arg_type = Int
        default = 100
end
parsed_args = parse_args(ARGS, s)
# load abalone data
df = DataFrame(CSV.File("../datasets/abalone.csv", header = 0))
data = convert(Matrix, df[:,2:8]) #length, diameter, height, whole weight, shucked weight, viscera weight, shell weight
target = convert(Array, df[:, 9]) #age
# shuffle data
randseed = 1234; rng = MersenneTwister(randseed)
ind_shuffle = randperm(rng, size(data, 1))
data = data[ind_shuffle, :]
target = target[ind_shuffle]
# training set
id_train = 1:parsed_args["ntrain"]; posx = 1:7; posc = 1:parsed_args["posc"]
n_train = length(id_train)
x = data[id_train, posx]
Fx = data[id_train, posc]
y = float(target[id_train])
ymax_train = maximum(y)
y ./= ymax_train
trainingData0 = trainingData(x, Fx, y)
d = getDimension(trainingData0); n = getNumPts(trainingData0); p = getCovDimension(trainingData0)

#parameter setting
# myquadtype = parsed_args["sparse"] ? ["SparseCarlo", "SparseCarlo"] : ["QuasiMonteCarlo", "QuasiMonteCarlo"]
myquadtype = ["Gaussian", "Gaussian"]
rangeλ = [-1.5 1.]
rangeθs = [500. 1000]
rangeθm = repeat(rangeθs, d, 1)
# rangeθ = parsed_args["single"] ? rangeθs : rangeθm
rangeθ = rangeθs
# build btg model
btg0 = btg(trainingData0, rangeθ, rangeλ; quadtype = myquadtype)
(pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw) = solve(btg0);

####################################
############### Test ###############
####################################
if parsed_args["test"]
    @info "Start Test..."
    before = Dates.now()
    id_test = 1001:(parsed_args["ntest"]+1000)
    n_test = length(id_test)
    id_fail = []
    id_nonproper = []
    x_test = data[id_test, posx]
    Fx_test = data[id_test, posc]
    y_test_true = target[id_test]
    count_test = 0
    error_abs = 0.
    error_sq = 0.
    nlpd = 0.
    error_abs_set = zeros(n_test)
    nlpd_set = zeros(n_test)
    for i in 1:n_test
        global error_abs, error_sq, nlpd, count_test
        # mod(i, 20) == 0 ? (@info i) : nothing
        # @info "i" i
        x_test_i = reshape(x_test[i, :], 1, length(posx))
        Fx_test_i = reshape(Fx_test[i, :], 1, length(posc))
        try
            pdf_test_i, cdf_test_i, dpdf_test_i, quantbound_test_i, support_test_i = pre_process(x_test_i, Fx_test_i, pdf0_raw, cdf0_raw, dpdf0_raw, quantInfo0_raw)
            y_test_i_true = y_test_true[i]
            median_test_i = ymax_train * quantile(cdf_test_i, quantbound_test_i, support_test_i)[1]
            # @info "True, median " y_test_i_true, median_test_i
            try
                CI_test_i = ymax_train .* credible_interval(cdf_test_i, quantbound_test_i, support_test_i; mode=:equal, wp=.95)[1]
                count_test += (y_test_i_true >= CI_test_i[1])&&(y_test_i_true <= CI_test_i[2]) ? 1 : 0
                # @info "CI" CI_test_i
            catch err
                append!(id_fail, i)
            end
            error_abs += abs(y_test_i_true - median_test_i)
            error_sq += (y_test_i_true - median_test_i)^2
            nlpd += log(pdf_test_i(y_test_i_true/ymax_train))
            error_abs_set[i] =  abs(y_test_i_true - median_test_i)
            nlpd_set[i] = log(pdf_test_i(y_test_i_true/ymax_train))
        # @info "Count, id_fail" count_test, id_fail
        catch err
            append!(id_nonproper, i)
        end
    end
    count_test /= n_test - length(id_fail) - length(id_nonproper)
    error_abs  /= n_test - length(id_nonproper)
    error_sq   /= n_test - length(id_nonproper)
    nlpd       /= -n_test - length(id_nonproper)
    after = Dates.now()
    elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)

    io = open("profile_abalone_errorhist.txt", "a")
    write(io, "\n$(Dates.now()), randseed: $randseed \n")
    write(io, "Data set: Abalone
    id_train:  $id_train;  id_test:  $id_test;   posx: $posx;   posc: $posc\n")
    write(io, "BTG model:
            $myquadtype  ;  rangeλ: $rangeλ;   rangeθ: $rangeθs (single length-scale: $(parsed_args["single"])) \n")
    write(io, "Absolute error history \n $error_abs_set \n")
    write(io, "Negative log predictive density history \n $nlpd_set \n")
    close(io)
    if parsed_args["GP"]
        global error_abs_GP, error_sq_GP, CI_test_GP, count_test_GP, nlpd_GP
        # training set
        x = data[id_train, posx]'
        y = float(target[id_train])
        # build and fit a GP
        mymean = MeanLin(zeros(d))
        # mymean = MeanZero()
        kern = SE(zeros(d),0.0)
        gp = GP(x, y, mymean, kern)
        optimize!(gp)
        # predict
        x_test = data[id_test, posx]'
        y_test_true = target[id_test]
        μ, σ² = predict_y(gp, x_test); stdv = sqrt.(σ²)
        error_GP = abs.(μ .- y_test_true)
        error_abs_GP = mean(error_GP)
        error_sq_GP = mean(error_GP.^2)
        CI_test_GP = hcat(-1.96.*stdv .+ μ, 1.96.*stdv .+ μ)
        count_test_GP = sum((y_test_true .>= CI_test_GP[:, 1]) .* (y_test_true .<= CI_test_GP[:,2]))/n_test
        nlpd_GP = -mean(log.(pdf.(Normal(), (y_test_true.-μ)./stdv)./stdv))
    end

    if parsed_args["logGP"]
        global error_abs_logGP, error_sq_logGP, CI_test_logGP, count_test_logGP, nlpd_logGP
        x = data[id_train, posx]'
        y = float(target[id_train])
        trans = BoxCox()
        g_fixed(x) = trans(x, 0.); dg(x) = partialx(trans, x, 0.)
        invg(x) = inverse(trans, x, 0.)
        gy = g_fixed.(y)
        # build and fit a GP
        mymean = MeanLin(zeros(d))
        # mymean = MeanZero()
        kern = SE(zeros(d),0.0)
        loggp = GP(x, gy, mymean, kern)
        optimize!(loggp)
        # predict
        x_test = data[id_test, posx]'
        y_test_true = target[id_test]
        μ, σ² = predict_y(loggp, x_test); stdv = sqrt.(σ²)
        CI_test_logGP = invg.(hcat(-1.96.*stdv .+ μ, 1.96.*stdv .+ μ))
        count_test_logGP = sum((y_test_true .>= CI_test_logGP[:, 1]) .* (y_test_true .<= CI_test_logGP[:,2]))/n_test
        y_pred = invg.(μ)
        error_logGP = abs.(y_pred .- y_test_true)
        error_abs_logGP = mean(error_logGP)
        error_sq_logGP = mean(error_logGP.^2)
        nlpd_logGP = -mean(log.( dg.(y_test_true) .* pdf.(Normal(), (g_fixed.(y_test_true).-μ)./stdv) ./stdv ))
    end

    @info "Average time for single test: " elapsedmin/n_test
    io1 = open("Profile_abalone_results.txt", "a")
    write(io1, "\n$(Dates.now()), randseed: $randseed \n" )
    write(io1, "Data set: Abalone
        id_train:  $id_train;  id_test:  $id_test;   posx: $posx;   posc: $posc\n")
    write(io1, "BTG model:
        $myquadtype  ;  rangeλ: $rangeλ;   rangeθ: $rangeθs (single length-scale: $(parsed_args["single"])) \n")
    if parsed_args["GP"] && parsed_args["logGP"]
        write(io1, "Compare test results: ")
        write(io1, "                               BTG               GP               logGP
        credible intervel accuracy percentage:   $(@sprintf("%11.8f", count_test))       $(@sprintf("%11.8f", count_test_GP))       $(@sprintf("%11.8f", count_test_logGP))
        mean absolute error:                     $(@sprintf("%11.8f", error_abs))       $(@sprintf("%11.8f", error_abs_GP))       $(@sprintf("%11.8f", error_abs_logGP))
        mean squared error:                      $(@sprintf("%11.8f", error_sq))       $(@sprintf("%11.8f", error_sq_GP))       $(@sprintf("%11.8f", error_sq_logGP))
        mean negative log predictive density:    $(@sprintf("%11.8f", nlpd))       $(@sprintf("%11.8f", nlpd_GP))       $(@sprintf("%11.8f", nlpd_logGP))
        Time cost by prediction: $elapsedmin
        BTG: Failed index in credible intervel:   $id_fail
        BTG: Failed index in pdf computation:     $id_nonproper\n")
    else
        write(io1, "BTG test results:
        credible intervel accuracy percentage:   $(@sprintf("%11.8f", count_test))
        mean absolute error:                     $(@sprintf("%11.8f", error_abs))
        mean squared error:                      $(@sprintf("%11.8f", error_sq))
        mean negative log predictive density:    $(@sprintf("%11.8f", nlpd))
        Time cost by prediction: $elapsedmin
        Failed index in credible intervel:       $id_fail
        BTG: Failed index in pdf computation:     $id_nonproper\n")
    end
    close(io1)

end

####################################
############ Validation ############
####################################
function lootd(td::AbstractTrainingData, i::Int64)
    x = getPosition(td)
    Fx = getCovariates(td)
    z = getLabel(td)
    x_minus_i = x[[1:i-1;i+1:end], :]
    Fx_minus_i = Fx[[1:i-1;i+1:end], :]
    z_minus_i = z[[1:i-1;i+1:end]]
    x_i = x[i:i, :]
    Fx_i = Fx[i:i, :]
    z_i = z[i:i, :]
    return trainingData(x_minus_i, Fx_minus_i, z_minus_i), x_i, Fx_i, z_i
end

mutable struct validation_plot_buffer
    ygrid::Array{Float64, 1}
    median::Float64
    CI::Array{Float64, 1}
    pdf_ytrue::Float64
    pdf_median_i::Float64
    pdf_CI::Array{Float64, 1}
    function validation_plot_buffer(ygrid_i, median_i, CI_i, pdf_ytrue_i, pdf_median_i, pdf_CI_i)
        return new(ygrid_i, median_i, CI_i, pdf_ytrue_i, pdf_median_i, pdf_CI_i)
    end
end

function unpack(a::validation_plot_buffer)
    return (a.ygrid, a.median, a.CI, a.pdf_ytrue, a.pdf_median_i, a.pdf_CI)
end

if parsed_args["validate"]
    @info "Start validation"
    before = Dates.now()
    id_fail_val = []
    id_nonproper_val = []
    count_val = 0
    error_abs_val = 0.
    error_sq_val = 0.
    nlpd_val = 0.
    # prepare space for plotting
    n_row = 6; n_col = 8
    xgrid = range(.01, stop=1.3, length=100)
    validation_plot_buffer_dict = Dict{T where T<:Int64, validation_plot_buffer}()
    for i = 1:n_train
        global count_val, error_abs_val, error_sq_val, nlpd_val
        mod(i, 10) == 0 ? (@info i) : nothing
        if parsed_args["fast"]
            (pdf_val_i_raw, cdf_val_i_raw, dpdf_val_i_raw, quantInfo_val_i_raw) = solve(btg0, validate = i)
        else
            (trainingdata_minus_i, x_i, Fx_i, z_i) = lootd(trainingData0, i)
            btg_val_naive = btg(trainingdata_minus_i, rangeθ, rangeλ; quadtype = myquadtype)
            (pdf_val_i_raw, cdf_val_i_raw, dpdf_val_i_raw, quantInfo_val_i_raw) = solve(btg_val_naive)
        end
        try
            x_i = x[i:i, :]; Fx_i = Fx[i:i, :]
            pdf_val_i, cdf_val_i, dpdf_val_i, quantbound_val_i, support_val_i = pre_process(x_i, Fx_i, pdf_val_i_raw, cdf_val_i_raw, dpdf_val_i_raw, quantInfo_val_i_raw)
            median_val_i = quantile(cdf_val_i, quantbound_val_i, support_val_i)[1]
            y_val_i_true = getLabel(btg0.trainingData)[i]
            pdf_ytrue_i = pdf_val_i(y_val_i_true)
            try
                CI_val_i = credible_interval(cdf_val_i, quantbound_val_i, support_val_i; mode=:equal, wp=.95)[1]
                count_val += (y_val_i_true >= CI_val_i[1])&&(y_val_i_true <= CI_val_i[2]) ? 1 : 0
                if length(validation_plot_buffer_dict) < n_row*n_col
                    ygrid_i = pdf_val_i.(xgrid)
                    pdf_CI_i = pdf_val_i.(CI_val_i)
                    pdf_median_i = pdf_val_i.(median_val_i)
                    push!(validation_plot_buffer_dict, i => validation_plot_buffer(ygrid_i, median_i, CI_val_i, pdf_ytrue_i, pdf_median_i, pdf_CI_i))
                end
            catch err
                append!(id_fail_val, i)
            end
            error_abs_val += abs(y_val_i_true - median_val_i)
            error_sq_val += (y_val_i_true - median_val_i)^2
            nlpd_val -= log(pdf_ytrue_i)
        # @info "Count, id_fail" count_test, id_fail
        catch err
            append!(id_nonproper_val, i)
        end
    end
    count_val      /= n_train - length(id_fail_val) - length(id_nonproper_val)
    error_abs_val  /= n_train - length(id_nonproper_val)
    error_sq_val   /= n_train - length(id_nonproper_val)
    nlpd_val       /= n_train - length(id_nonproper_val)
    after = Dates.now()
    elapsedmin = round(((after - before) / Millisecond(1000))/60, digits=5)

    # save results
    io2 = open("Exp_abalone_LOOCV.txt", "a")
    write(io2, "\n$(Dates.now()), randseed: $randseed \n" )
    write(io2, "Data set: Abalone
        id_train:  $id_train;   posx: $posx;   posc: $posc\n")
    write(io2, "BTG model:
        $myquadtype  ;  rangeλ: $rangeλ;   rangeθ: $rangeθs (single length-scale: $(parsed_args["single"]))
        Fast LOOCV: $(parsed_args["fast"]) \n")
    write(io2, "LOOCV on training set results:
    credible intervel accuracy percentage:   $(@sprintf("%11.8f", count_val))
    mean absolute error:                     $(@sprintf("%11.8f", error_abs_val))
    mean squared error:                      $(@sprintf("%11.8f", error_sq_val))
    mean negative log predictive density:    $(@sprintf("%11.8f", nlpd_val))
    Time validation took: $elapsedmin
    Failed index in credible intervel:       $id_fail_val
    BTG: Failed index in pdf computation:     $id_nonproper_val\n")
    close(io2)

    # Plot the first 48
    PyPlot.close("all") #close existing windows
    plt, axs = PyPlot.subplots(n_row, n_col)
    PyPlot.suptitle("Partial Plots of LOOCV", fontsize=10)
    id = 0
    for i in keys(validation_plot_buffer_dict)
        id += 1
        global id
        ind1 = Int64(ceil(id/n_col))
        ind2 = Int64(id - n_col*(floor((id-.1)/n_col)))
        y_val_i_true = getLabel(btg0.trainingData)[i]
        ygrid_i, median_i, CI_val_i, pdf_ytrue_i, pdf_median_i, pdf_CI_i = unpack(validation_plot_buffer_dict[i])
        axs[ind1, ind2].plot(xgrid, ygrid_i, linewidth = 1.0, linestyle = "-")
        # axs[ind1, ind2].plot(xgrid, cdf1.(xgrid), color = "orange", linewidth = 1.0, linestyle = "-")
        axs[ind1, ind2].vlines(y_val_i_true, 0, pdf_ytrue_i, label = "true value")
        axs[ind1, ind2].vlines(median_val_i, 0, pdf_median_i, label = "median")
        CI_id = (xgrid .> CI_val_i[1]) .* (xgrid .< CI_val_i[2])
        CI_xrange = vcat(CI_val_i[1], xgrid[CI_id], CI_val_i[2])
        CI_yrange = vcat(pdf_CI_i[1], ygrid_i[CI_id], pdf_CI_i[2]) # utilize previous evaluation results
        axs[ind1, ind2].fill_between(CI_xrange, 0, CI_yrange, alpha = 0.3, label = "95% confidence interval")
    end
    # Hide x labels and tick labels for top plots and y ticks for right plots.
    for ax in axs
        ax.label_outer()
    end
    PyPlot.savefig("figure/exp_abalone_$(id_train[1])_$(id_train[2])_$(myquadtype[1])$(myquadtype[2])_posc$(posc[2])_rθ_$(Int(rangeθs[1]))_$(Int(rangeθs[2]))_rλ_$(Int(rangeλ[1]))_$(Int(rangeλ[2]))_$(parsed_args["fast"]).pdf")

end

# after_all = Dates.now()
# elapsedmin_all = round(((after_all - before_all) / Millisecond(1000))/60, digits=5)
# @info "Total time" elapsedmin_all

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg