text:
abstract type AbstractSpectralStyle <: Fields.AbstractFieldStyle end

struct SpectralStyle <: AbstractSpectralStyle end

struct CompositeSpectralStyle <: AbstractSpectralStyle end

"""
    SpectralElementOperator

Represents an operation that is applied to each element.

Subtypes `Op` of this should define the following:
- [`operator_return_eltype(::Op, ElTypes...)`](@ref)
- [`allocate_work(::Op, args...)`](@ref)
- [`apply_slab(::Op, work, args...)`](@ref)

Additionally, the result type `OpResult <: OperatorSlabResult` of `apply_slab` should define `get_node(::OpResult, i, j)`.
"""
abstract type SpectralElementOperator end

"""
    operator_axes(space)

Return a tuple of the axis indicies a given field operator works over.
"""
function operator_axes end

operator_axes(space::Spaces.AbstractSpace) = ()
operator_axes(space::Spaces.SpectralElementSpace1D) = (1,)
operator_axes(space::Spaces.SpectralElementSpace2D) = (1, 2)
operator_axes(space::Spaces.ExtrudedFiniteDifferenceSpace) =
    operator_axes(space.horizontal_space)

"""
    SpectralBroadcasted{Style}(op, args[,axes[, work]])

This is similar to a `Base.Broadcast.Broadcasted` object, except it contains space for an intermediate `work` storage.

This is returned by `Base.Broadcast.broadcasted(op::SpectralElementOperator)`.
"""
struct SpectralBroadcasted{Style, Op, Args, Axes, InputSpace} <:
       Base.AbstractBroadcasted
    op::Op
    args::Args
    axes::Axes
    input_space::InputSpace
end
SpectralBroadcasted{Style}(
    op::Op,
    args::Args,
    axes::Axes = nothing,
    input_space::InputSpace = nothing,
) where {Style, Op, Args, Axes, InputSpace} =
    SpectralBroadcasted{Style, Op, Args, Axes, InputSpace}(
        op,
        args,
        axes,
        input_space,
    )

input_space(arg) = axes(arg)
input_space(::SpectralElementOperator, space) = space

input_space(sbc::SpectralBroadcasted) =
    isnothing(sbc.input_space) ?
    input_space(sbc.op, tuplemap(axes, sbc.args)...) : sbc.input_space

return_space(::SpectralElementOperator, space) = space

Base.axes(sbc::SpectralBroadcasted) =
    isnothing(sbc.axes) ? return_space(sbc.op, tuplemap(axes, sbc.args)...) :
    sbc.axes

Base.Broadcast.broadcasted(op::SpectralElementOperator, args...) =
    Base.Broadcast.broadcasted(SpectralStyle(), op, args...)

Base.Broadcast.broadcasted(
    ::SpectralStyle,
    op::SpectralElementOperator,
    args...,
) = SpectralBroadcasted{SpectralStyle}(op, args)

Base.eltype(sbc::SpectralBroadcasted) =
    operator_return_eltype(sbc.op, tuplemap(eltype, sbc.args)...)

function Base.Broadcast.instantiate(
    sbc::SpectralBroadcasted{Style},
) where {Style}
    op = sbc.op
    # recursively instantiate the arguments to allocate intermediate work arrays
    args = tuplemap(Base.Broadcast.instantiate, sbc.args)
    # axes: same logic as Broadcasted
    if sbc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        axes = Base.axes(sbc)
    else
        axes = sbc.axes
        Base.Broadcast.check_broadcast_axes(axes, args...)
    end
    if sbc.input_space isa Nothing
        inspace = input_space(sbc)
    else
        inspace = sbc.input_space
    end
    op = typeof(op)(axes)
    return SpectralBroadcasted{Style}(op, args, axes, inspace)
end

function Base.Broadcast.instantiate(
    bc::Base.Broadcast.Broadcasted{Style},
) where {Style <: AbstractSpectralStyle}
    # recursively instantiate the arguments to allocate intermediate work arrays
    args = tuplemap(Base.Broadcast.instantiate, bc.args)
    # axes: same logic as Broadcasted
    if bc.axes isa Nothing # Not done via dispatch to make it easier to extend instantiate(::Broadcasted{Style})
        axes = Base.Broadcast.combine_axes(args...)
    else
        axes = bc.axes
        Base.Broadcast.check_broadcast_axes(axes, args...)
    end
    return Base.Broadcast.Broadcasted{Style}(bc.f, args, axes)
end

function Base.similar(sbc::SpectralBroadcasted, ::Type{Eltype}) where {Eltype}
    space = axes(sbc)
    return Field(Eltype, space)
end

function Base.similar(
    bc::Base.Broadcast.Broadcasted{<:AbstractSpectralStyle},
    ::Type{Eltype},
) where {Eltype}
    space = axes(bc)
    return Field(Eltype, space)
end

function Base.copy(sbc::SpectralBroadcasted)
    # figure out return type
    dest = similar(sbc, eltype(sbc))
    # allocate dest
    copyto!(dest, sbc)
end
Base.Broadcast.broadcastable(sbc::SpectralBroadcasted) = sbc

function Base.Broadcast.materialize(sbc::SpectralBroadcasted)
    copy(Base.Broadcast.instantiate(sbc))
end

function Base.copyto!(field_out::Field, sbc::SpectralBroadcasted)
    data_out = Fields.field_values(field_out)
    space = axes(field_out)
    Nv = Spaces.nlevels(space)
    Nh = Topologies.nlocalelems(Spaces.topology(space))
    @inbounds for h in 1:Nh, v in 1:Nv
        slab_out = slab(field_out, v, h)
        out_slab_space = slab(axes(sbc), v, h)
        in_slab_space = slab(input_space(sbc), v, h)
        _slab_args = _apply_slab_args(slab_args(sbc.args, v, h), v, h)
        copy_slab!(
            slab_out,
            apply_slab(sbc.op, out_slab_space, in_slab_space, _slab_args...),
        )
    end
    return field_out
end

function Base.Broadcast.materialize!(dest, sbc::SpectralBroadcasted)
    copyto!(dest, Base.Broadcast.instantiate(sbc))
end

function slab(
    sbc::SpectralBroadcasted{Style},
    inds...,
) where {Style <: SpectralStyle}
    _args = slab_args(sbc.args, inds...)
    _axes = slab(axes(sbc), inds...)
    SpectralBroadcasted{Style}(sbc.op, _args, _axes, sbc.input_space)
end

function slab(
    bc::Base.Broadcast.Broadcasted{Style},
    inds...,
) where {Style <: AbstractSpectralStyle}
    _args = slab_args(bc.args, inds...)
    _axes = slab(axes(bc), inds...)
    Base.Broadcast.Broadcasted{Style}(bc.f, _args, _axes)
end

function Base.copyto!(
    field_out::Field,
    bc::Base.Broadcast.Broadcasted{Style},
) where {Style <: AbstractSpectralStyle}
    space = axes(field_out)
    topology = Spaces.topology(space)
    Nv = Spaces.nlevels(space)::Int
    Nh = Topologies.nlocalelems(topology)::Int
    @inbounds for h in 1:Nh, v in 1:Nv
        slab_out = slab(field_out, v, h)
        _slab_args = _apply_slab_args(slab_args(bc.args, v, h), v, h)
        copy_slab!(
            slab_out,
            Base.Broadcast.Broadcasted{Style}(bc.f, _slab_args, axes(slab_out)),
        )
    end
    return field_out
end

function copy_slab!(slab_out::Fields.SlabField1D, res)
    space = axes(slab_out)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)::Int
    @inbounds for i in 1:Nq
        set_node!(slab_out, i, get_node(res, i))
    end
    return slab_out
end

function copy_slab!(slab_out::Fields.SlabField2D, res)
    space = axes(slab_out)
    Nq = Quadratures.degrees_of_freedom(space.quadrature_style)::Int
    @inbounds for j in 1:Nq, i in 1:Nq
        set_node!(slab_out, i, j, get_node(res, i, j))
    end
    return slab_out
end

# 1D get/set node

# 1D intermediate slab data types
Base.@propagate_inbounds function get_node(v::MVector, i)
    v[i]
end

Base.@propagate_inbounds function get_node(v::SVector, i)
    v[i]
end

Base.@propagate_inbounds function get_node(scalar, i)
    scalar[]
end

Base.@propagate_inbounds function get_node(field::Fields.SlabField1D, i)
    getindex(Fields.field_values(field), i)
end

Base.@propagate_inbounds function get_node(bc::Base.Broadcast.Broadcasted, i)
    args = tuplemap(bc.args) do arg
        Base.@_propagate_inbounds_meta
        get_node(arg, i)
    end
    bc.f(args...)
end

Base.@propagate_inbounds set_node!(field::Fields.SlabField1D, i, val) =
    setindex!(Fields.field_values(field), val, i)

# 2D get/set node
Base.@propagate_inbounds function get_node(v::MMatrix, i, j)
    v[i, j]
end

Base.@propagate_inbounds function get_node(v::SMatrix, i, j)
    v[i, j]
end

Base.@propagate_inbounds function get_node(scalar, i, j)
    scalar[]
end

Base.@propagate_inbounds function get_node(field::Fields.SlabField2D, i, j)
    getindex(Fields.field_values(field), i, j)
end

Base.@propagate_inbounds function get_node(bc::Base.Broadcast.Broadcasted, i, j)
    args = tuplemap(bc.args) do arg
        Base.@_propagate_inbounds_meta
        get_node(arg, i, j)
    end
    bc.f(args...)
end

Base.@propagate_inbounds function set_node!(
    field::Fields.SlabField2D,
    i,
    j,
    val,
)
    setindex!(Fields.field_values(field), val, i, j)
end

# Broadcast recursive machinery
@inline _apply_slab_args(args::Tuple, inds...) = (
    _apply_slab(args[1], inds...),
    _apply_slab_args(Base.tail(args), inds...)...,
)
@inline _apply_slab_args(args::Tuple{Any}, inds...) =
    (_apply_slab(args[1], inds...),)
@inline _apply_slab_args(args::Tuple{}, inds...) = ()

@inline _apply_slab(x, inds...) = x
@inline _apply_slab(sbc::SpectralBroadcasted, inds...) = apply_slab(
    sbc.op,
    slab(axes(sbc), inds...),
    slab(input_space(sbc), inds...),
    _apply_slab_args(sbc.args, inds...)...,
)
@inline _apply_slab(
    bc::Base.Broadcast.Broadcasted{CompositeSpectralStyle},
    inds...,
) = Base.Broadcast.Broadcasted{CompositeSpectralStyle}(
    bc.f,
    _apply_slab_args(bc.args, inds...),
    bc.axes,
)

function Base.Broadcast.BroadcastStyle(
    ::Type{SB},
) where {SB <: SpectralBroadcasted}
    CompositeSpectralStyle()
end

function Base.Broadcast.BroadcastStyle(
    ::CompositeSpectralStyle,
    ::Fields.AbstractFieldStyle,
)
    CompositeSpectralStyle()
end

"""
    div = Divergence()
    div.(u)

Computes the per-element spectral (strong) divergence of a vector field ``u``.

The divergence of a vector field ``u`` is defined as
```math
\\nabla \\cdot u = \\sum_i \\frac{1}{J} \\frac{\\partial (J u^i)}{\\partial \\xi^i}
```
where ``J`` is the Jacobian determinant, ``u^i`` is the ``i``th contravariant
component of ``u``.

This is discretized by
```math
\\sum_i I \\left\\{\\frac{1}{J} \\frac{\\partial (I\\{J u^i\\})}{\\partial \\xi^i} \\right\\}
```
where ``I\\{x\\}`` is the interpolation operator that projects to the
unique polynomial interpolating ``x`` at the quadrature points. In matrix
form, this can be written as
```math
J^{-1} \\sum_i D_i J u^i
```
where ``D_i`` is the derivative matrix along the ``i``th dimension

## References
- [Taylor2010](@cite), equation 15
"""
struct Divergence{I} <: SpectralElementOperator end
Divergence() = Divergence{()}()
Divergence{()}(space) = Divergence{operator_axes(space)}()

