text:
"""
    DerivativeCoefficients

The coefficients of a derivative operator on a nonperiodic grid.
"""
struct DerivativeCoefficients{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                              LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients} <: AbstractDerivativeCoefficients{T}
    # coefficients defining the operator and its action
    left_boundary::LeftBoundary
    right_boundary::RightBoundary
    left_boundary_derivatives::LeftBoundaryDerivatives
    right_boundary_derivatives::RightBoundaryDerivatives
    lower_coef::SVector{LowerOffset, T}
    central_coef::T
    upper_coef::SVector{UpperOffset, T}
    left_weights::SVector{LeftWidth, T}
    right_weights::SVector{RightWidth, T}
    mode::ExecutionMode
    # corresponding orders etc.
    derivative_order::Int
    accuracy_order  ::Int
    symmetric       ::Bool
    source_of_coefficients::SourceOfCoefficients

    function DerivativeCoefficients(left_boundary::LeftBoundary, right_boundary::RightBoundary,
                                    left_boundary_derivatives::LeftBoundaryDerivatives,
                                    right_boundary_derivatives::RightBoundaryDerivatives,
                                    lower_coef::SVector{LowerOffset, T}, central_coef::T, upper_coef::SVector{UpperOffset, T},
                                    left_weights::SVector{LeftWidth, T}, right_weights::SVector{RightWidth, T},
                                    mode::ExecutionMode, derivative_order::Int, accuracy_order::Int,
                                    source_of_coefficients::SourceOfCoefficients) where {T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                                                                                LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients}
        @argcheck length(left_boundary) == LeftWidth
        @argcheck length(right_boundary) == RightWidth
        @argcheck length(left_boundary_derivatives) == length(right_boundary_derivatives)
        symmetric = LowerOffset == UpperOffset
        if symmetric
            @inbounds for i in Base.OneTo(LowerOffset)
                symmetric = symmetric && lower_coef[i] == upper_coef[i]
            end
        end
        symmetric = symmetric && length(left_boundary[end]) == length(left_boundary)
        #TODO: check boundary coefficients
        if derivative_order-1 != length(left_boundary_derivatives)
            warn("Derivative coefficients of degree $derivative_order should provide $(derivative_order-1) boundary derivatives.")
        end
        new{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients}(
            left_boundary, right_boundary, left_boundary_derivatives, right_boundary_derivatives,
            lower_coef, central_coef, upper_coef, left_weights, right_weights,
            mode, derivative_order, accuracy_order, symmetric, source_of_coefficients)
    end
end

@inline source_of_coefficients(coefficients::AbstractDerivativeCoefficients) = coefficients.source_of_coefficients

function Base.show(io::IO, coefficients::DerivativeCoefficients)
    if derivative_order(coefficients) == 1
        print(io, "Coefficients of the first-derivative operator")
    elseif  derivative_order(coefficients) == 2
        print(io, "Coefficients of the second-derivative operator")
    elseif  derivative_order(coefficients) == 3
        print(io, "Coefficients of the third-derivative operator")
    else
        print(io, "Coefficients of the ", derivative_order(coefficients),
              "-derivative operator")
    end
    print(io, " of order ", accuracy_order(coefficients), " of ")
    print(io, source_of_coefficients(coefficients))
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α, β)
    @unpack left_boundary, right_boundary, lower_coef, central_coef, upper_coef, mode = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > length(lower_coef)+length(upper_coef) DimensionMismatch
        @argcheck length(u) > length(left_boundary) + length(right_boundary) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, left_boundary, right_boundary, u, α, β)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β, static_length(left_boundary), static_length(right_boundary), mode)
end

# Compute `α*D*u` and store the result in `dest`.
function mul!(dest::AbstractVector, coefficients::DerivativeCoefficients, u::AbstractVector, α)
    @unpack left_boundary, right_boundary, lower_coef, central_coef, upper_coef, mode = coefficients

    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > length(lower_coef)+length(upper_coef) DimensionMismatch
        @argcheck length(u) > length(coefficients.left_boundary) + length(coefficients.right_boundary) DimensionMismatch
    end

    convolve_boundary_coefficients!(dest, left_boundary, right_boundary, u, α)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, static_length(left_boundary), static_length(right_boundary), mode)
end

"""
    DerivativeCoefficientRow{T,Start,Length}

A struct representing a row in the boundary block of an SBP derivative operator
with scalar type `T`.
"""
struct DerivativeCoefficientRow{T,Start,Length}
    coef::SVector{Length, T}
end

function Base.length(::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Length}
    Length
end

function offset(::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Length}
    Start
end

function Base.:-(coef_row::DerivativeCoefficientRow{T,Start,Length}) where {T,Start,Length}
    DerivativeCoefficientRow{T,Start,Length}(-coef_row.coef)
end

function widening_plus(row1::SVector{Length1,T}, row2::SVector{Length2,T}) where {T,Length1,Length2}
    Length = max(Length1, Length2)
    row = SVector{Length,T}(ntuple(i -> zero(T), Length))
    for i in 1:Length1
        row = Base.setindex(row, row[i] + row1[i], i)
    end
    for i in 1:Length2
        row = Base.setindex(row, row[i] + row2[i], i)
    end
    row
end

function Base.:+(coef_row1::DerivativeCoefficientRow{T,Start1,Length1}, coef_row2::DerivativeCoefficientRow{T,Start2,Length2}) where {T,Start1,Length1,Start2,Length2}
    End = max(Start1+Length1, Start2+Length2)
    Start = min(Start1, Start2)
    Length = End - Start
    row = SVector{Length,T}(ntuple(i -> zero(T), Length))
    for i in 1:Length1
        j = i-Start1+Start
        row = Base.setindex(row, row[j] + coef_row1.coef[i], j)
    end
    for i in 1:Length2
        j = i-Start2+Start
        row = Base.setindex(row, row[j] + coef_row2.coef[i], j)
    end
    DerivativeCoefficientRow{T,Start,Length}(row)
end

function Base.:/(coef_row::DerivativeCoefficientRow{T,Start,Length}, α::Number) where {T,Start,Length}
    DerivativeCoefficientRow{T,Start,Length}(coef_row.coef / α)
end

@inline function convolve_left_row(coef_row::DerivativeCoefficientRow{T,Start,Length}, u) where {T,Start,Length}
    @inbounds tmp = coef_row.coef[1]*u[Start]
    @inbounds for i in 2:Length
        tmp += coef_row.coef[i]*u[Start+i-1]
    end
    tmp
end

@inline function convolve_right_row(coef_row::DerivativeCoefficientRow{T,Start,Length}, u) where {T,Start,Length}
    @inbounds tmp = coef_row.coef[1]*u[end-Start+1]
    @inbounds for i in 2:Length
        tmp += coef_row.coef[i]*u[end-(Start+i-2)]
    end
    tmp
end

@unroll function convolve_boundary_coefficients!(dest::AbstractVector, left_boundary, right_boundary, u::AbstractVector, α, β)
    @unroll for i in 1:length(left_boundary)
        tmp = convolve_left_row(left_boundary[i], u)
        @inbounds dest[i] = β*dest[i] + α*tmp
    end
    @unroll for i in 1:length(right_boundary)
        tmp = convolve_right_row(right_boundary[i], u)
        @inbounds dest[end-i+1] = β*dest[end-i+1] + α*tmp
    end
end

@unroll function convolve_boundary_coefficients!(dest::AbstractVector, left_boundary, right_boundary, u::AbstractVector, α)
    @unroll for i in 1:length(left_boundary)
        tmp = convolve_left_row(left_boundary[i], u)
        @inbounds dest[i] = α*tmp
    end
    @unroll for i in 1:length(right_boundary)
        tmp = convolve_right_row(right_boundary[i], u)
        @inbounds dest[end-i+1] = α*tmp
    end
end

@generated function convolve_interior_coefficients!(
        dest::AbstractVector,
        lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
        u::AbstractVector, α, β,
        ::StaticInt{left_boundary_width}, ::StaticInt{right_boundary_width}, mode) where {LowerOffset, UpperOffset, left_boundary_width, right_boundary_width}
    if LowerOffset > 0
        ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] )
        for j in LowerOffset-1:-1:1
            ex = :( $ex + lower_coef[$j]*u[i-$j] )
        end
        ex = :( $ex + central_coef*u[i] )
    else
        ex = :( central_coef*u[i] )
    end
    for j in 1:UpperOffset
        ex = :( $ex + upper_coef[$j]*u[i+$j] )
    end

    if mode <: ThreadedMode
        quote
            Base.@_inline_meta
            @tturbo for i in (left_boundary_width+1):(length(dest)-right_boundary_width)
                dest[i] = β*dest[i] + α*$ex
            end
        end
    elseif mode <: SafeMode
        quote
            Base.@_inline_meta
            @inbounds for i in (left_boundary_width+1):(length(dest)-right_boundary_width)
                dest[i] = β*dest[i] + α*$ex
            end
        end
    else
        quote
            Base.@_inline_meta
            @turbo for i in (left_boundary_width+1):(length(dest)-right_boundary_width)
                dest[i] = β*dest[i] + α*$ex
            end
        end
    end
end

@generated function convolve_interior_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset}, u::AbstractVector, α, ::StaticInt{left_boundary_width}, ::StaticInt{right_boundary_width}, mode) where {LowerOffset, UpperOffset, left_boundary_width, right_boundary_width}
    if LowerOffset > 0
        ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] )
        for j in LowerOffset-1:-1:1
            ex = :( $ex + lower_coef[$j]*u[i-$j] )
        end
        ex = :( $ex + central_coef*u[i] )
    else
        ex = :( central_coef*u[i] )
    end
    for j in 1:UpperOffset
        ex = :( $ex + upper_coef[$j]*u[i+$j] )
    end

    if mode <: ThreadedMode
        quote
            Base.@_inline_meta
            @tturbo for i in (firstindex(dest) + $left_boundary_width):(lastindex(dest) - $right_boundary_width)
                dest[i] = α*$ex
            end
        end
    elseif mode <: SafeMode
        quote
            Base.@_inline_meta
            @simd ivdep for i in (firstindex(dest) + $left_boundary_width):(lastindex(dest) - $right_boundary_width)
                @inbounds dest[i] = α*$ex
            end
        end
    else
        quote
            Base.@_inline_meta
            @turbo for i in (firstindex(dest) + $left_boundary_width):(lastindex(dest) - $right_boundary_width)
                dest[i] = α*$ex
            end
        end
    end
end

