text:
module HeartBeats

export detect_heartbeats, example_ecg

using DelimitedFiles, DSP, Statistics

"""
    example_ecg()

Return example ECG data as obtained from `scipy.misc.electrocardiogram()`.
"""
function example_ecg()
    vec(readdlm(joinpath(pkgdir(HeartBeats, "data"), "ecg.txt")))
end

"""
    detect_heartbeats(ecg, fs)

Detect heartbeats in ECG signal `ecg` (recorded with sampling frequency `fs`).
"""
function detect_heartbeats(ecg::AbstractVector, fs::Real)
    # Pan-Tompkins detector based on https://github.com/cbrnr/sleepecg
    # bandpass filter between 5-30Hz
    filtered = filtfilt(digitalfilter(Bandpass(5, 30, fs=fs), Butterworth(2)), ecg)

    # set samples before first zero crossing (within first two seconds) to zero
    start = findfirst(!iszero, diff(signbit.(filtered[1:2*fs])))
    !isnothing(start) && filtered[1:start] .= 0

    # 5-point derivative
    derivative = xcorr(filtered, [-1, -2, 0, 2, 1])[3:end-2]

    integrated = conv(derivative.^2, ones(trunc(Int, 0.15 * fs)))
    middle, half = length(integrated) ÷ 2, length(derivative) ÷ 2
    integrated = integrated[middle - half + 1 : middle + half]  # mode="same"

    signal_len = length(filtered)
    beat_mask = zeros(Bool, signal_len)
    refractory_samples = trunc(Int, 0.2 * fs)  # 200ms
    t_wave_samples = trunc(Int, 0.36 * fs)  # 360ms

    n = trunc(Int, 2 * fs)
    s_f = maximum(filtered[1:n])  # filtered signal
    n_f = mean(filtered[1:n])  # filtered noise
    s_i = maximum(integrated[1:n])  # integrated signal
    n_i = mean(integrated[1:n])  # integrated noise
    threshold_i1 = n_i + 0.25 * (s_i - n_i)
    threshold_f1 = n_f + 0.25 * (s_f - n_f)

    rr_intervals = zeros(signal_len ÷ refractory_samples)
    num_peaks_found = 0

    peak_index = -refractory_samples + 1
    previous_peak_index = -refractory_samples + 1
    RR_missed_limit = nothing
    RR_low_limit = nothing
    RR_high_limit = nothing

    # searchback
    do_searchback = true
    index = 2
    while index < signal_len
        PEAKF = nothing
        PEAKI = nothing
        signal_peak_found = false
        noise_peak_found = false
        if (num_peaks_found > 1 && index - previous_peak_index > RR_missed_limit && do_searchback) ||
           (num_peaks_found == 0 && index > fs) ||
           (num_peaks_found == 1 && index - previous_peak_index > 1.5 * fs)
            for i = 1:16
                found_a_candidate = false
                searchback_divisor = 2^i
                best_searchback_index = previous_peak_index + refractory_samples
                best_candidate_amplitude = -1
                searchback_index = best_searchback_index
                while searchback_index < index
                    PEAKF = filtered[searchback_index]
                    if PEAKF > filtered[searchback_index + 1]
                        if PEAKF > filtered[searchback_index - 1]
                            PEAKI = integrated[searchback_index]
                            if (threshold_f1 / searchback_divisor < PEAKF && PEAKF < threshold_f1 && threshold_i1 / searchback_divisor < PEAKI) ||
                               (threshold_i1 / searchback_divisor < PEAKI && PEAKI < threshold_i1 && threshold_f1 / searchback_divisor < PEAKF)
                                if PEAKF > best_candidate_amplitude
                                    best_searchback_index = searchback_index
                                    best_candidate_amplitude = filtered[searchback_index]
                                    found_a_candidate = true
                                end
                            end
                        end
                        searchback_index += 1
                    end
                    searchback_index += 1
                end
                if found_a_candidate
                    s_i = 0.25 * PEAKI + 0.75 * s_i
                    s_f = 0.25 * PEAKF + 0.75 * s_f
                    signal_peak_found = true
                    peak_index = best_searchback_index
                    do_searchback = false
                    break
                end
            end
        elseif filtered[index] > filtered[index + 1]
            if filtered[index] > filtered[index - 1]
                PEAKF = filtered[index]
                PEAKI = integrated[index]
                if PEAKF > threshold_f1 && PEAKI > threshold_i1
                    s_f = 0.125 * PEAKF + 0.875 * s_f
                    s_i = 0.125 * PEAKI + 0.875 * s_i

                    signal_peak_found = true
                    peak_index = index
                else
                    noise_peak_found = true
                end
            end
            index += 1
        end
        if signal_peak_found && num_peaks_found > 0
            RR = peak_index - previous_peak_index

            # t-wave identification
            if RR < t_wave_samples
                reverse_index = peak_index
                max_slope_in_this_peak = -1
                while reverse_index > 0
                    amplitude_here = filtered[reverse_index]
                    amplitude_before = filtered[reverse_index - 1]
                    if amplitude_before > amplitude_here
                        break
                    end
                    slope = amplitude_here - amplitude_before
                    if slope > max_slope_in_this_peak
                        max_slope_in_this_peak = slope
                    end
                    reverse_index -= 1
                end

                reverse_index = previous_peak_index
                max_slope_in_previous_peak = -1
                while reverse_index > 0
                    amplitude_here = filtered[reverse_index]
                    amplitude_before = filtered[reverse_index - 1]
                    if amplitude_before > amplitude_here
                        break
                    end
                    slope = amplitude_here - amplitude_before
                    if (slope > max_slope_in_previous_peak)
                        max_slope_in_previous_peak = slope
                    end
                    reverse_index -= 1
                end

                if max_slope_in_this_peak < max_slope_in_previous_peak / 2.0
                    signal_peak_found = false
                    noise_peak_found = true
                end
            end
        end

        if signal_peak_found
            num_peaks_found += 1
            beat_mask[peak_index] = true

            threshold_i1 = n_i + 0.25 * (s_i - n_i)
            threshold_f1 = n_f + 0.25 * (s_f - n_f)

            if num_peaks_found > 1
                rr_intervals[num_peaks_found] = peak_index - previous_peak_index

                # learning phase 2
                if num_peaks_found == 2
                    RR_low_limit = 0.92 * rr_intervals[num_peaks_found]
                    RR_high_limit = 1.16 * rr_intervals[num_peaks_found]
                end

                # RR average 1 / RR average 2
                RR_sum = 0
                RR_count = 0
                irregular = false
                for i = num_peaks_found:-1:1
                    RR_n = rr_intervals[i]
                    if RR_low_limit < RR_n && RR_n < RR_high_limit
                        RR_sum += RR_n
                        RR_count += 1
                        if RR_count >= 8
                            break
                        end
                    else
                        irregular = true
                    end
                end

                RR_average = RR_sum / RR_count
                RR_low_limit = 0.92 * RR_average
                RR_high_limit = 1.16 * RR_average
                RR_missed_limit = 1.66 * RR_average

                if irregular
                    threshold_f1 /= 2
                    threshold_i1 /= 2
                end
            end

            do_searchback = true
            previous_peak_index = peak_index
            index = peak_index + refractory_samples

        elseif noise_peak_found
            n_i = 0.125 * PEAKI + 0.875 * n_i
            n_f = 0.125 * PEAKF + 0.875 * n_f
            threshold_i1 = n_i + 0.25 * (s_i - n_i)
            threshold_f1 = n_f + 0.25 * (s_f - n_f)
        end
        index += 1
    end
    findall(beat_mask)
end

end


text:
# Author: Mathias Louboutin, mlouboutin3@gatech.edu
# Date: December 2020

using SlimOptim

# Rosenbrock function
f(m) = ((1-m[1])^2 + 100*(m[2]-m[1]^2)^2)
g!(g,m) = (g .= [-2*(1-m[1]) - 200*(m[2]-m[1]^2)*2*m[1], 200*(m[2]-m[1]^2)])
prj(x) = x

function obj(x)
    g = zeros(2)
    ϕ = f(x)
    g!(g,x)
    ϕ, g
end

function fg!(g,x)
    !isnothing(g) && g!(g,x)
    f(x)
end

function fg!(ff,g,x)
    !isnothing(g) && g!(g,x)
    !isnothing(ff) && return f(x)
end

niter = 1000
for (opt, algo)=zip([spg_options, pqn_options], [spg, pqn])
    x0 = zeros(2)
    options = opt(maxIter=niter, progTol=1f-30, optTol=0, iniStep=2, memory=5)
    algo(obj,x0,prj,options)
    @test x0 ≈ [1, 1]

    x1 = zeros(2)
    algo(f,g!,fg!,x1,prj,options)
    @show x1 ≈ [1, 1]
end

ext: jl


text:
# The geometry of the staggerd grid used by FiniteDifference in the vertical (sometimes called the C-grid)
# is (in one dimension) shown below

#         face   cell   face   cell   face
# left                                        right
#                 i-1            i
#                  ↓             ↓
#           |      ×      |      ×      |
#           ↑             ↑             ↑
#          i-1            i            i+1

# Gradient:
#  scalar => d₃ = D₃ x = Covariant3Vector(x[i+half] - x[i-half])
# ∂ξ³∂x * d₃

# Divergence: 1/J * D₃ (Jv³):
#  Jv³ = J * contravariant3(v)
#  div = 1/J * (Jv³[i+half] - Jv³[i-half])

#====

interpolate, multiply, gradient

  --->
ca  fa  fb  cb
    1 x F
  x       \
1           1    Boundary indices
  \       /
    2 - 2        --------------
  /       \
2           2    Interior indices
  \       /
    3 - 3
  /       \l
3           3

gradient with dirichlet

ca  fa
D - 1   modified stencil
  /
1
  \
    2   regular stencil
  /
2
  \
    3
  /
3

c2c laplacian, dirichlet

D - 1
  /   \
1       1
  \   /
    2   ---------------------
  /   \
2       2
  \   /
    3
  /   \
3

c2c laplacian, neumann
c   f   c
    N
      \
1       1
  \   /
    2  ----------------------
  /   \
2       2
  \   /
    n
  /   \
n       n
      /
   n+1
    =N

f2f laplacian, dirichlet
LaplacianF2F(left=Dirichlet(D))(x) = GradientC2F(left=Neumann(0))(GradientF2C(left=Dirichlet(D))(x))

dθ/dt = - ∇^2 θ
  --->
fa  ca  fb
D       _   effectively ignored: -dD/dt since is set by equations: use 2-point, 1-sided stencil?
  \
    1       boundary window
  /   \
2       2
  \   /
    2     ---------------
  /   \
3       3   interior
  \   /
    3

interior_indices1 = 2:...
interior_indices2 = 3:...

f2f laplacian, neumann
LaplacianF2F(left=Neumann(N))(x) = GradientC2F(left=Dirichlet(N))(GradientF2C(left=nothing)(x))

  --->
fa  ca  fb
1   N - 1   set by 1-sided stencil
  \   /
    1       boundary window
  /   \     ---------------
2       2
  \   /
    2
  /   \
3       3   interior
  \   /
    3

interior_indices1 = 1:n
interior_indices2 = 2:n-1
===#

"""
    PlusHalf(i)

Represents `i + 1/2`, but stored as internally as an integer value. Used for
indexing into staggered finite difference meshes: the convention "half" values
are indexed at cell faces, whereas centers are indexed at cell centers.

Supports `+`, `-` and inequalities.

See also [`half`](@ref).
"""
struct PlusHalf{I <: Integer} <: Real
    i::I
end
PlusHalf{I}(h::PlusHalf{I}) where {I <: Integer} = h

"""
    const half = PlusHalf(0)
"""
const half = PlusHalf(0)

Base.:+(h::PlusHalf) = h
Base.:-(h::PlusHalf) = PlusHalf(-h.i - one(h.i))
Base.:+(i::Integer, h::PlusHalf) = PlusHalf(i + h.i)
Base.:+(h::PlusHalf, i::Integer) = PlusHalf(h.i + i)
Base.:+(h1::PlusHalf, h2::PlusHalf) = h1.i + h2.i + one(h1.i)
Base.:-(i::Integer, h::PlusHalf) = PlusHalf(i - h.i - one(h.i))
Base.:-(h::PlusHalf, i::Integer) = PlusHalf(h.i - i)
Base.:-(h1::PlusHalf, h2::PlusHalf) = h1.i - h2.i

Base.:<=(h1::PlusHalf, h2::PlusHalf) = h1.i <= h2.i
Base.:<(h1::PlusHalf, h2::PlusHalf) = h1.i < h2.i
Base.max(h1::PlusHalf, h2::PlusHalf) = PlusHalf(max(h1.i, h2.i))
Base.min(h1::PlusHalf, h2::PlusHalf) = PlusHalf(min(h1.i, h2.i))

Base.convert(::Type{P}, i::Integer) where {P <: PlusHalf} =
    throw(InexactError(:convert, P, i))
Base.convert(::Type{I}, h::PlusHalf) where {I <: Integer} =
    throw(InexactError(:convert, I, h))

Base.step(::AbstractUnitRange{PlusHalf{I}}) where {I} = one(I)

left_idx(::Spaces.CenterFiniteDifferenceSpace) = 1
right_idx(space::Spaces.CenterFiniteDifferenceSpace) =
    length(space.center_local_geometry)
left_idx(::Spaces.FaceFiniteDifferenceSpace) = half
right_idx(space::Spaces.FaceFiniteDifferenceSpace) =
    PlusHalf(length(space.center_local_geometry))