operator_return_eltype(op::Divergence{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(Geometry.divergence_result_type, S)

function apply_slab(op::Divergence{(1,)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        Jv¹ =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i),
            )
        for ii in 1:Nq
            out[ii] = out[ii] ⊞ (D[ii, i] ⊠ Jv¹)
        end
    end
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        out[i] = RecursiveApply.rdiv(out[i], local_geometry.J)
    end
    return SVector(out)
end

function apply_slab(op::Divergence{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        Jv¹ =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for ii in 1:Nq
            out[ii, j] = out[ii, j] ⊞ (D[ii, i] ⊠ Jv¹)
        end
        Jv² =
            local_geometry.J ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant2(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for jj in 1:Nq
            out[i, jj] = out[i, jj] ⊞ (D[jj, j] ⊠ Jv²)
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
    end
    return SMatrix(out)
end

"""
    wdiv = WeakDivergence()
    wdiv.(u)

Computes the "weak divergence" of a vector field `u`.

This is defined as the scalar field ``\\theta \\in \\mathcal{V}_0`` such that
for all ``\\phi\\in \\mathcal{V}_0``
```math
\\int_\\Omega \\phi \\theta \\, d \\Omega
=
- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
```
where ``\\mathcal{V}_0`` is the space of ``u``.

This arises as the contribution of the volume integral after by applying
integration by parts to the weak form expression of the divergence
```math
\\int_\\Omega \\phi (\\nabla \\cdot u) \\, d \\Omega
=
- \\int_\\Omega (\\nabla \\phi) \\cdot u \\,d \\Omega
+ \\oint_{\\partial \\Omega} \\phi (u \\cdot n) \\,d \\sigma
```

It can be written in matrix form as
```math
ϕ^\\top WJ θ = - \\sum_i (D_i ϕ)^\\top WJ u^i
```
which reduces to
```math
θ = -(WJ)^{-1} \\sum_i D_i^\\top WJ u^i
```
where
 - ``J`` is the diagonal Jacobian matrix
 - ``W`` is the diagonal matrix of quadrature weights
 - ``D_i`` is the derivative matrix along the ``i``th dimension
"""
struct WeakDivergence{I} <: SpectralElementOperator end
WeakDivergence() = WeakDivergence{()}()
WeakDivergence{()}(space) = WeakDivergence{operator_axes(space)}()

operator_return_eltype(::WeakDivergence{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(Geometry.divergence_result_type, S)

function apply_slab(op::WeakDivergence{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    RT = operator_return_eltype(op, eltype(slab_data))
    # allocate temp output
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        WJv¹ =
            local_geometry.WJ ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant1(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for ii in 1:Nq
            out[ii, j] = out[ii, j] ⊞ (D[i, ii] ⊠ WJv¹)
        end
        WJv² =
            local_geometry.WJ ⊠ RecursiveApply.rmap(
                v -> Geometry.contravariant2(v, local_geometry),
                get_node(slab_data, i, j),
            )
        for jj in 1:Nq
            out[i, jj] = out[i, jj] ⊞ (D[j, jj] ⊠ WJv²)
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], ⊟(local_geometry.WJ))
    end
    return SMatrix(out)
end

"""
    grad = Gradient()
    grad.(f)

Compute the (strong) gradient of `f` on each element, returning a
`CovariantVector`-field.

The ``i``th covariant component of the gradient is the partial derivative with
respect to the reference element:
```math
(\\nabla f)_i = \\frac{\\partial f}{\\partial \\xi^i}
```

Discretely, this can be written in matrix form as
```math
D_i f
```
where ``D_i`` is the derivative matrix along the ``i``th dimension.

## References
- [Taylor2010](@cite), equation 16
"""
struct Gradient{I} <: SpectralElementOperator end
Gradient() = Gradient{()}()
Gradient{()}(space) = Gradient{operator_axes(space)}()

operator_return_eltype(::Gradient{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)

function apply_slab(op::Gradient{(1,)}, slab_space, _, slab_data)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        x = get_node(slab_data, i)
        for ii in 1:Nq
            ∂f∂ξ = Geometry.Covariant1Vector(D[ii, i]) ⊗ x
            out[ii] += ∂f∂ξ
        end
    end
    return SVector(out)
end

function apply_slab(op::Gradient{(1, 2)}, slab_space, _, slab_data)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        x = get_node(slab_data, i, j)
        for ii in 1:Nq
            ∂f∂ξ₁ = Geometry.Covariant12Vector(D[ii, i], zero(eltype(D))) ⊗ x
            out[ii, j] = out[ii, j] ⊞ ∂f∂ξ₁
        end
        for jj in 1:Nq
            ∂f∂ξ₂ = Geometry.Covariant12Vector(zero(eltype(D)), D[jj, j]) ⊗ x
            out[i, jj] = out[i, jj] ⊞ ∂f∂ξ₂
        end
    end
    return SMatrix(out)
end

"""
    wgrad = WeakGradient()
    wgrad.(f)

Compute the "weak gradient" of `f` on each element.

This is defined as the the vector field ``\\theta \\in \\mathcal{V}_0`` such
that for all ``\\phi \\in \\mathcal{V}_0``
```math
\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
=
- \\int_\\Omega (\\nabla \\cdot \\phi) f \\, d\\Omega
```
where ``\\mathcal{V}_0`` is the space of ``f``.

This arises from the contribution of the volume integral after by applying
integration by parts to the weak form expression of the gradient
```math
\\int_\\Omega \\phi \\cdot (\\nabla f) \\, d \\Omega
=
- \\int_\\Omega f (\\nabla \\cdot \\phi) \\, d\\Omega
+ \\oint_{\\partial \\Omega} f (\\phi \\cdot n) \\, d \\sigma
```

In matrix form, this becomes
```math
{\\phi^i}^\\top W J \\theta_i = - ( J^{-1} D_i J \\phi^i )^\\top W J f
```
which reduces to
```math
\\theta_i = -W^{-1} D_i^\\top W f
```
where ``D_i`` is the derivative matrix along the ``i``th dimension.
"""
struct WeakGradient{I} <: SpectralElementOperator end
WeakGradient() = WeakGradient{()}()
WeakGradient{()}(space) = WeakGradient{operator_axes(space)}()

operator_return_eltype(::WeakGradient{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(T -> Geometry.gradient_result_type(Val(I), T), S)

function apply_slab(op::WeakGradient{(1,)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        W = local_geometry.WJ / local_geometry.J
        Wx = W ⊠ get_node(slab_data, i)
        for ii in 1:Nq
            Dᵀ₁Wf = Geometry.Covariant1Vector(D[i, ii]) ⊗ Wx
            out[ii] = out[ii] ⊟ Dᵀ₁Wf
        end
    end
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        W = local_geometry.WJ / local_geometry.J
        out[i] = RecursiveApply.rdiv(out[i], W)
    end
    return SVector(out)
end

function apply_slab(op::WeakGradient{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        W = local_geometry.WJ / local_geometry.J
        Wx = W ⊠ get_node(slab_data, i, j)
        for ii in 1:Nq
            Dᵀ₁Wf = Geometry.Covariant12Vector(D[i, ii], zero(eltype(D))) ⊗ Wx
            out[ii, j] = out[ii, j] ⊟ Dᵀ₁Wf
        end
        for jj in 1:Nq
            Dᵀ₂Wf = Geometry.Covariant12Vector(zero(eltype(D)), D[j, jj]) ⊗ Wx
            out[i, jj] = out[i, jj] ⊟ Dᵀ₂Wf
        end
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        W = local_geometry.WJ / local_geometry.J
        out[i, j] = RecursiveApply.rdiv(out[i, j], W)
    end
    return SMatrix(out)
end

abstract type CurlSpectralElementOperator <: SpectralElementOperator end

"""
    curl = Curl()
    curl.(u)

Computes the per-element spectral (strong) curl of a vector field ``u``.

The curl of a vector field ``u`` is a vector field with contravariant components
```math
(\\nabla \\times u)^i = \\frac{1}{J} \\sum_{jk} \\epsilon^{ijk} \\frac{\\partial u_k}{\\partial \\xi^j}
```
where ``J`` is the Jacobian determinant, ``u_k`` is the ``k``th covariant
component of ``u``, and ``\\epsilon^{ijk}`` are the [Levi-Civita
symbols](https://en.wikipedia.org/wiki/Levi-Civita_symbol#Three_dimensions_2).
In other words
```math
\\begin{bmatrix}
  (\\nabla \\times u)^1 \\\\
  (\\nabla \\times u)^2 \\\\
  (\\nabla \\times u)^3
\\end{bmatrix}
=
\\frac{1}{J} \\begin{bmatrix}
  \\frac{\\partial u_3}{\\partial \\xi^2} - \\frac{\\partial u_2}{\\partial \\xi^3} \\\\
  \\frac{\\partial u_1}{\\partial \\xi^3} - \\frac{\\partial u_3}{\\partial \\xi^1} \\\\
  \\frac{\\partial u_2}{\\partial \\xi^1} - \\frac{\\partial u_1}{\\partial \\xi^2}
\\end{bmatrix}
```

In matrix form, this becomes
```math
\\epsilon^{ijk} J^{-1} D_j u_k
```
Note that unused dimensions will be dropped: e.g. the 2D curl of a
`Covariant12Vector`-field will return a `Contravariant3Vector`.

## References
- [Taylor2010](@cite), equation 17
"""
struct Curl{I} <: CurlSpectralElementOperator end
Curl() = Curl{()}()
Curl{()}(space) = Curl{operator_axes(space)}()

operator_return_eltype(::Curl{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(T -> Geometry.curl_result_type(Val(I), T), S)

function apply_slab(op::Curl{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    if RT <: Geometry.Contravariant3Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            v₁ = Geometry.covariant1(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for jj in 1:Nq
                D₂v₁ = D[jj, j] ⊠ v₁
                out[i, jj] =
                    out[i, jj] ⊞ Geometry.Contravariant3Vector(
                        ⊟(D₂v₁),
                        slab_local_geometry[i, jj],
                    )
            end
            v₂ = Geometry.covariant2(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for ii in 1:Nq
                D₁v₂ = D[ii, i] ⊠ v₂
                out[ii, j] =
                    out[ii, j] ⊞ Geometry.Contravariant3Vector(
                        D₁v₂,
                        slab_local_geometry[ii, j],
                    )
            end
        end
    elseif RT <: Geometry.Contravariant12Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            v₃ = Geometry.covariant3(
                get_node(slab_data, i, j),
                slab_local_geometry[i, j],
            )
            for ii in 1:Nq
                D₁v₃ = D[ii, i] ⊠ v₃
                out[ii, j] =
                    out[ii, j] ⊞
                    Geometry.Contravariant12Vector(zero(D₁v₃), ⊟(D₁v₃))
            end
            for jj in 1:Nq
                D₂v₃ = D[jj, j] ⊠ v₃
                out[i, jj] =
                    out[i, jj] ⊞
                    Geometry.Contravariant12Vector(D₂v₃, zero(D₂v₃))
            end
        end
    else
        error("invalid return type: $RT")
    end
    @inbounds for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.J)
    end
    return SMatrix(out)
end

function apply_slab(op::Curl{(1,)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MVector{Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    if RT <: Geometry.Contravariant2Vector
        @inbounds for i in 1:Nq
            v₃ = Geometry.covariant3(
                get_node(slab_data, i),
                slab_local_geometry[i],
            )
            for ii in 1:Nq
                D₁v₃ = D[ii, i] ⊠ v₃
                out[ii] = out[ii] ⊞ Geometry.Contravariant2Vector(⊟(D₁v₃))
            end
        end
    else
        error("invalid return type: $RT")
    end
    @inbounds for i in 1:Nq
        local_geometry = slab_local_geometry[i]
        out[i] = RecursiveApply.rdiv(out[i], local_geometry.J)
    end
    return SVector(out)
end

"""
    wcurl = WeakCurl()
    wcurl.(u)

Computes the "weak curl" on each element of a vector field `u`.

This is defined as the vector field ``\\theta \\in \\mathcal{V}_0`` such that
for all ``\\phi \\in \\mathcal{V}_0``
```math
\\int_\\Omega \\phi \\cdot \\theta \\, d \\Omega
=
\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
```
where ``\\mathcal{V}_0`` is the space of ``f``.

This arises from the contribution of the volume integral after by applying
integration by parts to the weak form expression of the curl
```math
\\int_\\Omega \\phi \\cdot (\\nabla \\times u) \\,d\\Omega
=
\\int_\\Omega (\\nabla \\times \\phi) \\cdot u \\,d \\Omega
- \\oint_{\\partial \\Omega} (\\phi \\times u) \\cdot n \\,d\\sigma
```

In matrix form, this becomes
```math
{\\phi_i}^\\top W J \\theta^i = (J^{-1} \\epsilon^{kji} D_j \\phi_i)^\\top W J u_k
```
which, by using the anti-symmetry of the Levi-Civita symbol, reduces to
```math
\\theta^i = - \\epsilon^{ijk} (WJ)^{-1} D_j^\\top W u_k
```
"""
struct WeakCurl{I} <: CurlSpectralElementOperator end
WeakCurl() = WeakCurl{()}()
WeakCurl{()}(space) = WeakCurl{operator_axes(space)}()

operator_return_eltype(::WeakCurl{I}, ::Type{S}) where {I, S} =
    RecursiveApply.rmaptype(T -> Geometry.curl_result_type(Val(I), T), S)

function apply_slab(op::WeakCurl{(1, 2)}, slab_space, _, slab_data)
    slab_local_geometry = Spaces.local_geometry_data(slab_space)
    FT = Spaces.undertype(slab_space)
    QS = Spaces.quadrature_style(slab_space)
    Nq = Quadratures.degrees_of_freedom(QS)
    D = Quadratures.differentiation_matrix(FT, QS)
    # allocate temp output
    RT = operator_return_eltype(op, eltype(slab_data))
    out = StaticArrays.MMatrix{Nq, Nq, RT}(undef)
    DataLayouts._mzero!(out, FT)
    if RT <: Geometry.Contravariant3Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            local_geometry = slab_local_geometry[i, j]
            W = local_geometry.WJ / local_geometry.J
            Wv₁ =
                W ⊠ Geometry.covariant1(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for jj in 1:Nq
                Dᵀ₂Wv₁ = D[j, jj] ⊠ Wv₁
                out[i, jj] =
                    out[i, jj] ⊞ Geometry.Contravariant3Vector(
                        Dᵀ₂Wv₁,
                        slab_local_geometry[i, jj],
                    )
            end
            Wv₂ =
                W ⊠ Geometry.covariant2(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for ii in 1:Nq
                Dᵀ₁Wv₂ = D[i, ii] ⊠ Wv₂
                out[ii, j] =
                    out[ii, j] ⊞ Geometry.Contravariant3Vector(
                        ⊟(Dᵀ₁Wv₂),
                        slab_local_geometry[ii, j],
                    )
            end
        end
    elseif RT <: Geometry.Contravariant12Vector
        @inbounds for j in 1:Nq, i in 1:Nq
            local_geometry = slab_local_geometry[i, j]
            W = local_geometry.WJ / local_geometry.J
            Wv₃ =
                W ⊠ Geometry.covariant3(
                    get_node(slab_data, i, j),
                    slab_local_geometry[i, j],
                )
            for ii in 1:Nq
                Dᵀ₁Wv₃ = D[i, ii] ⊠ Wv₃
                out[ii, j] =
                    out[ii, j] ⊞
                    Geometry.Contravariant12Vector(zero(Dᵀ₁Wv₃), Dᵀ₁Wv₃)
            end
            for jj in 1:Nq
                Dᵀ₂Wv₃ = D[j, jj] ⊠ Wv₃
                out[i, jj] =
                    out[i, jj] ⊞
                    Geometry.Contravariant12Vector(⊟(Dᵀ₂Wv₃), zero(Dᵀ₂Wv₃))
            end
        end
    else
        error("invalid return type: $RT")
    end
    for j in 1:Nq, i in 1:Nq
        local_geometry = slab_local_geometry[i, j]
        out[i, j] = RecursiveApply.rdiv(out[i, j], local_geometry.WJ)
    end
    return SMatrix(out)
end

# interplation / restriction
abstract type TensorOperator <: SpectralElementOperator end

input_space(op::TensorOperator, inspace) = inspace
return_space(op::TensorOperator, inspace) = op.space
operator_return_eltype(::TensorOperator, ::Type{S}) where {S} = S

"""
    i = Interpolate(space)
    i.(f)

Interpolates `f` to the `space`. If `space` has equal or higher polynomial
degree as the space of `f`, this is exact, otherwise it will be lossy.

In matrix form, it is the linear operator
```math
I = \\bigotimes_i I_i
```
where ``I_i`` is the barycentric interpolation matrix in the ``i``th dimension.

See also [`Restrict`](@ref).
"""
struct Interpolate{I, S} <: TensorOperator
    space::S
end
Interpolate(space) = Interpolate{operator_axes(space), typeof(space)}(space)

function apply_slab(
    op::Interpolate{(1,)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    Imat = Quadratures.interpolation_matrix(FT, QS_out, QS_in)
    S = eltype(slab_data)
    slab_data_out = MVector{Nq_out, S}(undef)
    @inbounds for i in 1:Nq_out
        # manually inlined rmatmul with slab_getnode
        r = Imat[i, 1] ⊠ get_node(slab_data, 1)
        for ii in 2:Nq_in
            r = RecursiveApply.rmuladd(Imat[i, ii], get_node(slab_data, ii), r)
        end
        slab_data_out[i] = r
    end
    return slab_data_out
end

function apply_slab(
    op::Interpolate{(1, 2)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    Imat = Quadratures.interpolation_matrix(FT, QS_out, QS_in)
    S = eltype(slab_data)
    # temporary storage
    temp = MArray{Tuple{Nq_out, Nq_in}, S, 2, Nq_out * Nq_in}(undef)
    slab_data_out = MArray{Tuple{Nq_out, Nq_out}, S, 2, Nq_out * Nq_out}(undef)
    @inbounds for j in 1:Nq_in, i in 1:Nq_out
        # manually inlined rmatmul1 with slab get_node
        # we do this to remove one allocated intermediate array
        r = Imat[i, 1] ⊠ get_node(slab_data, 1, j)
        for ii in 2:Nq_in
            r = RecursiveApply.rmuladd(
                Imat[i, ii],
                get_node(slab_data, ii, j),
                r,
            )
        end
        temp[i, j] = r
    end
    @inbounds for j in 1:Nq_out, i in 1:Nq_out
        slab_data_out[i, j] = RecursiveApply.rmatmul2(Imat, temp, i, j)
    end
    return SMatrix(slab_data_out)
end

"""
    r = Restrict(space)
    r.(f)

Computes the projection of a field `f` on ``\\mathcal{V}_0`` to a lower degree
polynomial space `space` (``\\mathcal{V}_0^*``). `space` must be on the same
topology as the space of `f`, but have a lower polynomial degree.

It is defined as the field ``\\theta \\in \\mathcal{V}_0^*`` such that for all ``\\phi \\in \\mathcal{V}_0^*``
```math
\\int_\\Omega \\phi \\theta \\,d\\Omega = \\int_\\Omega \\phi f \\,d\\Omega
```
In matrix form, this is
```math
\\phi^\\top W^* J^* \\theta = (I \\phi)^\\top WJ f
```
where ``W^*`` and ``J^*`` are the quadrature weights and Jacobian determinant of
``\\mathcal{V}_0^*``, and ``I`` is the interpolation operator (see [`Interpolate`](@ref))
from ``\\mathcal{V}_0^*`` to ``\\mathcal{V}_0``. This reduces to
```math
\\theta = (W^* J^*)^{-1} I^\\top WJ f
```
"""
struct Restrict{I, S} <: TensorOperator
    space::S
end
Restrict(space) = Restrict{operator_axes(space), typeof(space)}(space)

function apply_slab(
    op::Restrict{(1,)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    ImatT = Quadratures.interpolation_matrix(FT, QS_in, QS_out)' # transpose
    S = eltype(slab_data)
    slab_data_out = MVector{Nq_out, S}(undef)
    slab_local_geometry_in = Spaces.local_geometry_data(slab_space_in)
    slab_local_geometry_out = Spaces.local_geometry_data(slab_space_out)
    WJ_in = slab_local_geometry_in.WJ
    WJ_out = slab_local_geometry_out.WJ
    @inbounds for i in 1:Nq_out
        # manually inlined rmatmul with slab get_node
        r = ImatT[i, 1] ⊠ (WJ_in[1] ⊠ get_node(slab_data, 1))
        for ii in 2:Nq_in
            WJ_node = WJ_in[ii] ⊠ get_node(slab_data, ii)
            r = RecursiveApply.rmuladd(ImatT[i, ii], WJ_node, r)
        end
        slab_data_out[i] = RecursiveApply.rdiv(r, WJ_out[i])
    end
    return slab_data_out
end

function apply_slab(
    op::Restrict{(1, 2)},
    slab_space_out,
    slab_space_in,
    slab_data,
)
    FT = Spaces.undertype(slab_space_out)
    QS_in = Spaces.quadrature_style(slab_space_in)
    QS_out = Spaces.quadrature_style(slab_space_out)
    Nq_in = Quadratures.degrees_of_freedom(QS_in)
    Nq_out = Quadratures.degrees_of_freedom(QS_out)
    ImatT = Quadratures.interpolation_matrix(FT, QS_in, QS_out)' # transpose
    S = eltype(slab_data)
    # temporary storage
    temp = MArray{Tuple{Nq_out, Nq_in}, S, 2, Nq_out * Nq_in}(undef)
    slab_data_out = MArray{Tuple{Nq_out, Nq_out}, S, 2, Nq_out * Nq_out}(undef)
    slab_local_geometry_in = Spaces.local_geometry_data(slab_space_in)
    slab_local_geometry_out = Spaces.local_geometry_data(slab_space_out)
    WJ_in = slab_local_geometry_in.WJ
    WJ_out = slab_local_geometry_out.WJ
    @inbounds for j in 1:Nq_in, i in 1:Nq_out
        # manually inlined rmatmul1 with slab get_node
        r = ImatT[i, 1] ⊠ (WJ_in[1, j] ⊠ get_node(slab_data, 1, j))
        for ii in 2:Nq_in
            WJ_node = WJ_in[ii, j] ⊠ get_node(slab_data, ii, j)
            r = RecursiveApply.rmuladd(ImatT[i, ii], WJ_node, r)
        end
        temp[i, j] = r
    end
    @inbounds for j in 1:Nq_out, i in 1:Nq_out
        slab_data_out[i, j] = RecursiveApply.rdiv(
            RecursiveApply.rmatmul2(ImatT, temp, i, j),
            WJ_out[i, j],
        )
    end
    return SMatrix(slab_data_out)
end

"""
    tensor_product!(out, in, M)
    tensor_product!(inout, M)

Computes the tensor product `out = (M ⊗ M) * in` on each element.
"""
function tensor_product! end

function tensor_product!(
    out::DataLayouts.Data1DX{S, Ni_out},
    in::DataLayouts.Data1DX{S, Ni_in},
    M::SMatrix{Ni_out, Ni_in},
) where {S, Ni_out, Ni_in}
    (_, _, _, Nv_in, Nh_in) = size(in)
    (_, _, _, Nv_out, Nh_out) = size(out)
    # TODO: assumes the same number of levels (horizontal only)
    @assert Nv_in == Nv_out
    @assert Nh_in == Nh_out
    for h in 1:Nh_out, v in 1:Nv_out
        in_slab = slab(in, v, h)
        out_slab = slab(out, v, h)
        @inbounds for i in 1:Ni_out
            r = M[i, 1] ⊠ in_slab[1]
            for ii in 2:Ni_in
                r = RecursiveApply.rmuladd(M[i, ii], in_slab[ii], r)
            end
            out_slab[i] = r
        end
    end
    return out
end

function tensor_product!(
    out::DataLayouts.Data2D{S, Nij_out},
    in::DataLayouts.Data2D{S, Nij_in},
    M::SMatrix{Nij_out, Nij_in},
) where {S, Nij_out, Nij_in}

    Nh = length(in)
    @assert Nh == length(out)

    # temporary storage
    temp = MArray{Tuple{Nij_out, Nij_in}, S, 2, Nij_out * Nij_in}(undef)

    for h in 1:Nh
        in_slab = slab(in, h)
        out_slab = slab(out, h)
        for j in 1:Nij_in, i in 1:Nij_out
            temp[i, j] = RecursiveApply.rmatmul1(M, in_slab, i, j)
        end
        for j in 1:Nij_out, i in 1:Nij_out
            out_slab[i, j] = RecursiveApply.rmatmul2(M, temp, i, j)
        end
    end
    return out
end

function tensor_product!(
    out_slab::DataLayouts.DataSlab2D{S, Nij_out},
    in_slab::DataLayouts.DataSlab2D{S, Nij_in},
    M::SMatrix{Nij_out, Nij_in},
) where {S, Nij_out, Nij_in}
    # temporary storage
    temp = MArray{Tuple{Nij_out, Nij_in}, S, 2, Nij_out * Nij_in}(undef)
    for j in 1:Nij_in, i in 1:Nij_out
        temp[i, j] = RecursiveApply.rmatmul1(M, in_slab, i, j)
    end
    for j in 1:Nij_out, i in 1:Nij_out
        out_slab[i, j] = RecursiveApply.rmatmul2(M, temp, i, j)
    end
    return out_slab
end

function tensor_product!(
    inout::Data2D{S, Nij},
    M::SMatrix{Nij, Nij},
) where {S, Nij}
    tensor_product!(inout, inout, M)
end

"""
    matrix_interpolate(field, quadrature)

Computes the tensor product given a uniform quadrature `out = (M ⊗ M) * in` on each element.
Returns a 2D Matrix for plotting / visualizing 2D Fields.
"""
function matrix_interpolate end

function matrix_interpolate(
    field::Fields.SpectralElementField2D,
    Q_interp::Quadratures.Uniform{Nu},
) where {Nu}
    S = eltype(field)
    space = axes(field)
    topology = Spaces.topology(space)
    quadrature_style = Spaces.quadrature_style(space)
    mesh = topology.mesh
    n1, n2 = size(Meshes.elements(mesh))
    interp_data =
        DataLayouts.IH1JH2{S, Nu}(Matrix{S}(undef, (Nu * n1, Nu * n2)))
    M = Quadratures.interpolation_matrix(Float64, Q_interp, quadrature_style)
    Operators.tensor_product!(interp_data, Fields.field_values(field), M)
    return parent(interp_data)
end

function matrix_interpolate(
    field::Fields.ExtrudedFiniteDifferenceField,
    Q_interp::Union{Quadratures.Uniform{Nu}, Quadratures.ClosedUniform{Nu}},
) where {Nu}
    S = eltype(field)
    space = axes(field)
    quadrature_style = Spaces.quadrature_style(space)
    nl = Spaces.nlevels(space)
    n1 = Topologies.nlocalelems(Spaces.topology(space))
    interp_data = DataLayouts.IV1JH2{S, Nu}(Matrix{S}(undef, (nl, Nu * n1)))
    M = Quadratures.interpolation_matrix(Float64, Q_interp, quadrature_style)
    Operators.tensor_product!(interp_data, Fields.field_values(field), M)
    return parent(interp_data)
end

"""
    matrix_interpolate(field, Nu::Integer)

Computes the tensor product given a uniform quadrature degree of Nu on each element.
Returns a 2D Matrix for plotting / visualizing 2D Fields.
"""
matrix_interpolate(field::Field, Nu::Integer) =
    matrix_interpolate(field, Quadratures.Uniform{Nu}())


text: register_amplifier(SimpleAmplifier("30DV300", 0.0V, 10.0V, 15.0, -20.0V, 0.15A)) #used on OCPI2. Note that there are really two current channels (one per piezo stack) with a max amperage of 150mA each.  Therefore the true current calculation is more complicated than we acknowledge, but it should be about the same.

ext: jl


text:
using ADCME
using PyCall
using LinearAlgebra
using PyPlot
using Random
Random.seed!(233)

Sys.isapple() && matplotlib.use("macosx")

function advection_jl(m, n, h, v, u)
    # v: velocity, (m*n) * 2
    # u: temperature, (m*n)
    a = zeros(m*n)
    for j=1:n
        for i=1:m
            vx = v[(j-1)*m+i, 1]
            vy = v[(j-1)*m+i, 2]
            if vx >= 0
                if i==1
                    a[(j-1)*m+i] += vx * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vx * h * (u[(j-1)*m+i] - u[(j-1)*m+i-1])
                end
            else
                if i==m
                    a[(j-1)*m+i] += - vx * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vx * h * (u[(j-1)*m+i+1] - u[(j-1)*m+i])
                end
            end

            if vy >= 0
                if j==1
                    a[(j-1)*m+i] += vy * 2 * h * u[(j-1)*m+i]
                else
                    a[(j-1)*m+i] += vy * h * (u[(j-1)*m+i] - u[(j-2)*m+i])
                end
            else
                if j==n
                    a[(j-1)*m+i] += - vy * 2 * h * u[(j-1)*m+i]

                else
                    a[(j-1)*m+i] += vy * h * (u[j*m+i] - u[(j-1)*m+i])
                end
            end
        end
    end
    return a
end

function advection(v,u,m,n,h)
    advection_ = load_op_and_grad("./build/libAdvection","advection")
    v,u,m,n,h = convert_to_tensor(Any[v,u,m,n,h], [Float64,Float64,Int64,Int64,Float64])
    advection_(v,u,m,n,h)
end

# TODO: specify your input parameters
m = 10
n = 10
h = 0.1
v = rand(m*n,2)
u = rand(m*n)
jl_u = advection_jl(m, n, h, v, u)
u = advection(v,u,m,n,h)
sess = Session(); init(sess)
@show run(sess, u)-jl_u

# uncomment it for testing gradients
# error()

# TODO: change your test parameter to `m`
#       in the case of `multiple=true`, you also need to specify which component you are testings
# gradient check -- v
function scalar_function(x)
    return sum(advection(x,u,m,n,h)^2)
end

# TODO: change `m_` and `v_` to appropriate values
# m_ = constant(rand(m*n))
# v_ = rand(m*n)

m_ = constant(rand(m*n,2))
v_ = rand(m*n,2)
y_ = scalar_function(m_)
dy_ = gradients(y_, m_)
ms_ = Array{Any}(undef, 5)
ys_ = Array{Any}(undef, 5)
s_ = Array{Any}(undef, 5)
w_ = Array{Any}(undef, 5)
gs_ =  @. 1 / 10^(1:5)

for i = 1:5
    g_ = gs_[i]
    ms_[i] = m_ + g_*v_
    ys_[i] = scalar_function(ms_[i])
    s_[i] = ys_[i] - y_
    w_[i] = s_[i] - g_*sum(v_.*dy_)
end

sess = Session(); init(sess)
sval_ = run(sess, s_)
wval_ = run(sess, w_)
close("all")
loglog(gs_, abs.(sval_), "*-", label="finite difference")
loglog(gs_, abs.(wval_), "+-", label="automatic differentiation")
loglog(gs_, gs_.^2 * 0.5*abs(wval_[1])/gs_[1]^2, "--",label="\$\\mathcal{O}(\\gamma^2)\$")
loglog(gs_, gs_ * 0.5*abs(sval_[1])/gs_[1], "--",label="\$\\mathcal{O}(\\gamma)\$")

plt.gca().invert_xaxis()
legend()
xlabel("\$\\gamma\$")
ylabel("Error")
savefig("gradv.png")


text:
mutable struct CoefficientSum <: AbstractVector{Coefficient}
    #=
    Symbolic sum of symbolic coefficients.
    =#

    coeffs::Vector{Coefficient}

end

# interface with abstract array
size(cs::CoefficientSum) = size(cs.coeffs)

getindex(cs::CoefficientSum, i) = getindex(cs.coeffs, i)

setindex!(cs::CoefficientSum, c::Coefficient, i) = setindex!(cs.coeffs, c, i)

empty!(cs::CoefficientSum) = empty!(cs.coeffs)

push!(cs::CoefficientSum, cfs...) = push!(cs.coeffs, cfs...)

function show(io::IO, cs::CoefficientSum)
    #=
    Formatted printing.
    =#

    print(io, "(", cs[1])
    for c in cs[2:end]
        print(io, " + ", c)
    end
    print(io, ")")

end

# a zero element
zero(T::Type{CoefficientSum}) = CoefficientSum([zero(Coefficient)])

function zero!(cs::CoefficientSum)
    #=
    Set to canonical zero.
    =#

    empty!(cs)
    push!(cs, zero(Coefficient))

    nothing

end

function iszero!(cs::CoefficientSum)
    #=
    Is the coefficient sum a zero? Simplifies coefficient sum in the process.
    =#

    # first simplify in case coefficients cancel, then check if coefficient is zero
    simplify!(cs)
    all(iszero.(cs))

end

function simplify!(cs::CoefficientSum)
    #=
    Simplify the coefficeint sum.
    =#

    # put all coefficents into canonical form
    canonicalize!.(cs)

    # if coefficients can be added together, do so
    newcoeffs = _compresssum(cs)

    # find any zeros
    zs = iszero.(newcoeffs)

    if all(zs)
        # set to canonical zero
        zero!(cs)
    else
        # remove the zero coefficents
        deleteat!(newcoeffs, zs)

        # replace coefficients with simplified version
        cs.coeffs = newcoeffs
    end

    nothing

end

function *(cs1::CoefficientSum, cs2::CoefficientSum)
    #=
    Product of coefficient sums.
    =#

    newcoeffs = [c1*c2 for c1 in cs1, c2 in cs2 if true] # a 1d array
    newcs = CoefficientSum(newcoeffs)

    # simplify the result
    simplify!(newcs)

    newcs

end

function +(cs1::CoefficientSum, cs2::CoefficientSum)
    #=
    Add coefficient sums.
    =#

    newcs = CoefficientSum([cs1.coeffs; cs2.coeffs])

    # simplify the result
    simplify!(newcs)

    newcs

end

max_stars_repo_name: aormorningstar/SymbolicMajoranas

max_issues_repo_name: aormorningstar/SymbolicMajoranas

max_forks_repo_name: aormorningstar/SymbolicMajoranas


text:
"""
Calculates the negative log of the likelihood of our data given the paramters, using the Kalman filter. It uses the
predicted_observation_distributions from the kalman_filter function. The entries of this array in the second and
third columns represent the probability of the future observation of mRNA and Protein respectively, given our current knowledge.

Parameters
----------

protein_at_observations : numpy array.
    Observed protein. The dimension is m x n x 2, where m is the number of data sets, n is the
    number of observation time points. For each data set, the first column is the time,
    and the second column is the observed protein copy number at that time.

model_parameters : numpy array.
    An array containing the model parameters in the following order:
    repression_threshold, hill_coefficient, mRNA_degradation_rate,
    protein_degradation_rate, basal_transcription_rate, translation_rate,
    transcription_delay.

measurement_variance : float.
    The variance in our measurement. This is given by Sigma_e in Calderazzo et. al. (2018).

Returns
-------

log_likelihood : float.
    The log of the likelihood of the data.
"""
function calculate_log_likelihood_at_parameter_point(protein_at_observations,model_parameters,measurement_variance)
    if any(model_parameters .< 0)
        return -Inf
    end

    log_likelihood = 0
    _, _, _, _, predicted_observation_distributions, _, _ = kalman_filter(protein_at_observations,
                                                                          model_parameters,
                                                                          measurement_variance)
    observations = protein_at_observations[:,2]
    mean = predicted_observation_distributions[:,2]
    sd = sqrt.(predicted_observation_distributions[:,3])

    for observation_index in 1:length(observations)
        log_likelihood += logpdf(Normal(mean[observation_index],sd[observation_index]),observations[observation_index])
    end #for
    return log_likelihood
end # function

"""
Calculates the log of the likelihood, and the derivative of the negative log likelihood wrt each parameter, of our data given the
paramters, using the Kalman filter. It uses the predicted_observation_distributions, predicted_observation_mean_derivatives, and
predicted_observation_variance_derivatives from the kalman_filter function. It returns the log likelihood as in the
calculate_log_likelihood_at_parameter_point function, and also returns an array of the derivative wrt each parameter.

Parameters
----------

protein_at_observations : numpy array.
    Observed protein. The dimension is n x 2, where n is the number of observation time points.
    The first column is the time, and the second column is the observed protein copy number at
    that time.

model_parameters : numpy array.
    An array containing the moderowl parameters in the following order:
    repression_threshold, hill_coefficient, mRNA_degradation_rate,
    protein_degradation_rate, basal_transcription_rate, translation_rate,
    transcription_delay.

measurement_variance : float.
    The variance in our measurement. This is given by Sigma_e in Calderazzo et. al. (2018).

Returns
-------

log_likelihood : float.
    The log of the likelihood of the data.

log_likelihood_derivative : numpy array.
    The derivative of the log likelihood of the data, wrt each model parameter
"""
function calculate_log_likelihood_and_derivative_at_parameter_point(protein_at_observations,model_parameters,measurement_variance)
    number_of_parameters = size(model_parameters,1)

    _, _, _, _, predicted_observation_distributions, predicted_observation_mean_derivatives, predicted_observation_variance_derivatives = kalman_filter(protein_at_observations,
                                                                                model_parameters,
                                                                                measurement_variance)
    # calculate log likelihood as before
    log_likelihood = 0
    number_of_observations = size(protein_at_observations,1)
    observations = protein_at_observations[:,2]
    mean = predicted_observation_distributions[:,2]
    sd = sqrt.(predicted_observation_distributions[:,3])

    for observation_index in 1:length(observations)
        log_likelihood += logpdf(Normal(mean[observation_index],sd[observation_index]),observations[observation_index])
    end #for
    # now for the computation of the derivative of the negative log likelihood. An expression of this can be found
    # at equation (28) in Mbalawata, Särkkä, Haario (2013)
    observation_transform = [0.0 1.0]
    helper_inverse = 1.0./predicted_observation_distributions[:,3]

    log_likelihood_derivative = zeros(number_of_parameters)

    for parameter_index in 1:number_of_parameters
        for time_index in 1:number_of_observations
            log_likelihood_derivative[parameter_index] -= 0.5*((helper_inverse[time_index]*tr(observation_transform*predicted_observation_variance_derivatives[time_index,parameter_index,:,:]*
                                                               transpose(observation_transform))) -
                                                               (dot(helper_inverse[time_index]*transpose(observation_transform*predicted_observation_mean_derivatives[time_index,parameter_index,:,:]),
                                                               (observations[time_index] - mean[time_index]))) -
                                                               (helper_inverse[time_index]^2*(observations[time_index] - mean[time_index])^2*
                                                               dot(observation_transform,predicted_observation_variance_derivatives[time_index,parameter_index,:,:]*
                                                               transpose(observation_transform))) -
                                                               (dot(helper_inverse[time_index]*(observations[time_index] - mean[time_index]),
                                                               observation_transform*predicted_observation_mean_derivatives[time_index,parameter_index,:,:])))
         end #for
     end #for
    return log_likelihood, log_likelihood_derivative
end #function

"""
todo
"""
function log_likelihood_and_derivative_with_prior_and_transformation(protein_at_observations,proposed_position,mean_protein,measurement_variance)
    number_of_parameters = size(proposed_position,1)

    # reparameterise
    reparameterised_proposed_position = copy(proposed_position)
    reparameterised_proposed_position[[3,4,5,6]] .= exp.(reparameterised_proposed_position[[3,4,5,6]])

    # priors
    if any([pdf(Uniform(50,2*mean_protein),reparameterised_proposed_position[1]) == 0,
            pdf(Uniform(2,6),reparameterised_proposed_position[2]) == 0,
            pdf(Uniform(log(2)/150,log(2)/10),reparameterised_proposed_position[3]) == 0,
            pdf(Uniform(log(2)/150,log(2)/10),reparameterised_proposed_position[4]) == 0,
            pdf(Uniform(0.01,120),reparameterised_proposed_position[5]) == 0,
            pdf(Uniform(0.01,40),reparameterised_proposed_position[6]) == 0,
            pdf(Uniform(1,40),reparameterised_proposed_position[7]) == 0])
        return -Inf, zeros(number_of_parameters)
    end

    # quick and dirty clean up (make data start at '0')
    protein_at_observations[:,1] .-= protein_at_observations[1,1]
    log_likelihood, log_likelihood_derivative = calculate_log_likelihood_and_derivative_at_parameter_point(protein_at_observations,
                                                                                reparameterised_proposed_position,
                                                                                mean_protein,
                                                                                measurement_variance)
    for i in [3,4,5,6]
        log_likelihood_derivative[i] = reparameterised_proposed_position[i]*log_likelihood_derivative[i]
    end
    return log_likelihood, log_likelihood_derivative
end


text:
module MCMC
using Distributions

export gibbs, TuningParam, metropolis, metropolisAdaptive, logpdfLogX
export metLogAdaptive, metLogitAdaptive

include("gibbs.jl")
include("TuningParam.jl")
include("Fitness.jl")
include("Util.jl")

function metropolisBase(curr::Float64, log_prob::Function, stepSD::Float64)
  cand = curr + randn() * stepSD
  logU = log(rand())
  logP = log_prob(cand) - log_prob(curr)
  accept = logP > logU
  draw = accept ? cand : curr
  return (draw, accept)
end

function metropolis(curr::Float64, log_prob::Function, stepSD::Float64)
  return metropolisBase(curr, log_prob, stepSD)[1]
end

function metropolis(curr::Vector{Float64}, log_prob::Function,
                    stepSD::Matrix{Float64})
  cand = rand(MvNormal(curr, stepSD))
  logU = log(rand())
  p = log_prob(cand) - log_prob(curr)
  out = p > logU ? cand : curr
  return out
end

"""
Adaptive metropolis (within Gibbs). See section 3 of the paper below:
  http://probability.ca/jeff/ftpdir/adaptex.pdf

Another useful website:
  https://m-clark.github.io/docs/ld_mcmc/index_onepage.html
"""
function metropolisAdaptive(curr::Float64, log_prob::Function,
                            tuner::TuningParam;
                            update::Function=update_tuning_param_default)
  draw, accept = metropolisBase(curr, log_prob, tuner.value)
  update(tuner, accept)
  return draw
end

function logit(p::Float64; a::Float64=0.0, b::Float64=1.0)
  return log(p - a) - log(b - p)
end

function sigmoid(x::Float64; a::Float64=0.0, b::Float64=1.0)
  out = 0.0

  if a == 0 && b == 1
    out = 1 / (1 + exp(-x))
  else
    ex = exp(x)
    out = (b * ex + a) / (1 + ex)
  end

  return out
end

"""
Log absolute jacobian term to be added when taking the log of a postiviely
supported random variable (x), so as to produce a random variable with
support on the real line. When added to the target log density, the resulting
expression can be used in a metropolis step conveniently.
"""
function logabsjacobian_logx(real_x::Float64)::Float64
  return real_x
end

"""
Log absolute jacobian term to be added when taking the logit of a
lower-and-upper-bounded random variable (x), so as to produce a random variable
with support on the real line. When added to the target log density, the
resulting expression can be used in a metropolis step conveniently.

# Arguments
- `real_x::Float64`: the unconstrained value of the positively-supported
                     parameter x (i.e. log x)
- `a::Float64`: lower bound for unconstrained parameter (default: 0)
- `b::Float64`: upper bound for unconstrained parameter (default: 1)
"""
function logabsjacobian_logitx(real_x::Float64;
                               a::Float64=0.0, b::Float64=1.0)::Float64
  return logpdf(Logistic(), real_x) + log(b - a)
end

function metLogAdaptive(curr::Float64, log_prob::Function, tuner::TuningParam;
                        update::Function=update_tuning_param_default)

  function log_prob_plus_logabsjacobian(log_x::Float64)
    x = exp(log_x)
    return log_prob(x) + logabsjacobian_logx(log_x)
  end

  log_x = metropolisAdaptive(log(curr), log_prob_plus_logabsjacobian, tuner,
                             update=update)

  return exp(log_x)
end

function metLogitAdaptive(curr::Float64, log_prob::Function,
                          tuner::TuningParam; a::Float64=0.0, b::Float64=1.0,
                          update::Function=update_tuning_param_default)

  function log_prob_plus_logabsjacobian(logit_x::Float64)
    x = sigmoid(logit_x, a=a, b=b)
    return log_prob(x) + logabsjacobian_logitx(logit_x)
  end

  logit_x = metropolisAdaptive(logit(curr,a=a,b=b),
                               log_prob_plus_logabsjacobian,
                               tuner, update=update)

  return sigmoid(logit_x, a=a, b=b)
end

function normalize(x::Vector{T})::Vector{T} where T
  if isprobvec(x)
    return x
  elseif length(x) == 1
    @assert x[1] > 0
    return [1.0]
  else
    out = x / sum(x)
    out[1] = 1.0 - sum(x[2:end])
    return out
  end
end

"""
weighted sampling: takes (unnormalized) log probs and returns index
"""
function wsample_logprob(logProbs::Vector{T}) where {T <: Number}
  log_p_max = maximum(logProbs)
  p = exp.(logProbs .- log_p_max)
  return Distributions.wsample(p)
end

function logsumexp(logx::Vector{T}) where {T <: Number}
  mx = maximum(logx)
  return log(sum(exp.(logx .- mx))) + mx
end

function logsumexp(logx::T...) where {T <: Number}
  mx = maximum(logx)
  return log(sum(exp.(logx .- mx))) + mx
end

"""
log 1 minus. For example, `log1m(.3) == log(1 - .3) == log1p(-.3) == log(.7)`.
"""
log1m(x::T) where {T <: Number} = log1p(-x)
end # MCMC

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
using JuLIP
using Base.Test

println("===================================================")
println("          TEST SOLVE ")
println("===================================================")

println("-----------------------------------------------------------------")
println("Testing `minimise!` with equilibration with LJ calculator to lattice")
println("-----------------------------------------------------------------")
calc = lennardjones(r0=rnn(:Al))
at = bulk(:Al, cubic=true) * 10
X0 = positions(at) |> mat
at = rattle!(at, 0.02)
set_calculator!(at, calc)
set_constraint!(at, FixedCell(at))
minimise!(at, precond=:id, verbose=2)
X1 = positions(at) |> mat
X0 .-= X0[:, 1]
X1 .-= X1[:, 1]
F = X1 / X0
println("check that the optimiser really converged to a lattice")
@show vecnorm(F'*F - eye(3), Inf)
@show vecnorm(F*X0 - X1, Inf)
@test vecnorm(F*X0 - X1, Inf) < 1e-4

println("-------------------------------------------------")
println("same test but large and with Exp preconditioner")
println("-------------------------------------------------")

at = bulk(:Al, cubic=true) * (20,20,2)
at = rattle!(at, 0.02)
set_calculator!(at, calc)
set_constraint!(at, FixedCell(at))
minimise!(at, precond = :exp, method = :lbfgs,
          robust_energy_difference = true, verbose=2)

println("-------------------------------------------------")
println("   Variable Cell Test")
println("-------------------------------------------------")
calc = lennardjones(r0=rnn(:Al))
at = set_pbc!(bulk(:Al, cubic=true), true)
set_calculator!(at, calc)
set_constraint!(at, VariableCell(at))
minimise!(at, verbose = 2)

println("-------------------------------------------------")
println(" FF preconditioner for StillingerWeber ")
println("-------------------------------------------------")

at = bulk(:Si, cubic=true) * (10,10,2)
at = set_pbc!(at, true)
at = rattle!(at, 0.02)
set_calculator!(at, StillingerWeber())
set_constraint!(at, FixedCell(at))
P = FF(at, StillingerWeber())
minimise!(at, precond = P, method = :lbfgs, robust_energy_difference = true, verbose=2)

println("-------------------------------------------------")
println(" FF preconditioner for EAM ")
println("-------------------------------------------------")

at = bulk(:W, cubic=true) * (10,10,2)
at = set_pbc!(at, true)
at = rattle!(at, 0.02)
X0 = positions(at)
set_calculator!(at, eam_W)
set_constraint!(at, FixedCell(at))
P = FF(at, eam_W)
minimise!(at, precond = P, method = :lbfgs, robust_energy_difference = true, verbose=2)

println("for comparison now with Exp")
set_positions!(at, X0)
minimise!(at, precond = :exp, method = :lbfgs, robust_energy_difference = true, verbose=2)


text:
"""
    module Trees

Module defining some commonly used tree structures, as well as an interface for
talking about trees inside the `WaveProp` organization.
"""
module Trees

import Base: ==, length, in, iterate, getindex, lastindex, isempty, eltype, keys
import Base: union, setdiff, intersect, issubset

using StaticArrays
using LinearAlgebra
using AbstractTrees
using Printf
using Statistics: median

using WavePropBase
using WavePropBase.Utils
using WavePropBase.Geometry
using WavePropBase.Interpolation

WavePropBase.@import_interface

include("abstracttree.jl")
include("clustertree.jl")
include("splitter.jl")

export
    # abstract types
    AbstractTree,
    # types
    ClusterTree,
    GeometricSplitter,
    GeometricMinimalSplitter,
    PrincipalComponentSplitter,
    DyadicSplitter,
    CardinalitySplitter,
    # re-exported from AbstractTrees
    Leaves,
    PreOrderDFS,
    PostOrderDFS,
    print_tree,
    # functions
    children,
    parent,
    isroot,
    isleaf,
    depth,
    container,
    loc2glob,
    root_elements,
    index_range,
    filter_tree,
    filter_tree!

end # module

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
@enum Trit False Maybe True
const trits = (False, Maybe, True)

Base.:!(a::Trit) = a == False ? True : a == Maybe ? Maybe : False
∧(a::Trit, b::Trit) = a == b == True ? True : (a, b) ∋ False ? False : Maybe
∨(a::Trit, b::Trit) = a == b == False ? False : (a, b) ∋ True ? True : Maybe
⊃(a::Trit, b::Trit) = a == False || b == True ? True : (a, b) ∋ Maybe ? Maybe : False
≡(a::Trit, b::Trit) = (a, b) ∋ Maybe ? Maybe : a == b ? True : False

println("Not (!):")
println(join(@sprintf("%10s%s is %5s", "!", t, !t) for t in trits))
println("And (∧):")
for a in trits
    println(join(@sprintf("%10s ∧ %5s is %5s", a, b, a ∧ b) for b in trits))
end
println("Or (∨):")
for a in trits
    println(join(@sprintf("%10s ∨ %5s is %5s", a, b, a ∨ b) for b in trits))
end
println("If Then (⊃):")
for a in trits
    println(join(@sprintf("%10s ⊃ %5s is %5s", a, b, a ⊃ b) for b in trits))
end
println("Equivalent (≡):")
for a in trits
    println(join(@sprintf("%10s ≡ %5s is %5s", a, b, a ≡ b) for b in trits))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function compute_gs_hs(::Type{CoefficientsFactory}, arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement, arg3::RealFieldElement, arg4::jint, arg5::Field)
    return jcall(CoefficientsFactory, "computeGsHs", Vector{Vector{RealFieldElement}}, (RealFieldElement, RealFieldElement, RealFieldElement, RealFieldElement, jint, Field), arg0, arg1, arg2, arg3, arg4, arg5)
end

function compute_gs_hs(::Type{CoefficientsFactory}, arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble, arg4::jint)
    return jcall(CoefficientsFactory, "computeGsHs", Vector{Vector{jdouble}}, (jdouble, jdouble, jdouble, jdouble, jint), arg0, arg1, arg2, arg3, arg4)
end

function compute_qns(::Type{CoefficientsFactory}, arg0::RealFieldElement, arg1::jint, arg2::jint)
    return jcall(CoefficientsFactory, "computeQns", Vector{Vector{RealFieldElement}}, (RealFieldElement, jint, jint), arg0, arg1, arg2)
end

function compute_qns(::Type{CoefficientsFactory}, arg0::jdouble, arg1::jint, arg2::jint)
    return jcall(CoefficientsFactory, "computeQns", Vector{Vector{jdouble}}, (jdouble, jint, jint), arg0, arg1, arg2)
end

function compute_vns(::Type{CoefficientsFactory}, arg0::jint)
    return jcall(CoefficientsFactory, "computeVns", TreeMap, (jint,), arg0)
end

function get_vmns(::Type{CoefficientsFactory}, arg0::jint, arg1::jint, arg2::jint)
    return jcall(CoefficientsFactory, "getVmns", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end


text:
"""
    sat_tracks(; tiles::Bool=false, position::Bool=false, kwargs...)

Compute satellite tracks using the TLE, or Two Line Elements set, a data format that contains
information about the orbit at a specific epoch of an Earth-orbiting object. It can also
calculate polygons arround the scene extents of AQUA and TERRA satellites as well as create
the scene names, which provides a mean to direct download that data.
  - `start`: A DateTime object or a string convertable to a DateTime with `DateTime(start)`
     specifying the start of the orbit calculation. If omited, current time in UTC will be used.
  - `duration`: Length of time for which the orbit is calculated. Accepts duration in days, hours,
     minutes or seconds. The default is minutes (100 minutes). To use other units use a string with
     the value appended with 'D', 'h', 'm' or 's'. _e.g._ `duration="55m"` to compute orbit 55 minutes from `start`
  - `step` or `inc` or `dt`: The time interval at which to compute locations along the orbit. The default
     unit here is seconds (30 sec) but minutes can be used as well by appending 'm'. _e.g._ `step="1m"`
  - `stop`: As alternative to `duration` provide the end date for the orbit. Same conditions as `start`
  - `position`: Computes only first location af the `start` time. Boolean, use `position=true`
  - `tle` or `TLE`: a file name with the TLE data for a specific satellite and period. It can also be a
    two elements string vector with the first and second lines of the TLE file.
  - `tiles`: Compute the scene limits and file names for some satellites. Currently AQUA only.
  - `sat`, `SAT` or `satellite`: Name of the satellite to use; choose from (string or symbols)
     :TERRA, :AQUA. Use only with the `tiles` option.

### Returns
A GMTdataset with the orbit or the scene polygons

## Example:
Compute ~one orbit of the AQUA satellite starting at current local time. Note, this
will be accurate for the month of September 2021. For other dates it needs an updated TLE.

    tle1 = "1 27424U 02022A   21245.83760660  .00000135  00000-0  39999-4 0  9997";
    tle2 = "2 27424  98.2123 186.0654 0002229  67.6025 313.3829 14.57107527 28342";
    orb = sat_tracks(tle=[tle1; tle2], duration=100);

and the orbit track can be visualized with

    imshow(orb,  proj=:Robinson, region=:global, coast=true)
"""
function sat_tracks(; geocentric::Bool=false, tiles::Bool=false, position::Bool=false, kwargs...)
	# ...
	(position && tiles) && error("Cannot require tiles and a single position. Makes no sense.")
	d = KW(kwargs)

	start = ((val = find_in_dict(d, [:start])[1]) === nothing) ? now(Dates.UTC) : getitDTime(val)

	(tiles) && (sat_name = get_sat_name(d))
	(tiles) && (start = round(start, Dates.Minute(5)))	# This is for MODIS only

	if ((val = find_in_dict(d, [:duration])[1]) !== nothing)
		if (isa(val, String))		# Accept duration in D(ays), h(ours), m(inutes) or s(econds)
			if     (endswith(val,"D"))  dur = Day(parse(Int, val[1:end-1]))
			elseif (endswith(val,"h"))  dur = Hour(parse(Int, val[1:end-1]))
			elseif (endswith(val,"m"))  dur = Minute(parse(Int, val[1:end-1]))
			elseif (endswith(val,"s"))  dur = Second(parse(Int, val[1:end-1]))
			else    error("Only 'D', 'h', 'm' or 's' are accepted in duration")
			end
		elseif (isa(val, Real))		# Assume duration was given in hours
			dur = Hour(trunc(Int, val))
		end
		stop = start + dur
	else
		if ((val = find_in_dict(d, [:stop])[1]) !== nothing)
			stop = getitDTime(val)
		else
			stop = start + Minute(100)		# Default is ~Terra period
		end
	end
	if ((val = find_in_dict(d, [:step :inc :dt])[1]) !== nothing)	# Steps are in seconds
		if (isa(val, String))
			if     (endswith(val,"m"))  dt = parse(Int, val[1:end-1]) * 60
			elseif (endswith(val,"s"))  dt = parse(Int, val[1:end-1])
			else    error("Only 's' or 'm' are accepted in increment")
			end
		else
			dt = trunc(Int, val)
		end
	else
		dt = 30
	end

	if ((val_tle = find_in_dict(d, [:tle_obj])[1]) !== nothing)  tle = val_tle	# Some other fun already got it
	else                                                         tle = loadTLE(d)
	end
	orbp = SatelliteToolbox.init_orbit_propagator(Val(:sgp4), tle[1])

	startmfe = (datetime2julian(DateTime(start)) - tle[1].epoch) * 24 * 3600
	stopmfe  = (datetime2julian(DateTime(stop))  - tle[1].epoch) * 24 * 3600
	(tiles) && (dt = 60)			# Arbitrary choice that works well for MODIS but may need revision for others
	t = startmfe:dt:stopmfe

	(position) && (t = [t[1]])		# Single position. Doing it here wastes work above but code is way cleaner

	out = Matrix{Float64}(undef, length(t), 4)
	r, = SatelliteToolbox.propagate!(orbp, t)

	for n = 1:length(t)
		jd = tle[1].epoch + t[n] / (24 * 3600)
		tt = SatelliteToolbox.r_eci_to_ecef(SatelliteToolbox.TEME(), SatelliteToolbox.PEF(), jd) * r[n]
		out[n,1], out[n,2], out[n,3], out[n, 4] = tt[1], tt[2], tt[3], jd
	end

	if (tiles)
		out = mapproject(out, E=true, I=true)
		return make_sat_tiles(out[1].data, SCENE_HALFW[sat_name], sat_name)
	end

	return (geocentric) ? out : mapproject(out, E=true, I=true)[1]
end

# --------------------------------------------------------------------------------------------
function getitDTime(val)
	if     isa(val, String) || isa(val, Tuple)  ret::DateTime = DateTime(val)
	elseif (isa(val, DateTime))  ret = val
	else   error("Bad input type $(typeof(val)). Must be a DateTime, a String or a Tuple(Int)")
	end
	return ret
end

# --------------------------------------------------------------------------------------------
function loadTLE(d::Dict)
	# Load a TLE or use a default one. In a function because it's used by two functions
	if ((val = find_in_dict(d, [:tle :TLE])[1]) !== nothing)
		if (isa(val, String))  tle = SatelliteToolbox.read_tle(val)
		elseif (isa(val, Vector{String}) && length(val) == 2)
			tle = SatelliteToolbox.read_tle_from_string(val[1], val[2])
		else
			error("BAD input TLE data")
		end
	else
		#tle = SatelliteToolbox.read_tle("C:\\v\\Landsat8.tle")
		tle = SatelliteToolbox.read_tle("C:\\v\\AQUA.tle")
	end
	tle
end

# --------------------------------------------------------------------------------------------
function get_sat_name(d::Dict)::String
	# Get the satellite name from kwargs (encoded in 'd'). Used by at least two functions
	((val = RemoteS.find_in_dict(d, [:sat :SAT :satellite])[1]) === nothing) &&
		error("Must provide the satellite name. Pick one of :TERRA, :AQUA")
	!isa(val, String) && (val = string(val))
	sat = uppercase(val)
	(sat != "TERRA" && sat != "AQUA" && sat != "LANDSAT8") &&
		error("Unknown satellite name $sat. Must be one of :TERRA, :AQUA")
	return sat
end

# --------------------------------------------------------------------------------------------
"""
sat_scenes(track, sat_name::String)

Compute polygons delimiting AQUA and TERRA scenes.

- `trac`: Is an orbit computed with `sat_tracks` at steps of 1 minute (crucial)
- `sat_name`: The satellite name. At this time only AQUA and TERRA are allowed.

Returns a GMTdataset vector with the polygons and the scene names in the dataset `header` field.

# Example
Imagine that `orb` was obtained with

orb = sat_tracks(tle=[tle1; tle2], start=DateTime("2021-09-02T13:30:00"),
	stop=DateTime("2021-09-02T13:40:00"), step="1m");

The scenes limits (two) are computed with:

```
Dscenes = sat_scenes(orb, "AQUA");
```

"""
function sat_scenes(track, sat_name::String)
	# Compute polygons with the scene limits. 'halfwidth' is half the scene width.
	halfwidth = get(RemoteS.SCENE_HALFW, sat_name, 0)
	(halfwidth == 0) && error("Currently only \"AQUA\" and \"TERRA\" satellite names are available")
	_, azim, = invgeod(track[1:end-1, 1:2], track[2:end, 1:2])	# distances and azimuths along the tracks
	append!(azim, azim[end])		# To make it same size of tracks
	D = Vector{GMTdataset}(undef, trunc(Int, (length(azim)-1)/5))
	n = 0
	for k = 1:5:length(azim)-5
		ll14 = geod(track[k,1:2],   [azim[k]+90, azim[k]-90], halfwidth)[1]
		ll23 = geod(track[k+5,1:2], [azim[k+5]+90, azim[k+5]-90], halfwidth)[1]
		sc = get_MODIS_scene_name(track[k,4], sat_name)
		D[n+=1] = GMTdataset([ll14[1:1,:]; ll23[1:1,:]; ll23[2:2,:]; ll14[2:2,:]; ll14[1:1,:]], String[], sc, String[], "", "", GMT.Gdal.wkbPolygon)
	end
	D[1].proj4 = GMT.prj4WGS84
	D
end

# --------------------------------------------------------------------------------------------
function get_MODIS_scene_name(jd::Float64, prefix::String, sst::Bool=true)
	# Build the scene name after it's acquisition time coded in the JD argin.
	DT = julian2datetime(jd)
	if (sst)
		p = (prefix[1] == 'A') ? "AQUA_MODIS." : "TERRA_MODIS."
		@sprintf("%s%.4d%.2d%.2dT%.2d%.2d01.L2.SST.NRT.nc", p, year(DT), month(DT), day(DT), hour(DT), minute(DT))
	else
		@sprintf("%s%.4d%.3d%.2d%.2d00.L2_LAC_OC.nc", prefix[1], year(DT), dayofyear(DT), hour(DT), minute(DT))
	end
end

# --------------------------------------------------------------------------------------------
"""
    clip_orbits(track, BoundingBox::Vector{<:Real})

Clips the orbits that are contained inside a rectangular geographical region

- `track`: A GMTdataset or a Mx2 matrix with the orbits [lon, lat] position. This is normally calculated
   with the `sat_tracks` function.
- `BoundingBox`: A vector the region limits made up with [lon_min, lon_max, lat_min, lat_max]

Returns a GMTdataset vector with the chunks of `tracks` that cross inside the `BoundingBox` region.

#Example
Suppose `orb` holds orbits computed with sat_tracks() during 2 days, clip them inside the 20W-10E, 30N-45N window

```
D = clip_orbits(orb, [-20, 10, 30, 50]);
```
"""
function clip_orbits(track, bb::Vector{<:Real})
	# Select the parts of the track that are contained inside the BB.
	(length(bb) != 4) && error("The BoundingBox vector must have 4 elements")
	segments = GMT.gmtselect(track, region=bb, f=:g)
	azim = invgeod(segments[1].data[1:end-1, 1:2], segments[1].data[2:end, 1:2])[2]
	inds = findall(abs.(diff(azim)) .> 10)			# Detect line breaks

	begin_seg = [1; inds[2:2:length(inds)].+1]		# Not easy to explain why those indices give seg boundaries
	end_seg   = [inds[2:2:length(inds)]; length(azim)+1]
	D = Vector{GMTdataset}(undef, length(begin_seg))
	for k = 1:length(begin_seg)
		D[k] = GMTdataset(segments[1].data[begin_seg[k]:end_seg[k], :], String[], "", String[], "", "", GMT.Gdal.wkbLineString)
	end
	D[1].proj4 = GMT.prj4WGS84
	D
end

# --------------------------------------------------------------------------------------------
"""
    findscenes(lon::Real, lat::Real; kwargs...)

Find the names of the scenes that cover the location point `lon`, `lat` in the period determined
by the dates and satellite set via kwargs.
  - `day`: Search only on the day time part of the orbits.
  - `night`: Search only on the night time part of the orbits.
  - `oc`: For the AQUA or TERRA satellites pick only the chlorophyl content scenes.
  - `sst`: For the AQUA or TERRA satellites pick only the Sae Surface Temperature content scenes.
  - `sat`, `SAT` or `satellite`: Name of the satellite to use; choose from (string or symbols)
     :TERRA, :AQUA
  - `start`: A DateTime object or a string convertable to a DateTime with `DateTime(start)`
     specifying the start of the looking period. If omited, current time in UTC will be used.
  - `duration`: Length of time for which the scenes are searched. The duration is expected in days
     and can be a negative number, meaning we'll look that span days from `start`.
  - `stop`: As alternative to `duration` provide the end date for the serch. Same conditions as `start`
  - `tle` or `TLE`: a file name with the TLE data for a specific satellite and period. It can also be a
    two elements string vector with the first and second lines of the TLE file.

### Returns
A string vector with the scene names

## Example:
Find the AQUA scenes with chlorophyl-a (oceancolor) that cover the point (-8, 36) in the two days
before "2021-09-07T17:00:00"
Note, this will be accurate for the month of September 2021. For other dates it needs an updated TLE.

    tle1 = "1 27424U 02022A   21245.83760660  .00000135  00000-0  39999-4 0  9997";
    tle2 = "2 27424  98.2123 186.0654 0002229  67.6025 313.3829 14.57107527 28342";
    findscenes(-8,36, start="2021-09-07T17:00:00", sat=:aqua, day=true, duration=-2, oc=1, tle=[tle1, tle2])

	2-element Vector{String}:
	"A2021251125500.L2_LAC_OC.nc"
	"A2021252134000.L2_LAC_OC.nc"
"""
function findscenes(lon::Real, lat::Real; kwargs...)
	# ...
	d = KW(kwargs)
	sat = get_sat_name(d)

	day   = (haskey(d, :day))   ? true : false
	night = (haskey(d, :night)) ? true : false
	sst   = (haskey(d, :sst))   ? true : false
	oc    = (haskey(d, :oc))    ? true : false

	start = ((val = find_in_dict(d, [:start])[1]) === nothing) ? now(UTC) - Day(2) : getitDTime(val)
	if ((val = find_in_dict(d, [:duration])[1]) !== nothing)
		period = round(Int, val)
		(period == 0) && (period = -1; @warn("Duration too short. Using -1 days"))
		if (period < 0)  stop, start  = start+Day(2), start+Day(2) - Day(-period)
		else             stop = start + Day(period)
		end
	else
		stop  = ((val = find_in_dict(d, [:stop])[1])  === nothing) ? start + Day(2) : getitDTime(val)
	end
	(start >= stop) && error("Start date ($start) is greater or equal to stop ($stop)")
	(haskey(d, :Vd)) && println("start = ",start, " stop = ", stop)

	tle = loadTLE(d)
	tracks = sat_tracks(start=start, stop=stop, tle_obj=tle)
	BB = [lon-10, lon+10, lat-10, lat+10]
	D = clip_orbits(tracks, BB)
	(day || night) && (D = day_night_orbits(D, day=day, night=night))
	isempty(D) && (@warn("No tracks left after this choice."); return nothing)

	scenes = Vector{String}(undef, 0)
	dists = mapproject(D, G=(lon,lat))
	for k = 1:length(dists)
		d, ind = findmin(view(dists[k].data, :,5))
		if (d <= SCENE_HALFW[sat])
			jd = datetime2julian(floor(julian2datetime(dists[k][ind,4]), Dates.Minute(5)))
			sc = get_MODIS_scene_name(jd, sat, sst == true)
			append!(scenes, [sc])
		end
	end
	scenes
end

# --------------------------------------------------------------------------------------------
function day_night_orbits(D; day::Bool=false, night::Bool=false)
	# Pick only the day or night orbits in D. D is normally the output of the clip_orbits() fun
	(!day && !night) && return D		# No selection requested

	pass = zeros(Bool, length(D))
	for k = 1:length(D)
		lon, lat = D[k][1,1:2]
		jd = D[k][1,4]
		raise, set = solar(@sprintf("-I%.4f/%.4f+d%s -C", lon, lat, string(julian2datetime(jd))))[1][5:6]
		# Subtract 0.5 because julian2datetime(0.0) = -4713-11-24T12:00:00. That is JDs are shifted of 0.5
		t1, t2 = jd - 0.5, D[k][end,4] - 0.5
		dday_first, dday_last = t1 - trunc(t1), t2 - trunc(t2)
		if (day)
			if ((raise < dday_first < set) && (raise < dday_last < set))
				pass[k] = true
			elseif (dday_first < raise && dday_last > raise) || (dday_first > raise && dday_last > set)
				@warn("Track $k contains day and night parts. Accepting it as day but be warned of this.")
				pass[k] = true
			end
		else		# night
			if !((raise < dday_first < set) && (raise < dday_last < set))
				pass[k] = true
			elseif (dday_first < raise && dday_last > set) || (dday_first > raise && dday_last > set)
				@warn("Track $k contains day and night parts. Accepting it as night but be warned of this.")
				pass[k] = true
			end
		end
	end
	D[pass]
end

max_stars_repo_name: GenericMappingTools/RemoteS

max_issues_repo_name: GenericMappingTools/RemoteS

max_forks_repo_name: GenericMappingTools/RemoteS


text:
using OpenCV
const cv = OpenCV

img = rand(UInt8, 3, 500, 500)
filter = rand(Float32, 1, 5, 5)/25

out = OpenCV.filter2D(img, Int32(-1), filter)

cv.namedWindow("orig")
cv.namedWindow("out")

cv.imshow("orig", img)
cv.imshow("out", out)

cv.waitKey(Int32(0))


text:
module DocServer

import Documenter
import Documenter: Selectors, Utilities, Documents, Expanders
import Documenter.Expanders: ExpanderPipeline
import Documenter.Builder: DocumentPipeline

using Mux
using WebIO

"""
"""
abstract type LiveBlocks <: ExpanderPipeline end
Selectors.order(::Type{LiveBlocks})     = 8.5
Selectors.matcher(::Type{LiveBlocks},  node, page, doc) = Expanders.iscode(node, r"^@live")

const live_blocks = Dict{String, Any}()

function Selectors.runner(::Type{LiveBlocks}, x, page, doc)
    matched = match(r"^@live[ ]?(.*)$", x.language)
    matched === nothing && error("invalid '@live' syntax: $(x.language)")
    # The sandboxed module -- either a new one or a cached one from this page.

    key = replace(page.source, "/", ":")
    name = matched[1]
    sym  = isempty(name) ? gensym("live-") : Symbol("live-", name)
    mod  = get!(page.globals.meta, sym, Module(sym))::Module
    content = []
    input   = Expanders.droplines(x.code)

    if !haskey(live_blocks, key)
        live_blocks[key] = []
    end
    push!(live_blocks[key], (page, doc, mod, input))

    id = "live-$(length(live_blocks[key]))"
    scr = """

live output $id

    """
    output = Documents.RawHTML(scr)

    # Only add content when there's actually something to add.
    isempty(input)  || push!(content, Markdown.Code("julia", input))
    isempty(output.code) || push!(content, output)
    # ... and finally map the original code block to the newly generated ones.
    page.mapping[x] = Markdown.MD(content)
end

"""
Adds the document to the Mux app
"""
abstract type ServeDocument <: DocumentPipeline end
Selectors.order(::Type{ServeDocument}) = typemax(Float64) - π

function htmlrender(m::Union{MIME"text/html",
                            MIME"text/plain",
                            MIME"text/svg"}, x)
    stringmime(m, x)
end

function htmlrender(m::Union{MIME"image/png", MIME"image/jpeg"})
    ""
end

const bestmimes = MIME.(["text/html", "text/svg", "image/png", "image/jpeg"])

function htmlrender(x)
    for m in bestmimes
        if mimewritable(m, x)
            return htmlrender(m, x)
       end
    end
    return string(x)
end

function server(req)
    key = req[:params][:key]
    @show keys(live_blocks) |> collect
    if haskey(live_blocks, key)
        blocks = []
        for (page, doc, mod, input) in live_blocks[key]
            result = nothing
            for (ex, str) in Utilities.parseblock(input, doc, page)
                (value, success, backtrace, text) = Utilities.withoutput() do
                    cd(dirname(page.build)) do
                        eval(mod, ex)
                    end
                end
                result = value
                if !success
                    Utilities.warn(page.source, "failed to run code block.\n\n$(value)")
                    return
                end
            end
            push!(blocks, htmlrender(result))
        end

        return Dict(
                    :headers => Dict("Content-Type" => "text/json"),
                    :body => JSON.json(blocks)
                   )
    end
    return Dict(
                :headers => Dict("Content-Type"=> "text/json"),
                :body => ""
               )
end

include(Pkg.dir("WebIO", "src", "providers", "mux_setup.jl"))
function Selectors.runner(::Type{ServeDocument}, doc::Documents.Document)
    Utilities.log(doc, "serving files...")
    dir = joinpath(doc.user.root, doc.user.build)
    fileserver = Mux.stack(Mux.files(dir, false),
                           branch(req -> Mux.validpath(dir, joinpath(req[:path]..., "index.html")),
                                  req -> Mux.fresp(joinpath(dir, req[:path]..., "index.html"))))
    cd(dir) do
        ware = route("", fileserver, Mux.notfound())
        webio_serve(Mux.stack(page("/liveblocks/:key", server), ware))
    end
end

end # module


text:
# Integrate solution by repeatedly calling the rhs! method on the solver solution.
# function timestep_XYZ!(solver::AbstractSolver, t, dt)

"""
    timestep_2N!(solver, t, dt, a, b, c)

Perform one timestep using an explicit Runge-Kutta method of the low-storage
class 2N of Williamson.
"""
@inline function timestep_2N!(solver, t, dt, a, b, c)
  for stage in eachindex(c)
    t_stage = t + dt * c[stage]
    @timeit timer() "rhs" rhs!(solver, t_stage)

    a_stage    = a[stage]
    b_stage_dt = b[stage] * dt
    @timeit timer() "Runge-Kutta step" begin
      Threads.@threads for i in eachindex(solver.elements.u)
        solver.elements.u_tmp2[i] = solver.elements.u_t[i] - solver.elements.u_tmp2[i] * a_stage
        solver.elements.u[i] += solver.elements.u_tmp2[i] * b_stage_dt
      end
    end
  end
  apply_positivity_preserving_limiter!(solver)
end

"""
    timestep_carpenter_kennedy_erk54_2N!(solver::AbstractSolver, t, dt)

Carpenter, Kennedy (1994) Fourth order 2N storage RK schemes, Solution 3
"""
function timestep_carpenter_kennedy_erk54_2N!(solver::AbstractSolver, t, dt)
  a = @SVector [0.0, 567301805773.0 / 1357537059087.0,2404267990393.0 / 2016746695238.0,
       3550918686646.0 / 2091501179385.0, 1275806237668.0 / 842570457699.0]
  b = @SVector [1432997174477.0 / 9575080441755.0, 5161836677717.0 / 13612068292357.0,
       1720146321549.0 / 2090206949498.0, 3134564353537.0 / 4481467310338.0,
       2277821191437.0 / 14882151754819.0]
  c = @SVector [0.0, 1432997174477.0 / 9575080441755.0, 2526269341429.0 / 6820363962896.0,
       2006345519317.0 / 3224310063776.0, 2802321613138.0 / 2924317926251.0]

  timestep_2N!(solver, t, dt, a, b, c)
end

"""
    timestep_carpenter_kennedy_erk43_2N!(solver::AbstractSolver, t, dt)

Carpenter, Kennedy (1994) Third order 2N storage RK schemes with error control
"""
function timestep_carpenter_kennedy_erk43_2N!(solver::AbstractSolver, t, dt)
  a = @SVector [0, 756391 / 934407, 36441873 / 15625000, 1953125 / 1085297]
  b = @SVector [8 / 141, 6627 / 2000, 609375 / 1085297, 198961 / 526383]
  c = @SVector [0, 8 / 141, 86 / 125, 1]

  timestep_2N!(solver, t, dt, a, b, c)
end

"""
    timestep_3Sstar!(solver, t, dt, gamma1, gamma2, gamma3, beta, delta, c)

Perform one timestep using an explicit Runge-Kutta method of the low-storage
class 3Sstar of Ketcheson.
"""
@inline function timestep_3Sstar!(solver, t, dt, gamma1, gamma2, gamma3, beta, delta, c)
  solver.elements.u_tmp2 .= zero(eltype(solver.elements.u_tmp2))
  solver.elements.u_tmp3 .= solver.elements.u
  for stage in eachindex(c)
    t_stage = t + dt * c[stage]
    @timeit timer() "rhs" rhs!(solver, t_stage)

    delta_stage   = delta[stage]
    gamma1_stage  = gamma1[stage]
    gamma2_stage  = gamma2[stage]
    gamma3_stage  = gamma3[stage]
    beta_stage_dt = beta[stage] * dt
    @timeit timer() "Runge-Kutta step" begin
      Threads.@threads for i in eachindex(solver.elements.u)
        solver.elements.u_tmp2[i] += delta_stage * solver.elements.u[i]
        solver.elements.u[i]       = (gamma1_stage * solver.elements.u[i] +
                                      gamma2_stage * solver.elements.u_tmp2[i] +
                                      gamma3_stage * solver.elements.u_tmp3[i] +
                                      beta_stage_dt * solver.elements.u_t[i])
      end
    end
  end
end

"""
    timestep_hyp_diff_N3_erk52_3Sstar!(solver::AbstractSolver, t, dt)

Five stage, second-order acurate explicit Runge-Kutta scheme with stability region optimized for
the hyperbolic diffusion equation with LLF flux and polynomials of degree polydeg=3.
"""
function timestep_hyp_diff_N3_erk52_3Sstar!(solver::AbstractSolver, t, dt)
  # # New 3Sstar coefficients optimized for polynomials of degree polydeg=3
  # # and examples/parameters_hyp_diff_llf.toml
  # # 5 stages, order 1
  # gamma1 = @SVector [0.0000000000000000E+00, 5.2910412316555866E-01, 2.8433964362349406E-01, -1.4467571130907027E+00, 7.5592215948661057E-02]
  # gamma2 = @SVector [1.0000000000000000E+00, 2.6366970460864109E-01, 3.7423646095836322E-01, 7.8786901832431289E-01, 3.7754129043053775E-01]
  # gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 8.0043329115077388E-01, 1.3550099149374278E-01]
  # beta   = @SVector [1.9189497208340553E-01, 5.4506406707700059E-02, 1.2103893164085415E-01, 6.8582252490550921E-01, 8.7914657211972225E-01]
  # delta  = @SVector [1.0000000000000000E+00, 7.8593091509463076E-01, 1.2639038717454840E-01, 1.7726945920209813E-01, 0.0000000000000000E+00]
  # c      = @SVector [0.0000000000000000E+00, 1.9189497208340553E-01, 1.9580448818599061E-01, 2.4241635859769023E-01, 5.0728347557552977E-01]
  # 5 stages, order 2
  gamma1 = @SVector [0.0000000000000000E+00, 5.2656474556752575E-01, 1.0385212774098265E+00, 3.6859755007388034E-01, -6.3350615190506088E-01]
  gamma2 = @SVector [1.0000000000000000E+00, 4.1892580153419307E-01, -2.7595818152587825E-02, 9.1271323651988631E-02, 6.8495995159465062E-01]
  gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 4.1301005663300466E-01, -5.4537881202277507E-03]
  beta   = @SVector [4.5158640252832094E-01, 7.5974836561844006E-01, 3.7561630338850771E-01, 2.9356700007428856E-02, 2.5205285143494666E-01]
  delta  = @SVector [1.0000000000000000E+00, 1.3011720142005145E-01, 2.6579275844515687E-01, 9.9687218193685878E-01, 0.0000000000000000E+00]
  c      = @SVector [0.0000000000000000E+00, 4.5158640252832094E-01, 1.0221535725056414E+00, 1.4280257701954349E+00, 7.1581334196229851E-01]
  # # 5 stages, order 3
  # gamma1 = @SVector [0.0000000000000000E+00, 6.9362208054011210E-01, 9.1364483229179472E-01, 1.3129305757628569E+00, -1.4615811339132949E+00]
  # gamma2 = @SVector [1.0000000000000000E+00, 1.3224582239681788E+00, 2.4213162353103135E-01, -3.8532017293685838E-01, 1.5603355704723714E+00]
  # gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 3.8306787039991996E-01, -3.5683121201711010E-01]
  # beta   = @SVector [8.4476964977404881E-02, 3.0834660698015803E-01, 3.2131664733089232E-01, 2.8783574345390539E-01, 8.2199204703236073E-01]
  # delta  = @SVector [1.0000000000000000E+00, -7.6832695815481578E-01, 1.2497251501714818E-01, 1.4496404749796306E+00, 0.0000000000000000E+00]
  # c      = @SVector [0.0000000000000000E+00, 8.4476964977404881E-02, 2.8110631488732202E-01, 5.7093842145029405E-01, 7.2999896418559662E-01]

  timestep_3Sstar!(solver, t, dt, gamma1, gamma2, gamma3, beta, delta, c)
end

"""
    timestep_parsani_ketcheson_deconinck_erk94_3Sstar!(solver::AbstractSolver, t, dt)

Parsani, Ketcheson, Deconinck (2013)
  Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems
[DOI: 10.1137/120885899](https://doi.org/10.1137/120885899)
"""
function timestep_parsani_ketcheson_deconinck_erk94_3Sstar!(solver::AbstractSolver, t, dt)
  gamma1 = @SVector [0.0000000000000000E+00, -4.6556413837561301E+00, -7.7202649689034453E-01, -4.0244202720632174E+00, -2.1296873883702272E-02, -2.4350219407769953E+00, 1.9856336960249132E-02, -2.8107894116913812E-01, 1.6894354373677900E-01]
  gamma2 = @SVector [1.0000000000000000E+00, 2.4992627683300688E+00, 5.8668202764174726E-01, 1.2051419816240785E+00, 3.4747937498564541E-01, 1.3213458736302766E+00, 3.1196363453264964E-01, 4.3514189245414447E-01, 2.3596980658341213E-01]
  gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00, 7.6209857891449362E-01, -1.9811817832965520E-01, -6.2289587091629484E-01, -3.7522475499063573E-01, -3.3554373281046146E-01, -4.5609629702116454E-02]
  beta   = @SVector [2.8363432481011769E-01, 9.7364980747486463E-01, 3.3823592364196498E-01, -3.5849518935750763E-01, -4.1139587569859462E-03, 1.4279689871485013E+00, 1.8084680519536503E-02, 1.6057708856060501E-01, 2.9522267863254809E-01]
  delta  = @SVector [1.0000000000000000E+00, 1.2629238731608268E+00, 7.5749675232391733E-01, 5.1635907196195419E-01, -2.7463346616574083E-02, -4.3826743572318672E-01, 1.2735870231839268E+00, -6.2947382217730230E-01, 0.0000000000000000E+00]
  c      = @SVector [0.0000000000000000E+00, 2.8363432481011769E-01, 5.4840742446661772E-01, 3.6872298094969475E-01, -6.8061183026103156E-01, 3.5185265855105619E-01, 1.6659419385562171E+00, 9.7152778807463247E-01, 9.0515694340066954E-01]

  timestep_3Sstar!(solver, t, dt, gamma1, gamma2, gamma3, beta, delta, c)
end

"""
    timestep_parsani_ketcheson_deconinck_erk32_3Sstar!(solver::AbstractSolver, t, dt)

Parsani, Ketcheson, Deconinck (2013)
  Optimized explicit RK schemes for the spectral difference method applied to wave propagation problems
[DOI: 10.1137/120885899](https://doi.org/10.1137/120885899)
"""
function timestep_parsani_ketcheson_deconinck_erk32_3Sstar!(solver::AbstractSolver, t, dt)
  gamma1 = @SVector [0.0000000000000000E+00, -1.2664395576322218E-01, 1.1426980685848858E+00]
  gamma2 = @SVector [1.0000000000000000E+00, 6.5427782599406470E-01, -8.2869287683723744E-02]
  gamma3 = @SVector [0.0000000000000000E+00, 0.0000000000000000E+00, 0.0000000000000000E+00]
  beta   = @SVector [7.2366074728360086E-01, 3.4217876502651023E-01, 3.6640216242653251E-01]
  delta  = @SVector [1.0000000000000000E+00, 7.2196567116037724E-01, 0.0000000000000000E+00]
  c      = @SVector [0.0000000000000000E+00, 7.2366074728360086E-01, 5.9236433182015646E-01]

  timestep_3Sstar!(solver, t, dt, gamma1, gamma2, gamma3, beta, delta, c)
end

# Add implemenations for coupled Euler-gravity simulations
include("timedisc_euler_gravity.jl")


text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots, LaTeXStrings
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.20"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
begin
	alpha, beta = 8, 2
	prior(lam) = pdf(Gamma(alpha, 1/beta), lam)
	data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

	like(lam) = *([pdf(Poisson(lam),x) for x in data]...)
	posteriorUpToK(lam) = like(lam)*prior(lam)

	sig = 0.5
	foldedNormalPDF(x,mu) = (1/sqrt(2*pi*sig^2))*(exp(-(x-mu)^2/2sig^2)
													+ exp(-(x+mu)^2/2sig^2))
	foldedNormalRV(mu) = abs(rand(Normal(mu,sig)))
end

# ╔═╡ f9da18ea-3afb-11eb-1c03-eb97ca5c6918
function sampler(piProb,qProp,rvProp)
    lam = 1
    warmN, N = 10^5, 10^6
    samples = zeros(N-warmN)

    for t in 1:N
        while true
            lamTry = rvProp(lam)
            L = piProb(lamTry)/piProb(lam)
            H = min(1,L*qProp(lam,lamTry)/qProp(lamTry,lam))
            if rand() < H
                lam = lamTry
                if t > warmN
                    samples[t-warmN] = lam
                end
                break
            end
        end
    end
    return samples
end

# ╔═╡ f9da4af4-3afb-11eb-3be6-c136b5d61c81
begin
	mcmcSamples = sampler(posteriorUpToK,foldedNormalPDF,foldedNormalRV)
	Text("MCMC Bayes Estimate: $(mean(mcmcSamples))")
end

# ╔═╡ f9dae8a6-3afb-11eb-3b01-15befda98071
begin
	stephist(mcmcSamples, bins=100,
		c=:black, normed=true, label="Histogram of MCMC samples")

	lamRange = 0:0.01:10
	plot!(lamRange, prior.(lamRange),
		c=:blue, label="Prior distribution")

	closedFormPosterior(lam)=pdf(Gamma(alpha + sum(data),1/(beta+length(data))),lam)
	plot!(lamRange, closedFormPosterior.(lamRange),
		c=:red, label="Posterior distribution",
		xlims=(0, 10), ylims=(0, 1.2),
		xlabel=L"\lambda",ylabel="Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.20"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═f9da18ea-3afb-11eb-1c03-eb97ca5c6918
# ╠═f9da4af4-3afb-11eb-3be6-c136b5d61c81
# ╠═f9dae8a6-3afb-11eb-3b01-15befda98071
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12