@static if VERSION >= v"1.6" # `reinterpret(reshape, ...)` was introduced in Julia v1.6
    # Specialized for vectors of `StaticVector`s
    @generated function convolve_interior_coefficients!(_dest::AbstractVector{<:StaticVector{N,T1}}, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset}, _u::AbstractVector{<:StaticVector{N,T2}}, α, β, ::StaticInt{left_boundary_width}, ::StaticInt{right_boundary_width}, mode) where {LowerOffset, UpperOffset, N, T1, T2, left_boundary_width, right_boundary_width}
        if LowerOffset > 0
            ex = :( lower_coef[$LowerOffset]*u[v, i-$LowerOffset] )
            for j in LowerOffset-1:-1:1
                ex = :( $ex + lower_coef[$j]*u[v, i-$j] )
            end
            ex = :( $ex + central_coef*u[v, i] )
        else
            ex = :( central_coef*u[v, i] )
        end
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[v, i+$j] )
        end

        if N == 1
            ex_dest = :( reshape(reinterpret($T1, _dest), 1, length(_dest)) )
            ex_u    = :( reshape(reinterpret($T2, _u),    1, length(_u)) )
        else
            ex_dest = :( reinterpret(reshape, $T1, _dest) )
            ex_u    = :( reinterpret(reshape, $T2, _u) )
        end

        if mode <: ThreadedMode
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @tturbo for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = β*dest[v, i] + α*$ex
                    end
                end
            end
        elseif mode <: SafeMode
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @inbounds for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = β*dest[v, i] + α*$ex
                    end
                end
            end
        else
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @turbo for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = β*dest[v, i] + α*$ex
                    end
                end
            end
        end
    end

    @generated function convolve_interior_coefficients!(_dest::AbstractVector{<:StaticVector{N,T1}}, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset}, _u::AbstractVector{<:StaticVector{N,T2}}, α, ::StaticInt{left_boundary_width}, ::StaticInt{right_boundary_width}, mode) where {LowerOffset, UpperOffset, N, T1, T2, left_boundary_width, right_boundary_width}
        if LowerOffset > 0
            ex = :( lower_coef[$LowerOffset]*u[v, i-$LowerOffset] )
            for j in LowerOffset-1:-1:1
                ex = :( $ex + lower_coef[$j]*u[v, i-$j] )
            end
            ex = :( $ex + central_coef*u[v, i] )
        else
            ex = :( central_coef*u[v, i] )
        end
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[v, i+$j] )
        end

        if N == 1
            ex_dest = :( reshape(reinterpret($T1, _dest), 1, length(_dest)) )
            ex_u    = :( reshape(reinterpret($T2, _u),    1, length(_u)) )
        else
            ex_dest = :( reinterpret(reshape, $T1, _dest) )
            ex_u    = :( reinterpret(reshape, $T2, _u) )
        end

        if mode <: ThreadedMode
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @tturbo for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = α*$ex
                    end
                end
            end
        elseif mode <: SafeMode
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @inbounds for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = α*$ex
                    end
                end
            end
        else
            quote
                Base.@_inline_meta
                dest = $ex_dest
                u    = $ex_u
                @turbo for i in (firstindex(dest, 2) + left_boundary_width):(lastindex(dest, 2) - right_boundary_width)
                    for v in LoopVectorization.indices((dest, u), (1, 1))
                        dest[v, i] = α*$ex
                    end
                end
            end
        end
    end
end

"""
    DerivativeOperator

A derivative operator on a nonperiodic finite difference grid.
See [`derivative_operator`](@ref).
"""
struct DerivativeOperator{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                          LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients,Grid} <: AbstractNonperiodicDerivativeOperator{T}
    coefficients::DerivativeCoefficients{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                                         LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients}
    grid::Grid
    Δx::T
    factor::T

    function DerivativeOperator(coefficients::DerivativeCoefficients{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                                                                     LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients},
                                grid::Grid) where {T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,
                                                   LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients,Grid}
        @argcheck length(grid) > LowerOffset+UpperOffset DimensionMismatch
        @argcheck length(grid) > length(coefficients.left_boundary) + length(coefficients.right_boundary) DimensionMismatch

        Δx = step(grid)
        factor = inv(Δx^derivative_order(coefficients))
        new{T,LeftBoundary,RightBoundary,LeftBoundaryDerivatives,RightBoundaryDerivatives,LowerOffset,UpperOffset,LeftWidth,RightWidth,ExecutionMode,SourceOfCoefficients,Grid}(
            coefficients, grid, Δx, factor)
    end
end

@inline source_of_coefficients(D::DerivativeOperator) = source_of_coefficients(D.coefficients)

function Base.show(io::IO, D::DerivativeOperator)
    if derivative_order(D) == 1
        print(io, "SBP first-derivative operator")
    elseif  derivative_order(D) == 2
        print(io, "SBP second-derivative operator")
    elseif  derivative_order(D) == 3
        print(io, "SBP third-derivative operator")
    else
        print(io, "SBP ", derivative_order(D),
              "-derivative operator")
    end
    print(io, " of order ", accuracy_order(D))
    if get(io, :compact, false) == false
        print(io, " on a grid in [", first(grid(D)), ", ", last(grid(D)),
                "] using ", length(grid(D)), " nodes \n")
        print(io, "and coefficients")
    end
    print(io, " of ", source_of_coefficients(D))
end

"""
    derivative_left(D::DerivativeOperator, u, der_order::Val{N})

Compute the `N`-th derivative of the function given by the coefficients `u` at
the left boundary of the grid.
"""
@inline function derivative_left(D::DerivativeOperator, u, der_order::Val{N}) where {N}
    convolve_left_row(D.coefficients.left_boundary_derivatives[N], u) / D.Δx^N
end

@inline function derivative_left(D::AbstractNonperiodicDerivativeOperator, u, der_order::Val{0})
    u[begin]
end
@inline function derivative_left(D::DerivativeOperator, u, der_order::Val{0})
    u[begin]
end

"""
    derivative_right(D::DerivativeOperator, u, der_order::Val{N})

Compute the `N`-th derivative of the function given by the coefficients `u` at
the right boundary of the grid.
"""
@inline function derivative_right(D::DerivativeOperator, u, der_order::Val{N}) where {N}
    convolve_right_row(D.coefficients.right_boundary_derivatives[N], u) / D.Δx^N
end

@inline function derivative_right(D::AbstractNonperiodicDerivativeOperator, u, der_order::Val{0})
    u[end]
end
@inline function derivative_right(D::DerivativeOperator, u, der_order::Val{0})
    u[end]
end

"""
    mul!(du, D::DerivativeOperator, u, α=true, β=false)

Efficient in-place version of `du = α * D * u + β * du`. Note that `du` must not
be aliased with `u`.
"""
function mul! end

"""
    mul_transpose_derivative_left!(u, D::DerivativeOperator, der_order::Val{N}, α=true, β=false)

Set the grid function `u` to `α` times the transposed `N`-th derivative functional
applied to `u` plus `β` times `u` in the domain of the `N`-th derivative functional
at the left boundary of the grid.
Thus, the coefficients `α, β` have the same meaning as in [`mul!`](@ref).
"""
@inline function mul_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α=true, β=false) where {N}
    factor = α / D.Δx^N
    coef = D.coefficients.left_boundary_derivatives[N].coef
    for i in eachindex(coef)
        u[i] = factor * coef[i] + β * u[i]
    end
end

@inline function mul_transpose_derivative_left!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α=true, β=false)
    @inbounds u[begin] = α * u[begin] + β * u[begin]
    return nothing
end

"""
    mul_transpose_derivative_right!(u, D::DerivativeOperator, der_order::Val{N}, α=true, β=false)

Set the grid function `u` to `α` times the transposed `N`-th derivative functional
applied to `u` plus `β` times `u` in the domain of the `N`-th derivative functional
at the right boundary of the grid.
Thus, the coefficients `α, β` have the same meaning as in [`mul!`](@ref).
"""
@inline function mul_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{N}, α=true, β=false) where {N}
    factor = α / D.Δx^N
    coef = D.coefficients.right_boundary_derivatives[N].coef
    for i in eachindex(coef)
        u[end-i+1] += factor * coef[i]
    end
end

@inline function mul_transpose_derivative_right!(u::AbstractVector, D::DerivativeOperator, der_order::Val{0}, α=true, β=false)
    @inbounds u[end] = α * u[end] + β * u[end]
    return nothing
end

struct OneHot <: AbstractVector{Bool}
    i::Int
    length::Int
end

Base.getindex(v::OneHot, i::Integer) = i == v.i
Base.size(v::OneHot) = (v.length,)

"""
    derivative_left(D::AbstractNonperiodicDerivativeOperator, der_order)

Get a representation of the linear functional evaluation the `N`th derivative at
the left boundary node as (dense) vector.
"""
function derivative_left(D::AbstractNonperiodicDerivativeOperator, der_order)
    [derivative_left(D, OneHot(i, size(D, 2)), der_order) for i in 1:size(D, 2)]
end

"""
    derivative_right(D::AbstractNonperiodicDerivativeOperator, der_order)

Get a representation of the linear functional evaluation the `N`th derivative at
the right boundary node as (dense) vector.
"""
function derivative_right(D::AbstractNonperiodicDerivativeOperator, der_order)
    [derivative_right(D, OneHot(i, size(D, 2)), der_order) for i in 1:size(D, 2)]
end

# Compute `α*D*u + β*dest` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor, β)
end

# Compute `α*D*u` and store the result in `dest`.
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::DerivativeOperator, u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor)
end

"""
    integrate(func, u, D::DerivativeOperator)

Map the function `func` to the coefficients `u` and integrate with respect to
the quadrature rule associated with the SBP derivative operator `D`.
"""
function integrate(func, u::AbstractVector, D::DerivativeOperator)
    @boundscheck begin
        length(u) == length(grid(D))
    end
    @unpack Δx = D
    @unpack left_weights, right_weights = D.coefficients

    @inbounds res = sum(func, view(u,1+length(left_weights):length(u)-length(right_weights)))
    @inbounds for i in Base.OneTo(length(left_weights))
        res += left_weights[i]*func(u[i])
    end
    @inbounds for i in Base.OneTo(length(right_weights))
        res += right_weights[i]*func(u[end-i+1])
    end

    Δx * res
end

"""
    derivative_operator(source_of_coefficients,
                        derivative_order, accuracy_order,
                        xmin, xmax, N, mode=FastMode())
    derivative_operator(source_of_coefficients;
                        derivative_order, accuracy_order,
                        xmin, xmax, N, mode=FastMode())

Create a [`DerivativeOperator`](@ref) approximating the `derivative_order`-th derivative
on a grid between `xmin` and `xmax` with `N` grid points up to order of accuracy
`accuracy_order`. with coefficients given by `source_of_coefficients`.
The evaluation of the derivative can be parallized using threads by chosing
`mode=ThreadedMode()`.
"""
function derivative_operator(source_of_coefficients, derivative_order, accuracy_order,
                            xmin, xmax, N, mode=FastMode())
    if mode === Val(:serial) || mode === Val(:threads)
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the argument `parallel` is deprecated." *
                     "Use `mode` instead.", :derivative_operator)
        mode = _parallel_to_mode(mode)
    end
    grid = construct_grid(source_of_coefficients, accuracy_order, xmin, xmax, N)
    if derivative_order == 1
        coefficients = first_derivative_coefficients(source_of_coefficients, accuracy_order, eltype(grid), mode)
    elseif derivative_order == 2
        coefficients = second_derivative_coefficients(source_of_coefficients, accuracy_order, eltype(grid), mode)
    elseif derivative_order == 3
        coefficients = third_derivative_coefficients(source_of_coefficients, accuracy_order, eltype(grid), mode)
    elseif derivative_order == 4
        coefficients = fourth_derivative_coefficients(source_of_coefficients, accuracy_order, eltype(grid), mode)
    else
        throw(ArgumentError("Derivative order $derivative_order not implemented."))
    end
    DerivativeOperator(coefficients, grid)
end

function derivative_operator(source_of_coefficients;
                             derivative_order, accuracy_order,
                             xmin, xmax, N, mode=FastMode(), parallel=nothing)
    if parallel !== nothing
        # TODO: deprecated in v0.5
        Base.depwarn("Providing the keyword argument `parallel` is deprecated." *
                     "Use `mode` instead.", :dissipation_operator)
        mode = _parallel_to_mode(parallel)
    end
    derivative_operator(source_of_coefficients, derivative_order, accuracy_order,
                        xmin, xmax, N, mode)