left_face_boundary_idx(space::Spaces.FiniteDifferenceSpace) =
    left_idx(Spaces.FaceFiniteDifferenceSpace(space))
right_face_boundary_idx(space::Spaces.FiniteDifferenceSpace) =
    right_idx(Spaces.FaceFiniteDifferenceSpace(space))
left_center_boundary_idx(space::Spaces.FiniteDifferenceSpace) =
    left_idx(Spaces.CenterFiniteDifferenceSpace(space))
right_center_boundary_idx(space::Spaces.FiniteDifferenceSpace) =
    right_idx(Spaces.CenterFiniteDifferenceSpace(space))

left_face_boundary_idx(arg) = left_face_boundary_idx(axes(arg))
right_face_boundary_idx(arg) = right_face_boundary_idx(axes(arg))
left_center_boundary_idx(arg) = left_center_boundary_idx(axes(arg))
right_center_boundary_idx(arg) = right_center_boundary_idx(axes(arg))

function Geometry.LocalGeometry(
    space::Spaces.FiniteDifferenceSpace,
    idx::Integer,
)
    space.center_local_geometry[idx]
end
function Geometry.LocalGeometry(
    space::Spaces.FiniteDifferenceSpace,
    idx::PlusHalf,
)
    space.face_local_geometry[idx.i + 1]
end

abstract type BoundaryCondition end

"""
    SetValue(val)

Set the value at the boundary to be `val`. In the case of gradient operators,
this will set the input value from which the gradient is computed.
"""
struct SetValue{S} <: BoundaryCondition
    val::S
end

"""
    SetGradient(val)

Set the gradient at the boundary to be `val`. In the case of gradient operators
this will set the output value of the gradient.
"""
struct SetGradient{S} <: BoundaryCondition
    val::S
end

"""
    SetDivergence(val)

Set the divergence at the boundary to be `val`.
"""
struct SetDivergence{S} <: BoundaryCondition
    val::S
end

"""
    Extrapolate()

Set the value at the boundary to be the same as the closest interior point.
"""
struct Extrapolate <: BoundaryCondition end

abstract type Location end
abstract type Boundary <: Location end
abstract type BoundaryWindow <: Location end

struct Interior <: Location end
struct LeftBoundaryWindow{name} <: BoundaryWindow end
struct RightBoundaryWindow{name} <: BoundaryWindow end

abstract type FiniteDifferenceOperator end

return_eltype(::FiniteDifferenceOperator, arg) = eltype(arg)

boundary_width(op::FiniteDifferenceOperator, bc) =
    error("Boundary $(typeof(bc)) is not supported for operator $(typeof(op))")

get_boundary(
    op::FiniteDifferenceOperator,
    ::LeftBoundaryWindow{name},
) where {name} = getproperty(op.bcs, name)

get_boundary(
    op::FiniteDifferenceOperator,
    ::RightBoundaryWindow{name},
) where {name} = getproperty(op.bcs, name)

has_boundary(
    op::FiniteDifferenceOperator,
    ::LeftBoundaryWindow{name},
) where {name} = hasproperty(op.bcs, name)

has_boundary(
    op::FiniteDifferenceOperator,
    ::RightBoundaryWindow{name},
) where {name} = hasproperty(op.bcs, name)

abstract type AbstractStencilStyle <: Fields.AbstractFieldStyle end

# the .f field is an operator
struct StencilStyle <: AbstractStencilStyle end

# the .f field is a function, but one of the args (or their args) is a StencilStyle
struct CompositeStencilStyle <: AbstractStencilStyle end

"""
    stencil_interior_width(op::Op)

The width of the interior stencil for operator of type `Op`. This should return a tuple
of 2-tuples: each 2-tuple should be the lower and upper bounds of the index offsets of
the stencil for each argument in the stencil.
"""
function stencil_interior_width end

abstract type InterpolationOperator <: FiniteDifferenceOperator end

"""
    InterpolateF2C()

Interpolate from face to center mesh. No boundary conditions are required (or supported).
"""
struct InterpolateF2C{BCS} <: InterpolationOperator
    bcs::BCS
end
InterpolateF2C(; kwargs...) = InterpolateF2C(NamedTuple(kwargs))

return_space(::InterpolateF2C, space::Spaces.FaceFiniteDifferenceSpace) =
    Spaces.CenterFiniteDifferenceSpace(space)

return_space(
    ::InterpolateF2C,
    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
) = Spaces.CenterExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::InterpolateF2C) = ((-half, half),)

function stencil_interior(::InterpolateF2C, loc, idx, arg)
    a⁺ = getidx(arg, loc, idx + half)
    a⁻ = getidx(arg, loc, idx - half)
    RecursiveApply.rdiv(a⁺ ⊞ a⁻, 2)
end

"""
    I = InterpolateC2F(;boundaries..)
    I.(x)

Interpolate a center-valued field `x` to faces, using the stencil
```math
I(x)[i] = \\frac{1}{2} (x[i+\\tfrac{1}{2}] + x[i-\\tfrac{1}{2}])
```

Supported boundary conditions are:

- [`SetValue(x₀)`](@ref): set the value at the boundary face to be `x₀`. On the
  left boundary the stencil is
```math
I(x)[\\tfrac{1}{2}] = x₀
```
- [`SetGradient(v)`](@ref): set the value at the boundary such that the gradient
  is `v`. At the left boundary the stencil is
```math
I(x)[\\tfrac{1}{2}] = x[1] - \\frac{1}{2} v³
```
- [`Extrapolate`](@ref): use the closest interior point as the boundary value.
  At the left boundary the stencil is
```math
I(x)[\\tfrac{1}{2}] = x[1]
```
"""
struct InterpolateC2F{BCS} <: InterpolationOperator
    bcs::BCS
end
InterpolateC2F(; kwargs...) = InterpolateC2F(NamedTuple(kwargs))

return_space(::InterpolateC2F, space::Spaces.CenterFiniteDifferenceSpace) =
    Spaces.FaceFiniteDifferenceSpace(space)

return_space(
    ::InterpolateC2F,
    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::InterpolateC2F) = ((-half, half),)

function stencil_interior(::InterpolateC2F, loc, idx, arg)
    a⁺ = getidx(arg, loc, idx + half)
    a⁻ = getidx(arg, loc, idx - half)
    RecursiveApply.rdiv(a⁺ ⊞ a⁻, 2)
end

boundary_width(op::InterpolateC2F, ::SetValue) = 1

function stencil_left_boundary(::InterpolateC2F, bc::SetValue, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    bc.val
end

function stencil_right_boundary(::InterpolateC2F, bc::SetValue, loc, idx, arg)
    @assert idx == right_face_boundary_idx(arg)
    bc.val
end

boundary_width(op::InterpolateC2F, ::SetGradient) = 1

function stencil_left_boundary(::InterpolateC2F, bc::SetGradient, loc, idx, arg)
    space = axes(arg)
    @assert idx == left_face_boundary_idx(space)
    a⁺ = getidx(arg, loc, idx + half)
    v₃ = Geometry.covariant3(bc.val, Geometry.LocalGeometry(space, idx))
    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
end

function stencil_right_boundary(
    ::InterpolateC2F,
    bc::SetGradient,
    loc,
    idx,
    arg,
)
    space = axes(arg)
    @assert idx == right_face_boundary_idx(space)
    a⁻ = getidx(arg, loc, idx - half)
    v₃ = Geometry.covariant3(bc.val, Geometry.LocalGeometry(space, idx))
    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
end

boundary_width(op::InterpolateC2F, ::Extrapolate) = 1

function stencil_left_boundary(::InterpolateC2F, bc::Extrapolate, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    a⁺ = getidx(arg, loc, idx + half)
    a⁺
end

function stencil_right_boundary(
    ::InterpolateC2F,
    bc::Extrapolate,
    loc,
    idx,
    arg,
)
    @assert idx == right_face_boundary_idx(arg)
    a⁻ = getidx(arg, loc, idx - half)
    a⁻
end

"""
    WI = WeightedInterpolateF2C(; boundaries)
    WI.(w, x)

Interpolate a face-valued field `x` to centers, weighted by a face-valued field
`w`, using the stencil
```math
WI(w, x)[i] = \\frac{
        w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
    }{
        2 (w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}])
    }
```

No boundary conditions are required (or supported)
"""
struct WeightedInterpolateF2C{BCS} <: InterpolationOperator
    bcs::BCS
end
WeightedInterpolateF2C(; kwargs...) = InterpolateF2C(NamedTuple(kwargs))

function return_space(
    ::WeightedInterpolateF2C,
    weight_space::Spaces.FaceFiniteDifferenceSpace,
    value_space::Spaces.FaceFiniteDifferenceSpace,
)
    Spaces.CenterFiniteDifferenceSpace(value_space)
end

function return_space(
    ::WeightedInterpolateF2C,
    weight_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
    value_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
)
    Spaces.CenterExtrudedFiniteDifferenceSpace(value_space)
end

stencil_interior_width(::WeightedInterpolateF2C) = ((0, 0), (-half, half))

function stencil_interior(
    ::WeightedInterpolateF2C,
    loc,
    idx,
    weight_field,
    value_field,
)
    w⁺ = getidx(weight_field, loc, idx + half)
    w⁻ = getidx(weight_field, loc, idx - half)
    a⁺ = getidx(value_field, loc, idx + half)
    a⁻ = getidx(value_field, loc, idx - half)
    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (2 ⊠ (w⁺ ⊞ w⁻)))
end

"""
    WI = WeightedInterpolateC2F(; boundaries)
    WI.(w, x)

Interpolate a center-valued field `x` to faces, weighted by a center-valued field
`w`, using the stencil
```math
WI(w, x)[i] = \\frac{
    w[i+\\tfrac{1}{2}] x[i+\\tfrac{1}{2}] +  w[i-\\tfrac{1}{2}] x[i-\\tfrac{1}{2}])
}{
    2 (w[i+\\tfrac{1}{2}] + w[i-\\tfrac{1}{2}])
}
```

Supported boundary conditions are:

- [`SetValue(val)`](@ref): set the value at the boundary face to be `val`.
- [`SetGradient`](@ref): set the value at the boundary such that the gradient is `val`.
- [`Extrapolate`](@ref): use the closest interior point as the boundary value
"""
struct WeightedInterpolateC2F{BCS} <: InterpolationOperator
    bcs::BCS
end
WeightedInterpolateC2F(; kwargs...) = WeightedInterpolateC2F(NamedTuple(kwargs))

function return_space(
    ::WeightedInterpolateC2F,
    weight_space::Spaces.CenterFiniteDifferenceSpace,
    value_space::Spaces.CenterFiniteDifferenceSpace,
)
    Spaces.FaceFiniteDifferenceSpace(value_space)
end

function return_space(
    ::WeightedInterpolateC2F,
    weight_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
    value_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
)
    Spaces.FaceExtrudedFiniteDifferenceSpace(value_space)
end

stencil_interior_width(::WeightedInterpolateC2F) = ((0, 0), (0, 1))

function stencil_interior(
    ::WeightedInterpolateC2F,
    loc,
    idx,
    weight_field,
    value_field,
)
    w⁺ = getidx(weight_field, loc, idx + half)
    w⁻ = getidx(weight_field, loc, idx - half)
    a⁺ = getidx(value_field, loc, idx + half)
    a⁻ = getidx(value_field, loc, idx - half)
    RecursiveApply.rdiv((w⁺ ⊠ a⁺) ⊞ (w⁻ ⊠ a⁻), (2 ⊠ (w⁺ ⊞ w⁻)))
end

boundary_width(op::WeightedInterpolateC2F, ::SetValue) = 1

function stencil_left_boundary(
    ::WeightedInterpolateC2F,
    bc::SetValue,
    loc,
    idx,
    weight_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_face_boundary_idx(space)
    bc.val
end

function stencil_right_boundary(
    ::WeightedInterpolateC2F,
    bc::SetValue,
    loc,
    idx,
    weight_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_face_boundary_idx(space)
    bc.val
end

boundary_width(op::WeightedInterpolateC2F, ::SetGradient) = 1

function stencil_left_boundary(
    ::WeightedInterpolateC2F,
    bc::SetGradient,
    loc,
    idx,
    weight_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_face_boundary_idx(space)
    a⁺ = getidx(arg, loc, idx + half)
    v₃ = Geometry.covariant3(bc.val, Geometry.LocalGeometry(space, idx))
    a⁺ ⊟ RecursiveApply.rdiv(v₃, 2)
end

function stencil_right_boundary(
    ::WeightedInterpolateC2F,
    bc::SetGradient,
    loc,
    idx,
    weight_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_face_boundary_idx(space)
    a⁻ = getidx(value_field, loc, idx - half)
    v₃ = Geometry.covariant3(bc.val, Geometry.LocalGeometry(space, idx))
    a⁻ ⊞ RecursiveApply.rdiv(v₃, 2)
end

boundary_width(op::WeightedInterpolateC2F, ::Extrapolate) = 1

function stencil_left_boundary(
    ::WeightedInterpolateC2F,
    bc::Extrapolate,
    loc,
    idx,
    weight_field,
    value_field,
)
    @assert idx == left_face_boundary_idx(space)
    a⁺ = getidx(arg, loc, idx + half)
    a⁺
end

function stencil_right_boundary(
    ::WeightedInterpolateC2F,
    bc::Extrapolate,
    loc,
    idx,
    weight_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_face_boundary_idx(space)
    a⁻ = getidx(value_field, loc, idx - half)
    a⁻
end

"""
    L = LeftBiasedC2F(;boundaries)
    L.(x)

Interpolate a left-value field to a face-valued field from the left.
```math
L(x)[i] = x[i-\\tfrac{1}{2}]
```

Only the left boundary condition should be set. Currently supported is:
- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
```math
L(x)[\\tfrac{1}{2}] = x_0
```
"""
struct LeftBiasedC2F{BCS} <: InterpolationOperator
    bcs::BCS
end
LeftBiasedC2F(; kwargs...) = LeftBiasedC2F(NamedTuple(kwargs))

return_space(::LeftBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) =
    Spaces.FaceFiniteDifferenceSpace(space)

return_space(
    ::LeftBiasedC2F,
    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::LeftBiasedC2F) = ((-half, -half),)

stencil_interior(::LeftBiasedC2F, loc, idx, arg) = getidx(arg, loc, idx - half)

boundary_width(op::LeftBiasedC2F, ::SetValue) = 1

function stencil_left_boundary(::LeftBiasedC2F, bc::SetValue, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    bc.val
end

"""
    R = RightBiasedC2F(;boundaries)
    R.(x)

Interpolate a center-valued field to a face-valued field from the right.
```math
R(x)[i] = x[i+\\tfrac{1}{2}]
```

Only the right boundary condition should be set. Currently supported is:
- [`SetValue(x₀)`](@ref): set the value to be `x₀` on the boundary.
```math
R(x)[n+\\tfrac{1}{2}] = x_0
```
"""
struct RightBiasedC2F{BCS} <: InterpolationOperator
    bcs::BCS
end
RightBiasedC2F(; kwargs...) = RightBiasedC2F(NamedTuple(kwargs))

return_space(::RightBiasedC2F, space::Spaces.CenterFiniteDifferenceSpace) =
    Spaces.FaceFiniteDifferenceSpace(space)

return_space(
    ::RightBiasedC2F,
    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::RightBiasedC2F) = ((half, half),)

stencil_interior(::RightBiasedC2F, loc, idx, arg) = getidx(arg, loc, idx + half)

boundary_width(op::RightBiasedC2F, ::SetValue) = 1

function stencil_right_boundary(::RightBiasedC2F, bc::SetValue, loc, idx, arg)
    @assert idx == right_face_boundary_idx(arg)
    bc.val
end

abstract type AdvectionOperator <: FiniteDifferenceOperator end

return_eltype(::AdvectionOperator, velocity, arg) = eltype(arg)

"""
    U = UpwindBiasedProductC2F(;boundaries)
    U.(v, x)

Compute the product of a face-valued vector field `v` and a center-valued field
`x` at cell faces by upwinding `x` according to `v`
```math
U(v,x)[i] = \\begin{cases}
  v[i] x[i-\\tfrac{1}{2}] \\quad v^3[i] > 0 \\\\
  v[i] x[i+\\tfrac{1}{2}] \\quad v^3[i] < 0
  \\end{cases}
```

Supported boundary conditions are:
- [`SetValue(x₀)`](@ref): set the value of `x` to be `x₀` on the boundary. On
  the left boundary the stencil is
  ```math
  U(v,x)[\\tfrac{1}{2}] = \\begin{cases}
    v[\\tfrac{1}{2}] x_0 \\quad v^3[\\tfrac{1}{2}] > 0 \\\\
    v[\\tfrac{1}{2}] x[1] \\quad v^3[\\tfrac{1}{2}] < 0
    \\end{cases}
  ```
- [`Extrapolate()`](@ref): set the value of `x` to be the same as the closest
  interior point. On the left boundary, the stencil is
  U(v,x)[\\tfrac{1}{2}] = U(v,x)[1 + \\tfrac{1}{2}]
"""
struct UpwindBiasedProductC2F{BCS} <: AdvectionOperator
    bcs::BCS
end
UpwindBiasedProductC2F(; kwargs...) = UpwindBiasedProductC2F(NamedTuple(kwargs))

return_eltype(::UpwindBiasedProductC2F, V, A) =
    Geometry.Contravariant3Vector{eltype(eltype(V))}

function return_space(
    ::UpwindBiasedProductC2F,
    velocity_space::Spaces.FaceFiniteDifferenceSpace,
    field_space::Spaces.CenterFiniteDifferenceSpace,
)
    velocity_space
end

function return_space(
    ::UpwindBiasedProductC2F,
    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
    field_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
)
    velocity_space
end

stencil_interior_width(::UpwindBiasedProductC2F) = ((0, 0), (-half, half))

function upwind_biased_product(v, a⁻, a⁺)
    RecursiveApply.rdiv(
        ((v ⊞ RecursiveApply.rmap(abs, v)) ⊠ a⁻) ⊞
        ((v ⊟ RecursiveApply.rmap(abs, v)) ⊠ a⁺),
        2,
    )
end

function stencil_interior(
    ::UpwindBiasedProductC2F,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    a⁻ = stencil_interior(LeftBiasedC2F(), loc, idx, value_field)
    a⁺ = stencil_interior(RightBiasedC2F(), loc, idx, value_field)
    vᶠ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx),
        Geometry.LocalGeometry(space, idx),
    )
    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, a⁺))
end

boundary_width(op::UpwindBiasedProductC2F, ::SetValue) = 1

function stencil_left_boundary(
    ::UpwindBiasedProductC2F,
    bc::SetValue,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_face_boundary_idx(space)
    aᴸᴮ = bc.val
    a⁺ = stencil_interior(RightBiasedC2F(), loc, idx, value_field)
    vᶠ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx),
        Geometry.LocalGeometry(space, idx),
    )
    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, aᴸᴮ, a⁺))
end

function stencil_right_boundary(
    ::UpwindBiasedProductC2F,
    bc::SetValue,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_face_boundary_idx(space)
    a⁻ = stencil_interior(LeftBiasedC2F(), loc, idx, value_field)
    aᴿᴮ = bc.val
    vᶠ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx),
        Geometry.LocalGeometry(space, idx),
    )
    return Geometry.Contravariant3Vector(upwind_biased_product(vᶠ, a⁻, aᴿᴮ))