end

@inline construct_grid(source_of_coefficients, accuracy_order, xmin, xmax, N) = range(xmin, stop=xmax, length=N)

function lower_bandwidth(D::DerivativeOperator)
    @unpack left_boundary, right_boundary, lower_coef = D.coefficients

    l = length(lower_coef)
    for (i,coef_row) in enumerate(left_boundary)
        l = max(l, i-offset(coef_row))
    end
    for (i,coef_row) in enumerate(right_boundary)
        l = max(l, 1+length(coef_row)-offset(coef_row)-i)
    end
    l
end

function upper_bandwidth(D::DerivativeOperator)
    @unpack left_boundary, right_boundary, upper_coef = D.coefficients

    u = length(upper_coef)
    for (i,coef_row) in enumerate(left_boundary)
        u = max(u, 1+length(coef_row)-offset(coef_row)-i)
    end
    for (i,coef_row) in enumerate(right_boundary)
        u = max(u, i-offset(coef_row))
    end
    u
end

function left_boundary_weight(D::DerivativeOperator)
    @inbounds retval = D.Δx * D.coefficients.left_weights[1]
    retval
end

function right_boundary_weight(D::DerivativeOperator)
    @inbounds retval = D.Δx * D.coefficients.right_weights[1]
    retval
end


text:
# """Extract the dimension of the space that this physical property lives in"""
# dim(p::AbstractSymmetry{Dim}) where {Dim} = Dim

"""
    Specie

Represents a set of particles which are all the same. The type of particle is given by `Specie.particle` and the volume fraction this specie occupies is given by `Specie.volume_fraction`.

We can use `Specie.numberofparticles` to specify the number of particles, otherwise for an infinite `Specie.numberofparticles = Inf`.

The minimum distance between any two particles will equal `outer_radius(Specie) * Specie.exclusion_distance`.
"""
struct Specie{T<:AbstractFloat,Dim,P<:AbstractParticle{Dim}}
    particle::P
    volume_fraction::T
    exclusion_distance::T
end

# Convenience constructor which does not require explicit types/parameters
function Specie(p::AbstractParticle{Dim}; number_density::T = 0.0, volume_fraction::T = number_density * volume(p), exclusion_distance::T = 1.005) where {Dim,T<:AbstractFloat}
    Specie{T,Dim,typeof(p)}(p,volume_fraction,exclusion_distance)
end

function Specie(medium::P,s::S; kws...) where {Dim,T,P<:PhysicalMedium{Dim},S<:Shape{Dim}}
    Specie(Particle(medium, s); kws...)
end

function Specie(medium::P, radius::T; kws...) where {T,P<:PhysicalMedium}
    Specie(Particle(medium, radius); kws...)
end

# Shorthand for all Vectors of species
Species{T<:AbstractFloat,Dim,P} = Vector{S} where S<:Specie{T,Dim,P}

"Returns the volume fraction of the specie."
volume_fraction(s::Specie) = s.volume_fraction
volume_fraction(ss::Species) = sum(volume_fraction.(ss))

import MultipleScattering.volume
volume(s::Specie) = volume(s.particle)

import MultipleScattering.outer_radius

"""
    number_density(s::Specie)

Gives the number of particles per unit volume. Note this is given exactly by `N / V` where `V` is the volume of the region containing the origins of all particles. For consistency, [`volume_fraction`](@ref) is given by `N * volume(s) / V`.
"""
number_density(s::Specie) = s.volume_fraction / volume(s)
# number_density(s::Specie{T,2}) where {T} = s.volume_fraction / (outer_radius(s.particle)^2 * pi)
# number_density(s::Specie{T,3}) where {T} = s.volume_fraction / (T(4/3) * outer_radius(s.particle)^3 * pi)
number_density(ss::Species) = sum(number_density.(ss))

outer_radius(s::Specie) = outer_radius(s.particle)

import MultipleScattering.get_t_matrices
import MultipleScattering.t_matrix

get_t_matrices(medium::PhysicalMedium, species::Vector{S}, ω::AbstractFloat, Nh::Integer) where S<:Specie = get_t_matrices(medium, [s.particle for s in species], ω, Nh)

t_matrix(s::Specie, medium::PhysicalMedium, ω::AbstractFloat, order::Integer) = t_matrix(s.particle, medium, ω, order)

"""
    Material(region::Shape, species::Species [, numberofparticles = Inf])

Creates a material filled with [`Specie`](@ref)'s inside `region`.
"""
struct Material{Dim,S<:Shape,Sps<:Species}
    shape::S
    species::Sps
    numberofparticles::Number
    # Enforce that the Dims and Types are all the same
    function Material{Dim,S,Sps}(shape::S,species::Sps,numberofparticles::Number = Inf) where {T,Dim,S<:Shape{Dim},Sps<:Species{T,Dim}}
        new{Dim,S,Sps}(shape,species,numberofparticles)
    end
end

# Convenience constructor which does not require explicit types/parameters
function Material(shape::S,species::Sps) where {T,Dim,S<:Shape{Dim},Sps<:Species{T,Dim}}

    V = volume(shape)
    numberofparticles = round(sum(
        volume_fraction(s) * V / volume(s)
    for s in species))

    Material{Dim,S,Sps}(shape,species,numberofparticles)
end

function Material(shape::S,specie::Sp) where {T,Dim,S<:Shape{Dim},Sp<:Specie{T,Dim}}
    Material(shape,[specie])
end

import MultipleScattering.PhysicalMedium
import MultipleScattering.Symmetry

PhysicalMedium(s::Specie) = typeof(s.particle.medium)
PhysicalMedium(m::Material) = PhysicalMedium(m.species[1])

Symmetry(m::Material) = Symmetry(m.shape)

# """
#     setupsymmetry(source::AbstractSource, material::Material)
#
# Returns the shared symmetries between the `source` and `materail`.
# """
# setupsymmetry(source::AbstractSource, material::Material{Dim}) where Dim = WithoutSymmetry{Dim}()

# function setupsymmetry(source::PlaneSource{T,3,1}, material::Material{3,Sphere{T,3}};
#         basis_order::Int = 4, ω::T = 0.9) where T
#
#     ls, ms = spherical_harmonics_indices(basis_order)
#
#     gs = regular_spherical_coefficients(source)(basis_order,origin(material.shape),ω)
#     azimuthal_symmetry = norm((ms[n] != 0) ? gs[n] : zero(T) for n in eachindex(gs)) ≈ zero(T)
#
#     return if azimuthal_symmetry
#         AzimuthalSymmetry{3}()
#     else
#         WithoutSymmetry{3}()
#     end
# end
#
# function setupsymmetry(psource::PlaneSource{T,Dim}, material::Material{Dim,S}) where {T<:AbstractFloat, Dim, S<:Union{Halfspace{T,Dim},Plate{T,Dim}}}
#
#     hv = material.shape.normal
#     kv = psource.direction
#
#     if abs(dot(hv, kv)^2) ≈ abs(dot(hv, hv) * dot(kv, kv))
#         # for direct incidence
#         return PlanarAzimuthalSymmetry{Dim}()
#     else
#         return PlanarSymmetry{Dim}()
#     end
# end

ext: jl

max_stars_repo_name: arturgower/EffectiveWaves


text:
"""
Recursively descends an expression's AST to find all of the symbols contained
in it. Inserts any unquoted symbols that are found into the set argument,
`s`. Note that this function is not designed to handle assignment-like
expressions: it is intended for application to value expressions only.

Arguments:

* s::Set{Symbol}: A set of symbols that will be mutated whenever any new
    symbols are found.
* e::Any: An expression-like object that will be descended through to find new
    symbols.

Returns:

* Void: This function is used exclusively to mutate the argument `s`.
"""
function find_symbols!(
    ds::Dict{Symbol, Set{Symbol}}, srcs::Set{Symbol}, e::Any,
    index::Dict{Symbol, Symbol},
)::Void
    if isa(e, Expr)
        if e.head == :call
            # ignore e.args[1], which cannot be a token-field ref
            for i in 2:length(e.args)
                find_symbols!(ds, srcs, e.args[i], index)
            end
        elseif e.head in (:(||), :(&&), :if)
            for i in 1:length(e.args)
                find_symbols!(ds, srcs, e.args[i], index)
            end
        elseif e.head == :.
            _token, _field = e.args[1], e.args[2]
            field = _field.args[1] # field comes wrapped in an Expr with head :quote
            if haskey(index, _token)
                s = get!(ds, _token, Set{Symbol}()) # Create a set s to hold the fields
                push!(s, field)
                push!(srcs, index[_token]) # remember sources we need
            end
        else
            # NOTE: Do not descend when e is a QuoteNode.
        end
    end
    return
end

"""
Recursively descends an expression's AST to find all of the symbols contained
in it. Returns all found symbols in a `Set{Symbol}` object.

Arguments:

* e::Any: An expression-like object that will be descended through to find new
    symbols.

Returns:

* s::Set{Symbol}: A set containing all of the symbols found by descending
    through the expression-like object's AST.
"""
function find_symbols(e, index)::Dict{Symbol, Set{Symbol}}
    ds = Dict{Symbol, Set{Symbol}}()
    find_symbols!(ds, e)
    return ds
end

"""
Produce a mapping from symbols to numeric indices and a reverse mapping from
numeric indices to symbols.

Arguments:

* s::Set{Symbol}: A set of symbols that should be assigned numeric indices.

Returns:

* mapping::Dict{Symbol, Int}: A mapping from symbols to indices.
* reverse_mapping::Vector{Symbol}: A mapping from indices to symbols.
"""
function map_symbols(
        ds::Dict{Symbol, Set{Symbol}}
)::Tuple{Dict{Symbol, Dict{Symbol, Int}}, Dict{Symbol, Vector{Symbol}}}

    smaps = Dict{Symbol, Dict{Symbol, Int}}()
    reverse_smaps = Dict{Symbol, Vector{Symbol}}()

    for token in collect(keys(ds))
        mapping = Dict{Symbol, Int}()
        s = ds[token]
        reverse_mapping = Array(Symbol, length(s))
        for (i, sym) in enumerate(s)
            mapping[sym] = i
            reverse_mapping[i] = sym
        end
        smaps[token] = mapping
        reverse_smaps[token] = reverse_mapping
    end
    return smaps, reverse_smaps
end

"""
Traverse an AST-like object and replace a fixed set of symbols with
tuple-indexing expressions.

Arguments:

* e::Any: An AST-like object.
* mapping:Dict{Symbol, Int}: A mapping from symbols to numeric indices.
* tuple_name::Symbol: The name of the tuple that will be indexed into.

Returns:

* new_e::Any: A new AST-like object with all symbols replaced with
    tuple-indexing operations.
"""
function replace_symbols!(ai::ArgsIndex, e::Any, index, primary)::Any
    # smaps::Dict{Symbol, Dict{Symbol, Int}}
# )::Any
    if isa(e, Expr)
        # To ensure purity, we copy any Expr objects rather than mutate them.
        new_e = copy(e)
        if e.head == :call
            # Do two things:
            #   1)  Escape function `f` being called so its not rooted to the
            #       module StructuredQueries
            #   2)  replace `f(xs...)` with `lift(f, xs...)`
            args_copy = copy(e.args)
            lifted_e = Expr(:call)
            push!(lifted_e.args, Expr(:., :StructuredQueries, QuoteNode(:lift)))
            push!(lifted_e.args, esc(args_copy[1]))
            for i in 2:length(args_copy)
                push!(lifted_e.args,
                      replace_symbols!(ai, args_copy[i], index, primary))
            end
            return lifted_e
        elseif new_e.head in (:(||), :(&&), :if)
            for i in 1:length(new_e.args)
                new_e.args[i] = replace_symbols!(ai, new_e.args[i], index, primary)
            end
        elseif e.head == :. # token.attr
            _token = new_e.args[1]
            _arg = new_e.args[2]
            if haskey(index, _token)
                arg = _arg.args[1] #grab the attribute name from inside the :quote Expr
                s = get!(ai, _token)
                push!(s, arg)
            end
            return new_e
            # token, _field = e.args[1], e.args[2]
            # field = _field.args[1] # field comes wrapped in Expr with head :quote
            # if haskey(smaps, token)
            #     return Expr(:ref, token, smaps[token][field])
            # else
            #     return e
            # end
        else
            # TODO: Handle other kinds of Expr heads.
            error(
                @sprintf("Unknown Expr head type %s in %s", new_e.head, new_e)
            )
        end
        return new_e
    elseif isa(e, Symbol)
        if haskey(index, e)
            s = get!(ai, e)
        end
        if !isnull(primary)
            token = index[get(primary)]
            s = get!(ai, token)
            push!(s, e)
            return Expr(:call, :getfield, token, QuoteNode(e))
        end
    # elseif isa(e, Symbol)
    #     # Replace unquoted symbols with tuple indexing expressions.
    #     return Expr(:ref, tuple_name, mapping[e])
    # elseif isa(e, QuoteNode)
    #     # Replace quoted symbols with raw symbols.
    #     return e.value
    else
        # Hopefully we have a literal here since we stop going down the AST
        # when we hit this branch.
        return e
    end
    return
end


text:
"""
PDE_Tools

This file has the operators I need for the PDE example. They
live in a separate file to make the CI easier for me to organize.
"""
# Famous sparse matrices
"""
Dx2d(n)

returns x partial on n x n grid.
Unit square, homogeneous Dirichlet BC
"""
function Dx2d(n)
    h = 1 / (n + 1)
    ssdiag = ones(n^2 - 1) / (2 * h)
    for iz = n:n:n^2-1
        ssdiag[iz] = 0.0
    end
    updiag = Pair(1, ssdiag)
    lowdiag = Pair(-1, -ssdiag)
    Dx = spdiagm(lowdiag, updiag)
    return Dx
end

"""
Dy2d(n)

returns y partial on n x n grid.
Unit square, homogeneous Dirichlet BC
"""
function Dy2d(n)
    h = 1 / (n + 1)
    ssdiag = ones(n^2 - n) / (2 * h)
    updiag = Pair(n, ssdiag)
    lowdiag = Pair(-n, -ssdiag)
    Dy = spdiagm(lowdiag, updiag)
    return Dy
end

"""
Lap2d(n)

returns the negative Laplacian in two space dimensions
on n x n grid.

Unit square, homogeneous Dirichlet BC
"""
function Lap2d(n)
    # hm2=1/h^2
    hm2 = (n + 1.0)^2
    maindiag = fill(4 * hm2, (n^2,))
    sxdiag = fill(-hm2, (n^2 - 1,))
    sydiag = fill(-hm2, (n^2 - n,))
    for iz = n:n:n^2-1
        sxdiag[iz] = 0.0
    end
    D2 = spdiagm(-n => sydiag, -1 => sxdiag, 0 => maindiag, 1 => sxdiag, n => sydiag)
    return D2
end

"""
u=fish2d(f, fdata)

Fast Poisson solver in two space dimensions.
Same as the Matlab code.
Unit squre + homogeneous Dirichlet BCs.

Grid is nx by nx

You give me f as a two-dimensional vector f(x,y).
I return the solution u.
"""
function fish2d(f, fdata)
    u = fdata.utmp
    v = fdata.uhat
    T = fdata.T
    ST = fdata.ST
    (nx, ny) = size(f)
    nx == ny || error("need a square grid in fish2d")
    u .= f
    u = ST * u
    u = u'
    u1 = reshape(u, (nx * nx,))
    v1 = reshape(v, (nx * nx,))
    v1 .= u1
    ldiv!(u1, T, v1)
    u = u'
    u .= ST * u
    u ./= (2 * nx + 2)
    return u
end

"""
fishinit(n)

Run FFTW.plan_r2r to set up the solver. Do not mess
with this function.
"""
function fishinit(n)
    #
    # Get the sine transform from FFTW. This is faster/better/cleaner
    # than what I did in the Matlab codes.
    #
    zstore = zeros(n, n)
    ST = FFTW.plan_r2r!(zstore, FFTW.RODFT00, 1)
    uhat = zeros(n, n)
    fishu = zeros(n, n)
    TD = newT(n)
    T = lu!(TD)
    fdata = (ST = ST, uhat = uhat, utmp = zstore, T = T, fishu = fishu)
    return fdata
end

"""
T = newT(n)

Builds the n^2 x n^2 sparse tridiagonal matrix for
the 2D fast Poisson solver.
"""
function newT(n)
    N = n * n
    h = 1 / (n + 1)
    x = h:h:1-h
    h2 = 1 / (h * h)
    LE = 2 * (2 .- cos.(pi * x)) * h2
    fn = ones(N - 1) * h2
    gn = ones(N - 1) * h2
    dx = zeros(N)
    for k = 1:n-1
        fn[k*n] = 0.0
        gn[k*n] = 0.0
        dx[(k-1)*n+1:n*k] = LE[k] * ones(n)
    end
    dx[(n-1)*n+1:n*n] = LE[n] * ones(n)
    T = Tridiagonal(-fn, dx, -gn)
    return T
end

"""
Use fish2d and reshape for preconditioning.
"""
function Pfish2d(v, fdata)
    n2 = length(v)
    n = Int(sqrt(n2))
    (n * n == n2) || error("input to Pfish2d not a square array")
    v2 = reshape(v, (n, n))
    u = fish2d(v2, fdata)
    u = reshape(u, (n2,))
    return u
end

"""
Pvec2d(v, u, pdata)

Returns inverse Laplacian * v

u is a dummy argument to make nsoli happy

Preconditioner for nsoli
"""
function Pvec2d(v, u, pdata)
    fdata = pdata.fdata
    p = Pfish2d(v, fdata)
    return p
end


text:
"""
    track_propagated_discontinuities!(integrator::DDEIntegrator)

Try to find a propagated discontinuity in the time interval `[integrator.t, integrator.t +
integrator.dt]` and add it to the set of discontinuities and grid points of the
`integrator`.
"""
function track_propagated_discontinuities!(integrator::DDEIntegrator)
  # calculate interpolation points
  interp_points = integrator.discontinuity_interp_points
  Θs = range(zero(integrator.t); stop = oneunit(integrator.t), length = interp_points)

  # for dependent lags and previous discontinuities
  for lag in integrator.sol.prob.dependent_lags, discontinuity in integrator.tracked_discontinuities
    # obtain time of previous discontinuity
    T = discontinuity.t

    # estimate subinterval of current integration step that contains a propagated
    # discontinuity induced by the lag and the previous discontinuity
    interval = discontinuity_interval(integrator, lag, T, Θs)

    # if a discontinuity exists in the current integration step
    if interval !== nothing
      # estimate time point of discontinuity
      t = discontinuity_time(integrator, lag, T, interval)

      # add new discontinuity of correct order at the estimated time point
      if integrator.sol.prob.neutral
        d = Discontinuity(t, discontinuity.order)
      else
        d = Discontinuity(t, discontinuity.order + 1)
      end
      push!(integrator.opts.d_discontinuities, d)
      push!(integrator.opts.tstops, t)

      # analogously to RADAR5 we do not strive for finding the first discontinuity
      break
    end
  end

  nothing
end

"""
    discontinuity_function(integrator::DDEIntegrator, lag, T, t)

Evaluate function ``f(x) = T + lag(u(x), p, x) - x`` at time point `t`, where `T` is time
point of a previous discontinuity and `lag` is a dependent delay.
"""
function discontinuity_function(integrator::DDEIntegrator, lag, T, t)
  tmp = get_tmp_cache(integrator)
  cache = tmp === nothing ? nothing : first(tmp)

  # estimate state at the given time
  if cache === nothing
    ut = integrator(t, Val{0})
  else
    integrator(cache, t, Val{0})
    ut = cache
  end

  T + lag(ut, integrator.p, t) - t
end

"""
    discontinuity_interval(integrator::DDEIntegrator, lag, T, Θs)

Return an estimated subinterval of the current integration step of the `integrator` that
contains a propagated discontinuity induced by the dependent delay `lag` and the
discontinuity at time point `T`, or `nothing`.

The interval is estimated by checking the signs of `T + lag(u(t), p, t) - t` for time points
`integrator.t .+ θs` in the interval `[integrator.t, integrator.t + integrator.dt]`.
"""
function discontinuity_interval(integrator::DDEIntegrator, lag, T, Θs)
  # use start and end point of last time interval to check for discontinuities
  previous_condition = discontinuity_function(integrator, lag, T, integrator.t)
  if isapprox(previous_condition, 0;
              rtol = integrator.discontinuity_reltol,
              atol = integrator.discontinuity_abstol)
    prev_sign = 0
  else
    prev_sign = cmp(previous_condition, zero(previous_condition))
  end
  new_condition = discontinuity_function(integrator, lag, T, integrator.t + integrator.dt)
  new_sign = cmp(new_condition, zero(new_condition))

  # if signs are different we already know that a discontinuity exists
  if prev_sign * new_sign < 0
    return (zero(eltype(Θs)), one(eltype(Θs)))
  end

  # recheck interpolation intervals if no discontinuity found yet
  prev_Θ = zero(eltype(Θs))

  for i in 2:length(Θs)
    # evaluate sign at next interpolation point
    new_Θ = Θs[i]
    new_t = integrator.t + new_Θ * integrator.dt
    new_condition = discontinuity_function(integrator, lag, T, new_t)
    new_sign = cmp(new_condition, zero(new_condition))

    # return estimated interval if we find a root or observe a switch of signs
    if new_sign == 0
      return (new_Θ, new_Θ)
    elseif prev_sign * new_sign < 0
      return (prev_Θ, new_Θ)
    else
      # otherwise update lower estimate of subinterval
      prev_sign = new_sign
      prev_Θ = new_Θ
    end
  end

  nothing
end

"""
    discontinuity_time(integrator::DDEIntegrator, lag, T, interval)

Estimate time point of the propagated discontinuity induced by the dependent delay
`lag` and the discontinuity at time point `T` inside the `interval` of the current
integration step of the `integrator`.
"""
function discontinuity_time(integrator::DDEIntegrator, lag, T, (bottom_Θ, top_Θ))
  if bottom_Θ == top_Θ
    # in that case we have already found the time point of a discontinuity
    Θ = top_Θ
  else
    # define function for root finding
    zero_func = let integrator = integrator, lag = lag, T = T, t = integrator.t, dt = integrator.dt
      (θ, p=nothing) -> discontinuity_function(integrator, lag, T, t + θ * dt)
    end

    Θ = NonlinearSolve.solve(
      NonlinearSolve.NonlinearProblem{false}(zero_func, (bottom_Θ, top_Θ)),
      NonlinearSolve.Falsi(),
    ).left
  end

  integrator.t + Θ * integrator.dt