end

boundary_width(op::UpwindBiasedProductC2F, ::Extrapolate) = 1

function stencil_left_boundary(
    op::UpwindBiasedProductC2F,
    ::Extrapolate,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_face_boundary_idx(space)
    stencil_interior(op, loc, idx + 1, velocity_field, value_field)
end

function stencil_right_boundary(
    op::UpwindBiasedProductC2F,
    ::Extrapolate,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_face_boundary_idx(space)
    stencil_interior(op, loc, idx - 1, velocity_field, value_field)
end

"""
    A = AdvectionF2F(;boundaries)
    A.(v, θ)

Vertical advection operator at cell faces, for a face-valued velocity field `v` and face-valued
variables `θ`, approximating ``v^3 \\partial_3 \\theta``.

It uses the following stencil
```math
A(v,θ)[i] = \\frac{1}{2} (θ[i+1] - θ[i-1]) v³[i]
```

No boundary conditions are currently supported.
```
"""
struct AdvectionF2F{BCS} <: AdvectionOperator
    bcs::BCS
end

AdvectionF2F(; kwargs...) = AdvectionF2F(NamedTuple(kwargs))

stencil_interior_width(::AdvectionF2F) = ((0, 0), (-1, 1))

function return_space(
    ::AdvectionF2F,
    velocity_space::Spaces.FaceFiniteDifferenceSpace,
    field_space::Spaces.FaceFiniteDifferenceSpace,
)
    field_space
end

function return_space(
    ::AdvectionF2F,
    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
    field_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
)
    field_space
end

function stencil_interior(::AdvectionF2F, loc, idx, velocity_field, value_field)
    space = axes(value_field)
    θ⁺ = getidx(value_field, loc, idx + 1)
    θ⁻ = getidx(value_field, loc, idx - 1)
    w³ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx),
        Geometry.LocalGeometry(space, idx),
    )
    ∂θ₃ = RecursiveApply.rdiv(θ⁺ ⊟ θ⁻, 2)
    return w³ ⊠ ∂θ₃
end

"""
    A = AdvectionC2C(;boundaries)
    A.(v, θ)

Vertical advection operator at cell centers, for cell face velocity field `v` cell center
variables `θ`, approximating ``v^3 \\partial_3 \\theta``.

It uses the following stencil
```math
A(v,θ)[i] = \\frac{1}{2} \\{ (θ[i+1] - θ[i]) v³[i+\\tfrac{1}{2}] + (θ[i] - θ[i-1])v³[i-\\tfrac{1}{2}]\\}
```

Supported boundary conditions:

- [`SetValue(θ₀)`](@ref): set the value of `θ` at the boundary face to be `θ₀`.
  At the lower boundary, this is:
```math
A(v,θ)[1] = \\frac{1}{2} \\{ (θ[2] - θ[1]) v³[1 + \\tfrac{1}{2}] + (θ[1] - θ₀)v³[\\tfrac{1}{2}]\\}
```
- [`Extrapolate`](@ref): use the closest interior point as the boundary value.
  At the lower boundary, this is:
```math
A(v,θ)[1] = (θ[2] - θ[1]) v³[1 + \\tfrac{1}{2}] \\}
```
"""
struct AdvectionC2C{BCS} <: AdvectionOperator
    bcs::BCS
end

AdvectionC2C(; kwargs...) = AdvectionC2C(NamedTuple(kwargs))

stencil_interior_width(::AdvectionC2C) = ((-half, +half), (-1, 1))

function return_space(
    ::AdvectionC2C,
    velocity_space::Spaces.FaceFiniteDifferenceSpace,
    field_space::Spaces.CenterFiniteDifferenceSpace,
)
    field_space
end

function return_space(
    ::AdvectionC2C,
    velocity_space::Spaces.FaceExtrudedFiniteDifferenceSpace,
    field_space::Spaces.CenterExtrudedFiniteDifferenceSpace,
)
    field_space
end

function stencil_interior(::AdvectionC2C, loc, idx, velocity_field, value_field)
    space = axes(value_field)
    θ⁺ = getidx(value_field, loc, idx + 1)
    θ = getidx(value_field, loc, idx)
    θ⁻ = getidx(value_field, loc, idx - 1)
    w³⁺ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    w³⁻ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    ∂θ₃⁺ = θ⁺ ⊟ θ
    ∂θ₃⁻ = θ ⊟ θ⁻
    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
end

boundary_width(op::AdvectionC2C, ::SetValue) = 1

function stencil_left_boundary(
    ::AdvectionC2C,
    bc::SetValue,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_center_boundary_idx(space)
    θ⁺ = getidx(value_field, loc, idx + 1)
    θ = getidx(value_field, loc, idx)
    θ⁻ = bc.val # defined at face, not the center
    w³⁺ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    w³⁻ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    ∂θ₃⁺ = θ⁺ ⊟ θ
    ∂θ₃⁻ = 2 ⊠ (θ ⊟ θ⁻)
    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
end

function stencil_right_boundary(
    ::AdvectionC2C,
    bc::SetValue,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_center_boundary_idx(space)
    θ⁺ = bc.val # value at the face
    θ = getidx(value_field, loc, idx)
    θ⁻ = getidx(value_field, loc, idx - 1)
    w³⁺ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    w³⁻ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    ∂θ₃⁺ = 2 ⊠ (θ⁺ ⊟ θ)
    ∂θ₃⁻ = θ ⊟ θ⁻
    return RecursiveApply.rdiv((w³⁺ ⊠ ∂θ₃⁺) ⊞ (w³⁻ ⊠ ∂θ₃⁻), 2)
end

boundary_width(op::AdvectionC2C, ::Extrapolate) = 1

function stencil_left_boundary(
    ::AdvectionC2C,
    ::Extrapolate,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == left_center_boundary_idx(space)
    θ⁺ = getidx(value_field, loc, idx + 1)
    θ = getidx(value_field, loc, idx)
    w³⁺ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    ∂θ₃⁺ = θ⁺ ⊟ θ
    return (w³⁺ ⊠ ∂θ₃⁺)
end

function stencil_right_boundary(
    ::AdvectionC2C,
    ::Extrapolate,
    loc,
    idx,
    velocity_field,
    value_field,
)
    space = axes(value_field)
    @assert idx == right_center_boundary_idx(space)
    θ = getidx(value_field, loc, idx)
    θ⁻ = getidx(value_field, loc, idx - 1)
    w³⁻ = Geometry.contravariant3(
        getidx(velocity_field, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    ∂θ₃⁻ = θ ⊟ θ⁻
    return (w³⁻ ⊠ ∂θ₃⁻)
end

abstract type BoundaryOperator <: FiniteDifferenceOperator end

"""
    SetBoundaryOperator(;boundaries...)

This operator only modifies the values at the boundary:
 - [`SetValue(val)`](@ref): set the value to be `val` on the boundary.
"""
struct SetBoundaryOperator{BCS} <: BoundaryOperator
    bcs::BCS
end
SetBoundaryOperator(; kwargs...) = SetBoundaryOperator(NamedTuple(kwargs))

return_space(::SetBoundaryOperator, space::Spaces.FaceFiniteDifferenceSpace) =
    space

return_space(
    ::SetBoundaryOperator,
    space::Spaces.FaceExtrudedFiniteDifferenceSpace,
) = space

stencil_interior_width(::SetBoundaryOperator) = ((0, 0),)

function stencil_interior(::SetBoundaryOperator, loc, idx, arg)
    getidx(arg, loc, idx)
end

boundary_width(op::SetBoundaryOperator, ::SetValue) = 1

function stencil_left_boundary(
    ::SetBoundaryOperator,
    bc::SetValue,
    loc,
    idx,
    arg,
)
    @assert idx == left_face_boundary_idx(arg)
    bc.val
end

function stencil_right_boundary(
    ::SetBoundaryOperator,
    bc::SetValue,
    loc,
    idx,
    arg,
)
    @assert idx == right_face_boundary_idx(arg)
    bc.val
end

abstract type GradientOperator <: FiniteDifferenceOperator end
# TODO: we should probably make the axis the operator is working over as part of the operator type
# similar to the spectral operators, hardcoded to vertical only `(3,)` for now
return_eltype(::GradientOperator, arg) =
    Geometry.gradient_result_type(Val((3,)), eltype(arg))
#return_eltype(::GradientOperator, arg) = Geometry.Covariant3Vector{eltype(arg)}

"""
    G = GradientF2C(;boundaryname=boundarycondition...)
    G.(x)

Compute the gradient of a face-valued field `x`, returning a center-valued
`Covariant3` vector field, using the stencil:
```math
G(x)[i]^3 = x[i+\\tfrac{1}{2}] - x[i-\\tfrac{1}{2}]
```

The following boundary conditions are supported:
 - by default, the value of `x` at the boundary face will be used.
 - [`SetValue(x₀)`](@ref): calculate the gradient assuming the value at the
   boundary is `x₀`. For the left boundary, this becomes:
```math
G(x)[1]³ = x[1+\\tfrac{1}{2}] - x₀
```
- [`Extrapolate()`](@ref): set the value at the center closest to the boundary
to be the same as the neighbouring interior value. For the left boundary, this becomes:
```math
G(x)[1]³ = G(x)[2]³
```
"""
struct GradientF2C{BCS} <: GradientOperator
    bcs::BCS
end
GradientF2C(; kwargs...) = GradientF2C(NamedTuple(kwargs))

return_space(::GradientF2C, space::Spaces.FaceFiniteDifferenceSpace) =
    Spaces.CenterFiniteDifferenceSpace(space)

return_space(::GradientF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) =
    Spaces.CenterExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::GradientF2C) = ((-half, half),)

function stencil_interior(::GradientF2C, loc, idx, arg)
    Geometry.Covariant3Vector(1) ⊗
    (getidx(arg, loc, idx + half) ⊟ getidx(arg, loc, idx - half))
end

boundary_width(op::GradientF2C, ::SetValue) = 1

function stencil_left_boundary(::GradientF2C, bc::SetValue, loc, idx, arg)
    @assert idx == left_center_boundary_idx(arg)
    Geometry.Covariant3Vector(1) ⊗ (getidx(arg, loc, idx + half) ⊟ bc.val)
end

function stencil_right_boundary(::GradientF2C, bc::SetValue, loc, idx, arg)
    @assert idx == right_center_boundary_idx(arg)
    Geometry.Covariant3Vector(1) ⊗ (bc.val ⊟ getidx(arg, loc, idx - half))
end

boundary_width(op::GradientF2C, ::Extrapolate) = 1

function stencil_left_boundary(op::GradientF2C, ::Extrapolate, loc, idx, arg)
    space = axes(arg)
    @assert idx == left_center_boundary_idx(arg)
    Geometry.transform(
        Geometery.Covariant3Axis(),
        stencil_interior(op, loc, idx + 1, arg),
        Geometry.LocalGeometry(space, idx),
    )
end
function stencil_right_boundary(op::GradientF2C, ::Extrapolate, loc, idx, arg)
    space = axes(arg)
    @assert idx == right_center_boundary_idx(arg)
    Geometry.transform(
        Geometry.Covariant3Axis(),
        stencil_interior(op, loc, idx - 1, arg),
        Geometry.LocalGeometry(space, idx),
    )
end

"""
    G = GradientC2F(;boundaryname=boundarycondition...)
    G.(x)

Compute the gradient of a center-valued field `x`, returning a face-valued
`Covariant3` vector field, using the stencil:
```math
G(x)[i]^3 = x[i+\\tfrac{1}{2}] - x[i-\\tfrac{1}{2}]
```

The following boundary conditions are supported:
- [`SetValue(x₀)`](@ref): calculate the gradient assuming the value at the
  boundary is `x₀`. For the left boundary, this becomes:
  ```math
  G(x)[\\tfrac{1}{2}]³ = 2 (x[1] - x₀)
  ```
- [`SetGradient(v₀)`](@ref): set the value of the gradient at the boundary to be
  `v₀`. For the left boundary, this becomes:
  ```math
  G(x)[\\tfrac{1}{2}] = v₀
  ```
"""
struct GradientC2F{BC} <: GradientOperator
    bcs::BC
end
GradientC2F(; kwargs...) = GradientC2F(NamedTuple(kwargs))

return_space(::GradientC2F, space::Spaces.CenterFiniteDifferenceSpace) =
    Spaces.FaceFiniteDifferenceSpace(space)

return_space(::GradientC2F, space::Spaces.CenterExtrudedFiniteDifferenceSpace) =
    Spaces.FaceExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::GradientC2F) = ((-half, half),)

function stencil_interior(::GradientC2F, loc, idx, arg)
    Geometry.Covariant3Vector(1) ⊗
    (getidx(arg, loc, idx + half) ⊟ getidx(arg, loc, idx - half))
end

boundary_width(op::GradientC2F, ::SetValue) = 1
function stencil_left_boundary(::GradientC2F, bc::SetValue, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    # ∂x[i] = 2(∂x[i + half] - val)
    Geometry.Covariant3Vector(2) ⊗ (getidx(arg, loc, idx + half) ⊟ bc.val)
end
function stencil_right_boundary(::GradientC2F, bc::SetValue, loc, idx, arg)
    @assert idx == right_face_boundary_idx(arg)
    Geometry.Covariant3Vector(2) ⊗ (bc.val ⊟ getidx(arg, loc, idx - half))
end

# left / right SetGradient boundary conditions
boundary_width(op::GradientC2F, ::SetGradient) = 1
function stencil_left_boundary(::GradientC2F, bc::SetGradient, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    space = axes(arg)
    # imposed flux boundary condition at left most face
    Geometry.transform(
        Geometry.Covariant3Axis(),
        bc.val,
        Geometry.LocalGeometry(space, idx),
    )
end
function stencil_right_boundary(::GradientC2F, bc::SetGradient, loc, idx, arg)
    @assert idx == right_face_boundary_idx(arg)
    space = axes(arg)
    # imposed flux boundary condition at right most face
    Geometry.transform(
        Geometry.Covariant3Axis(),
        bc.val,
        Geometry.LocalGeometry(space, idx),
    )
end

abstract type DivergenceOperator <: FiniteDifferenceOperator end
return_eltype(::DivergenceOperator, arg) =
    Geometry.divergence_result_type(eltype(arg))

"""
    D = DivergenceF2C(;boundaryname=boundarycondition...)
    D.(v)

Compute the vertical contribution to the divergence of a face-valued field
vector `v`, returning a center-valued scalar field, using the stencil
```math
D(v)[i] = (Jv³[i+\\tfrac{1}{2}] - Jv³[i-\\tfrac{1}{2}]) / J[i]
```
where `Jv³` is the Jacobian multiplied by the third contravariant component of
`v`.

The following boundary conditions are supported:
 - by default, the value of `v` at the boundary face will be used.
 - [`SetValue(v₀)`](@ref): calculate the divergence assuming the value at the
   boundary is `v₀`. For the left boundary, this becomes:
```math
D(v)[1] = (Jv³[1+\\tfrac{1}{2}] - Jv³₀) / J[i]
```
- [`Extrapolate()`](@ref): set the value at the center closest to the boundary
  to be the same as the neighbouring interior value. For the left boundary, this
  becomes:
```math
D(v)[1]³ = D(v)[2]³
```
"""
struct DivergenceF2C{BCS} <: DivergenceOperator
    bcs::BCS
end
DivergenceF2C(; kwargs...) = DivergenceF2C(NamedTuple(kwargs))

return_space(::DivergenceF2C, space::Spaces.FaceFiniteDifferenceSpace) =
    Spaces.CenterFiniteDifferenceSpace(space)

return_space(::DivergenceF2C, space::Spaces.FaceExtrudedFiniteDifferenceSpace) =
    Spaces.CenterExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::DivergenceF2C) = ((-half, half),)

function stencil_interior(::DivergenceF2C, loc, idx, arg)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³₊ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    Ju³₋ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    RecursiveApply.rdiv(Ju³₊ ⊟ Ju³₋, local_geometry.J)
end

boundary_width(op::DivergenceF2C, ::SetValue) = 1
function stencil_left_boundary(::DivergenceF2C, bc::SetValue, loc, idx, arg)
    @assert idx == left_center_boundary_idx(arg)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³₊ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    Ju³₋ = Geometry.Jcontravariant3(
        bc.val,
        Geometry.LocalGeometry(space, idx - half),
    )
    RecursiveApply.rdiv(Ju³₊ ⊟ Ju³₋, local_geometry.J)
end
function stencil_right_boundary(::DivergenceF2C, bc::SetValue, loc, idx, arg)
    @assert idx == right_center_boundary_idx(arg)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³₊ = Geometry.Jcontravariant3(
        bc.val,
        Geometry.LocalGeometry(space, idx + half),
    )
    Ju³₋ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    RecursiveApply.rdiv(Ju³₊ ⊟ Ju³₋, local_geometry.J)
end

boundary_width(op::DivergenceF2C, ::Extrapolate) = 1
function stencil_left_boundary(op::DivergenceF2C, ::Extrapolate, loc, idx, arg)
    @assert idx == left_center_boundary_idx(arg)
    stencil_interior(op, loc, idx + 1, arg)
end
function stencil_right_boundary(op::DivergenceF2C, ::Extrapolate, loc, idx, arg)
    @assert idx == right_center_boundary_idx(arg)
    stencil_interior(op, loc, idx - 1, arg)
end

"""
    D = DivergenceC2F(;boundaryname=boundarycondition...)
    D.(v)

Compute the vertical contribution to the divergence of a center-valued field
vector `v`, returning a face-valued scalar field, using the stencil
```math
D(v)[i] = (Jv³[i+\\tfrac{1}{2}] - Jv³[i-\\tfrac{1}{2}]) / J[i]
```
where `Jv³` is the Jacobian multiplied by the third contravariant component of
`v`.

The following boundary conditions are supported:
- [`SetValue(v₀)`](@ref): calculate the divergence assuming the value at the
   boundary is `v₀`. For the left boundary, this becomes:
  ```math
  D(v)[\\tfrac{1}{2}] = \\frac{1}{2} (Jv³[1] - Jv³₀) / J[i]
  ```
- [`SetDivergence(x)`](@ref): set the value of the divergence at the boundary to be `x`.
  ```math
  D(v)[\\tfrac{1}{2}] = x
  ```
"""
struct DivergenceC2F{BC} <: DivergenceOperator
    bcs::BC
end
DivergenceC2F(; kwargs...) = DivergenceC2F(NamedTuple(kwargs))

return_space(::DivergenceC2F, space::Spaces.CenterFiniteDifferenceSpace) =
    Spaces.FaceFiniteDifferenceSpace(space)

return_space(
    ::DivergenceC2F,
    space::Spaces.CenterExtrudedFiniteDifferenceSpace,
) = Spaces.FaceExtrudedFiniteDifferenceSpace(space)

stencil_interior_width(::DivergenceC2F) = ((-half, half),)

function stencil_interior(::DivergenceC2F, loc, idx, arg)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³₊ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    Ju³₋ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    RecursiveApply.rdiv(Ju³₊ ⊟ Ju³₋, local_geometry.J)
end

boundary_width(op::DivergenceC2F, ::SetValue) = 1
function stencil_left_boundary(::DivergenceC2F, bc::SetValue, loc, idx, arg)
    @assert idx == left_face_boundary_idx(arg)
    # ∂x[i] = 2(∂x[i + half] - val)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³₊ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx + half),
        Geometry.LocalGeometry(space, idx + half),
    )
    Ju³ = Geometry.Jcontravariant3(bc.val, local_geometry)
    RecursiveApply.rdiv(Ju³₊ ⊟ Ju³, local_geometry.J / 2)
end
function stencil_right_boundary(::DivergenceC2F, bc::SetValue, loc, idx, arg)
    @assert idx == right_face_boundary_idx(arg)
    space = axes(arg)
    local_geometry = Geometry.LocalGeometry(space, idx)
    Ju³ = Geometry.Jcontravariant3(bc.val, local_geometry)
    Ju³₋ = Geometry.Jcontravariant3(
        getidx(arg, loc, idx - half),
        Geometry.LocalGeometry(space, idx - half),
    )
    RecursiveApply.rdiv(Ju³ ⊟ Ju³₋, local_geometry.J / 2)
end

# left / right SetDivergence boundary conditions
boundary_width(op::DivergenceC2F, ::SetDivergence) = 1
function stencil_left_boundary(
    ::DivergenceC2F,
    bc::SetDivergence,
    loc,
    idx,
    arg,
)
    @assert idx == left_face_boundary_idx(arg)
    # imposed flux boundary condition at left most face
    bc.val
end
function stencil_right_boundary(
    ::DivergenceC2F,
    bc::SetDivergence,
    loc,
    idx,
    arg,
)
    @assert idx == right_face_boundary_idx(arg)
    # imposed flux boundary condition at right most face
    bc.val
end

boundary_width(obj, loc) = 0

boundary_width(bc::Base.Broadcast.Broadcasted{StencilStyle}, loc) =
    has_boundary(bc.f, loc) ? boundary_width(bc.f, get_boundary(bc.f, loc)) : 0

function left_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
    _,
    loc::LeftBoundaryWindow,
)
    space = axes(bc)
    arg_idxs = map(bc.args, stencil_interior_width(bc.f)) do arg, width
        left_interor_window_idx(arg, space, loc) - width[1]
    end
    return max(maximum(arg_idxs), left_idx(space) + boundary_width(bc, loc))
end

function right_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
    _,
    loc::RightBoundaryWindow,
)
    space = axes(bc)
    arg_idxs = map(bc.args, stencil_interior_width(bc.f)) do arg, width
        right_interor_window_idx(arg, space, loc) - width[2]
    end
    return min(minimum(arg_idxs), right_idx(space) - boundary_width(bc, loc))
end

function left_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{CompositeStencilStyle},
    _,
    loc::LeftBoundaryWindow,
)
    space = axes(bc)
    maximum(arg -> left_interor_window_idx(arg, space, loc), bc.args)