end

max_stars_repo_name: JuliaDiffEq/DelayDiffEq

max_issues_repo_name: JuliaDiffEq/DelayDiffEq

max_forks_repo_name: JuliaDiffEq/DelayDiffEq


text:
@time @testset "PSSE Parsing " begin
    files = readdir(PSSE_RAW_DIR)
    if length(files) == 0
        error("No test files in the folder")
    end

    for f in files[1:1]
        @info "Parsing $f ..."
        pm_dict = PowerSystems.parse_file(joinpath(PSSE_RAW_DIR, f))
        @info "Successfully parsed $f to PowerModels dict"
        PowerSystems.make_mixed_units(pm_dict)
        @info "Successfully converted $f to mixed_units"
        ps_dict = PowerSystems.pm2ps_dict(pm_dict)
        @info "Successfully parsed $f to PowerSystems dict"
        Buses, Generators, Storage, Branches, Loads, LoadZones, Shunts =
            PowerSystems.ps_dict2ps_struct(ps_dict)
        @info "Successfully parsed $f to PowerSystems devices"
        sys_test = System(Buses, Generators, Loads, Branches, Storage,
            float(ps_dict["baseMVA"])) # TODO: Add DClines, Shunts
        @info "Successfully parsed $f to System struct"
    end

    # Test bad input
    pm_dict = PowerSystems.parse_file(joinpath(PSSE_RAW_DIR, files[1]))
    pm_dict["bus"] = Dict{String, Any}()
    @test_throws PowerSystems.DataFormatError PowerSystems.pm2ps_dict(pm_dict)
end


text:
export TableOfContents

"""Generate Table of Contents using Markdown cells. Headers h1-h6 are used.

# Keyword arguments:
`title` header to this element, defaults to "Table of Contents"

`indent` flag indicating whether to vertically align elements by hierarchy

`depth` value to limit the header elements, should be in range 1 to 6 (default = 3)

`aside` fix the element to right of page, defaults to true

# Examples:

```julia
TableOfContents()

TableOfContents(title="Experiments 🔬")

TableOfContents(title="📚 Table of Contents", indent=true, depth=4, aside=true)
```
"""
Base.@kwdef struct TableOfContents
    title::AbstractString="Table of Contents"
    indent::Bool=true
    depth::Integer=3
    aside::Bool=true
end

function Base.show(io::IO, ::MIME"text/html", toc::TableOfContents)

    withtag(io, :script) do
        print(io, """
            const getParentCellId = el => {
                // Traverse up the DOM tree until you reach a pluto-cell
                while (el.nodeName != 'PLUTO-CELL') {
                    el = el.parentNode;
                    if (!el) return null;
                }
                return el.id;
            }

            const getElementsByNodename = nodeName => Array.from(
                document.querySelectorAll(
                    "pluto-notebook pluto-output " + nodeName
                )
            ).map(el => {
                return {
                    "nodeName" : el.nodeName,
                    "parentCellId": getParentCellId(el),
                    "innerText": el.innerText
                }
            })

            const getPlutoCellIds = () => Array.from(
                document.querySelectorAll(
                    "pluto-notebook pluto-cell"
                )
            ).map(el => el.id)

            const getHeaders = () => {
                const depth = Math.max(1, Math.min(6, $(toc.depth))) // should be in range 1:6
                const range = Array.from({length: depth}, (x, i) => i+1) // [1, ... depth]
                let headers = [].concat.apply([], range.map(i => getElementsByNodename("h"+i))); // flatten [[h1s...], [h2s...], ...]
                const plutoCellIds = getPlutoCellIds()
                headers.sort((a,b) => plutoCellIds.indexOf(a.parentCellId) - plutoCellIds.indexOf(b.parentCellId)); // sort in the order of appearance
                return headers
            }

            const tocIndentClass = '$(toc.indent ? "-indent" : "")'
            const render = (el) => `\${el.map(h => `

                                            {document.getElementById('\${h.parentCellId}').firstElementChild.classList.add('highlight-pluto-cell-shoulder')})()"
                                                onmouseout="(()=>{document.getElementById('\${h.parentCellId}').firstElementChild.classList.remove('highlight-pluto-cell-shoulder')})()"
                                                onclick="((e)=>{
                                                    e.preventDefault();
                                                    document.getElementById('\${h.parentCellId}').scrollIntoView({
                                                        behavior: 'smooth',
                                                        block: 'center'
                                                    });
                                                })(event)"
                                                > \${h.innerText}

`).join('')}`

            const tocClass = '$(toc.aside ? "toc-aside" : "")'

            const tocContentNode = html`<div class="toc-content" id="toc-content"></div>`
            const tocNode = html`

$(toc.title)

                    \${tocContentNode}

`
            const updateCallback = e => {
                tocContentNode.innerHTML = render(getHeaders())
            }
            updateCallback()

            const observers = {
                current: [],
            }

            const notebook = document.querySelector("pluto-notebook")
            const createCellObservers = () => {
                observers.current.forEach((o) => o.disconnect())
                observers.current = Array.from(notebook.querySelectorAll("pluto-cell")).map(el => {
                    const o = new MutationObserver(updateCallback)
                    o.observe(el, {attributeFilter: ["class"]})
                    return o
                })
            }
            const notebookObserver = new MutationObserver(createCellObservers)
            notebookObserver.observe(notebook, {childList: true})

            createCellObservers()

            invalidation.then(() => {
                notebookObserver.disconnect()
                observers.current.forEach((o) => o.disconnect())
            })

            return tocNode
        """)
    end

    withtag(io, :style) do
        print(io, """
            @media screen and (min-width: 1081px) {
                .toc-aside {
                    position:fixed;
                    right: 1rem;
                    top: 5rem;
                    width:25%;
                    padding: 10px;
                    border: 3px solid rgba(0, 0, 0, 0.15);
                    border-radius: 10px;
                    box-shadow: 0 0 11px 0px #00000010;
                    max-height: 500px;
                    overflow: auto;
                }
            }

            .toc-title{
                display: block;
                font-size: 1.5em;
                margin-top: 0.67em;
                margin-bottom: 0.67em;
                margin-left: 0;
                margin-right: 0;
                font-weight: bold;
                border-bottom: 2px solid rgba(0, 0, 0, 0.15);
            }

            .toc-row {
                white-space: nowrap;
                overflow: hidden;
                text-overflow: ellipsis;
                padding-bottom: 2px;
            }

            .highlight-pluto-cell-shoulder {
                background: rgba(0, 0, 0, 0.05);
                background-clip: padding-box;
            }

            a {
                text-decoration: none;
                font-weight: normal;
                color: gray;
            }
            a:hover {
                color: black;
            }
            a.H1-indent {
                padding: 0px 0px;
            }
            a.H2-indent {
                padding: 0px 10px;
            }
            a.H3-indent {
                padding: 0px 20px;
            }
            a.H4-indent {
                padding: 0px 30px;
            }
            a.H5-indent {
                padding: 0px 40px;
            }
            a.H6-indent {
                padding: 0px 50px;
            }
            """)
    end
end


text:
using Test

function run(nums, target_turn=2020)
    last_said_on_turn = Dict()

    turn = 1
    # initialize
    for n in nums[1:length(nums) - 1]
        last_said_on_turn[n] = turn
        turn += 1
        # println("Turn = ", turn, " Val = ", n)
    end

    prev = last(nums)
    while turn < target_turn
        # compute next number to say
        # If that was the first time the number has been spoken, the current player
        # says 0. (i.e. default to 0)
        next = 0
        if haskey(last_said_on_turn, prev)
            # Otherwise, the number had been spoken before; the current player
            # announces how many turns apart the number is from when it was previously
            # spoken.
            next = turn - last_said_on_turn[prev]
        end

        # speak!
        last_said_on_turn[prev] = turn
        # println("Turn = ", turn, " Val = ", prev)

        prev = next
        turn += 1
    end

    println("Turn = ", turn, " Val = ", prev)
    return prev
end

println("\n[Part 1] Examples:")
@test run([0,3,6], 4) == 0
@test run([0,3,6], 5) == 3
@test run([0,3,6], 6) == 3
@test run([0,3,6], 7) == 1
@test run([0,3,6], 8) == 0
@test run([0,3,6], 9) == 4
@test run([0,3,6], 10) == 0
@test run([0,3,6]) == 436

println("\n[Part 1] Result:")
run([19,0,5,1,10,13])

println("\n[Part 2] Examples:")
@test run([0,3,6], 30000000) == 175594.
# Given 1,3,2, the 30000000th number spoken is 2578.
# Given 2,1,3, the 30000000th number spoken is 3544142.
# Given 1,2,3, the 30000000th number spoken is 261214.
# Given 2,3,1, the 30000000th number spoken is 6895259.
# Given 3,2,1, the 30000000th number spoken is 18.
# Given 3,1,2, the 30000000th number spoken is 362.

println("\n[Part 2] Result:")
run([19,0,5,1,10,13], 30000000)


text:
module Parsso

include("CoreDataStructure/CoreDataStructure.jl")
include("BasicIO/BasicIO.jl")

end # module


text:
using Test
using Random

using IndexFunArrays
using BenchmarkTools

BenchmarkTools.DEFAULT_PARAMETERS.samples = 2
BenchmarkTools.DEFAULT_PARAMETERS.evals = 2

include("other_ifas.jl")
include("test_center_scaling.jl")
include("constructors.jl")
include("concrete_generators.jl")
include("utils.jl")

@testset "Test defined array operation" begin
    f(ind) = ind[1]
    arr = IndexFunArray(Int, f, (10, 10))
    @test arr[2] == arr[12] == arr[2, 2] == 2

    arr2 = similar(arr, (11, 11))
    @test size(arr2) == (11,11)
    @test eltype(arr2) == eltype(arr)

    for s in [(10,), (10,1,2), (100, 100), (20,2,2,2,2)]
        @test s == size(IndexFunArray(ComplexF32, x -> zero(ComplexF32), s))
    end

end

@testset "Check errors" begin
    @test_throws ArgumentError IndexFunArray(Int, x -> zero(Float32), (10, 10))

    a = IndexFunArray(Int, x -> zero(Int), (10, 10))
    @test_throws ErrorException (a[1] = 1)
end

include("performance_tests.jl")

return


text:
module MeshManager

using CoreExtended
using FileManager

using ..StorageManager

type Mesh
	typ::Symbol
	vertices::StorageGroup
	normals::StorageGroup
	uvs::StorageGroup
	indicies::StorageGroup

	function Mesh()
		this=new(:_)
		this.vertices = StorageManager.StorageGroup(:DATA, :VERTEX)
		this.uvs = StorageManager.StorageGroup(:DATA, :UV)
		this.normals = StorageManager.StorageGroup(:DATA, :NORMAL)
		this.indicies = StorageManager.StorageGroup(:DATA, :INDEX)
		this
	end
end

# preset for manager: includes code here
presetManager(MeshManager,Mesh)

isEmpty(mesh::Mesh) = isEmptyVertices(mesh) && isEmptyUVs(mesh) && isEmptyNormals(mesh) && isEmptyIndicies(mesh)
isEmptyVertices(mesh::Mesh) = length(mesh.vertices.data.values) == 0
isEmptyUVs(mesh::Mesh) = length(mesh.uvs.data.values) == 0
isEmptyNormals(mesh::Mesh) = length(mesh.normals.data.values) == 0
isEmptyIndicies(mesh::Mesh) = length(mesh.indicies.data.values) == 0