end
function right_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{CompositeStencilStyle},
    _,
    loc::RightBoundaryWindow,
)
    space = axes(bc)
    minimum(arg -> right_interor_window_idx(arg, space, loc), bc.args)
end

function left_interor_window_idx(field::Field, _, loc::LeftBoundaryWindow)
    left_idx(axes(field))
end

function right_interor_window_idx(field::Field, _, loc::RightBoundaryWindow)
    right_idx(axes(field))
end

function left_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{Style},
    _,
    loc::LeftBoundaryWindow,
) where {Style <: Fields.AbstractFieldStyle}
    left_idx(axes(bc))
end

function right_interor_window_idx(
    bc::Base.Broadcast.Broadcasted{Style},
    _,
    loc::RightBoundaryWindow,
) where {Style <: Fields.AbstractFieldStyle}
    right_idx(axes(bc))
end

function left_interor_window_idx(_, space, loc::LeftBoundaryWindow)
    left_idx(space)
end

function right_interor_window_idx(_, space, loc::RightBoundaryWindow)
    right_idx(space)
end

function getidx(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
    loc::Interior,
    idx,
)
    stencil_interior(bc.f, loc, idx, bc.args...)
end

function getidx(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
    loc::LeftBoundaryWindow,
    idx,
)
    op = bc.f
    space = axes(bc)
    if idx < left_idx(space) + boundary_width(bc, loc)
        stencil_left_boundary(op, get_boundary(op, loc), loc, idx, bc.args...)
    else
        # fallback to interior stencil
        stencil_interior(op, loc, idx, bc.args...)
    end
end

function getidx(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
    loc::RightBoundaryWindow,
    idx,
)
    op = bc.f
    space = axes(bc)
    if idx > (right_idx(space) - boundary_width(bc, loc))
        stencil_right_boundary(op, get_boundary(op, loc), loc, idx, bc.args...)
    else
        # fallback to interior stencil
        stencil_interior(op, loc, idx, bc.args...)
    end
end

# broadcasting a StencilStyle gives a CompositeStencilStyle
Base.Broadcast.BroadcastStyle(
    ::Type{<:Base.Broadcast.Broadcasted{StencilStyle}},
) = CompositeStencilStyle()

Base.Broadcast.BroadcastStyle(
    ::AbstractStencilStyle,
    ::Fields.AbstractFieldStyle,
) = CompositeStencilStyle()

Base.eltype(bc::Base.Broadcast.Broadcasted{StencilStyle}) =
    return_eltype(bc.f, bc.args...)

function getidx(
    bc::Fields.CenterFiniteDifferenceField,
    ::Location,
    idx::Integer,
)
    Fields.field_values(bc)[idx]
end

function getidx(bc::Fields.FaceFiniteDifferenceField, ::Location, idx::PlusHalf)
    Fields.field_values(bc)[idx.i + 1]
end

getidx(field::Fields.Field, ::Location, idx) = error(
    "Invalid index type $(typeof(idx)) for field on space $(summary(axes(field)))",
)
getidx(field::Base.Broadcast.Broadcasted{StencilStyle}, ::Location, idx) =
    error(
        "Invalid index type $(typeof(idx)) for field on space $(summary(axes(field)))",
    )

getidx(scalar, ::Location, idx) = scalar

# unwap boxed scalars
getidx(scalar::Ref, loc::Location, idx) = getidx(scalar[], loc, idx)

function getidx(bc::Base.Broadcast.Broadcasted, loc::Location, idx)
    args = map(arg -> getidx(arg, loc, idx), bc.args)
    bc.f(args...)
end

function setidx!(bc::Fields.CenterFiniteDifferenceField, idx::Integer, val)
    Fields.field_values(bc)[idx] = val
end

function setidx!(bc::Fields.FaceFiniteDifferenceField, idx::PlusHalf, val)
    Fields.field_values(bc)[idx.i + 1] = val
end

function Base.Broadcast.broadcasted(op::FiniteDifferenceOperator, args...)
    Base.Broadcast.broadcasted(StencilStyle(), op, args...)
end

function Base.Broadcast.broadcasted(
    ::StencilStyle,
    op::FiniteDifferenceOperator,
    args...,
)
    ax = return_space(op, map(axes, args)...)
    Base.Broadcast.Broadcasted{StencilStyle}(op, args, ax)
end

Base.Broadcast.instantiate(bc::Base.Broadcast.Broadcasted{StencilStyle}) = bc
Base.Broadcast._broadcast_getindex_eltype(
    bc::Base.Broadcast.Broadcasted{StencilStyle},
) = eltype(bc)

# check that inferred output field space is equal to dest field space
@inline function Base.Broadcast.materialize!(
    ::Base.Broadcast.BroadcastStyle,
    dest::Fields.Field,
    bc::Base.Broadcast.Broadcasted{Style},
) where {Style <: AbstractStencilStyle}
    dest_space, result_space = axes(dest), axes(bc)
    if result_space !== dest_space
        error(
            "dest space `$(summary(dest_space))` is not equal to the inferred broadcasted result space `$(summary(result_space))`",
        )
    end
    # the default Base behavior is to instantiate a Broadcasted object with the same axes as the dest
    return copyto!(
        dest,
        Base.Broadcast.instantiate(
            Base.Broadcast.Broadcasted{Style}(bc.f, bc.args, dest_space),
        ),
    )
end

function column(
    bc::Base.Broadcast.Broadcasted{Style},
    inds...,
) where {Style <: AbstractStencilStyle}
    _args = map(a -> column(a, inds...), bc.args)
    _axes = column(axes(bc), inds...)
    Base.Broadcast.Broadcasted{Style}(bc.f, _args, _axes)
end

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

function Base.copyto!(
    field_out::Field,
    bc::Base.Broadcast.Broadcasted{S},
) where {S <: AbstractStencilStyle}
    space = axes(bc)
    local_geometry = Spaces.local_geometry_data(space)
    (Ni, Nj, _, _, Nh) = size(local_geometry)
    for h in 1:Nh, j in 1:Nj, i in 1:Ni
        column_field_out = column(field_out, i, j, h)
        column_bc = column(bc, i, j, h)
        apply_stencil!(column_field_out, column_bc)
    end
    return field_out
end

function apply_stencil!(field_out, bc)
    space = axes(bc)
    lbw = LeftBoundaryWindow{Spaces.left_boundary_name(space)}()
    rbw = RightBoundaryWindow{Spaces.right_boundary_name(space)}()
    li = left_idx(space)
    lw = left_interor_window_idx(bc, space, lbw)
    ri = right_idx(space)
    rw = right_interor_window_idx(bc, space, rbw)
    @assert li <= lw <= rw <= ri
    for idx in li:(lw - 1)
        setidx!(field_out, idx, getidx(bc, lbw, idx))
    end
    for idx in lw:rw
        setidx!(field_out, idx, getidx(bc, Interior(), idx))
    end
    for idx in (rw + 1):ri
        setidx!(field_out, idx, getidx(bc, rbw, idx))
    end
    return field_out
end


text:
mutable struct DiscreteIterator
    prob
    u0
    u1

    function DiscreteIterator(prob)
        u0 = copy(prob.u0)
        u1 = similar(prob.u0)
        new(prob, u0, u1)
    end
end

get_integrator(prob::DiscreteProblem) = DiscreteIterator(prob)

DiffEqBase.set_u!(diter::DiscreteIterator, u0) = diter.u0 = u0

function step!(diter::DiscreteIterator, dt::Int = 1, _ignored::Bool = false)
    f = diter.prob.f
    p = diter.prob.p
    u0 = diter.u0
    u1 = diter.u1
    for t in 0:dt-1
        f(u1, u0, p, t)
        u0, u1 = u1, u0
    end
    diter.u0 = u0
    diter.u1 = u1
end
# See:
# - [[../../OrdinaryDiffEq/src/perform_step/fixed_timestep_perform_step.jl::perform_step!.*DiscreteConstantCache]]
# - [[../../OrdinaryDiffEq/src/solve.jl::solve!.*::ODEIntegrator]]
# - http://devdocs.juliadiffeq.org/latest/contributing/diffeq_internals.html

@inline assert_success(integrator::DiscreteIterator) = nothing

const DiscreteLEProblem = LEProblem{DiscreteProblem}

"""
    DiscreteLEProblem(phase_dynamics, u0 [, p];
                      t_attr=<number>, <keyword arguments>)

This is a short-hand notation for:

```julia
LEProblem(DiscreteProblem(phase_dynamics, u0 [, p]), t_attr)
```

For the list of usable keyword arguments, see [`LEProblem`](@ref).
"""
DiscreteLEProblem(phase_dynamics, u0, p=nothing;
                  tspan=(0, 100), kwargs...) =
    LEProblem(DiscreteProblem(phase_dynamics, u0, tspan, p);
              kwargs...)

init_phase_state(integrator::DiscreteIterator) = integrator.u0[:, 1]
init_tangent_state(integrator::DiscreteIterator) = integrator.u0[:, 2:end]

current_state(integrator::DiscreteIterator) = integrator.u0


text:
import Base.-
import Base.getindex
import Base.setindex!
import Base.*
import Base.+
import Base./
import Base.copy
import Base.zero
using Distributions
mutable struct MvGaussianExponential
    n::Tuple{Vector{Float64},Matrix{Float64}} # Expectation Parameter
    gradG::Tuple{Vector{Float64},Matrix{Float64}}
    G::Float64