function update(mesh::Mesh, typ::Symbol, values::AbstractArray, elems=0)
	mesh.typ = typ
	if length(values)>0 StorageManager.setValues(mesh.vertices, values, elems, :TRIANGLES) end
	mesh
end

function update(mesh::Mesh, typ::Symbol, f::Function)
	data = f()
	mesh.typ = typ

	drawtyp	= data[1]
	vertices = data[2]
	normals	= data[3]
	uvs = data[4]
	indicies = data[5]

	if length(vertices)>0 StorageManager.setValues(mesh.vertices, vertices, 3, drawtyp) end
	if length(uvs)>0 StorageManager.setValues(mesh.uvs, uvs, 2, drawtyp) end
	if length(normals)>0 StorageManager.setValues(mesh.normals, normals, 3, drawtyp) end
	if length(indicies)>0 StorageManager.setValues(mesh.indicies, indicies, 1, drawtyp) end

	mesh
end

end #MeshManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
# Configuration for nearly neutral simuation
const nn_simtype = 1    # nearly neutral infinite alleles
const type_str = "N3"
const mu_list_flag = false
@everywhere const N_list     = [50,100]     # popsize
const N_mu_list= [1.0]   # Mutation rate as a multiple of 1.0/N
const L = 1           # number of loci
const ngens = 10000      # Generations after burn-in
const fix_minimum = 0.45  # Fraction of population needed to "fix"
const burn_in= 2.0    # generations of burn_in as a multiple of 1/mu
dfe = dfe_neutral
dfe_str = "dfe_neutral"
const use_poplist=true   # if true, save a list of populations as a check of correctness

max_stars_repo_name: um-tech-evolution/nearly-neutral

max_issues_repo_name: um-tech-evolution/nearly-neutral

max_forks_repo_name: um-tech-evolution/nearly-neutral


text:
using PyPlot
using Random
using Distributions

# Some info on graphs
function get_graph_info()
    transition_matrix_g1 =  [3 7; 5 6; 4 2; 7 1; 4 6; 1 3; 5 2]
    state_txt_g1 = ["G", "i1","s1","bG1","i2","bG2","s2"]
    stepsToGoal_g1 = [0 0; 3 2; 2 3; 4 1; 2 2; 1 3; 3 3]
    correctAction_g1 = [0, 2, 1, 2, 0, 1, 0]
    states_of_interest_g1 = [2, 3, 4, 6] # States where one action is correct and the other wrong
    # minimal_dist_g1 = [0 3 2 4 2 1 3; 0 2 3 1 2 3 3] # same as stepsToGoal_g1
    minimal_dist_g1= [0 2 2 1 2 1 3; 0 3 3 4 2 3 3]
    Vvalues_g1 = [1. 0.6769 0.6559 0.7807 0.7091 0.7952 0.6237]
    state_order_g1 = [4, 6, 2, 3, 5, 7]

    transition_matrix_g3 = [ 7 6; 1 3; 1 4; 6 2; 7 2; 5 3; 4 5]
    state_txt_g3 = ["G", "bG1","bG2","i1","i2","s1","s2"]
    stepsToGoal_g3 = [0 0; 1 2; 1 3; 3 2; 4 2; 3 2; 3 3]
    correctAction_g3 = [0, 1, 1, 2, 2, 2, 0]
    states_of_interest_g3 = [2, 3, 4, 5, 6] # States where one action is correct and the other wrong
    # minimal_dist_g3 = [0 1 1 3 4 3 3; 0 2 3 2 2 2 3] # same as stepsToGoal_g3
    minimal_dist_g3 = [0 1 1 2 2 2 3; 0 2 3 3 4 3 3]
    Vvalues_g3 = [1. 0.8641 0.8090 0.6867 0.6619 0.6619 0.6069]
    state_order_g3 = [3, 2, 5, 4, 6, 7]

    graph_info = Dict("transition_matrix_g1" => transition_matrix_g1,
                    "state_txt_g1" => state_txt_g1,
                    "stepsToGoal_g1" => stepsToGoal_g1,
                    "correctAction_g1" => correctAction_g1,
                    "states_of_interest_g1" => states_of_interest_g1,
                    "minimal_dist_g1" => minimal_dist_g1,
                    "Vvalues_g1" => Vvalues_g1,
                    "state_order_g1" => state_order_g1,
                    #
                    "transition_matrix_g3" => transition_matrix_g3,
                    "state_txt_g3" => state_txt_g3,
                    "stepsToGoal_g3" => stepsToGoal_g3,
                    "correctAction_g3" => correctAction_g3,
                    "states_of_interest_g3" => states_of_interest_g3,
                    "minimal_dist_g3" => minimal_dist_g3,
                    "Vvalues_g3" => Vvalues_g3,
                    "state_order_g3" => state_order_g3
                    )
end
function get_of_interest(graph_info, graph_id)
    if graph_id == 1
        states_of_interest = graph_info["states_of_interest_g1"]
        state_txt = graph_info["state_txt_g1"]
        stepsToGoal = graph_info["stepsToGoal_g1"]
    elseif graph_id == 3
        states_of_interest = graph_info["states_of_interest_g3"]
        state_txt = graph_info["state_txt_g3"]
        stepsToGoal = graph_info["stepsToGoal_g3"]
    end
    states_of_interest, state_txt, stepsToGoal
end


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

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

"""

Lists all available virtual machine sizes for a subscription in a location.
Param: location::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualMachineSizeListResult
"""
function virtualMachineSizesList(_api::VirtualMachineSizesApi, location::String, api_version::String, subscriptionId::String; _mediaType=nothing)

    _ctx = Swagger.Ctx(_api.client, "GET", VirtualMachineSizeListResult, "/subscriptions/{subscriptionId}/providers/Microsoft.Compute/locations/{location}/vmSizes", ["azure_auth"])
    Swagger.set_param(_ctx.path, "location", location)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualMachineSizesList


text:
module HigherKindedPolymorphisms
export inj, prj, @def_hkt, App

const HKP = HigherKindedPolymorphisms
using CanonicalTraits
using MLStyle
@use UppercaseCapturing

function type_constructor_from_hkt end
function type_argument_from_hkt end
function type_app end

# type app representation
struct App{Cons, K₀}
    injected :: Any
end

@trait Higher{Cons, K₀, K₁} where {
    Cons=type_constructor_from_hkt(K₁),
    K₀=type_argument_from_hkt(K₁),
    K₁=type_app(Cons, K₀)
} begin
    inj :: K₁ => App{Cons, K₀}
    inj(data::K₁) = App{Cons, K₀}(data)
    prj :: App{Cons, K₀} => K₁
    prj(data::App{Cons, K₀})::K₁ = data.injected
end

macro def_hkt(expr)
    cg = @when :($Cons{$K₀} = $K₁) = expr begin
        quote
        $HKP.@implement $HKP.Higher{$Cons, $K₀, $K₁}
        $Base.@pure $HKP.type_constructor_from_hkt(::Type{$K₁}) = $Cons
        $Base.@pure $HKP.type_argument_from_hkt(::Type{$K₁}) = $K₀
        $Base.@pure $HKP.type_app(::Type{$Cons}, ::Type{$K₀}) = $K₁
        end
    # NOTE: tparams should occurred in K₁
    @when :($Cons{$K₀} where {$(tparams...)} = $K₁) = expr
        quote
        $HKP.@implement $HKP.Higher{$Cons, $K₀, $K₁} where {$(tparams...)}
        $Base.@pure $HKP.type_constructor_from_hkt(::Type{$K₁}) where {$(tparams...)} = $Cons
        $Base.@pure $HKP.type_argument_from_hkt(::Type{$K₁}) where {$(tparams...)} = $K₀
        $Base.@pure $HKP.type_app(::Type{$Cons}, ::Type{$K₀}) where {$(tparams...)} = $K₁
        end
    end
    esc(cg)
end

end # module


text:
module TestAccessors

using Test
using ..TestParsing: load_sample

function test_default()
    perf = load_sample("default")

    @test perf.context_switches == 124
    @test perf.cpu_migrations == 0
    @test perf.page_faults == 92541
    @test perf.cycles == 3410157225
    @test perf.stalled_cycles_frontend == 74012542
    @test perf.stalled_cycles_backend == 1498764411
    @test perf.instructions == 5107038786
    @test perf.instructions_per_cycle ≈ 1.5 rtol = 0.005
    @test perf.stalled_cycles_per_instructions ≈ 0.29 atol = 0.005

    @test perf.ratio.stalled_cycles_frontend == 74012542 / 3410157225
    @test perf.ratio.stalled_cycles_backend == 1498764411 / 3410157225
end

function test_d1_with_llc()
    perf = load_sample("d1-with-llc")

    @test perf.branches == 402774448
    @test perf.branch_misses == 6375125
    @test perf.l1_dcache_loads == 650213241
    @test perf.l1_dcache_load_misses == 19264237
    @test perf.l1_icache_loads === nothing
    @test perf.l1_icache_load_misses === nothing
    @test perf.llc_loads == 2029190
    @test perf.llc_load_misses == 1158158

    @test perf.ratio.branch_misses == 6375125 / 402774448
    @test perf.ratio.l1_dcache_load_misses == 19264237 / 650213241
    @test perf.ratio.l1_icache_load_misses === nothing
    @test perf.ratio.llc_load_misses == 1158158 / 2029190

    @test perf.percent.branch_misses == 6375125 / 402774448 * 100
    @test perf.percent.l1_dcache_load_misses == 19264237 / 650213241 * 100
    @test perf.percent.l1_icache_load_misses === nothing
    @test perf.percent.llc_load_misses == 1158158 / 2029190 * 100
end

end  # module


text:
using WordTokenizers

@testset "NER" begin
    ner = NERTagger()

    @testset "Basic" begin
        str = "Mr. Foo Bar works in Google, California."
        @test ner(str) == ["O", "PER", "PER", "O", "O", "ORG", "O", "LOC", "O"]

        str = "If the Irish win the World Cup this year, it will be their 3rd time in a row."
        @test ner(str) == [ "O", "O", "MISC", "O", "O", "MISC", "MISC", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O", "O"]
    end

    @testset "Unknown Unicode characters" begin
        # Making sure that the NER model handles for unknown unicode characters
        str = "आ β⬰ 5¥ "
        @test length(ner(str)) == length(WordTokenizers.tokenize(str))

        str = "You owe John Doe 5¥."
        @test ner(str) ==  [ "O", "O", "PER", "PER", "O", "O", "O"]
    end

    @testset "Documents and Corpus" begin
        text1 = "We aRE vErY ClOSE tO ThE HEaDQuarTeRS."
        text2 = "The World Health Organization (WHO) is a specialized agency of the United Nations that is concerned with international public health."

        sd = StringDocument(text1)
        td = TokenDocument(text2)

        tags = ner(sd)
        @test length(tags) == length(WordTokenizers.split_sentences(text1))
        @test length(tags[1]) == length(WordTokenizers.tokenize(text1))
        @test unique(vcat(tags...)) == ["O"]

        tags = ner(td)
        @test length(tags) == length(WordTokenizers.split_sentences(text2))
        @test length(tags[1]) == length(WordTokenizers.tokenize(text2))
        u =  unique(vcat(tags...))
        @test "O" ∈ u && "ORG" ∈ u

        crps = Corpus([sd, td])
        tags = ner(crps)

        @test length(tags) == length(crps.documents)
        @test tags[1] == ner(crps.documents[1])
        @test tags[2] == ner(crps.documents[2])
    end
end


text:
@doc raw"""
Holds the tableau of a stochastic implicit partitioned Runge-Kutta method.
qdrift, pdrift hold the RK coefficients for the drift part,
and qdiff, pdiff hold the RK coefficients for the diffusion part of the SDE.

Order of the tableau is not included, because unlike in the deterministic
setting, it depends on the properties of the noise (e.g., the dimension of
the Wiener process and the commutativity properties of the diffusion matrix).

Orders stored in qdrift and qdiff are understood as the classical orders of these methods.
"""
struct TableauSIPRK{T} <: AbstractTableauIRK{T}
    name::Symbol
    s::Int
    qdrift::CoefficientsRK{T}
    qdiff::CoefficientsRK{T}
    pdrift::CoefficientsRK{T}
    pdiff::CoefficientsRK{T}

    function TableauSIPRK{T}(name, s, qdrift, qdiff, pdrift, pdiff) where {T}
        @assert s == qdrift.s == qdiff.s == pdrift.s == pdiff.s
        new(name, s, qdrift, qdiff, pdrift, pdiff)
    end
end

function TableauSIPRK(name::Symbol, qdrift::CoefficientsRK{T}, qdiff::CoefficientsRK{T}, pdrift::CoefficientsRK{T}, pdiff::CoefficientsRK{T}) where {T}
    TableauSIPRK{T}(name, qdrift.s, qdrift, qdiff, pdrift, pdiff)
end

function TableauSIPRK(name::Symbol, qorder_drift::Int, qa_drift::Matrix{T}, qb_drift::Vector{T}, qc_drift::Vector{T}, qorder_diff::Int, qa_diff::Matrix{T}, qb_diff::Vector{T}, qc_diff::Vector{T},
                                    porder_drift::Int, pa_drift::Matrix{T}, pb_drift::Vector{T}, pc_drift::Vector{T}, porder_diff::Int, pa_diff::Matrix{T}, pb_diff::Vector{T}, pc_diff::Vector{T}) where {T}
    TableauSIPRK{T}(name, length(qc_drift), CoefficientsRK(name, qorder_drift, qa_drift, qb_drift, qc_drift), CoefficientsRK(name, qorder_diff, qa_diff, qb_diff, qc_diff), CoefficientsRK(name, porder_drift, pa_drift, pb_drift, pc_drift), CoefficientsRK(name, porder_diff, pa_diff, pb_diff, pc_diff))
end

# TODO function readTableauSFIRKFromFile(dir::AbstractString, name::AbstractString)

@doc raw"""
Parameters for right-hand side function of implicit Runge-Kutta methods.

`A` - if positive, the upper bound of the Wiener process increments; if A=0.0, no truncation
"""
mutable struct ParametersSIPRK{DT, TT, D, M, S, ET <: NamedTuple} <: Parameters{DT,TT}
    equ::ET
    tab::TableauSIPRK{TT}
    Δt::TT
    ΔW::Vector{DT}
    ΔZ::Vector{DT}
    A::DT

    t::TT
    q::Vector{DT}
    p::Vector{DT}

    function ParametersSIPRK{DT,D,M}(equ::ET, tab::TableauSIPRK{TT}, Δt::TT, ΔW::Vector{DT}, ΔZ::Vector{DT}, A::DT) where {DT, TT, D, M, ET <: NamedTuple}
        @assert M == length(ΔW) == length(ΔZ)
        new{DT, TT, D, M, tab.s, ET}(equ, tab, Δt, ΔW, ΔZ, A, zero(TT), zeros(DT,D), zeros(DT,D))
    end
end

@doc raw"""
Structure for holding the internal stages Q, the values of the drift vector
and the diffusion matrix evaluated at the internal stages `V=v(Q)`, `B=B(Q)`,
and the increments `Y = Δt*a_drift*v(Q) + a_diff*B(Q)*ΔW`.
"""
struct IntegratorCacheSIPRK{DT,D,M,S} <: PSDEIntegratorCache{DT,D,M}
    Q::Vector{Vector{DT}}
    P::Vector{Vector{DT}}
    V::Vector{Vector{DT}}
    F::Vector{Vector{DT}}
    B::Vector{Matrix{DT}}
    G::Vector{Matrix{DT}}
    Y::Vector{Vector{DT}}
    Z::Vector{Vector{DT}}

    v::Vector{DT}
    f::Vector{DT}
    b::Matrix{DT}
    g::Matrix{DT}
    y::Vector{DT}
    z::Vector{DT}

    V1::Vector{DT}
    V2::Vector{DT}
    F1::Vector{DT}
    F2::Vector{DT}
    B1::Matrix{DT}
    B2::Matrix{DT}
    G1::Matrix{DT}
    G2::Matrix{DT}
    Δw::Vector{DT}
    ΔQ::Vector{DT}
    Y1::Vector{DT}
    Y2::Vector{DT}
    ΔP::Vector{DT}
    Z1::Vector{DT}
    Z2::Vector{DT}
    Δy::Vector{DT}
    Δz::Vector{DT}

    function IntegratorCacheSIPRK{DT,D,M,S}() where {DT,D,M,S}

        # create internal stage vectors
        Q = create_internal_stage_vector(DT, D, S)
        P = create_internal_stage_vector(DT, D, S)
        V = create_internal_stage_vector(DT, D, S)
        F = create_internal_stage_vector(DT, D, S)
        B = create_internal_stage_matrix(DT, D, M, S)
        G = create_internal_stage_matrix(DT, D, M, S)
        Y = create_internal_stage_vector(DT, D, S)
        Z = create_internal_stage_vector(DT, D, S)

        # create velocity and update vector
        v = zeros(DT,D)
        f = zeros(DT,D)
        b = zeros(DT,D,M)
        g = zeros(DT,D,M)
        y = zeros(DT,D)
        z = zeros(DT,D)

        # create temporary arrays
        V1 = zeros(DT,D)
        V2 = zeros(DT,D)
        F1 = zeros(DT,D)
        F2 = zeros(DT,D)
        B1 = zeros(DT,D,M)
        B2 = zeros(DT,D,M)
        G1 = zeros(DT,D,M)
        G2 = zeros(DT,D,M)
        ΔQ = zeros(DT,D)
        Y1 = zeros(DT,D)
        Y2 = zeros(DT,D)
        ΔP = zeros(DT,D)
        Z1 = zeros(DT,D)
        Z2 = zeros(DT,D)
        Δw = zeros(DT,M)
        Δy = zeros(DT,M)
        Δz = zeros(DT,M)

        new(Q, P, V, F, B, G, Y, Z, v, f, b, g, y, z, V1, V2, F1, F2, B1, B2, G1, G2, Δw, ΔQ, Y1, Y2, ΔP, Z1, Z2, Δy, Δz)
    end
end

function Integrators.IntegratorCache{ST}(params::ParametersSIPRK{DT,TT,D,M,S}; kwargs...) where {ST,DT,TT,D,M,S}
    IntegratorCacheSIPRK{ST,D,M,S}(; kwargs...)
end

@inline Integrators.CacheType(ST, params::ParametersSIPRK{DT,TT,D,M,S}) where {DT,TT,D,M,S} = IntegratorCacheSIPRK{ST,D,M,S}

"Stochastic implicit partitioned Runge-Kutta integrator."
struct IntegratorSIPRK{DT, TT, D, M, S,
                PT <: ParametersSIPRK{DT,TT},
                ST <: NonlinearSolver{DT}} <: StochasticIntegratorPRK{DT,TT,D,M,S}
    params::PT
    solver::ST
    caches::CacheDict{PT}

    function IntegratorSIPRK(params::ParametersSIPRK{DT,TT,D,M,S}, solver::ST, caches) where {DT,TT,D,M,S,ST}
        new{DT, TT, D, M, S, typeof(params), ST}(params, solver, caches)
    end

    function IntegratorSIPRK{DT,D,M}(equations::NamedTuple, tableau::TableauSIPRK{TT}, Δt::TT; K::Int=0) where {DT,TT,D,M}
        # get number of stages
        S = tableau.s

        # K - the integer in the bound A = √(2 K Δt |log Δt|) due to Milstein & Tretyakov; K=0 no truncation
        K==0 ? A = 0.0 : A = sqrt( 2*K*Δt*abs(log(Δt)) )

        # create params
        params = ParametersSIPRK{DT,D,M}(equations, tableau, Δt, zeros(DT,M), zeros(DT,M), A)

        # create cache dict
        caches = CacheDict(params)

        # create solver
        solver = create_nonlinear_solver(DT, 2*D*S, params, caches)

        # create integrator
        IntegratorSIPRK(params, solver, caches)
    end

    function IntegratorSIPRK(equation::PSDE{DT,TT}, tableau::TableauSIPRK{TT}, Δt::TT; kwargs...) where {DT,TT}
        IntegratorSIPRK{DT, ndims(equation), equation.m}(get_function_tuple(equation), tableau, Δt; kwargs...)
    end
end

@doc raw"""
This function computes initial guesses for `Y`, `Z` and assigns them to int.solver.x
The prediction is calculated using an explicit integrator.

SIMPLE SOLUTION
The simplest initial guess for `Y`, `Z` is 0:
`int.solver.x .= zeros(eltype(int), 2*tableau(int).s*ndims(int))`

USING AN EXPLICIT INTEGRATOR TO COMPUTE AN INITIAL GUESS
Below we use the R2 method of Burrage & Burrage to calculate
the internal stages at the times c[1]...c[s].
This approach seems to give very good approximations if the time step
and magnitude of noise are not too large. If the noise intensity is too big,
one may have to perform a few iterations of the explicit method with a smaller
time step, use a higher-order explicit method (e.g. CL or G5), or use
the simple solution above.

When calling this function, int.params should contain the data:
`int.params.q` - the q solution at the previous time step
`int.params.p` - the p solution at the previous time step
`int.params.t` - the time of the previous step
`int.params.ΔW`- the increment of the Brownian motion for the current step
"""
function initial_guess!(int::IntegratorSIPRK{DT,TT}, sol::AtomicSolutionPSDE{DT,TT},
                        cache::IntegratorCacheSIPRK{DT}=int.caches[DT]) where {DT,TT}

    local t2::TT
    local Δt_local::TT

    # Evaluating the functions v and B at t,q - same for all stages
    int.params.equ[:v](int.params.t, int.params.q, int.params.p, cache.V1)
    int.params.equ[:B](int.params.t, int.params.q, int.params.p, cache.B1)
    int.params.equ[:f](int.params.t, int.params.q, int.params.p, cache.F1)
    int.params.equ[:G](int.params.t, int.params.q, int.params.p, cache.G1)

    # Calculating the positions q at the points qdrift.c[i]
    # if qdrift.c==pdrift.c, then also calculating the momenta
    for i in eachstage(int)
        # Taking the c[i] from the qdrift tableau.
        Δt_local = tableau(int).qdrift.c[i]  * timestep(int)
        cache.Δw .= tableau(int).qdrift.c[i] .* int.params.ΔW

        mul!(cache.ΔQ, cache.B1, cache.Δw)
        mul!(cache.ΔP, cache.G1, cache.Δw)
        @. cache.Q[i] = int.params.q + 2. / 3. * Δt_local * cache.V1 + 2. / 3. * cache.ΔQ
        @. cache.P[i] = int.params.p + 2. / 3. * Δt_local * cache.F1 + 2. / 3. * cache.ΔP

        t2 = int.params.t + 2. / 3. * Δt_local

        int.params.equ.v(t2, cache.Q[i], cache.P[i], cache.V2)
        int.params.equ.B(t2, cache.Q[i], cache.P[i], cache.B2)
        int.params.equ.f(t2, cache.Q[i], cache.P[i], cache.F2)
        int.params.equ.G(t2, cache.Q[i], cache.P[i], cache.G2)

        #Calculating the Y's and assigning them to the array int.solver.x as initial guesses
        mul!(cache.Y1, cache.B1, cache.Δw)
        mul!(cache.Y2, cache.B2, cache.Δw)
        for j in eachdim(int)
            int.solver.x[(i-1)*ndims(int)+j] =  Δt_local*(1. / 4. * cache.V1[j] + 3. / 4. * cache.V2[j]) + 1. / 4. * cache.Y1[j] + 3. / 4. * cache.Y2[j]
        end

        # if the collocation points are the same for both q and p parts
        mul!(cache.Z1, cache.G1, cache.Δw)
        mul!(cache.Z2, cache.G2, cache.Δw)
        if tableau(int).qdrift.c==tableau(int).pdrift.c
            for j in 1:ndims(int)
                int.solver.x[(tableau(int).s+i-1)*ndims(int)+j] =  Δt_local*(1. / 4. * cache.F1[j] + 3. / 4. * cache.F2[j]) + 1. / 4. * cache.Z1[j] + 3. / 4. * cache.Z2[j]
            end
        end
    end

    # If qdrift.c != pdrift.c, then calculating the momenta p at the points pdrift.c[i]
    if tableau(int).qdrift.c != tableau(int).pdrift.c
        for i in eachstage(int)
            # Taking the c[i] from the pdrift tableau.
            Δt_local = tableau(int).pdrift.c[i]  * timestep(int)
            cache.Δw .= tableau(int).pdrift.c[i] .* int.params.ΔW

            mul!(cache.ΔQ, cache.B1, cache.Δw)
            mul!(cache.ΔP, cache.G1, cache.Δw)
            @. cache.Q[i] = int.params.q + 2. / 3. * Δt_local * cache.V1 + 2. / 3. * cache.ΔQ
            @. cache.P[i] = int.params.p + 2. / 3. * Δt_local * cache.F1 + 2. / 3. * cache.ΔP

            t2 = int.params.t + 2. / 3. * Δt_local

            int.params.equ[:v](t2, cache.Q[i], cache.P[i], cache.V2)
            int.params.equ[:B](t2, cache.Q[i], cache.P[i], cache.B2)
            int.params.equ[:f](t2, cache.Q[i], cache.P[i], cache.F2)
            int.params.equ[:G](t2, cache.Q[i], cache.P[i], cache.G2)

            # Calculating the Z's and assigning them to the array int.solver.x as initial guesses
            # The guesses for the Y's have already been written to x above
            mul!(cache.Z1, cache.G1, cache.Δw)
            mul!(cache.Z2, cache.G2, cache.Δw)
            for j in eachdim(int)
                int.solver.x[(tableau(int).s+i-1)*ndims(int)+j] =  Δt_local*(1. / 4. * cache.F1[j] + 3. / 4. * cache.F2[j]) + 1. / 4. * cache.Z1[j] + 3. / 4. * cache.Z2[j]
            end
        end
    end
end

@doc raw"""
Unpacks the data stored in `x = (Y[1][1], Y[1][2], ... Y[1][D], Y[2][1], ..., Z[1][1], Z[1][2], ... Z[1][D], Z[2][1], ...)`
into `Y`, `Z`, calculates the internal stages `Q`, `P`, the values of the RHS
of the SDE ( `v(Q,P)`, `f(Q,P)`, `B(Q,P)` and `G(Q,P)` ), and assigns them to `V`, `F`, `B` and `G`.
Unlike for FIRK, here
`Y = Δt a_drift v(Q,P) + a_diff B(Q,P) ΔW`,
`Z = Δt â_drift v(Q,P) + â_diff B(Q,P) ΔW`.
"""
function compute_stages!(x::Vector{ST}, Q::Vector{Vector{ST}}, P::Vector{Vector{ST}},
                                        V::Vector{Vector{ST}}, F::Vector{Vector{ST}},
                                        B::Vector{Matrix{ST}}, G::Vector{Matrix{ST}},
                                        Y::Vector{Vector{ST}}, Z::Vector{Vector{ST}},
                                        params::ParametersSIPRK{DT,TT,D,M,S}) where {ST,DT,TT,D,M,S}

    local tqᵢ::TT       #times for the q internal stages
    local tpᵢ::TT       #times for the p internal stages

    @assert S == length(Q) == length(V) == length(B) == length(P) == length(F) == length(G)

    # copy x to Y, Z and calculate Q, P
    for i in 1:S
        @assert D == length(Q[i]) == length(V[i]) == size(B[i],1) == length(P[i]) == length(F[i]) == size(G[i],1)
        @assert M == size(B[i],2) == size(G[i],2)

        for k in 1:D
            Y[i][k] = x[D*(  i-1)+k]
            Z[i][k] = x[D*(S+i-1)+k]
        end

        Q[i] .= params.q .+ Y[i]
        P[i] .= params.p .+ Z[i]
    end

    # compute V=v(t,Q,P), F=f(t,Q,P), B=B(t,Q,P) and G=G(t,Q,P)
    for i in 1:S
        tqᵢ = params.t + params.Δt * params.tab.qdrift.c[i]
        tpᵢ = params.t + params.Δt * params.tab.pdrift.c[i]
        # calculates v(t,Q,P) and f(t,Q,P) and assigns to the i-th column of V and F
        params.equ[:v](tqᵢ, Q[i], P[i], V[i])
        params.equ[:f](tpᵢ, Q[i], P[i], F[i])
        # calculates B(t,Q,P) and assigns to the matrices B[i] and G[i]
        params.equ[:B](tqᵢ, Q[i], P[i], B[i])
        params.equ[:G](tpᵢ, Q[i], P[i], G[i])
    end
end

"Compute stages of implicit Runge-Kutta methods."
function Integrators.function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersSIPRK{DT,TT,D,M,S},
                caches::CacheDict) where {ST,DT,TT,D,M,S}

    # get cache for internal stages
    cache = caches[ST]

    compute_stages!(x, cache.Q, cache.P, cache.V, cache.F, cache.B, cache.G, cache.Y, cache.Z, params)

    local y1::ST
    local y2::ST
    local z1::ST
    local z2::ST

    # compute b = - (Y-AV)
    for i in 1:S
        for k in 1:D
            y1 = 0
            y2 = 0
            z1 = 0
            z2 = 0
            for j in 1:S
                y1 += params.tab.qdrift.a[i,j] * cache.V[j][k] * params.Δt
                y2 += params.tab.qdrift.â[i,j] * cache.V[j][k] * params.Δt
                z1 += params.tab.pdrift.a[i,j] * cache.F[j][k] * params.Δt
                z2 += params.tab.pdrift.â[i,j] * cache.F[j][k] * params.Δt
                for l in 1:M
                    y1 += params.tab.qdiff.a[i,j] * cache.B[j][k,l] * params.ΔW[l]
                    y2 += params.tab.qdiff.â[i,j] * cache.B[j][k,l] * params.ΔW[l]
                    z1 += params.tab.pdiff.a[i,j] * cache.G[j][k,l] * params.ΔW[l]
                    z2 += params.tab.pdiff.â[i,j] * cache.G[j][k,l] * params.ΔW[l]
                end
            end
            b[D*(  i-1)+k] = - cache.Y[i][k] + (y1 + y2)
            b[D*(S+i-1)+k] = - cache.Z[i][k] + (z1 + z2)
        end
    end