end
function MvGaussianExponential(data::Matrix)
    Σ = cov(data, dims=2)
    μ = vec(mean(data, dims=2))
    n = (μ, -(Σ + μ*μ'))
    (G, gradG) = compute_G(n[1], n[2])
    return MvGaussianExponential(
        n,
        gradG,
        G)
end
function sufficient_statistic(::Type{MvGaussianExponential}, x::Vector)
    x1 = zeros(1,size(x,1))
    x2 = zeros(1,1,size(x,1))
    for i=1:size(x,2)
        v = x[i]
        x2[1,1,i] = -v*v'
        x1[1,i] = v
    end
    return (x1, x2)
end

function sufficient_statistic(::Type{MvGaussianExponential}, x::AbstractArray{<:Number,2})
    x1 = zeros(size(x,1), size(x,2))
    x2 = zeros(size(x,1), size(x,1), size(x,2))
    for i=1:size(x,2)
        v = vec(x[:,i])
        x2[:,:,i] = -v*v'
        x1[:,i]  = v
    end
    return (x1, x2)
end
function dimension(g::MvGaussianExponential)
    c = length(g.n)
    c= round(Integer,(-1+sqrt(1+4*c))/2)
    return c
end
function get_parameters(v::Vector{Float64}, d::Integer)
     n = v[1:d]
     H = reshape(v[d+1:end],(d,d))
    return (n,H)
end

function compute_G(n::Vector, H::Matrix)
    d = size(n, 1)
    invH = inv(H)
    G =  -0.5 * real(log(Complex(1.0001+ n'*invH*n))) -
        0.5*log(det(-H)+eps()) -
        (d/2)*log(2*pi*2.71)
    #Satistical exponential families:  A digest with flash cards
    #Frank Nielsen,  and Vincent Garcia
    tmp = inv(H+ n*n')
    gradG = (tmp*(-n), -0.5*tmp)
    return (G, gradG)
end
"""
From expectation parameters to natural parameters
"""
function update_parameters!(g::MvGaussianExponential)
    d = dimension(g)
    (G, gradG) = compute_G(g.n[1], g.n[2])
    g.G = G
    g.gradG = gradG
end
function bregman_divergence(g::MvGaussianExponential, x1::Vector)
    #Simplification and hierarchical representations of mixtures    of exponential families.    V. Garcia, F. Nielsen
    n = g.n
    x = get_parameters(x1, dimension(g))
    return  g.G +  dot(x[1]-n[1], g.gradG[1]) + trace((x[2]-n[2])*g.gradG[2]')
end
function trace_axb(a, b)
    return dot(view(a,1,:),view(b,:,1)) + dot(view(a,2,:),view(b,:,2)) + dot(view(a,3,:),view(b,:,3))
end
function bregman_divergence(g::MvGaussianExponential, x1::Matrix, x2::Array{<:AbstractFloat,3})
    n = g.n
    D = size(x1,2)
    X1 = x1.-n[1]
    X2 = (x2.-n[2])
    trace_term = zeros(D)
    for i=1:D
	trace_term[i] = trace_axb(view(X2,:,:,i),g.gradG[2])
    end
    res =  g.G .+ (vec(g.gradG[1]'*X1)) + trace_term
    return res
end

function expectation_parameters(g::Type{MvGaussianExponential}, data::Matrix)
    μ = vec(mean(data, dims=2))
    Σ = cov(data, dims=2)
    return (μ, -(Σ + μ*μ'))
end


text:
#ubound-division.jl

#division on ubounds.

################################################################################
## division

@universal function udiv(a::Ubound, b::Unum)
  aln = is_negative(a.lower)
  ahn = is_negative(a.upper)
  bn = is_negative(b)

  is_zero(b) && return nan(U)

  if (aln != ahn)
    if bn
      return B(resolve_lower(a.upper / b), resolve_upper(a.lower / b))
    else
      return B(resolve_lower(a.lower / b), resolve_upper(a.upper / b))
    end
  end

  outer_result = a.upper / b
  inner_result = a.lower / b

  if is_negative(resolve_lower(outer_result))
    outer_result = resolve_lower(outer_result)
    inner_result = resolve_upper(inner_result)
  else
    outer_result = resolve_upper(outer_result)
    inner_result = resolve_lower(inner_result)
  end

  if is_ulp(inner_result) && is_ulp(outer_result)
    bn ? resolve_as_utype!(outer_result, inner_result) : resolve_as_utype!(inner_result, outer_result)
  else
    bn ? B(outer_result, inner_result) : B(inner_result, outer_result)
  end
end

@universal function udiv(a::Unum, b::Ubound)
  #throw nans for division by zero on either side.
  is_zero(b.lower) && return nan(U)
  is_zero(b.upper) && return nan(U)

  blp = is_pos_def(b.lower)
  bup = is_pos_def(b.upper)
  #if the dividend straddles, then we have nan.
  (blp != bup) && return nan(U)
  #check for zero values in the upper part.
  is_zero(a) && return zero(U)

  if is_inf(a)
    #double infs are no fun.
    is_inf(b.lower) && return nan(U)
    is_inf(b.upper) && return nan(U)

    #match the inf with the sign of the product.
    return inf(U, @signof(a) $ @signof(b.lower))
  end

  (b_inner, b_outer) = blp ? (b.lower, b.upper) : (b.upper, b.lower)

  if (is_negative(a) != blp)
    #result should resolve to a positive number.
    lower_result = resolve_lower(a / b_outer)
    upper_result = resolve_upper(a / b_inner)
    (is_ulp(lower_result) && is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  else
    #result should resolve to a negative number.
    lower_result = resolve_lower(a / b_inner)
    upper_result = resolve_upper(a / b_outer)
    (is_ulp(lower_result) && is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  end
end

@universal function udiv(a::Ubound, b::Ubound)
  signcode::UInt16 = 0

  #check some divisions by exact zero, which gives exact infinite bounds.
  is_zero(b.lower) && return nan(U)
  is_zero(b.upper) && return nan(U)

  is_negative(a.lower) && (signcode += 1)
  is_negative(a.upper) && (signcode += 2)
  is_negative(b.lower) && (signcode += 4)
  is_negative(b.upper) && (signcode += 8)

  if (signcode == 0) #everything is positive
    lower_result = resolve_lower(a.lower / b.upper)
    upper_result = resolve_upper(a.upper / b.lower)

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 1) #only a.lowbound is negative
    B(resolve_lower(a.lower / b.lower), resolve_upper(a.upper / b.lower))
  #signcode 2 is not possible
  elseif (signcode == 3) #a is negative and b is positive

    lower_result = resolve_lower(a.lower / b.lower)
    upper_result = resolve_upper(a.upper / b.upper)

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 4) #only b.lowbound is negative
    #b straddles zero so we'll output NaN
    return nan(U)
  elseif (signcode == 5) #a.lowbound and b.lowbound are negative
    #b straddles zero so we'll output NaN
    return nan(U)
  #signcode 6 is not possible
  elseif (signcode == 7) #only b.highbound is positive
    #b straddles zero so we'll output NaN
    return nan(U)
  #signcode 8, 9, 10, 11 are not possible
  elseif (signcode == 12) #b is negative, a is positive
    lower_result = resolve_lower(a.upper / b.lower)
    upper_result = resolve_upper(a.lower / b.upper)

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  elseif (signcode == 13) #b is negative, a straddles
    B(resolve_lower(a.upper / b.upper), resolve_upper(a.lower / b.upper))
  #signcode 14 is not possible
  elseif (signcode == 15) #everything is negative
    lower_result = resolve_lower(a.upper / b.lower)
    upper_result = resolve_upper(a.lower / b.upper)

    (is_ulp(lower_result) & is_ulp(upper_result)) ? resolve_as_utype!(lower_result, upper_result) : B(lower_result, upper_result)
  else
    throw(ArgumentError("error dividing ubounds $a and $b, throws invalid signcode $signcode."))
  end
end

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
""
function constraint_model_voltage_active(pm::_PM.AbstractPowerModel; nw::Int=nw_id_default, cnd::Int=pm.ccnd)
    constraint_model_voltage_active(pm, nw, cnd)
end

""
function constraint_generation_active(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    gen = _PM.ref(pm, nw, :gen, i)

    _PM.constraint_gen_power_on_off(pm, nw, i, gen["pmin"], gen["pmax"], gen["qmin"], gen["qmax"])
    _PMR.constraint_gen_bus_connection(pm, nw, i, gen["gen_bus"])

end

""
function constraint_load_active(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    load = _PM.ref(pm, nw, :load, i)
    _PMR.constraint_load_bus_connection(pm, nw, i, load["load_bus"])
end

""
function constraint_shunt_active(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    shunt = _PM.ref(pm, nw, :shunt, i)
    _PMR.constraint_shunt_bus_connection(pm, nw, i, shunt["shunt_bus"])
end

""
function constraint_branch_active(pm::_PM.AbstractPowerModel, i::Int; nw::Int=nw_id_default)
    branch = _PM.ref(pm, nw, :branch, i)

    bus_fr = haskey(_PM.ref(pm, nw, :bus), branch["f_bus"])
    bus_to = haskey(_PM.ref(pm, nw, :bus), branch["t_bus"])

    _PMR.constraint_branch_bus_connection(pm, nw, i, branch["f_bus"])
    _PMR.constraint_branch_bus_connection(pm, nw, i, branch["t_bus"])
end


text:
using Revise, BinaryTraits

# See https://docs.julialang.org/en/v1/manual/interfaces

# -----------------------------------------------------------------------------
# Iteration interface
# -----------------------------------------------------------------------------
import Base: iterate
@trait Iterable prefix Is,Not
@implement IsIterable by iterate()::Any
@implement IsIterable by iterate(state::Any)::Any

# Example from https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration-1
struct Squares
    count::Int
end

Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1)

# Let's assign the Squares type to Iterable
@assign Squares with IsIterable
@check(Squares)
#=
julia> @check(Squares)
✅ Squares has implemented:
1. IterableTrait: IsIterable ⇢ iterate(::<Type>)::Any
2. IterableTrait: IsIterable ⇢ iterate(::<Type>, ::Any)::Any
=#

# -----------------------------------------------------------------------------
# Indexing interface
# -----------------------------------------------------------------------------
import Base: getindex, setindex!, firstindex, lastindex

@trait Indexable prefix Is,Not
@implement IsIndexable by getindex(i::Any)

@trait IndexableFromBeginning prefix Is,Not
@implement IsIndexableFromBeginning by firstindex()

@trait IndexableAtTheEnd prefix Is,Not
@implement IsIndexableAtTheEnd by lastindex()

# Make sure that `i` is untyped (i.e. `Any`) to adhere to the contract
function Base.getindex(S::Squares, i)
    1 <= i <= S.count || throw(BoundsError(S, i))
    return i*i
end

@assign Squares with IsIndexable
@check(Squares)
#=
✅ Squares has implemented:
1. IterableTrait: IsIterable ⇢ iterate(::<Type>)::Any
2. IterableTrait: IsIterable ⇢ iterate(::<Type>, ::Any)::Any
3. IndexableTrait: IsIndexable ⇢ getindex(::<Type>, ::Any)::Any
=#

# We want to have the traits for indexing from beginning and at the end
@assign Squares with IsIndexableFromBeginning, IsIndexableAtTheEnd
@check(Squares)
#=
julia> @check(Squares)
┌ Warning: Missing implementation: IndexableAtTheEndTrait: IsIndexableAtTheEnd ⇢ lastindex(::Squares)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:200
┌ Warning: Missing implementation: IndexableFromBeginningTrait: IsIndexableFromBeginning ⇢ firstindex(::Squares)::Any
└ @ BinaryTraits ~/.julia/dev/BinaryTraits/src/interface.jl:200
✅ Squares has implemented:
1. IterableTrait: IsIterable ⇢ iterate(::<Type>)::Any
2. IterableTrait: IsIterable ⇢ iterate(::<Type>, ::Any)::Any
3. IndexableTrait: IsIndexable ⇢ getindex(::<Type>, ::Any)::Any
❌ Squares is missing these implementations:
1. IndexableAtTheEndTrait: IsIndexableAtTheEnd ⇢ lastindex(::<Type>)::Any
2. IndexableFromBeginningTrait: IsIndexableFromBeginning ⇢ firstindex(::<Type>)::Any
=#

# Let's implement them now.
Base.firstindex(S::Squares) = 1
Base.lastindex(S::Squares) = length(S)
@check(Squares)
#=
julia> @check(Squares)
✅ Squares has implemented:
1. IterableTrait: IsIterable ⇢ iterate(::<Type>)::Any
2. IterableTrait: IsIterable ⇢ iterate(::<Type>, ::Any)::Any
3. IndexableAtTheEndTrait: IsIndexableAtTheEnd ⇢ lastindex(::<Type>)::Any
4. IndexableFromBeginningTrait: IsIndexableFromBeginning ⇢ firstindex(::<Type>)::Any
5. IndexableTrait: IsIndexable ⇢ getindex(::<Type>, ::Any)::Any
=#


text:
using MCMCDiagnosticTools

using Distributions
using MLJBase
using MLJLIBSVMInterface
using MLJXGBoostInterface

using Test

const xgboost_deterministic = Pipeline(XGBoostClassifier(); operation=predict_mode)

@testset "rstar.jl" begin
    classifiers = (XGBoostClassifier(), xgboost_deterministic, SVC())
    N = 1_000

    @testset "examples (classifier = $classifier)" for classifier in classifiers
        # Compute R⋆ statistic for a mixed chain.
        samples = randn(N, 2)
        dist = rstar(classifier, randn(N, 2), rand(1:3, N))

        # Mean of the statistic should be focused around 1, i.e., the classifier does not
        # perform better than random guessing.
        if classifier isa MLJBase.Deterministic
            @test dist isa Float64
        else
            @test dist isa LocationScale
            @test dist.ρ isa PoissonBinomial
            @test minimum(dist) == 0
            @test maximum(dist) == 3
        end
        @test mean(dist) ≈ 1 rtol = 0.2

        # Compute R⋆ statistic for a mixed chain.
        samples = randn(4 * N, 8)
        chain_indices = repeat(1:4, N)
        dist = rstar(classifier, samples, chain_indices)

        # Mean of the statistic should be closte to 1, i.e., the classifier does not perform
        # better than random guessing.
        if classifier isa MLJBase.Deterministic
            @test dist isa Float64
        else
            @test dist isa LocationScale
            @test dist.ρ isa PoissonBinomial
            @test minimum(dist) == 0
            @test maximum(dist) == 4
        end
        @test mean(dist) ≈ 1 rtol = 0.15

        # Compute the R⋆ statistic for a non-mixed chain.
        samples = [
            sin.(1:N) cos.(1:N)
            100 .* cos.(1:N) 100 .* sin.(1:N)
        ]
        chain_indices = repeat(1:2; inner=N)
        dist = rstar(classifier, samples, chain_indices)

        # Mean of the statistic should be close to 2, i.e., the classifier should be able to
        # learn an almost perfect decision boundary between chains.
        if classifier isa MLJBase.Deterministic
            @test dist isa Float64
        else
            @test dist isa LocationScale
            @test dist.ρ isa PoissonBinomial
            @test minimum(dist) == 0
            @test maximum(dist) == 2
        end
        @test mean(dist) ≈ 2 rtol = 0.15
    end

    @testset "exceptions (classifier = $classifier)" for classifier in classifiers
        @test_throws DimensionMismatch rstar(classifier, randn(N - 1, 2), rand(1:3, N))
        for subset in (-0.3, 0, 1 / (3 * N), 1 - 1 / (3 * N), 1, 1.9)
            @test_throws ArgumentError rstar(
                classifier, randn(N, 2), rand(1:3, N); subset=subset
            )
        end
    end
end


text:
module Huffman

include("binarytree.jl")
include("paths.jl")

export HuffmanTree, allpaths

end


text:
const _default_light = Vec3f0[Vec3f0(1.0,1.0,1.0), Vec3f0(0.1,0.1,0.1), Vec3f0(0.9,0.9,0.9), Vec3f0(20,20,20)]
function default(main, s, data)
    data = _default(main, s, copy(data))
    @gen_defaults! data begin # make sure every object has these!
        model      	     = eye(Mat4f0)
        light      	     = _default_light
        preferred_camera = :perspective
        is_transparent_pass = Cint(false)
    end
end

"""
Creates a default visualization for any value.
The defaults can be customized via the key word arguments and the style parameter.
The style can change the the look completely (e.g points displayed as lines, or particles),
while the key word arguments just alter the parameters of one visualization.
Always returns a context, which can be displayed on a window via view(::Context, [display]).
"""
visualize(main, s::Symbol=:default; kw_args...) = visualize(main, Style{s}(), Dict{Symbol, Any}(kw_args))::Context
visualize(main, s::Style, data::Dict) = assemble_shader(default(main, s, data))::Context
visualize(c::Composable, s::Symbol=:default; kw_args...) = Context(c)
visualize(c::Context, s::Symbol=:default; kw_args...) = c

function Base.push!{Pre}(screen::Screen, robj::RenderObject{Pre})
    # find renderlist specialized to current prerender function
    index = findfirst(screen.renderlist) do renderlist
        prerendertype(eltype(renderlist)) == Pre
    end
    if index == 0
        # add new specialised renderlist, if none found
        screen.renderlist = (screen.renderlist..., RenderObject{Pre}[])
        index = length(screen.renderlist)
    end
    # only add to renderlist if not already in there
    in(robj, screen.renderlist[index]) || push!(screen.renderlist[index], robj)
    nothing
end

function _view(
		robj::RenderObject, screen=current_screen();
		camera = robj.uniforms[:preferred_camera],
		position = Vec3f0(2), lookat=Vec3f0(0)
	)
    if isa(camera, Camera)
    	real_camera = camera
    elseif haskey(screen.cameras, camera)
        real_camera = screen.cameras[camera]
    elseif camera == :perspective
        inside = screen.inputs[:mouseinside]
		real_camera = PerspectiveCamera(screen.inputs, position, lookat, keep=inside)
	elseif camera == :fixed_pixel
		real_camera = DummyCamera(window_size=screen.area)
	elseif camera == :orthographic_pixel
        inside = screen.inputs[:mouseinside]
        real_camera = OrthographicPixelCamera(screen.inputs, keep=inside)
	elseif camera == :nothing
        push!(screen, robj, :nothing)
		return nothing
	else
         error("Method $camera not a known camera type")
	end
    camsym = Symbol(string(camera))
    screen.cameras[camsym] = real_camera
	merge!(robj.uniforms, collect(real_camera), Dict( # add display dependant values
		:resolution => get!(screen.inputs, :resolution, const_lift(Vec2f0, const_lift(x->Vec2f0(x.w,x.h), screen.area))),
		:fixed_projectionview => get!(screen.cameras, :fixed_pixel, DummyCamera(window_size=screen.area)).projectionview
	))
    push!(screen, robj)
	nothing
end

_view(robjs::Vector, screen=current_screen(); kw_args...) = for robj in robjs
	_view(robj, screen; kw_args...)
end
_view(c::Composable, screen=current_screen(); kw_args...) = _view(extract_renderable(c), screen; kw_args...)


text:
using MyExample
using Test

@testset "MyExample.jl" begin
    @test my_f(2, 1) == 7
    @test my_f(2, 3) == 13
    @test my_f(1, 3) == 11
end

@testset "Derivative" begin
    @test derivative_of_my_f(2, 3) == 2
    @test derivative_of_my_f(1, 5) == 2
end


text:
abstract Queryable{T,Provider}

abstract QueryProvider


text:
using Documenter, HTTP

makedocs(
    modules = [HTTP],
    sitename = "HTTP.jl",
    pages = [
             "Home" => "index.md",
             "public_interface.md",
             "internal_architecture.md",
             "internal_interface.md",
             ],
)

deploydocs(
    repo = "github.com/JuliaWeb/HTTP.jl.git",
)


text:
using Revise, Parameters, Logging, Test
using CUDA
includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")
using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..BasicStructs
using Shuffle,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates
using ..MainOverlap, ..RandIndex , ..ProbabilisticMetrics , ..VolumeMetric ,..InformationTheorhetic
using ..CUDAAtomicUtils, ..TpfpfnKernel, ..InterClassCorrKernel

dimx = 529
dimy = 556
dimz = 339
arrGold = CUDA.ones(UInt8, (dimx,dimy,dimz) )
arrAlgo = CUDA.ones(UInt8, (dimx,dimy,dimz) )
numberToLooFor = UInt8(1)
argsMain, threadsMain,  blocksMain,threadsMean,blocksMean,argsMean, totalNumbOfVoxels=InterClassCorrKernel.prepareInterClassCorrKernel(arrGold ,arrAlgo,numberToLooFor)

globalICC=calculateInterclassCorr(arrGold,arrAlgo,argsMain, threadsMain,  blocksMain,threadsMean,blocksMean,argsMean, totalNumbOfVoxels)::Float64

Int64(argsMain[1][1]-dimx*dimy*dimz)
argsMain[2][1]==dimx*dimy*dimz


text:
"""
    deleteAnomalies!(df::DataFrame, ColOverflow::String, ColPct::String, threshold::Float64)

Deletes rows with precipitations under a certain threshold but with overflow.

Data provided must not contain any missing values for overflows' duration or occurrence.

# Arguments

- `df::DataFrame`: The dataframe containing the data.
- `ColOverflow::String`: The name of the column of `df` containing the binary vector of overflows.
- `ColPct::String`: The name of the column of `df` that allows to know if it has rained or not during the time division used.
- `threshold::Float64`: The precipitation threshold under which an anomalie is considered.
"""
function deleteAnomalies!(df::DataFrame, ColOverflow::String, ColPct::String, threshold::Float64)
    filter!([Symbol(ColOverflow), Symbol(ColPct)] => (x, y) -> (x == 1 && y > threshold) || (x == 0), df)
end

max_stars_repo_name: AlexLsn/CSO

max_issues_repo_name: AlexLsn/CSO

max_forks_repo_name: AlexLsn/CSO


text:
# Test the embedding table backprop against Zygote
# @testset "Testing Embedding Backprop" begin
#     weights_size = (20, 10)
#     indices_size = 5
#
#     bias_size = (first(weights_size), indices_size)
#     bias = randn(Float32, bias_size)
#
#     # Construct a function returning a scalar
#     f = (indices, weights) -> sum(bias .* embedding(indices, weights))
#
#     # Make the weights and indices
#     #weights = randn(Float32, weights_size)
#     weights = randn(Float32, weights_size)
#
#     # Zygote has a bug with how it deals with repeated elements - here we just make sure
#     # all are unique to get around it.
#     indices = shuffle(Int32(1):Int32(last(weights_size)))[1:indices_size]
#     backend = nGraph.Backend("CPU")
#
#     # Make sure embedding lookup works
#     F = nGraph.compile(backend, embedding, indices, weights)
#     @test embedding(indices, weights) == read(F())
#
#     # Make sure round trip of "f" works
#     F = nGraph.compile(backend, f, indices, weights)
#     @test isapprox(f(indices, weights), read(F())[])
#
#     # Now for the tricky bit - check if the gradients are the same.
#     #
#     # First - we take the Zygote gradient for the weights - which will be the second return
#     # value from Zygote.gradient
#     zygote_gradient = Zygote.gradient(f, indices, weights)[2]
#
#     # Next - we compile the nGraph function with the `Gradient` optimizer
#     #
#     # For now, we have to wrap `weights` in a `param` so the nGraph converter will actually
#     # take a gradient with respect to it.
#     P = Flux.param(weights)
#     F = nGraph.compile(backend, x -> f(x, P), indices; optimizer = nGraph.Gradient)
#
#     display(zygote_gradient)
#     println()
#
#     # Run the function - extract the gradient from F
#     F()
#     ngraph_gradient = read(F.optimizer._id[P][2]::nGraph.Tensor)
#     display(ngraph_gradient)
#     println()
#
#     @test isapprox(zygote_gradient, ngraph_gradient)
#
#     #####
#     ##### Test 3 dimensional stuff
#     #####
#
#     # weights = randn(Float32, weights_size)
#
#     # index_size = 2
#     # indices_total_size = index_size^2
#     # indices = shuffle(Int32(1):Int32(last(weights_size)))[1:indices_total_size]
#     # indices = reshape(indices, index_size, index_size)
#
#     # bias = randn(Float32, (first(weights_size), index_size, index_size))
#
#     # f = (indices, weights) -> sum(bias .+ embedding(indices, weights))
#
#     # zygote_gradient = Zygote.gradient(f, indices, weights)[2]
#     # @info "Zygote Gradient"
#     # println()
#     # display(zygote_gradient)
#     # println()
# end


text:
using LyceumDevTools
using Test

@testset "LyceumDevTools.jl" begin
    # Write your own tests here.
end


text:
# Fallbacks!
NLᵁ(m, i) = 0
NLⱽ(m, i) = 0
NLᵀ(m, i) = 0
NLˢ(m, i) = 0
NLᵉ(m, i) = 0

∂z_NLᵁ(m, i) = 0
∂z_NLⱽ(m, i) = 0
∂z_NLᵀ(m, i) = 0
∂z_NLˢ(m, i) = 0
∂z_NLᵉ(m, i) = 0

#####
##### Counter gradient flux model proposed by Large et al (1994)
#####

Base.@kwdef struct CounterGradientFlux{T} <: AbstractParameters
    Cᴺ :: T = 1.0 # Mass flux proportionality constant
end

const CGModel = Model{L, K, W, <:CounterGradientFlux} where {L, K, W}

#∂z_NLᵀ(m::CGModel, i) where K = ∂NL∂z(m.nonlocal_flux.Cᴺ, m.state.Qθ, m)
#∂z_NLˢ(m::CGModel, i) where K = ∂NL∂z(m.nonlocal_flux.Cᴺ, m.state.Qs, m)

#####
##### Diagnostic plume model
#####

abstract type AbstractDiagnosticPlumeModel <: AbstractParameters end

Base.@kwdef struct WitekDiagnosticPlumeModel{T} <: AbstractDiagnosticPlumeModel
     Ca :: T = 0.1
    Cbw :: T = 2.86
     Ce :: T = 0.1
    Cew :: T = 0.572
     CQ :: T = 1.0
end

const ModelWithPlumes = Model{L, K, W, <:AbstractDiagnosticPlumeModel} where {L, K, W}

instantiate_plume(::AbstractDiagnosticPlumeModel, grid) =
    (T=CellField(grid), S=CellField(grid), W²=FaceField(grid), e=nothing)

#####
##### Entrainment
#####

"""
    dynamic_entrainment(Ce, ΔB̆, W̆²)

Returns

``C^e * min(0, ΔB̆ / W̆²) ``,

which is always less than zero.
"""
@inline dynamic_entrainment(Ce, ΔB̆, W̆²::T) where T =
    ifelse(W̆²==0, zero(T), Ce * min(zero(T), ΔB̆ / W̆²))

@inline function tracer_entrainment(m::ModelWithPlumes, i)
    ΔB̆ = plume_buoyancy_excess(i, m.grid, m)
    W̆² = oncell(m.state.plume.W², i)
    return dynamic_entrainment(m.nonlocal_flux.Ce, ΔB̆, W̆²)
end

@inline function momentum_entrainment(m::ModelWithPlumes, i)
    ΔB̆ = onface(i, m.grid, plume_buoyancy_excess, m)
    W̆² = @inbounds m.state.plume.W²[i]
    return dynamic_entrainment(m.nonlocal_flux.Ce, ΔB̆, W̆²)
end

#####
##### Plume boundary conditions
#####

function set_tracer_plume_bc!(ϕ̆, Φ, Qϕ, CQ, model)
    n = ϕ̆.grid.N

    v★ = maxsqrt(model.solution.e, n) # turbulent velocity scale at the surface

    # Surface layer model: √e Δϕ̆ = - C Qϕ, where Δϕ̆ is plume excess.
    # Also, plumes can't exist if v★=0. This makes no sense, but our best option atm.
    @inbounds ϕ̆[n] = ifelse(v★==0, Φ[n], Φ[n] - CQ * Qϕ / v★)

    return nothing
end

function set_vertical_momentum_plume_bc!(W)
    @inbounds W[W.grid.N+1] = 0
    return nothing
end

#####
##### Plume equations
#####

#####
##### Buoyancy
#####

@inline plume_buoyancy_excess(i, grid, T̆, S̆, T, S, α, β, g) =
    @inbounds g * (α*(T̆[i] - T[i]) - β*(S̆[i] - S[i]))

@inline plume_buoyancy_excess(i, grid, model) =
    plume_buoyancy_excess(i, grid, model.state.plume.T, model.state.plume.S,
                                   model.solution.T, model.solution.S,
                                   model.constants.α, model.constants.β, model.constants.g)

function clip_infinite!(ϕ)
    for i in eachindex(ϕ)
        @inbounds begin
            ϕᵢ = ϕ[i]
            if !(isfinite(ϕᵢ))
                 ϕ[i] = 0
            end
        end
    end
    return nothing
end

function clip_positive!(ϕ)
    for i in eachindex(ϕ)
        @inbounds begin
            ϕᵢ = ϕ[i]
            if ϕᵢ > 0
                ϕ[i] = 0
            end
        end
    end
    return nothing
end

function clip_negative!(ϕ)
    for i in eachindex(ϕ)
        @inbounds begin
            ϕᵢ = ϕ[i]
            if ϕᵢ < 0
                ϕ[i] = 0
            end
        end
    end
    return nothing
end

function update_nonlocal_flux!(model::ModelWithPlumes) where K

    set_tracer_plume_bc!(model.state.plume.T, model.solution.T,
                         model.state.Qθ, model.nonlocal_flux.CQ, model)

    set_tracer_plume_bc!(model.state.plume.S, model.solution.S,
                         model.state.Qs, model.nonlocal_flux.CQ, model)

    set_vertical_momentum_plume_bc!(model.state.plume.W²)

    integrate_plume_equations!(model.state.plume.T, model.state.plume.S, model.state.plume.W²,
                               model.solution.T, model.solution.S, model.grid, model)

    clip_infinite!(model.state.plume.W²)
    clip_negative!(model.state.plume.W²)
    clip_infinite!(model.state.plume.T)

    return nothing
end

function integrate_plume_equations!(T̆, S̆, W̆², T, S, grid, model)

    n = grid.N

    # Vertical momentum at the nᵗʰ cell interface, approximating excess buoyancy at
    # interface with excess at top cell center:
    ΔB̆ᵢ₊₁ = plume_buoyancy_excess(n, grid, model)
    @inbounds W̆²[n] = -model.nonlocal_flux.Cbw * Δf(grid, n) * ΔB̆ᵢ₊₁

    # Integrate from surface cell `N-1` downwards
    for i in n-1 : -1 : 1
        if W̆²[i+1] <= 0 # plume is stopped
            @inbounds W̆²[i] = 0
            @inbounds T̆[i] = T[i]
            @inbounds S̆[i] = S[i]
        else # plume still lives

            # Integrate temperature and salinity
            @inbounds T̆[i] = T̆[i+1] + Δc(grid, i+1) * tracer_entrainment(model, i) * (T̆[i+1] - T[i+1])
            @inbounds S̆[i] = S̆[i+1] + Δc(grid, i+1) * tracer_entrainment(model, i) * (S̆[i+1] - S[i+1])

            # Integrate vertical momentum
            ΔB̆ᵢ₊₁ = onface(i+1, grid, plume_buoyancy_excess, model)

            @inbounds W̆²[i] = W̆²[i+1] - Δf(grid, i+1) * (
                                model.nonlocal_flux.Cbw * ΔB̆ᵢ₊₁
                              - model.nonlocal_flux.Cew * momentum_entrainment(model, i) * W̆²[i+1])
        end
    end

    return nothing
end

#####
##### Mass flux contribution to tracer budgets
#####

maxzero(ϕ::T) where T = max(zero(T), ϕ)

@inline M(m, i) = @inbounds -m.nonlocal_flux.Ca * maxsqrt(m.state.plume.W², i)

"Returns the non-local flux of temperature at cell centers."
@inline function NLᵀ(m::ModelWithPlumes, i)
    Mᵢ = oncell(M, m, i)
    ΔTᵢ = @inbounds m.state.plume.T[i] - m.solution.T[i]
    return Mᵢ * ΔTᵢ
end

"Returns the non-local flux of salinity at cell centers."
@inline function NLˢ(m::ModelWithPlumes, i)
    Mᵢ = oncell(M, m, i)
    ΔSᵢ = @inbounds m.state.plume.S[i] - m.solution.S[i]
    return Mᵢ * ΔSᵢ
end

"Returns the temperature transport for a plume model, with mass flux on cell interfaces."
@inline ∂z_NLᵀ(m::ModelWithPlumes, i) = 1 / (2 * Δc(m.grid, i+1)) * (NLᵀ(m, i+1) - NLᵀ(m, i))

"Returns the salinity transport for a plume model, with mass flux on cell interfaces."
@inline ∂z_NLˢ(m::ModelWithPlumes, i) = 1 / (2 * Δc(m.grid, i+1)) * (NLˢ(m, i+1) - NLˢ(m, i))


text:
# Generate voicings exercises.

import JSON

# This ID is incremented for every problem generated.  Barring changes
# to the generation, this should generate the same ids on subsequent
# runs.
cur_id = 1

function emit(question, answer)
    global cur_id
    println("INSERT INTO \"probs\" VALUES($cur_id, '$question', '$answer');")
    cur_id += 1
end

# The ii-V-I progressions through the cycle.  The DeGregg puts some
# work into keeping these in range.  However, since we will be mostly
# doing these out of order, the exact order here doesn't matter all
# that much, and the simple Circle of V is fine.  F♯ and G♭ and both
# printed (making for 13 of each voicing), since both are common.
major_iivi = [
    ["Dm7", "G7", "CΔ"],
    ["Am7", "D7", "GΔ"],
    ["Em7", "A7", "DΔ"],
    ["Bm7", "E7", "AΔ"],
    ["F♯m7", "B7", "EΔ"],
    ["C♯m7", "F♯7", "BΔ"],
    ["G♯m7", "C♯7", "F♯Δ"],
    ["A♭m7", "D♭7", "G♭Δ"],
    ["E♭m7", "A♭7", "D♭Δ"],
    ["B♭m7", "E♭7", "A♭Δ"],
    ["Fm7", "B♭7", "E♭Δ"],
    ["Cm7", "F7", "B♭Δ"],
    ["Gm7", "C7", "FΔ"] ]

# Diminished tri-tone sub.  The idea of this pattern is to alternate
# the inversion of the chords so that the pattern only moves slightly
# on the piano.  See p39 in the deGreg book for an example.
min_tritone_sub = [
    ["CΔ", "C♯°7", "Dm7", "D♯°7", "Em7", "E♭7", "Dm7", "D♭7", "CΔ"],
    ["FΔ", "F♯°7", "Gm7", "G♯°7", "Am7", "A♭7", "Gm7", "G♭7", "FΔ"],
    ["B♭Δ", "B°7", "Cm7", "C♯°7", "Dm7", "D♭7", "Cm7", "B7", "B♭Δ"],
    ["E♭Δ", "E°7", "Fm7", "F♯°7", "Gm7", "G♭7", "Fm7", "E7", "E♭Δ"],

    ["A♭Δ", "A°7", "B♭m7", "B°7", "Cm7", "B7", "B♭m7", "A7", "A♭Δ"],
    ["D♭Δ", "D°7", "E♭m7", "E°7", "Fm7", "E7", "E♭m7", "D7", "D♭Δ"],
    ["G♭Δ", "G°7", "A♭m7", "A°7", "B♭m7", "A7", "A♭m7", "G7", "G♭Δ"],
    ["BΔ", "C°7", "C♯m7", "D°7", "D♯m7", "D7", "C♯m7", "C7", "BΔ"],

    ["EΔ", "F°7", "F♯m7", "G°7", "G♯m7", "G7", "F♯m7", "F7", "EΔ"],
    ["AΔ", "A♯°7", "Bm7", "C°7", "C♯m7", "C7", "Bm7", "B♭7", "AΔ"],
    ["DΔ", "D♯°7", "Em7", "F°7", "F♯m7", "F7", "Em7", "E♭7", "DΔ"],
    ["GΔ", "G♯°7", "Am7", "A♯°7", "Bm7", "B♭7", "Am7", "A♭7", "GΔ"] ]

# Extract the root of a chord, the characters before any modifier.
function root_name(chord)
    m = match(r"^[A-G][♭♯]?", chord)
    if m === nothing
        error("Invalid chord name")
    else
        m.match
    end
end

base_notes = Dict("C" => 36, "C♯" => 37,
    "D♭" => 37, "D" => 38, "D♯" => 39,
    "E♭" => 39, "E" => 40, "E♯" => 41,
    "F♭" => 40, "F" => 41, "F♯" => 42,
    "G♭" => 42, "G" => 43, "G♯" => 44,
    "A♭" => 44, "A" => 45, "A♯" => 46,
    "B♭" => 46, "B" => 47, "B♯" => 48)

# A Chord Voicing describes how to play a particular chord.
type ChordVoicing
    # Indicates movement of the root from the previous chord.
    root :: Int8
    # Intervals beyond the "base" to be played by the left hand.
    # e.g., starting with a '0' would indicate that the root should be
    # played.  All intervals are in half steps.
    left :: Vector{Int8}
    # Intervals to be played by the right hand, starting with the last
    # note in the left hand.
    right :: Vector{Int8}
end

function make_notes(base, v :: ChordVoicing, first=false)
    notes = Vector{Int8}()
    cur = base
    for lh in v.left
        cur += lh
        push!(notes, cur)
    end
    rh_start = cur + v.right[1]
    for rh in v.right
        cur += rh
        push!(notes, cur)
    end
    # TODO: Adjust octave.
    while first && rh_start < 54
        rh_start += 12
        notes .+= 12
        base += 12
    end
    (base, notes)
end

# A voicing is then a full pattern, for multiple chords.
type Voicing
    name :: String
    chords :: Vector{ChordVoicing}
end

function make_notes(chords, v :: Voicing)
    base = base_notes[root_name(chords)]
    chords = Vector{Vector{Int8}}()
    first = true
    for ch in v.chords
        base += ch.root
        (base, notes) = make_notes(base, ch, first)
        first = false
        push!(chords, notes)
    end
    Dict("type" => "voicing",
        "chords" => chords)
end

# Generate patterns for a particular voicing.  The name will be
# printed before, and the hint afterwards.  This should be enough for
# the reader to understand how to voice the chord.  The chord start
# can be used to skip part of the progression.
function gen_cycle(exercise, progression, voicing, hint; chord_start=1)
    name = voicing.name
    for prog in progression
        chord = join(prog[chord_start:end], "-")
        notes = make_notes(prog[1], voicing)
        notes["chords"] = notes["chords"][chord_start:end]
        emit("$name $chord $hint ($exercise)", JSON.json(notes))
    end
end

# The "shell" voicing, with a 1-7 played first in the left hand.
# However, we'll encode the 7 in the right hand, so we can
# automatically find the right octave.
shell17_voicing = let
    minor = ChordVoicing(0, [0], [10])
    dominant = ChordVoicing(5, [0], [4])
    major = ChordVoicing(-7, [0], [11])
    Voicing("shell", [minor, dominant, major])
end

shell13_voicing = let
    minor = ChordVoicing(0, [0], [3])
    dominant = ChordVoicing(-7, [0], [10])
    major = ChordVoicing(5, [0], [4])
    Voicing("shell", [minor, dominant, major])
end

guide173_voicing = let
    minor = ChordVoicing(0, [0], [10, 5])
    dominant = ChordVoicing(-7, [0], [16, 6])
    major = ChordVoicing(5, [0], [11, 5])
    Voicing("guide", [minor, dominant, major])
end

guide137_voicing = let
    minor = ChordVoicing(0, [0], [15, 7])
    dominant = ChordVoicing(5, [0], [10, 6])
    major = ChordVoicing(-7, [0], [16, 7])
    Voicing("guide", [minor, dominant, major])
end

guide73_voicing = let
    minor = ChordVoicing(0, [], [10, 5])
    dominant = ChordVoicing(-7, [], [16, 6])
    major = ChordVoicing(5, [], [11, 5])
    Voicing("guide", [minor, dominant, major])
end

guide37_voicing = let
    minor = ChordVoicing(0, [], [15, 7])
    dominant = ChordVoicing(5, [], [10, 6])
    major = ChordVoicing(-7, [], [16, 7])
    Voicing("guide", [minor, dominant, major])
end

function generate()
    println("DELETE FROM probs;")

    # Exercise 1-1A, ii-V-I, voiced with a 1-7 shell first.
    gen_cycle("1-1A", major_iivi, shell17_voicing, "1-7")

    # Exercise 1-1B, ii-V-I, voiced with a 1-3 shell first.
    gen_cycle("1-1B", major_iivi, shell13_voicing, "1-3")

    # Exercise 1-1A, ii-V-I, voiced with a 1-7 shell first.
    gen_cycle("1-1A", major_iivi, shell13_voicing, "1-7", chord_start=2)

    # Exercise 1-1B, ii-V-I, voiced with a 1-3 shell first.
    gen_cycle("1-1B", major_iivi, shell17_voicing, "1-3", chord_start=2)

    # Exercise 1-1A, ii-V-I, voiced with a 1-7 shell first.
    gen_cycle("1-1A", major_iivi, shell17_voicing, "1-7", chord_start=3)

    # Exercise 1-1B, ii-V-I, voiced with a 1-3 shell first.
    gen_cycle("1-1B", major_iivi, shell13_voicing, "1-3", chord_start=3)

    # Exercise 2-1A, ii-V-I, guide tones, 1-7-3 first.
    gen_cycle("2-1A", major_iivi, guide173_voicing, "1-7-3")

    # Exercise 2-1B, ii-V-I, guide tones, 1-3-7 first.
    gen_cycle("2-1B", major_iivi, guide137_voicing, "1-3-7")

    # Exercise 2-1A, ii-V-I, guide tones, 1-7-3 first.
    gen_cycle("2-1A", major_iivi, guide137_voicing, "1-7-3", chord_start=2)

    # Exercise 2-1B, ii-V-I, guide tones, 1-3-7 first.
    gen_cycle("2-1B", major_iivi, guide173_voicing, "1-3-7", chord_start=2)

    # Exercise 2-1A, ii-V-I, guide tones, 1-7-3 first.
    gen_cycle("2-1A", major_iivi, guide173_voicing, "1-7-3", chord_start=3)

    # Exercise 2-1B, ii-V-I, guide tones, 1-3-7 first.
    gen_cycle("2-1B", major_iivi, guide137_voicing, "1-3-7", chord_start=3)

    # Repeat of the above, but rootless.
    gen_cycle("2-1Ar", major_iivi, guide73_voicing, "7-3")
    gen_cycle("2-1Br", major_iivi, guide37_voicing, "3-7")
    gen_cycle("2-1Ar", major_iivi, guide37_voicing, "7-3", chord_start=2)
    gen_cycle("2-1Br", major_iivi, guide73_voicing, "3-7", chord_start=2)
    gen_cycle("2-1Ar", major_iivi, guide73_voicing, "7-3", chord_start=3)
    gen_cycle("2-1Br", major_iivi, guide37_voicing, "3-7", chord_start=3)

    if false
        # Exercise 2-4A, Diminshed / Tri-tone sub, guide tones, 1-7-3 first.
        gen_cycle("2-4A", min_tritone_sub, "guide", "1-7-3")

        # Exercise 2-4B, Diminshed / Tri-tone sub, guide tones, 1-3-7 first.
        gen_cycle("2-4B", min_tritone_sub, "guide", "1-3-7")
    end
end

generate()