end

function update_solution!(sol::AtomicSolutionPSDE{T}, tab::TableauSIPRK{T}, Δt::T, ΔW::Vector{T},
                          cache::IntegratorCacheSIPRK{T}) where {T}

    update_solution!(sol, cache.V, cache.F, cache.B, cache.G,
                     tab.qdrift.b, tab.qdiff.b,
                     tab.pdrift.b, tab.pdiff.b,
                     Δt, ΔW, cache.Δy, cache.Δz)

    update_solution!(sol, cache.V, cache.F, cache.B, cache.G,
                     tab.qdrift.b̂, tab.qdiff.b̂,
                     tab.pdrift.b̂, tab.pdiff.b̂,
                     Δt, ΔW, cache.Δy, cache.Δz)
end

"Integrate PSDE with a stochastic implicit partitioned Runge-Kutta integrator."
function Integrators.integrate_step!(int::IntegratorSIPRK{DT,TT}, sol::AtomicSolutionPSDE{DT,TT},
                                     cache::IntegratorCacheSIPRK{DT}=int.caches[DT]) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int, sol)

    # compute initial guess and assign to int.solver.x
    initial_guess!(int, sol, cache)

    # reset cache
    reset!(sol, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params)

    # compute the drift vector field and the diffusion matrix at internal stages
    compute_stages!(int.solver.x, cache.Q, cache.P, cache.V, cache.F, cache.B, cache.G, cache.Y, cache.Z, int.params)

    # compute final update
    update_solution!(sol, tableau(int), timestep(int), int.params.ΔW, cache)
end


text:
function metropolis_hastings(sel::K,
                             call::C) where {K <: UnconstrainedSelection, C <: CallSite}
    ret, cl, w, retdiff, d = regenerate(sel, call)
    log(rand()) < w && return (cl, true)
    return (call, false)
end

function metropolis_hastings(sel::K,
                             call::C,
                             proposal::Function,
                             proposal_args::Tuple) where {K <: UnconstrainedSelection, C <: CallSite}
    p_ret, p_cl, p_w = propose(proposal, call, proposal_args...)
    s = selection(p_cl)
    u_ret, u_cl, u_w, retdiff, d = update(s, call.fn, NoChange(), call.args...)
    d_s = selection(d)
    s_ret, s_w = score(d_s, proposal, u_cl, proposal_args...)
    ratio = u_w - p_w + s_w
    log(rand()) < ratio && return (u_cl, true)
    return (call, false)
end

function metropolis_hastings!(sel::K,
                              ps::Particles) where K <: UnconstrainedSelection
    num_particles = length(ps)
    Threads.@threads for i in 1 : num_particles
        old_w = ps.lws[i]
        old_score = get_score(ps.calls[i])
        ps.calls[i], _ = metropolis_hastings(sel, ps.calls[i])
        ps.lws[i] = old_w + get_score(ps.calls[i]) - old_score
    end
end

function metropolis_hastings!(sel::K,
                              ps::Particles,
                              proposal::Function,
                              proposal_args::Tuple) where K <: UnconstrainedSelection
    num_particles = length(ps)
    Threads.@threads for i in 1 : num_particles
        old_w = ps.lws[i]
        old_score = get_score(ps.calls[i])
        ps.calls[i], _ = metropolis_hastings(sel, ps.calls[i], proposal, proposal_args)
        ps.lws[i] = old_w + get_score(ps.calls[i]) - old_score
    end
end