text:
"""
    Euler{T<:Real,Dim}

The compressible Euler equations in `Dim` dimensions using `T` as scalar type.
"""
struct Euler{T<:Real,Dim} <: AbstractBalanceLaw{Dim}
    γ::T
end

function Euler(γ=1.4, Dim=1)
  Euler{typeof(γ), Dim}(γ)
end

function Base.show(io::IO, model::Euler{T,Dim}) where {T,Dim}
  print(io, "Compressible Euler equations with γ=", model.γ, " {T=", T, ", Dim=", Dim, "}")
end

"""
    EulerVar1D{T<:Real}

Conserved variables (ϱ, ϱv, ϱe) of the Euler equations in one space dimension.
"""
struct EulerVar1D{T<:Real} <: FieldVector{3,T}
  ϱ ::T
  ϱv::T
  ϱe::T
end

function (::Type{EulerVar1D{T}})(val::Real) where {T}
  EulerVar2D{T}(val, val, val)
end

function similar_type(::EulerVar1D{T}) where {T}
  EulerVar1D{T}
end

@inline variables(model::Euler{T,1}) where {T} = EulerVar1D{T}

"""
    EulerVar2D{T<:Real}

Conserved variables (ϱ, ϱvx, ϱvy, ϱe) of the Euler equations in two space dimensions.
"""
struct EulerVar2D{T<:Real} <: FieldVector{4,T}
  ϱ  ::T
  ϱvx::T
  ϱvy::T
  ϱe ::T
end

function (::Type{EulerVar2D{T}})(val::Real) where {T}
  EulerVar2D{T}(val, val, val, val)
end

function similar_type(::EulerVar2D{T}) where {T}
  EulerVar2D{T}
end

@inline variables(model::Euler{T,2}) where {T} = EulerVar2D{T}

"""
    EulerVar3D{T<:Real}

Conserved variables (ϱ, ϱvx, ϱvy, ϱvz, ϱe) of the Euler equations in three space dimensions.
"""
struct EulerVar3D{T<:Real} <: FieldVector{5,T}
  ϱ  ::T
  ϱvx::T
  ϱvy::T
  ϱvz::T
  ϱe ::T
end

function (::Type{EulerVar3D{T}})(val::Real) where {T}
  EulerVar3D{T}(val, val, val, val, val)
end

function similar_type(::EulerVar3D{T}) where {T}
  EulerVar3D{T}
end

@inline variables(model::Euler{T,3}) where {T} = EulerVar3D{T}

"""
    IntegralQuantitiesEuler{T<:Real}

Some integrated quantities of interest for the Euler equations. Can be used in
callbacks of DifferentialEquations.jl.
"""
struct IntegralQuantitiesEuler{T<:Real} <: FieldVector{2,T}
    kinetic_energy::T
    entropy::T
end

function IntegralQuantitiesEuler(u, model::Euler)
    IntegralQuantitiesEuler(
        kinetic_energy(u, model),
        entropy(u, model)
    )
end

@inline function primitive_variables(u::EulerVar1D, model::Euler)
    @unpack ϱ, ϱv, ϱe = u
    @unpack γ = model
    if ϱ ≈ 0
        v = zero(ϱ)
    else
        v = ϱv / ϱ
    end
    p = (γ-1) * (ϱe - ϱv*v/2)

    ϱ, v, p
end

@inline function primitive_variables(u::EulerVar2D, model::Euler)
    @unpack ϱ, ϱvx, ϱvy, ϱe = u
    @unpack γ = model
    if ϱ ≈ 0
        vx = zero(ϱ)
        vy = zero(ϱ)
    else
        vx = ϱvx / ϱ
        vy = ϱvy / ϱ
    end
    p = (γ-1) * (ϱe - ϱvx*vx/2 - ϱvy*vy/2)

    ϱ, vx, vy, p
end

@inline function primitive_variables(u::EulerVar3D, model::Euler)
    @unpack ϱ, ϱvx, ϱvy, ϱvz, ϱe = u
    @unpack γ = model
    if ϱ ≈ 0
        vx = zero(ϱ)
        vy = zero(ϱ)
        vz = zero(ϱ)
    else
        vx = ϱvx / ϱ
        vy = ϱvy / ϱ
        vz = ϱvz / ϱ
    end
    p = (γ-1) * (ϱe - ϱvx*vx/2 - ϱvy*vy/2 - ϱvz*vz/2)

    ϱ, vx, vy, vz, p
end

@inline function conserved_variables(ϱ, v, p, model::Euler)
    @unpack γ = model

    variables(model)(ϱ, ϱ*v, ϱ*v^2/2+p/(γ-1))
end

@inline function conserved_variables(ϱ, vx, vy, p, model::Euler)
    @unpack γ = model

    variables(model)(ϱ, ϱ*vx, ϱ*vy, (ϱ*vx^2+ϱ*vy^2)/2+p/(γ-1))
end

@inline function conserved_variables(ϱ, vx, vy, vz, p, model::Euler)
    @unpack γ = model

    variables(model)(ϱ, ϱ*vx, ϱ*vy, ϱ*vz, (ϱ*vx^2+ϱ*vy^2+ϱ*vz^2)/2+p/(γ-1))
end

@inline function entropy_variables(u::EulerVar1D, model::Euler)
    entropy_variables(primitive_variables(u, model)..., model)
end
@inline function entropy_variables(ϱ, v, p, model::Euler)
    @unpack γ = model

    SVector((γ-log(p/ϱ^γ))/(γ-1)-ϱ*v^2/2p, ϱ*v/p, -ϱ/p)
end

@inline function flux_potential(u::EulerVar1D, model::Euler)
    flux_potential(primitive_variables(u, model)..., model)
end
@inline function flux_potential(ϱ, v, p, model::Euler)
    ϱ*v
end

@inline function satisfies_physical_constraints(u::EulerVar2D, model::Euler)
    ϱ, vx, vy, p = primitive_variables(u, model)
    ϱ >= 0 && p >= 0
end

@inline function satisfies_physical_constraints(u::EulerVar3D, model::Euler)
    ϱ, vx, vy, vz, p = primitive_variables(u, model)
    ϱ >= 0 && p >= 0
end

Base.@pure function kinetic_energy(u::EulerVar2D, model::Euler)
    ϱ, vx, vy, p = primitive_variables(u, model)

    ϱ * (vx^2 + vy^2) / 2
end

Base.@pure function kinetic_energy(u::EulerVar3D, model::Euler)
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    ϱ * (vx^2 + vy^2 + vz^2) / 2
end

Base.@pure function entropy(u::EulerVar1D, model::Euler)
    @unpack γ = model
    ϱ, v, p = primitive_variables(u, model)

    -ϱ*log(p/ϱ^γ) / (γ-1)
end

Base.@pure function entropy(u::EulerVar2D, model::Euler)
    @unpack γ = model
    ϱ, vx, vy, p = primitive_variables(u, model)

    -ϱ*log(p/ϱ^γ) / (γ-1)
end

Base.@pure function entropy(u::EulerVar3D, model::Euler)
    @unpack γ = model
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    -ϱ*log(p/ϱ^γ) / (γ-1)
end

@inline function flux(u::EulerVar1D{T}, model::Euler) where {T}
    @unpack ϱv, ϱe = u
    ϱ, v, p = primitive_variables(u, model)

    SVector{3,T}(ϱv, ϱv*v + p, (ϱe+p)*v)
end

@inline function flux(u::EulerVar2D{T}, model::Euler, dir::Val{:x}) where {T}
    @unpack ϱvx, ϱvy, ϱe = u
    ϱ, vx, vy, p = primitive_variables(u, model)

    SVector{4,T}(ϱvx, ϱvx*vx + p, ϱvx*vy, (ϱe+p)*vx)
end

@inline function flux(u::EulerVar2D{T}, model::Euler, dir::Val{:y}) where {T}
    @unpack ϱvx, ϱvy, ϱe = u
    ϱ, vx, vy, p = primitive_variables(u, model)

    SVector{4,T}(ϱvy, ϱvy*vx, ϱvy*vy + p, (ϱe+p)*vy)
end

@inline function flux(u::EulerVar3D{T}, model::Euler, dir::Val{:x}) where {T}
    @unpack ϱvx, ϱvy, ϱvz, ϱe = u
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    SVector{5,T}(ϱvx, ϱvx*vx + p, ϱvy*vx, ϱvz*vx, (ϱe+p)*vx)
end

@inline function flux(u::EulerVar3D{T}, model::Euler, dir::Val{:y}) where {T}
    @unpack ϱvx, ϱvy, ϱvz, ϱe = u
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    SVector{5,T}(ϱvy, ϱvx*vy, ϱvy*vy + p, ϱvz*vy, (ϱe+p)*vy)
end

@inline function flux(u::EulerVar3D{T}, model::Euler, dir::Val{:z}) where {T}
    @unpack ϱvx, ϱvy, ϱvz, ϱe = u
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    SVector{5,T}(ϱvz, ϱvx*vz, ϱvy*vz, ϱvz*vz + p, (ϱe+p)*vz)
end

@inline function max_abs_speed(u::EulerVar1D, model::Euler)
    @unpack γ = model
    ϱ, v, p = primitive_variables(u, model)

    abs(v) + sqrt(γ * p / ϱ)
end

"""
    max_abs_speed(uₗ::EulerVar1D, uᵣ::EulerVar1D, model::Euler)

Estimate the maximal speed in the solution of the Riemann problem with left and
right states `uₗ`, `uᵣ` for `model` via the two rarefaction approximation, see
Guermond and Popov (2016) Fast Estimation from above of the maximum wave speed in
the Riemann roblem for the Euler equations, Journal of Computational Physics 321,
pp. 908-926, and
Chen and Shu (2017) Entropy stable high order discontinuous Galerkin methods with
suitable quadrature rules for hyperbolic conservation laws, Journal of Computational
Physics 345, pp. 427-461.
"""
@inline function max_abs_speed(uₗ::EulerVar1D, uᵣ::EulerVar1D, model::Euler)
    λ₋, λ₊ = min_max_speed(uₗ, uᵣ, model)

    max(abs(λ₋), abs(λ₊))
end

"""
    min_max_speed(uₗ::EulerVar1D, uᵣ::EulerVar1D, model::Euler)

Estimate the maximal speed in the solution of the Riemann problem with left and
right states `uₗ`, `uᵣ` for `model` via the two rarefaction approximation, see
Guermond and Popov (2016) Fast Estimation from above of the maximum wave speed in
the Riemann roblem for the Euler equations, Journal of Computational Physics 321,
pp. 908-926, and
Chen and Shu (2017) Entropy stable high order discontinuous Galerkin methods with
suitable quadrature rules for hyperbolic conservation laws, Journal of Computational
Physics 345, pp. 427-461.
"""
@inline function min_max_speed(uₗ::EulerVar1D, uᵣ::EulerVar1D, model::Euler)
    @unpack γ = model
    ϱₗ, vₗ, pₗ = primitive_variables(uₗ, model)
    aₗ = sqrt(γ * pₗ / ϱₗ)
    ϱᵣ, vᵣ, pᵣ = primitive_variables(uᵣ, model)
    aᵣ = sqrt(γ * pᵣ / ϱᵣ)

    pₘ = ( (aₗ + aᵣ + (γ-1)*(vₗ-vᵣ)/2) / (aₗ/pₗ^((γ-1)/2γ) + aᵣ/pᵣ^((γ-1)/2γ)) )^(2γ/(γ-1))
    if pₘ <= pₗ
        # approximated by rarefaction wave
        λ₋ = vₗ - aₗ
    else
        # approximated by shock wave
        λ₋ = vₗ - aₗ*sqrt(1 + (γ+1)/2γ * (pₘ/pₗ - 1))
    end
    if pₘ <= pᵣ
        # approximated by rarefaction wave
        λ₊ = vᵣ + aᵣ
    else
        # approximated by shock wave
        λ₊ = vᵣ + aᵣ*sqrt(1 + (γ+1)/2γ * (pₘ/pᵣ - 1))
    end

    λ₋, λ₊
end

@inline function max_abs_speed(u::EulerVar2D, model::Euler)
    @unpack γ = model
    ϱ, vx, vy, p = primitive_variables(u, model)

    max(abs(vx),abs(vy)) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar2D, model::Euler, dir::Val{:x})
    @unpack γ = model
    ϱ, vx, vy, p = primitive_variables(u, model)

    abs(vx) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar2D, model::Euler, dir::Val{:y})
    @unpack γ = model
    ϱ, vx, vy, p = primitive_variables(u, model)

    abs(vy) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar3D, model::Euler)
    @unpack γ = model
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    max(abs(vx),abs(vy),abs(vz)) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar3D, model::Euler, dir::Val{:x})
    @unpack γ = model
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    abs(vx) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar3D, model::Euler, dir::Val{:y})
    @unpack γ = model
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    abs(vy) + sqrt(γ * p / ϱ)
end

@inline function max_abs_speed(u::EulerVar3D, model::Euler, dir::Val{:z})
    @unpack γ = model
    ϱ, vx, vy, vz, p = primitive_variables(u, model)

    abs(vz) + sqrt(γ * p / ϱ)
end

@recipe function f{Ξ,T}(ξumodel::Tuple{Ξ,Vector{EulerVar1D{T}},Euler{T,1}})
  ξ, u, model = ξumodel

  ϱ  = mappedarray(u->u.ϱ, u)
  v  = mappedarray(u->primitive_variables(u,model)[2], u)
  p  = mappedarray(u->primitive_variables(u,model)[3], u)

  size   --> (1000,400)
  layout --> (1,3)
  legend --> false

  @series begin
      subplot := 1
      yguide --> L"\varrho"
      ξ, ϱ
  end

  @series begin
      subplot := 2
      yguide --> L"v"
      ξ, v
  end

  @series begin
      subplot := 3
      yguide --> L"p"
      ξ, p
  end
end

@recipe function f{X,Y,T}(xyumodel::Tuple{X,Y,AbstractArray{EulerVar2D{T}},Euler})
    x, y, u, model = xyumodel

    ϱ  = mappedarray(u->primitive_variables(u,model)[1], u)
    vx = mappedarray(u->primitive_variables(u,model)[2], u)
    vy = mappedarray(u->primitive_variables(u,model)[3], u)
    p  = mappedarray(u->primitive_variables(u,model)[4], u)

    layout --> (2,2)
    legend --> false

    @series begin
        subplot := 1
        title  --> L"\varrho"
        x, y, ϱ
    end

    @series begin
        subplot := 2
        title  --> L"p"
        x, y, p
    end

    @series begin
        subplot := 3
        title  --> L"v_x"
        x, y, vx
    end

    @series begin
        subplot := 4
        title  --> L"v_y"
        x, y, vy
    end
end

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

"""
    (fnum::SuliciuFlux)(uₗ::EulerVar1D{T}, uᵣ::EulerVar1D{T}, model::Euler{T,1}) where T

The Suliciu relaxation solver / numerical flux, see
Bouchut (2004) Nonlinear Stability of Finite Volume Methods for Hyperbolic
Conservation Laws and Well-Balanced Schemes for Sources.
"""
@inline function (fnum::SuliciuFlux)(uₗ::EulerVar1D{T}, uᵣ::EulerVar1D{T}, model::Euler{T,1}) where T
    @unpack γ = model
    ϱeₗ = uₗ.ϱe
    ϱeᵣ = uᵣ.ϱe
    ϱₗ, vₗ, pₗ = primitive_variables(uₗ, model)
    cₗ = sqrt(γ * pₗ / ϱₗ)
    ɛₗ = ϱeₗ/ϱₗ - vₗ*vₗ/2
    ϱᵣ, vᵣ, pᵣ = primitive_variables(uᵣ, model)
    cᵣ = sqrt(γ * pᵣ / ϱᵣ)
    ɛᵣ = ϱeᵣ/ϱᵣ - vᵣ*vᵣ/2
    α = (γ+1)/2

    # compute speeds
    cₗ_ϱₗ = zero(T)
    cᵣ_ϱᵣ = zero(T)
    if pₗ <= pᵣ && 0 < pᵣ
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(ϱᵣ*cᵣ) + vₗ - vᵣ )
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(cₗ_ϱₗ*ϱₗ) + vₗ - vᵣ )
    elseif pᵣ <= pₗ && 0 < pₗ
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(ϱₗ*cₗ) + vₗ - vᵣ )
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(cᵣ_ϱᵣ*ϱᵣ) + vₗ - vᵣ )
    end

    # compute intermediate values
    cₗ = cₗ_ϱₗ*ϱₗ
    cᵣ = cᵣ_ϱᵣ*ϱᵣ
    vs = ( cₗ*vₗ + cᵣ*vᵣ + pₗ - pᵣ ) / ( cₗ+cᵣ )
    vs = ifelse(isnan(vs), zero(T), vs)
    pₗs = pᵣs = ( cᵣ*pₗ + cₗ*pᵣ - cₗ*cᵣ*(vᵣ-vₗ) ) / (cₗ+cᵣ)
    pₗs = pᵣs = ifelse(isnan(pₗs), zero(T), pₗs)
    ϱₗs = 1 / ( 1/ϱₗ + (cᵣ*(vᵣ-vₗ) + pₗ - pᵣ ) / (cₗ*(cₗ+cᵣ)) )
    ϱₗs = ifelse(isnan(ϱₗs), zero(T), ϱₗs)
    ϱᵣs = 1 / ( 1/ϱᵣ + (cₗ*(vᵣ-vₗ) + pᵣ - pₗ ) / (cᵣ*(cₗ+cᵣ)) )
    ϱᵣs = ifelse(isnan(ϱᵣs), zero(T), ϱᵣs)
    ɛₗs = ɛₗ + (pₗs*pₗs - pₗ*pₗ)/(2*cₗ*cₗ)
    ɛᵣs = ɛᵣ + (pᵣs*pᵣs - pᵣ*pᵣ)/(2*cᵣ*cᵣ)

    # compute fluxes
    fϱ  = zero(T)
    fϱv = zero(T)
    fϱe = zero(T)
    if 0 <= vₗ-cₗ_ϱₗ
        fϱ  = ϱₗ*vₗ
        fϱv = ϱₗ*vₗ*vₗ + pₗ
        fϱe = (ϱₗ*vₗ*vₗ/2 + ϱₗ*ɛₗ + pₗ)*vₗ
    elseif 0 <= vs
        fϱ  = ϱₗs*vs
        fϱv = ϱₗs*vs*vs + pₗs
        fϱe = (ϱₗs*vs*vs/2 + ϱₗs*ɛₗs + pₗs)*vs
    elseif 0 <= vᵣ+cᵣ_ϱᵣ
        fϱ  = ϱᵣs*vs
        fϱv = ϱᵣs*vs*vs + pᵣs
        fϱe = (ϱᵣs*vs*vs/2 + ϱᵣs*ɛᵣs + pᵣs)*vs
    else
        fϱ  = ϱᵣ*vᵣ
        fϱv = ϱᵣ*vᵣ*vᵣ + pᵣ
        fϱe = (ϱᵣ*vᵣ*vᵣ/2 + ϱᵣ*ɛᵣ + pᵣ)*vᵣ
    end

    SVector(fϱ, fϱv, fϱe)
end

@inline function (fnum::SuliciuFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱeₗ = uₗ.ϱe
    ϱeᵣ = uᵣ.ϱe
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    cₗ = sqrt(γ * pₗ / ϱₗ)
    ɛₗ = ϱeₗ/ϱₗ - vxₗ*vxₗ/2 - vyₗ*vyₗ/2
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    cᵣ = sqrt(γ * pᵣ / ϱᵣ)
    ɛᵣ = ϱeᵣ/ϱᵣ - vxᵣ*vxᵣ/2 - vyᵣ*vyᵣ/2
    α = (γ+1)/2

    # compute speeds
    cₗ_ϱₗ = zero(T)
    cᵣ_ϱᵣ = zero(T)
    if pₗ <= pᵣ && 0 < pᵣ
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(ϱᵣ*cᵣ) + vxₗ - vxᵣ )
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(cₗ_ϱₗ*ϱₗ) + vxₗ - vxᵣ )
    elseif pᵣ <= pₗ && 0 < pₗ
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(ϱₗ*cₗ) + vxₗ - vxᵣ )
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(cᵣ_ϱᵣ*ϱᵣ) + vxₗ - vxᵣ )
    end

    # compute intermediate values
    cₗ = cₗ_ϱₗ*ϱₗ
    cᵣ = cᵣ_ϱᵣ*ϱᵣ
    vxs = ( cₗ*vxₗ + cᵣ*vxᵣ + pₗ - pᵣ ) / ( cₗ+cᵣ )
    vxs = ifelse(isnan(vxs), zero(T), vxs)
    pₗs = pᵣs = ( cᵣ*pₗ + cₗ*pᵣ - cₗ*cᵣ*(vxᵣ-vxₗ) ) / (cₗ+cᵣ)
    pₗs = pᵣs = ifelse(isnan(pₗs), zero(T), pₗs)
    ϱₗs = 1 / ( 1/ϱₗ + (cᵣ*(vxᵣ-vxₗ) + pₗ - pᵣ ) / (cₗ*(cₗ+cᵣ)) )
    ϱₗs = ifelse(isnan(ϱₗs), zero(T), ϱₗs)
    ϱᵣs = 1 / ( 1/ϱᵣ + (cₗ*(vxᵣ-vxₗ) + pᵣ - pₗ ) / (cᵣ*(cₗ+cᵣ)) )
    ϱᵣs = ifelse(isnan(ϱᵣs), zero(T), ϱᵣs)
    ɛₗs = ɛₗ + (pₗs*pₗs - pₗ*pₗ)/(2*cₗ*cₗ)
    ɛᵣs = ɛᵣ + (pᵣs*pᵣs - pᵣ*pᵣ)/(2*cᵣ*cᵣ)

    # compute fluxes
    fϱ   = zero(T)
    fϱvx = zero(T)
    fϱvy = zero(T)
    fϱe  = zero(T)
    if 0 <= vxₗ-cₗ_ϱₗ
        fϱ   = ϱₗ*vxₗ
        fϱvx = ϱₗ*vxₗ*vxₗ + pₗ
        fϱvy = ϱₗ*vxₗ*vyₗ
        fϱe  = (ϱₗ*vxₗ*vxₗ/2 + ϱₗ*vyₗ*vyₗ/2 + ϱₗ*ɛₗ + pₗ)*vxₗ
    elseif 0 <= vxs
        fϱ   = ϱₗs*vxs
        fϱvx = ϱₗs*vxs*vxs + pₗs
        fϱvy = ϱₗ*vxs*vyₗ
        fϱe  = (ϱₗs*vxs*vxs/2 + ϱₗs*vyₗ*vyₗ/2 + ϱₗs*ɛₗs + pₗs)*vxs
    elseif 0 <= vxᵣ+cᵣ_ϱᵣ
        fϱ   = ϱᵣs*vxs
        fϱvx = ϱᵣs*vxs*vxs + pᵣs
        fϱvy = ϱₗ*vxs*vyᵣ
        fϱe  = (ϱᵣs*vxs*vxs/2 + ϱᵣs*vyᵣ*vyᵣ/2 + ϱᵣs*ɛᵣs + pᵣs)*vxs
    else
        fϱ   = ϱᵣ*vxᵣ
        fϱvx = ϱᵣ*vxᵣ*vxᵣ + pᵣ
        fϱvy = ϱₗ*vxᵣ*vyᵣ
        fϱe  = (ϱᵣ*vxᵣ*vxᵣ/2 + ϱᵣ*vyᵣ*vyᵣ/2 + ϱᵣ*ɛᵣ + pᵣ)*vxᵣ
    end

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::SuliciuFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱeₗ = uₗ.ϱe
    ϱeᵣ = uᵣ.ϱe
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    cₗ = sqrt(γ * pₗ / ϱₗ)
    ɛₗ = ϱeₗ/ϱₗ - vxₗ*vxₗ/2 - vyₗ*vyₗ/2
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    cᵣ = sqrt(γ * pᵣ / ϱᵣ)
    ɛᵣ = ϱeᵣ/ϱᵣ - vxᵣ*vxᵣ/2 - vyᵣ*vyᵣ/2
    α = (γ+1)/2

    # compute speeds
    cₗ_ϱₗ = zero(T)
    cᵣ_ϱᵣ = zero(T)
    if pₗ <= pᵣ && 0 < pᵣ
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(ϱᵣ*cᵣ) + vyₗ - vyᵣ )
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(cₗ_ϱₗ*ϱₗ) + vyₗ - vyᵣ )
    elseif pᵣ <= pₗ && 0 < pₗ
      cᵣ_ϱᵣ = cᵣ + α*max( zero(T), (pₗ-pᵣ)/(ϱₗ*cₗ) + vyₗ - vyᵣ )
      cₗ_ϱₗ = cₗ + α*max( zero(T), (pᵣ-pₗ)/(cᵣ_ϱᵣ*ϱᵣ) + vyₗ - vyᵣ )
    end

    # compute intermediate values
    cₗ = cₗ_ϱₗ*ϱₗ
    cᵣ = cᵣ_ϱᵣ*ϱᵣ
    vys = ( cₗ*vyₗ + cᵣ*vyᵣ + pₗ - pᵣ ) / ( cₗ+cᵣ )
    vys = ifelse(isnan(vys), zero(T), vys)
    pₗs = pᵣs = ( cᵣ*pₗ + cₗ*pᵣ - cₗ*cᵣ*(vyᵣ-vyₗ) ) / (cₗ+cᵣ)
    pₗs = pᵣs = ifelse(isnan(pₗs), zero(T), pₗs)
    ϱₗs = 1 / ( 1/ϱₗ + (cᵣ*(vyᵣ-vyₗ) + pₗ - pᵣ ) / (cₗ*(cₗ+cᵣ)) )
    ϱₗs = ifelse(isnan(ϱₗs), zero(T), ϱₗs)
    ϱᵣs = 1 / ( 1/ϱᵣ + (cₗ*(vyᵣ-vyₗ) + pᵣ - pₗ ) / (cᵣ*(cₗ+cᵣ)) )
    ϱᵣs = ifelse(isnan(ϱᵣs), zero(T), ϱᵣs)
    ɛₗs = ɛₗ + (pₗs*pₗs - pₗ*pₗ)/(2*cₗ*cₗ)
    ɛᵣs = ɛᵣ + (pᵣs*pᵣs - pᵣ*pᵣ)/(2*cᵣ*cᵣ)

    # compute fluxes
    fϱ   = zero(T)
    fϱvx = zero(T)
    fϱvy = zero(T)
    fϱe  = zero(T)
    if 0 <= vyₗ-cₗ_ϱₗ
        fϱ   = ϱₗ*vyₗ
        fϱvx = ϱₗ*vxₗ*vyₗ
        fϱvy = ϱₗ*vyₗ*vyₗ + pₗ
        fϱe  = (ϱₗ*vxₗ*vxₗ/2 + ϱₗ*vyₗ*vyₗ/2 + ϱₗ*ɛₗ + pₗ)*vyₗ
    elseif 0 <= vys
        fϱ   = ϱₗs*vys
        fϱvx = ϱₗs*vxₗ*vys
        fϱvy = ϱₗ*vys*vys + pₗs
        fϱe  = (ϱₗs*vxₗ*vxₗ/2 + ϱₗs*vys*vys/2 + ϱₗs*ɛₗs + pₗs)*vys
    elseif 0 <= vyᵣ+cᵣ_ϱᵣ
        fϱ   = ϱᵣs*vys
        fϱvx = ϱᵣs*vxᵣ*vys
        fϱvy = ϱₗ*vys*vys + pᵣs
        fϱe  = (ϱᵣs*vxᵣ*vxᵣ/2 + ϱᵣs*vys*vys/2 + ϱᵣs*ɛᵣs + pᵣs)*vys
    else
        fϱ   = ϱᵣ*vyᵣ
        fϱvx = ϱᵣ*vxᵣ*vyᵣ
        fϱvy = ϱₗ*vyᵣ*vyᵣ + pᵣ
        fϱe  = (ϱᵣ*vxᵣ*vxᵣ/2 + ϱᵣ*vyᵣ*vyᵣ/2 + ϱᵣ*ɛᵣ + pᵣ)*vyᵣ
    end

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

# TODO: (fnum::SuliciuFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T}, model::Euler{T,3}, dir::Val{:x}) where T
# TODO: (fnum::SuliciuFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T}, model::Euler{T,3}, dir::Val{:y}) where T
# TODO: (fnum::SuliciuFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T}, model::Euler{T,3}, dir::Val{:z}) where T

@inline function (fnum::ChandrashekarFluxEC)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                                model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    βₗ = ϱₗ / 2pₗ
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    βᵣ = ϱᵣ / 2pᵣ

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2 + vxᵣ^2+vyᵣ^2) / 2
    β    = (βₗ + βᵣ) / 2
    βlog = logmean(βₗ, βᵣ)

    fϱ   = ϱlog*vx
    fϱvx = vx*fϱ + ϱ/2β
    fϱvy = vy*fϱ
    fϱe  = 1/(2γ-2)*fϱ/βlog - v2*fϱ/2 + vx*fϱvx + vy*fϱvy

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::ChandrashekarFluxEC)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                                model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    βₗ = ϱₗ / 2pₗ
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    βᵣ = ϱᵣ / 2pᵣ

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2 + vxᵣ^2+vyᵣ^2) / 2
    β    = (βₗ + βᵣ) / 2
    βlog = logmean(βₗ, βᵣ)

    fϱ   = ϱlog*vy
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ + ϱ/2β
    fϱe  = 1/(2γ-2)*fϱ/βlog - v2*fϱ/2 + vx*fϱvx + vy*fϱvy

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::ChandrashekarFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    βₗ = ϱₗ / 2pₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    βᵣ = ϱᵣ / 2pᵣ

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    β    = (βₗ + βᵣ) / 2
    βlog = logmean(βₗ, βᵣ)

    fϱ   = ϱlog*vx
    fϱvx = vx*fϱ + ϱ/2β
    fϱvy = vy*fϱ
    fϱvz = vz*fϱ
    fϱe  = 1/(2γ-2)*fϱ/βlog - v2*fϱ/2 + vx*fϱvx + vy*fϱvy + vz*fϱvz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::ChandrashekarFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    βₗ = ϱₗ / 2pₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    βᵣ = ϱᵣ / 2pᵣ

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    β    = (βₗ + βᵣ) / 2
    βlog = logmean(βₗ, βᵣ)

    fϱ   = ϱlog*vy
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ + ϱ/2β
    fϱvz = vz*fϱ
    fϱe  = 1/(2γ-2)*fϱ/βlog - v2*fϱ/2 + vx*fϱvx + vy*fϱvy + vz*fϱvz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::ChandrashekarFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    βₗ = ϱₗ / 2pₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    βᵣ = ϱᵣ / 2pᵣ

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    β    = (βₗ + βᵣ) / 2
    βlog = logmean(βₗ, βᵣ)

    fϱ   = ϱlog*vz
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ
    fϱvz = vz*fϱ + ϱ/2β
    fϱe  = 1/(2γ-2)*fϱ/βlog - v2*fϱ/2 + vx*fϱvx + vy*fϱvy + vz*fϱvz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

Base.@pure function roe_variables(u::EulerVar2D, model::Euler)
    ϱ, vx, vy, p = primitive_variables(u, model)
    z1 = sqrt(ϱ/p)
    z2 = z1*vx
    z3 = z1*vy
    z5 = p*z1

    z1, z2, z3, z5
end

Base.@pure function roe_variables(u::EulerVar3D, model::Euler)
    ϱ, vx, vy, vz, p = primitive_variables(u, model)
    z1 = sqrt(ϱ/p)
    z2 = z1*vx
    z3 = z1*vy
    z4 = z1*vz
    z5 = p*z1

    z1, z2, z3, z4, z5
end

@inline function (fnum::IsmailRoeFluxEC)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                            model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    z1ₗ, z2ₗ, z3ₗ, z5ₗ = roe_variables(uₗ, model)
    z1ᵣ, z2ᵣ, z3ᵣ, z5ᵣ = roe_variables(uᵣ, model)

    z1    = (z1ₗ + z1ᵣ) / 2
    z1log = logmean(z1ₗ, z1ᵣ)
    z2    = (z2ₗ + z2ᵣ) / 2
    z3    = (z3ₗ + z3ᵣ) / 2
    z5    = (z5ₗ + z5ᵣ) / 2
    z5log = logmean(z5ₗ, z5ᵣ)

    ϱ  = z1*z5log
    vx = z2/z1
    vy = z3/z1
    p1 = z5/z1
    p2 = ( (γ+1)*z5log/z1log + (γ-1)*z5/z1 ) / 2γ
    h  = γ/(γ-1) * p2/ϱ + (vx*vx + vy*vy)/2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p1
    fϱvy = ϱ*vx*vy
    fϱe  = ϱ*vx*h

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::IsmailRoeFluxEC)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                            model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    z1ₗ, z2ₗ, z3ₗ, z5ₗ = roe_variables(uₗ, model)
    z1ᵣ, z2ᵣ, z3ᵣ, z5ᵣ = roe_variables(uᵣ, model)

    z1    = (z1ₗ + z1ᵣ) / 2
    z1log = logmean(z1ₗ, z1ᵣ)
    z2    = (z2ₗ + z2ᵣ) / 2
    z3    = (z3ₗ + z3ᵣ) / 2
    z5    = (z5ₗ + z5ᵣ) / 2
    z5log = logmean(z5ₗ, z5ᵣ)

    ϱ  = z1*z5log
    vx = z2/z1
    vy = z3/z1
    p1 = z5/z1
    p2 = ( (γ+1)*z5log/z1log + (γ-1)*z5/z1 ) / 2γ
    h  = γ/(γ-1) * p2/ϱ + (vx*vx + vy*vy)/2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p1
    fϱe  = ϱ*vy*h

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::IsmailRoeFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    z1ₗ, z2ₗ, z3ₗ, z4ₗ, z5ₗ = roe_variables(uₗ, model)
    z1ᵣ, z2ᵣ, z3ᵣ, z4ᵣ, z5ᵣ = roe_variables(uᵣ, model)

    z1    = (z1ₗ + z1ᵣ) / 2
    z1log = logmean(z1ₗ, z1ᵣ)
    z2    = (z2ₗ + z2ᵣ) / 2
    z3    = (z3ₗ + z3ᵣ) / 2
    z4    = (z4ₗ + z4ᵣ) / 2
    z5    = (z5ₗ + z5ᵣ) / 2
    z5log = logmean(z5ₗ, z5ᵣ)

    ϱ  = z1*z5log
    vx = z2/z1
    vy = z3/z1
    vz = z4/z1
    p1 = z5/z1
    p2 = ( (γ+1)*z5log/z1log + (γ-1)*z5/z1 ) / 2γ
    h  = γ/(γ-1) * p2/ϱ + (vx^2+vy^2+vz^2)/2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p1
    fϱvy = ϱ*vx*vy
    fϱvz = ϱ*vx*vz
    fϱe  = ϱ*vx*h

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::IsmailRoeFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    z1ₗ, z2ₗ, z3ₗ, z4ₗ, z5ₗ = roe_variables(uₗ, model)
    z1ᵣ, z2ᵣ, z3ᵣ, z4ᵣ, z5ᵣ = roe_variables(uᵣ, model)

    z1    = (z1ₗ + z1ᵣ) / 2
    z1log = logmean(z1ₗ, z1ᵣ)
    z2    = (z2ₗ + z2ᵣ) / 2
    z3    = (z3ₗ + z3ᵣ) / 2
    z4    = (z4ₗ + z4ᵣ) / 2
    z5    = (z5ₗ + z5ᵣ) / 2
    z5log = logmean(z5ₗ, z5ᵣ)

    ϱ  = z1*z5log
    vx = z2/z1
    vy = z3/z1
    vz = z4/z1
    p1 = z5/z1
    p2 = ( (γ+1)*z5log/z1log + (γ-1)*z5/z1 ) / 2γ
    h  = γ/(γ-1) * p2/ϱ + (vx^2+vy^2+vz^2)/2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p1
    fϱvz = ϱ*vy*vz
    fϱe  = ϱ*vy*h

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::IsmailRoeFluxEC)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    z1ₗ, z2ₗ, z3ₗ, z4ₗ, z5ₗ = roe_variables(uₗ, model)
    z1ᵣ, z2ᵣ, z3ᵣ, z4ᵣ, z5ᵣ = roe_variables(uᵣ, model)

    z1    = (z1ₗ + z1ᵣ) / 2
    z1log = logmean(z1ₗ, z1ᵣ)
    z2    = (z2ₗ + z2ᵣ) / 2
    z3    = (z3ₗ + z3ᵣ) / 2
    z4    = (z4ₗ + z4ᵣ) / 2
    z5    = (z5ₗ + z5ᵣ) / 2
    z5log = logmean(z5ₗ, z5ᵣ)

    ϱ  = z1*z5log
    vx = z2/z1
    vy = z3/z1
    vz = z4/z1
    p1 = z5/z1
    p2 = ( (γ+1)*z5log/z1log + (γ-1)*z5/z1 ) / 2γ
    h  = γ/(γ-1) * p2/ϱ + (vx^2+vy^2+vz^2)/2

    fϱ   = ϱ*vz
    fϱvx = ϱ*vx*vz
    fϱvy = ϱ*vy*vz
    fϱvz = ϱ*vz^2 + p1
    fϱe  = ϱ*vz*h

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

@inline function (fnum::RanochaFluxECandKEP)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                                model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2 + vxᵣ^2+vyᵣ^2) / 2
    p    = (pₗ + pᵣ) / 2
    ϱ_p_log = logmean(ϱₗ/pₗ, ϱᵣ/pᵣ)

    fϱ   = ϱlog*vx
    fϱvx = vx*fϱ + p
    fϱvy = vy*fϱ
    fϱe  = (ϱlog/((γ-1)*ϱ_p_log) + p + ϱlog*(vx^2+vy^2-v2/2)) * vx - (pᵣ-pₗ)*(vxᵣ-vxₗ)/4

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::RanochaFluxECandKEP)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                                model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2 + vxᵣ^2+vyᵣ^2) / 2
    p    = (pₗ + pᵣ) / 2
    ϱ_p_log = logmean(ϱₗ/pₗ, ϱᵣ/pᵣ)

    fϱ   = ϱlog*vy
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ + p
    fϱe  = (ϱlog/((γ-1)*ϱ_p_log) + p + ϱlog*(vx^2+vy^2-v2/2)) * vy - (pᵣ-pₗ)*(vyᵣ-vyₗ)/4

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::RanochaFluxECandKEP)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    p    = (pₗ + pᵣ) / 2
    ϱ_p_log = logmean(ϱₗ/pₗ, ϱᵣ/pᵣ)

    fϱ   = ϱlog*vx
    fϱvx = vx*fϱ + p
    fϱvy = vy*fϱ
    fϱvz = vz*fϱ
    fϱe  = (ϱlog/((γ-1)*ϱ_p_log) + p + ϱlog*(vx^2+vy^2+vz^2-v2/2)) * vx - (pᵣ-pₗ)*(vxᵣ-vxₗ)/4

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::RanochaFluxECandKEP)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    p    = (pₗ + pᵣ) / 2
    ϱ_p_log = logmean(ϱₗ/pₗ, ϱᵣ/pᵣ)

    fϱ   = ϱlog*vy
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ + p
    fϱvz = vz*fϱ
    fϱe  = (ϱlog/((γ-1)*ϱ_p_log) + p + ϱlog*(vx^2+vy^2+vz^2-v2/2)) * vy - (pᵣ-pₗ)*(vyᵣ-vyₗ)/4

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::RanochaFluxECandKEP)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                                model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)

    ϱ    = (ϱₗ + ϱᵣ) / 2
    ϱlog = logmean(ϱₗ, ϱᵣ)
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    v2   = (vxₗ^2+vyₗ^2+vzₗ^2 + vxᵣ^2+vyᵣ^2+vzᵣ^2) / 2
    p    = (pₗ + pᵣ) / 2
    ϱ_p_log = logmean(ϱₗ/pₗ, ϱᵣ/pᵣ)

    fϱ   = ϱlog*vz
    fϱvx = vx*fϱ
    fϱvy = vy*fϱ
    fϱvz = vz*fϱ + p
    fϱe  = (ϱlog/((γ-1)*ϱ_p_log) + p + ϱlog*(vx^2+vy^2+vz^2-v2/2)) * vz - (pᵣ-pₗ)*(vzᵣ-vzₗ)/4

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

@inline function (fnum::MorinishiFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ = uₗ.ϱvx
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ = uᵣ.ϱvx

    ϱvx  = (ϱvxₗ + ϱvxᵣ) / 2
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    p    = (pₗ + pᵣ) / 2
    pvx  = (pₗ*vxₗ + pᵣ*vxᵣ) / 2
    ϱvxx = (ϱvxₗ*vxₗ + ϱvxᵣ*vxᵣ) / 2
    ϱvxy = (ϱvxₗ*vyₗ + ϱvxᵣ*vyᵣ) / 2
    ϱvxv2= (ϱvxₗ*(vxₗ^2+vyₗ^2) + ϱvxᵣ*(vxᵣ^2+vyᵣ^2)) / 2

    fϱ   = ϱvx
    fϱvx = ϱvx*vx + p
    fϱvy = ϱvx*vy
    fϱe  = γ/(γ-1)*pvx + ϱvxx*vx + ϱvxy*vy - ϱvxv2/2

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::MorinishiFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱvyₗ = uₗ.ϱvy
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvyᵣ = uᵣ.ϱvy

    ϱvy  = (ϱvyₗ + ϱvyᵣ) / 2
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    p    = (pₗ + pᵣ) / 2
    pvy  = (pₗ*vyₗ + pᵣ*vyᵣ) / 2
    ϱvyy = (ϱvyₗ*vyₗ + ϱvyᵣ*vyᵣ) / 2
    ϱvyx = (ϱvyₗ*vxₗ + ϱvyᵣ*vxᵣ) / 2
    ϱvyv2= (ϱvyₗ*(vxₗ^2+vyₗ^2) + ϱvyᵣ*(vxᵣ^2+vyᵣ^2)) / 2

    fϱ   = ϱvy
    fϱvx = ϱvy*vx
    fϱvy = ϱvy*vy + p
    fϱe  = γ/(γ-1)*pvy + ϱvyx*vx + ϱvyy*vy - ϱvyv2/2

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::MorinishiFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ = uₗ.ϱvx
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ = uᵣ.ϱvx

    ϱvx  = (ϱvxₗ + ϱvxᵣ) / 2
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    p    = (pₗ + pᵣ) / 2
    pvx  = (pₗ*vxₗ + pᵣ*vxᵣ) / 2
    ϱvxx = (ϱvxₗ*vxₗ + ϱvxᵣ*vxᵣ) / 2
    ϱvxy = (ϱvxₗ*vyₗ + ϱvxᵣ*vyᵣ) / 2
    ϱvxz = (ϱvxₗ*vzₗ + ϱvxᵣ*vzᵣ) / 2
    ϱvxv2= (ϱvxₗ*(vxₗ^2+vyₗ^2+vzₗ^2) + ϱvxᵣ*(vxᵣ^2+vyᵣ^2+vzᵣ^2)) / 2

    fϱ   = ϱvx
    fϱvx = ϱvx*vx + p
    fϱvy = ϱvx*vy
    fϱvz = ϱvx*vz
    fϱe  = γ/(γ-1)*pvx + ϱvxx*vx + ϱvxy*vy + ϱvxz*vz - ϱvxv2/2

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::MorinishiFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvyₗ = uₗ.ϱvy
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvyᵣ = uᵣ.ϱvy

    ϱvy  = (ϱvyₗ + ϱvyᵣ) / 2
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    p    = (pₗ + pᵣ) / 2
    pvy  = (pₗ*vyₗ + pᵣ*vyᵣ) / 2
    ϱvyx = (ϱvyₗ*vxₗ + ϱvyᵣ*vxᵣ) / 2
    ϱvyy = (ϱvyₗ*vyₗ + ϱvyᵣ*vyᵣ) / 2
    ϱvyz = (ϱvyₗ*vzₗ + ϱvyᵣ*vzᵣ) / 2
    ϱvyv2= (ϱvyₗ*(vxₗ^2+vyₗ^2+vzₗ^2) + ϱvyᵣ*(vxᵣ^2+vyᵣ^2+vzᵣ^2)) / 2

    fϱ   = ϱvy
    fϱvx = ϱvy*vx
    fϱvy = ϱvy*vy + p
    fϱvz = ϱvy*vz
    fϱe  = γ/(γ-1)*pvy + ϱvyx*vx + ϱvyy*vy + ϱvyz*vz - ϱvyv2/2

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::MorinishiFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvzₗ = uₗ.ϱvz
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvzᵣ = uᵣ.ϱvz

    ϱvz  = (ϱvzₗ + ϱvzᵣ) / 2
    vx   = (vxₗ + vxᵣ) / 2
    vy   = (vyₗ + vyᵣ) / 2
    vz   = (vzₗ + vzᵣ) / 2
    p    = (pₗ + pᵣ) / 2
    pvz  = (pₗ*vzₗ + pᵣ*vzᵣ) / 2
    ϱvzx = (ϱvzₗ*vxₗ + ϱvzᵣ*vxᵣ) / 2
    ϱvzy = (ϱvzₗ*vyₗ + ϱvzᵣ*vyᵣ) / 2
    ϱvzz = (ϱvzₗ*vzₗ + ϱvzᵣ*vzᵣ) / 2
    ϱvzv2= (ϱvzₗ*(vxₗ^2+vyₗ^2+vzₗ^2) + ϱvzᵣ*(vxᵣ^2+vyᵣ^2+vzᵣ^2)) / 2

    fϱ   = ϱvz
    fϱvx = ϱvz*vx
    fϱvy = ϱvz*vy
    fϱvz = ϱvz*vz + p
    fϱe  = γ/(γ-1)*pvz + ϱvzx*vx + ϱvzy*vy + ϱvzz*vz - ϱvzv2/2

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

@inline function (fnum::DucrosEtAlFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ, ϱvyₗ, ϱeₗ = uₗ.ϱvx, uₗ.ϱvy, uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ, ϱvyᵣ, ϱeᵣ = uᵣ.ϱvx, uᵣ.ϱvy, uᵣ.ϱe

    ϱ   = (ϱₗ + ϱᵣ) / 2
    ϱvx = (ϱvxₗ + ϱvxᵣ) / 2
    ϱvy = (ϱvyₗ + ϱvyᵣ) / 2
    ϱe  = (ϱeₗ + ϱeᵣ) / 2
    vx  = (vxₗ + vxᵣ) / 2
    vy  = (vyₗ + vyᵣ) / 2
    p   = (pₗ + pᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱvx*vx + p
    fϱvy = ϱvy*vx
    fϱe  = (ϱe+p)*vx

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::DucrosEtAlFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ, ϱvyₗ, ϱeₗ = uₗ.ϱvx, uₗ.ϱvy, uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ, ϱvyᵣ, ϱeᵣ = uᵣ.ϱvx, uᵣ.ϱvy, uᵣ.ϱe

    ϱ   = (ϱₗ + ϱᵣ) / 2
    ϱvx = (ϱvxₗ + ϱvxᵣ) / 2
    ϱvy = (ϱvyₗ + ϱvyᵣ) / 2
    ϱe  = (ϱeₗ + ϱeᵣ) / 2
    vx  = (vxₗ + vxᵣ) / 2
    vy  = (vyₗ + vyᵣ) / 2
    p   = (pₗ + pᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱvx*vy
    fϱvy = ϱvy*vy + p
    fϱe  = (ϱe+p)*vy

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::DucrosEtAlFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ, ϱvyₗ, ϱvzₗ, ϱeₗ = uₗ.ϱvx, uₗ.ϱvy, uₗ. ϱvz, uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ, ϱvyᵣ, ϱvzᵣ, ϱeᵣ = uᵣ.ϱvx, uᵣ.ϱvy, uᵣ.ϱvz, uᵣ.ϱe

    ϱ   = (ϱₗ + ϱᵣ) / 2
    ϱvx = (ϱvxₗ + ϱvxᵣ) / 2
    ϱvy = (ϱvyₗ + ϱvyᵣ) / 2
    ϱvz = (ϱvzₗ + ϱvzᵣ) / 2
    ϱe  = (ϱeₗ + ϱeᵣ) / 2
    vx  = (vxₗ + vxᵣ) / 2
    vy  = (vyₗ + vyᵣ) / 2
    vz  = (vzₗ + vzᵣ) / 2
    p   = (pₗ + pᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱvx*vx + p
    fϱvy = ϱvy*vx
    fϱvz = ϱvz*vx
    fϱe  = (ϱe+p)*vx

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::DucrosEtAlFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ, ϱvyₗ, ϱvzₗ, ϱeₗ = uₗ.ϱvx, uₗ.ϱvy, uₗ. ϱvz, uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ, ϱvyᵣ, ϱvzᵣ, ϱeᵣ = uᵣ.ϱvx, uᵣ.ϱvy, uᵣ.ϱvz, uᵣ.ϱe

    ϱ   = (ϱₗ + ϱᵣ) / 2
    ϱvx = (ϱvxₗ + ϱvxᵣ) / 2
    ϱvy = (ϱvyₗ + ϱvyᵣ) / 2
    ϱvz = (ϱvzₗ + ϱvzᵣ) / 2
    ϱe  = (ϱeₗ + ϱeᵣ) / 2
    vx  = (vxₗ + vxᵣ) / 2
    vy  = (vyₗ + vyᵣ) / 2
    vz  = (vzₗ + vzᵣ) / 2
    p   = (pₗ + pᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱvx*vy
    fϱvy = ϱvy*vy + p
    fϱvz = ϱvz*vy
    fϱe  = (ϱe+p)*vy

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::DucrosEtAlFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱvxₗ, ϱvyₗ, ϱvzₗ, ϱeₗ = uₗ.ϱvx, uₗ.ϱvy, uₗ. ϱvz, uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱvxᵣ, ϱvyᵣ, ϱvzᵣ, ϱeᵣ = uᵣ.ϱvx, uᵣ.ϱvy, uᵣ.ϱvz, uᵣ.ϱe

    ϱ   = (ϱₗ + ϱᵣ) / 2
    ϱvx = (ϱvxₗ + ϱvxᵣ) / 2
    ϱvy = (ϱvyₗ + ϱvyᵣ) / 2
    ϱvz = (ϱvzₗ + ϱvzᵣ) / 2
    ϱe  = (ϱeₗ + ϱeᵣ) / 2
    vx  = (vxₗ + vxᵣ) / 2
    vy  = (vyₗ + vyᵣ) / 2
    vz  = (vzₗ + vzᵣ) / 2
    p   = (pₗ + pᵣ) / 2

    fϱ   = ϱ*vz
    fϱvx = ϱvx*vz
    fϱvy = ϱvy*vz
    fϱvz = ϱvz*vz + p
    fϱe  = (ϱe+p)*vz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

@inline function (fnum::KennedyGruberFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                            model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    eₗ = uₗ.ϱe / ϱₗ
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    eᵣ = uᵣ.ϱe / ϱᵣ

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    e  = (eₗ + eᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p
    fϱvy = ϱ*vx*vy
    fϱe  = (ϱ*e+p)*vx

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::KennedyGruberFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                            model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    eₗ = uₗ.ϱe / ϱₗ
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    eᵣ = uᵣ.ϱe / ϱᵣ

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    e  = (eₗ + eᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p
    fϱe  = (ϱ*e+p)*vy

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::KennedyGruberFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    eₗ = uₗ.ϱe / ϱₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    eᵣ = uᵣ.ϱe / ϱᵣ

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    e  = (eₗ + eᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p
    fϱvy = ϱ*vx*vy
    fϱvz = ϱ*vx*vz
    fϱe  = (ϱ*e+p)*vx

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::KennedyGruberFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    eₗ = uₗ.ϱe / ϱₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    eᵣ = uᵣ.ϱe / ϱᵣ

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    e  = (eₗ + eᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p
    fϱvz = ϱ*vy*vz
    fϱe  = (ϱ*e+p)*vy

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::KennedyGruberFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                            model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    eₗ = uₗ.ϱe / ϱₗ
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    eᵣ = uᵣ.ϱe / ϱᵣ

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    e  = (eₗ + eᵣ) / 2

    fϱ   = ϱ*vz
    fϱvx = ϱ*vx*vz
    fϱvy = ϱ*vy*vz
    fϱvz = ϱ*vz^2 + p
    fϱe  = (ϱ*e+p)*vz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

@inline function (fnum::PirozzoliFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱeₗ = uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱeᵣ = uᵣ.ϱe

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    h  = ((ϱeₗ+pₗ)/ϱₗ + (ϱeᵣ+pᵣ)/ϱᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p
    fϱvy = ϱ*vx*vy
    fϱe  = ϱ*h*vx

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end
@inline function (fnum::PirozzoliFlux)(uₗ::EulerVar2D{T}, uᵣ::EulerVar2D{T},
                                        model::Euler{T,2}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, pₗ = primitive_variables(uₗ, model)
    ϱeₗ = uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱeᵣ = uᵣ.ϱe

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    h  = ((ϱeₗ+pₗ)/ϱₗ + (ϱeᵣ+pᵣ)/ϱᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p
    fϱe  = ϱ*h*vy

    SVector(fϱ, fϱvx, fϱvy, fϱe)
end

@inline function (fnum::PirozzoliFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:x}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱeₗ = uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱeᵣ = uᵣ.ϱe

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    h  = ((ϱeₗ+pₗ)/ϱₗ + (ϱeᵣ+pᵣ)/ϱᵣ) / 2

    fϱ   = ϱ*vx
    fϱvx = ϱ*vx^2 + p
    fϱvy = ϱ*vx*vy
    fϱvz = ϱ*vx*vz
    fϱe  = ϱ*h*vx

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::PirozzoliFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:y}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱeₗ = uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱeᵣ = uᵣ.ϱe

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    h  = ((ϱeₗ+pₗ)/ϱₗ + (ϱeᵣ+pᵣ)/ϱᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vy
    fϱvy = ϱ*vy^2 + p
    fϱvz = ϱ*vy*vz
    fϱe  = ϱ*h*vy

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end
@inline function (fnum::PirozzoliFlux)(uₗ::EulerVar3D{T}, uᵣ::EulerVar3D{T},
                                        model::Euler{T,3}, dir::Val{:z}) where T
    @unpack γ = model
    ϱₗ, vxₗ, vyₗ, vzₗ, pₗ = primitive_variables(uₗ, model)
    ϱeₗ = uₗ.ϱe
    ϱᵣ, vxᵣ, vyᵣ, vzᵣ, pᵣ = primitive_variables(uᵣ, model)
    ϱeᵣ = uᵣ.ϱe

    ϱ  = (ϱₗ + ϱᵣ) / 2
    vx = (vxₗ + vxᵣ) / 2
    vy = (vyₗ + vyᵣ) / 2
    vz = (vzₗ + vzᵣ) / 2
    p  = (pₗ + pᵣ) / 2
    h  = ((ϱeₗ+pₗ)/ϱₗ + (ϱeᵣ+pᵣ)/ϱᵣ) / 2

    fϱ   = ϱ*vy
    fϱvx = ϱ*vx*vz
    fϱvy = ϱ*vy*vz
    fϱvz = ϱ*vz^2 + p
    fϱe  = ϱ*h*vz

    SVector(fϱ, fϱvx, fϱvy, fϱvz, fϱe)
end

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

struct EulerRiemannSolution{T,T1} <: AbstractRiemannSolution
    prob::RiemannProblem{Euler{T,1},EulerVar1D{T},T1}
    uₘ⁻::EulerVar1D{T} # left middle state
    uₘ⁺::EulerVar1D{T} # left middle state
    σₗ⁻::T # slow speed of left wave
    σₗ⁺::T # fast speed of left wave
    σₘ::T  # speed of the middle (rarefaction) wave
    σᵣ⁻::T # slow speed of second family
    σᵣ⁺::T # fast speed of second family

    function EulerRiemannSolution{T,T1}(prob::RiemannProblem{Euler{T,1},EulerVar1D{T},T1}) where {T,T1}
        uₘ⁻, uₘ⁺, σₗ⁻, σₗ⁺, σₘ, σᵣ⁻, σᵣ⁺ = compute_state_and_speeds(prob.uₗ, prob.uᵣ, prob.model)
        new(prob, uₘ⁻, uₘ⁺, σₗ⁻, σₗ⁺, σₘ, σᵣ⁻, σᵣ⁺)
    end
end

function solve(prob::RiemannProblem{Euler{T,1},EulerVar1D{T},T1}) where {T,T1}
    EulerRiemannSolution{T,T1}(prob)
end

function minmax_speeds(sol::EulerRiemannSolution)
    sol.σₗ⁻, sol.σᵣ⁺
end

function (sol::EulerRiemannSolution)(ξ::Real)
    @unpack uₘ⁻, uₘ⁺, σₗ⁻, σₗ⁺, σₘ, σᵣ⁻, σᵣ⁺ = sol
    @unpack uₗ, uᵣ, model = sol.prob

    if ξ <= σₗ⁻
        uₗ
    elseif ξ <= σₗ⁺
        rarefaction_wave_1(ξ, uₗ, model)
    elseif ξ <= σₘ
        uₘ⁻
    elseif ξ <= σᵣ⁻
        uₘ⁺
    elseif ξ <= σᵣ⁺
        rarefaction_wave_2(ξ, uᵣ, model)
    else # σᵣ⁺ < ξ
        uᵣ
    end
end

function rarefaction_wave_1(ξ, uₗ::EulerVar1D, model::Euler{T,1}) where {T}
    @unpack γ = model
    ϱₗ, vₗ, pₗ = primitive_variables(uₗ, model)
    aₗ = sqrt(γ * pₗ / ϱₗ)
    # Toro (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics,
    # equation (4.56)
    ϱ = ϱₗ * ( 2/(γ+1) + (γ-1)/(aₗ*(γ+1))*(vₗ-ξ) )^(2/(γ-1))
    v = (aₗ + (γ-1)/2*vₗ + ξ) * 2/(γ+1)
    p = pₗ * ( 2/(γ+1) + (γ-1)/(aₗ*(γ+1))*(vₗ-ξ) )^(2γ/(γ-1))
    conserved_variables(ϱ, v, p, model)
end

function rarefaction_wave_2(ξ, uᵣ::EulerVar1D, model::Euler{T,1}) where {T}
    @unpack γ = model
    ϱᵣ, vᵣ, pᵣ = primitive_variables(uᵣ, model)
    aᵣ = sqrt(γ * pᵣ / ϱᵣ)
    # Toro (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics,
    # equation (4.63)
    ϱ = ϱᵣ * ( 2/(γ+1) - (γ-1)/(aᵣ*(γ+1))*(vᵣ-ξ) )^(2/(γ-1))
    v = (-aᵣ + (γ-1)/2*vᵣ + ξ) * 2/(γ+1)
    p = pᵣ * ( 2/(γ+1) - (γ-1)/(aᵣ*(γ+1))*(vᵣ-ξ) )^(2γ/(γ-1))
    conserved_variables(ϱ, v, p, model)
end

function (sol::EulerRiemannSolution)(t::Real, x::Real)
  @unpack x₀, t₀ = sol.prob

  sol((x-x₀)/(t-t₀))
end

function compute_state_and_speeds(uₗ::EulerVar1D{T}, uᵣ::EulerVar1D{T}, model::Euler) where {T}
    @unpack γ = model
    ϱₗ, vₗ, pₗ = primitive_variables(uₗ, model)
    aₗ = sqrt(γ * pₗ / ϱₗ)
    ϱᵣ, vᵣ, pᵣ = primitive_variables(uᵣ, model)
    aᵣ = sqrt(γ * pᵣ / ϱᵣ)

    pₘ, vₘ = compute_pₘ_vₘ(ϱₗ, vₗ, pₗ, aₗ, ϱᵣ, vᵣ, pᵣ, aᵣ, γ)

    if pₘ > pₗ
        # left shock wave
        ϱₘ⁻ = ϱₗ * (pₘ/pₗ + (γ-1)/(γ+1)) / (1 + (γ-1)/(γ+1)*pₘ/pₗ)
        σₗ⁻ = σₗ⁺ = vₗ - aₗ*sqrt((γ+1)/2γ * pₘ/pₗ + (γ-1)/2γ)
    else
        # left rarefaction wave
        ϱₘ⁻ = ϱₗ * (pₘ/pₗ)^(1/γ)
        aₘ⁻ = aₗ * (pₘ/pₗ)^((γ-1)/2γ)
        σₗ⁻ = vₗ - aₗ
        σₗ⁺ = vₘ - aₘ⁻
    end

    if pₘ > pᵣ
        # right shock wave
        ϱₘ⁺ = ϱᵣ * (pₘ/pᵣ + (γ-1)/(γ+1)) / (1 + (γ-1)/(γ+1)*pₘ/pᵣ)
        σᵣ⁻ = σᵣ⁺ = vᵣ + aᵣ*sqrt((γ+1)/2γ * pₘ/pᵣ + (γ-1)/2γ)
    else
        # right rarefaction wave
        ϱₘ⁺ = ϱᵣ * (pₘ/pᵣ)^(1/γ)
        aₘ⁺ = aᵣ * (pₘ/pᵣ)^((γ-1)/2γ)
        σᵣ⁻ = vₘ + aₘ⁺
        σᵣ⁺ = vₗ + aₗ
    end

    uₘ⁻ = conserved_variables(ϱₘ⁻, vₘ, pₘ, model)
    uₘ⁺ = conserved_variables(ϱₘ⁺, vₘ, pₘ, model)
    σₘ = vₘ

    uₘ⁻, uₘ⁺, σₗ⁻, σₗ⁺, σₘ, σᵣ⁻, σᵣ⁺
end

function compute_pₘ_vₘ(ϱₗ, vₗ, pₗ, aₗ, ϱᵣ, vᵣ, pᵣ, aᵣ, γ)
    Aₗ = 2/((γ+1)*ϱₗ)
    Bₗ = (γ-1)/(γ+1)*pₗ
    Aᵣ =  2/((γ+1)*ϱᵣ)
    Bᵣ = (γ-1)/(γ+1)*pᵣ

    fₗ = p -> p > pₗ ? (p-pₗ)*sqrt(Aₗ/(p+Bₗ)) : 2aₗ/(γ-1) * ((p/pₗ)^((γ-1)/2γ) - 1)
    fᵣ = p -> p > pᵣ ? (p-pᵣ)*sqrt(Aᵣ/(p+Bᵣ)) : 2aᵣ/(γ-1) * ((p/pᵣ)^((γ-1)/2γ) - 1)

    # two rarefaction approximation
    # Toro (2009), Riemann Solvers and Numerical Methods for Fluid Dynamics,
    # equation (4.46)
    p₀ = ( (aₗ + aᵣ - (γ-1)*(vᵣ-vₗ)/2) / (aₗ/pₗ^((γ-1)/2γ) + aᵣ/pᵣ^((γ-1)/2γ) ) )^(2γ/(γ-1))
    # TODO: Improve; Newton with explicit derivative?
    pₘ = Roots.fzero(p -> fₗ(p)+fᵣ(p)+vᵣ-vₗ, p₀)

    pₘ, (vₗ+vᵣ)/2 + (fᵣ(pₘ) - fₗ(pₘ))/2
end


text:
using LinearAlgebraicRepresentation, ViewerGL
Lar = LinearAlgebraicRepresentation; GL = ViewerGL
using QHull, LinearAlgebra, DataStructures

function makecycle(points,edges_list)
    edgedict = SortedDict(zip(edges_list[:,1],edges_list[:,2]))
    polygon = [1]
    while true
        h = polygon[end]
        push!(polygon,edgedict[h])
        if polygon[end] == polygon[1] break end
    end
    return polygon # as ordered sequence of vertex indices
end

function vectors(V,polygon)
    vect = Array{Array{Float64,1},1}(undef,length(polygon))
    vect[1] = [0.0,0.0]
    for k=2:length(polygon)
        vect[k] = (V[polygon[k],:] - V[polygon[k-1],:])
    end
    return vect
end

function paramsum(vect)
    out = Array{Float64}(undef,length(vect))
    out[1] = 0.0
    for k=2:length(vect)
        out[k] = out[k-1] + norm(vect[k])
    end
    return out ./ out[end]
end

# problematic dataset
points = [-0.709583 -0.632635; 0.520341 0.561664; 0.918992 1.04452; 0.571694 0.592953; -0.120997 -0.0175708; -0.595547 -0.530965; 0.277832 0.413903; -0.980521 -0.823759]
points_map = [1, 2, 3, 4, 5, 6, 7, 8]
edges_list = [1 6; 3 8; 7 2; 2 4; 4 3; 5 7; 6 5; 8 1]

# problematic polygon
npoints = size(points,1)
polygon = makecycle(points,edges_list)
vects = vectors(points,polygon)
params = paramsum(vects)
vdict = Dict(zip([points[k,:] for k=1:npoints], 1:npoints))

# convex hull polygon
ch = QHull.chull(points)
chverts = ch.vertices
outverts = setdiff(1:npoints, chverts)
# place the outverts on the convex hupp polygon

ext: jl


text:
"""

#### References

Two-stage extreme learning machine for regression.

Yuan Lan, Yeng Chai Soh, and Guang-Bin Huang.

Neurocomputing, 2010 vol. 73 (16-18) pp. 3028-3038.

http://linkinghub.elsevier.com/retrieve/pii/S0925231210003401
"""
mutable struct TSELM{TN,TV} <: AbstractSLFN where TN<:AbstractNodeInput where TV<:AbstractArray{Float64}
    p::Int  # Number of training points
    q::Int  # Dimensionality of function domai
    s::Int  # maximum number of neurons
    ngroup::Int  # number of groups
    npg::Int  # nodes per group
    neuron_type::TN
    μx::Vector{Float64}
    σx::Vector{Float64}
    Wt::Matrix{Float64}
    d::Vector{Float64}
    v::TV

    function TSELM(p::Int, q::Int, s::Int, ngroup::Int, npg::Int,
                   neuron_type::TN, μx, σx, y::TV) where TN where TV
        new{TN,TV}(p, q, s, ngroup, npg, neuron_type, μx, σx)
    end
end

function TSELM(x::AbstractArray, y::TV;
               neuron_type::TN=Linear(Tanh()),
               s::Int=size(x, 1),
               ngroup::Int=5,
               npg::Int=ceil(Int, size(x, 1)/10),
               reg::AbstractLinReg=LSSVD()) where TN<:AbstractNodeInput where TV<:AbstractArray
    q = size(x, 2)  # dimensionality of function domain
    p = size(x, 1)  # number of training points
    s = min(p, s)  # can't have more neurons than obs
    npg = min(ceil(Int, s/4), npg)  # ensure at least 4 groups
    xn, μx, σx = standardize(x[:, :])
    out = TSELM(p, q, s, ngroup, npg, neuron_type, μx, σx, y)
    fit!(out, xn, y, reg)
end

## helper methods
# take every other data point so that when a is something like
# a linspace we still get decent coverage of the whole domain
_split_data(a::AbstractVector) = (a[1:2:end], a[2:2:end])
_split_data(a::AbstractMatrix) = (a[1:2:end, :], a[2:2:end, :])

function forward_selection!(elm::TSELM, x::AbstractArray, y::AbstractVector,
                            reg::AbstractLinReg)
    # split data sets -- but only when we have enough data for that to make
    # sense
    should_split = size(x, 1) >= 20
    if should_split
        train_x, validate_x = _split_data(x)
        train_y, validate_y = _split_data(y)
    else
        train_x, validate_x = x, x
        train_y, validate_y = y, y
    end
    N = size(train_x, 1)
    N_validate = size(validate_x, 1)

    # initialize
    L = 0
    R = zeros(Float64, N, N) + I
    W = zeros(elm.q, elm.s + elm.npg)
    d = zeros(elm.s + elm.npg)

    while L < elm.s
        L += elm.npg
        inds = (L-elm.npg)+1:L

        # variables to hold max ΔJ and corresponding δH for this
        # batch of groups. Define empty matrices to get mutable struct stability
        ΔJ = -Inf
        ΔR = zeros(0, 0)
        Ak = zeros(0, 0)
        dk = zeros(0)

        for i in 1:elm.ngroup
            # Step 1: randomly generate hidden parameters for this group
            wiT = 2*rand(elm.q, elm.npg) .- 1  # uniform [-1, 1]
            di = rand(elm.npg)                # uniform [0, 1]

            # Step 2: generate the hidden output for this group
            δHi = hidden_out(elm, train_x, wiT, di)

            # Step 3: compute the contribution of this group to cost function
            # TODO: sometimes (δHi'R*δHi) is singular. I can loop over above
            #       until it works
            ΔRi = (R*δHi)*(pinv(δHi'R*δHi) * (δHi'R'))
            ΔJi = dot(train_y, ΔRi*train_y)

            # Step 4: keep this group if it is the best we've seen so far
            if ΔJi > ΔJ
                ΔR = ΔRi
                Ak = wiT
                dk = di
            end
        end

        # Step 5: Update hidden node parameters (A, b), R matrix, H
        W[:, inds] = Ak
        d[inds] = dk
        R -= ΔR
    end

    # Step 6: Find the optimal number of neurons pstar
    pstar = 0
    fpe_min = Inf
    for p in elm.npg:elm.npg:elm.s
        Hp = hidden_out(elm, validate_x, W[:, 1:p], d[1:p])
        βp = regress(reg, Hp, validate_y)
        SSEp = should_add_intercept(reg) ?
            norm(validate_y - (Hp*βp[2:end] .+ βp[1]), 2) :
            norm(validate_y - (Hp*βp), 2)
        fpe = SSEp/N_validate * (N_validate+p)/(N_validate-p)

        if fpe < fpe_min
            fpe_min = fpe
            pstar = p
        end
    end

    # Step 7: Rebuild the net with the selected neurons and fit
    #         with entire training set
    elm.Wt = W[:, 1:pstar]
    elm.d = d[1:pstar]
    H = hidden_out(elm, x, elm.Wt, elm.d)
    elm.v = regress(reg, H, y)

    elm, H
end

function backward_elimination!(elm::TSELM, x::AbstractArray, y::AbstractVector,
                               H::AbstractMatrix)
    # TODO: come back to this
    return
    pstar = size(elm.Wt, 2)
    Hr = copy(H)
    Hpstar = H[:, pstar]
    N = size(H, 1)

    y = H*elm.v

    # Step 1: Compute press1 for all nodes
    press1 = zeros(Float64, pstar -1)
    for k in pstar-1:-1:1
        copy!(Hr, H)
        Hk = view(H, :, k)
        Hr[:, pstar] = Hk
        Hr[:, k] = Hpstar

        ϵ = y - Hr*elm.v
        M = inv(Hr'Hr)

        for i in 1:N
            hri = Hr[i, :]
            ϵi = ϵ[i] / (1 - dot(hri, M*hri))
            press1[k] += ϵi^2
        end
        press1[k] /= N

    end

end

function fit!(elm::TSELM, x::AbstractArray, y::AbstractVector,
              reg::AbstractLinReg)
    elm, H = forward_selection!(elm, x, y, reg)
    backward_elimination!(elm, x, y, H)
    elm
end

function Base.show(io::IO, elm::TSELM{TA,TN}) where TA where TN
    s =
    """
    TSELM with
      - $(TA) Activation function
      - $(TN) Neuron type
      - $(elm.q) input dimension(s)
      - $(size(elm.Wt, 2)) neuron(s)
      - $(elm.p) training point(s)
      - Algorithm parameters:
          - $(elm.ngroup) trials per group
          - $(elm.npg) neurons per group
          - $(elm.s) max neurons
    """
    print(io, s)
end


text:
abstract QueryNode
abstract JoinNode <: QueryNode
typealias QueryArg Union{Symbol, Expr}

"""
    DataNode <: QueryNode

Represent a data source in a `QueryNode` manipulation graph.

Note: the "leaves" of any `QueryNode` graph must be `DataNode`s.
"""
immutable DataNode{T} <: QueryNode
    input::T
end

immutable Node{V} <: QueryNode
    inputs::Tuple{Vararg{QueryNode}}
    args::Vector{QueryArg}
    dos::Vector{V} # NOTE: the plural of "do" (as a noun...)
end


text:
__precompile__()

module Ensemblearn

include("util.jl")
include("decisiontree.jl")
include("randomforest.jl")

Ensemblearn

end # module


text:
using EconUtils
using Test
using DataFrames

df = DataFrames(x = 1:5)

# write your own tests here
@test df == DataFrames(x = 1:5)


text:
using Pkg
using Documenter

push!(LOAD_PATH,"../src/")

Pkg.activate(".")
using EcologicalNetworks

makedocs(
    sitename = "EcologicalNetworks.jl",
    authors = "Timothée Poisot",
    modules = [EcologicalNetworks],
    pages = [
        "Index" => "index.md",
        "Contents" => "toc.md",
        "Interface" => [
            "Types" => "interface/types.md",
            "Conversions" => "interface/conversions.md",
            "Core functions" => "interface/highlevel.md"
        ],
        "Network measures" => [
            "Links" => "properties/links.md",
            "Modularity" => "properties/modularity.md",
            "Nestedness" => "properties/nestedness.md",
            "Motifs" => "properties/motifs.md",
            "Centrality and paths" => "properties/paths.md",
            "Overlap and similarity" => "properties/overlap.md",
            "Null models" => "properties/nullmodels.md",
            "Beta-diversity" => "properties/betadiversity.md"
        ]
    ]
)

deploydocs(
    deps   = Deps.pip("pygments", "python-markdown-math"),
    repo   = "github.com/PoisotLab/EcologicalNetworks.jl.git",
    devbranch = "develop"
)


text:
using domainUtils

url="http://lp-sandbox.ekstep.org:8080/taxonomy-service/v2/analytics/domain/map"

cd("/Users/soma/github/ekStep/Learning-Platform-Analytics/platform-scripts/domain-model")

# try reading the domainModel
status,domainModelData = domainUtils.readDomainModel(url);

if(!status)
  error("Aborting: can not process domain model")
end

# read the concepts
concepts, conceptHash = domainUtils.getConceptsFromConceptMap(domainModelData);

# read the parent-of relationship graph
A = domainUtils.getConceptGraph(domainModelData,"isParentOf");

# from the parent-child graph, get concept-concept similarity matrix
C1 = domainUtils.getConceptSimilarityByRelation(cGraph)

# get concept-grade coverage matrix
cGraph,grades,gradeHash,concepts,conceptHash = domainUtils.getConceptGradeCoverageVec(domainModelData);

# try computing grade similarity
x1 = [1 0 0 1 200]
x2 = [1 1 1 1 0]
x3 = [1 0 1 1 1]

w = [0 1 2 3 4;1 0 1 2 3;2 1 0 1 2;3 2 1 0 1;4 3 2 1 0]
wext = [1 1 1 1 1]
println("Source: ",x)
println("Target: ",y)
status, d1, mdl, flow, cost = domainUtils.getMinCostDistance(x1,x2,w,wext,false);
println("dist: ",d1)
status, d2, mdl, flow, cost = domainUtils.getMinCostDistance(x2,x3,w,wext,false);
println("dist: ",d2)
status, d3, mdl, flow, cost = domainUtils.getMinCostDistance(x1,x3,w,wext,false);
println("dist: ",d3)
println("triangle: ",[d1+d2,d3])

# read content and get its attributes
url="http://lp-sandbox.ekstep.org:8080/taxonomy-service/v2/content/org.ekstep.money.worksheet"
using Requests
url="http://lp-sandbox.ekstep.org:8080/taxonomy-service/v2/content/org.ekstep.daysoftheweek"
resp = get(url)

max_stars_repo_name: harshavardhanc/sunbird-analytics

max_issues_repo_name: harshavardhanc/sunbird-analytics

max_forks_repo_name: harshavardhanc/sunbird-analytics


text:
module BeamformerRecipes

export BeamformerRecipe
export DimInfo
export TelInfo
export ObsInfo
export CalInfo
export BeamInfo
export DelayInfo
export to_hdf5

abstract type Info end

include("diminfo.jl")
include("telinfo.jl")
include("obsinfo.jl")
include("calinfo.jl")
include("beaminfo.jl")
include("delayinfo.jl")

Base.Base.@kwdef mutable struct BeamformerRecipe
    diminfo::Union{DimInfo, Nothing} = nothing
    telinfo::Union{TelInfo, Nothing} = nothing
    obsinfo::Union{ObsInfo, Nothing} = nothing
    calinfo::Union{CalInfo, Nothing} = nothing
    beaminfo::Union{BeamInfo, Nothing} = nothing
    delayinfo::Union{DelayInfo, Nothing} = nothing
end # mutable struct BeamformerRecipe

include("hdf5utils.jl")

function to_hdf5(h5::HDF5.File, recipe::BeamformerRecipe)
    for name in propertynames(recipe)
        to_hdf5(h5, getproperty(recipe, name))
    end
end

function to_hdf5(h5name::AbstractString, recipe::BeamformerRecipe)
    # "cw" = read-write, create file if not existing, preserve existing contents
    h5open(h5name, "cw") do h5
        to_hdf5(h5, recipe)
    end
end

end # module BeamformerRecipes


text:
# Flags represented as types to dispatch on
abstract type FunctionType end
struct Dynamics <: FunctionType end
struct Cost <: FunctionType end
struct Algebraic <: FunctionType end
struct AllFunctions <: FunctionType end

abstract type JacobianType end
struct State <: JacobianType end
struct Control <: JacobianType end
struct Time <: JacobianType end
struct Static <: JacobianType end
struct AllJacobians <: JacobianType end


text:
#some random matrices for testing
module Randm
export randsym, randposdef, randstable, randnormal, randorth, randunitary

#%  .. currentmodule:: Randm
#%

#%  Module Randm
#%  ------------
#%
#%
#%  Introduction
#%  ~~~~~~~~~~~~
#%

#%  Random matrices for testing purposes. I did not figure out the actual distributions
#%  the matrices are drawn from.
#%

#%  Reference
#%  ~~~~~~~~~
#%

##%  .. function:: randsym(d)
##%
##%  	Random symmetric matrix. Note: randsym
##%

#function randsym(d)
#	x = 0.5*randn(d, d)
#	x' + x
#end

#%  .. function:: randposdef(d)
#%
#%  	Random positive definite matrix of dimension ``d``.
#%

function randposdef(d)
	x = randn(d, d)
	x'* x
end

#%  .. function:: randstable(d)
#%
#%  	Random stable matrix (matrix with eigenvalues with negative real part) with
#%  	dimension ``d``.

function randstable(d)
	# positive definite matrix
	x = randn(d, d)
	a = x'*x

	# anti symmetric matrix
	above_diag =	[ i<j ? 1 : 0 for i in 1:d, j in 1:d] # nonzero pattern for above diagonal
	y = randn(d, d)
	b = above_diag.*y - above_diag'.*y'

	# return stable matrix
	b - a
end

#%  .. function:: randunitary(d)
#%
#%  	Random unitary matrix of dimension ``d``.
#%

randunitary(d) = expm(im* randposdef(d))

#%  .. function:: randorth(d)
#%
#%  	Orthogonal matrix drawn according to the Haar measure on the group of orthogonal matrices.
#%

randorth(d) = qr(randn(d,d))[1]

#%  .. function:: randnormal(d)
#%
#%  	Random normal matrix.
#%

function tridiag(di, u, l)
	assert(length(di) -1 == length(u) == length(l))
	M = diagm(di)
	for i in 1:length(u)
	  M[i,i+1] = l[i]
  	  M[i+1,i] = u[i]
	end
	M

end
tridiag(alpha, beta) = tridiag(alpha, beta, -beta)

function randnormal(d)
	Q = randorth(d)
	B = schur(randn(d,d))[1] #same eigenvalues as a randn(d,d)
	m = d
	alpha = zeros(d)
	beta = zeros(d-1)
	while m > 1
	   s = abs(B[m-1,m-1]) + abs(B[m,m])
	   if s + abs(B[m,m-1]) > s # if significant offdiagonal value: evaluate submatrix m-1:m;
		spur = B[m,m] + B[m-1,m-1]
		dis2 = abs2(B[m,m]- B[m-1,m-1]) + 4.0*B[m,m-1]*B[m-1,m] # = spur^2 - 4det
		alpha[m] = alpha[m-1] = 0.5*spur

		beta[m-1] = 0.5*sqrt(-dis2)
		m -= 2
	   else
	        alpha[m] = B[m,m]
	        beta[m-1] = 0
	        m -= 1
	   end

	end
	if (m == 1)
		alpha[1] = B[1,1]
	end

	Q' * tridiag(alpha, beta) * Q
	#tridiag(alpha, beta)

end

end #module


text:
function read_data(experiment)
    #= Read data from one of the experiments (exp={SPR,ET,EEG}) and do preprocessing:
    - Log-transform RTs
    - Add columns for previous-word covariates
    - Remove rejected data. Note that log_freq_prev is missing for words following punctuation and clitics, so these
      items will not be analysed even though they are not explicitly rejected.
    - Normalize quantitative covariates. Surprisal values will be normalized upon use
    =#
    print("Reading "* experiment * " data files\n")
    data = CSV.read("data_" * experiment * ".csv", delim='\t')

    if experiment=="SPR"
        data[:logRT]         = log.(data[:RT])
        data[:log_freq_prev] = [missing; data[1:end-1,:log_freq]]
        data[:nr_char_prev]  = [missing; data[1:end-1,:nr_char]]
        data[:logRT_prev]    = [missing; data[1:end-1,:logRT]]

        data[[:log_freq, :nr_char, :word_pos, :log_freq_prev, :nr_char_prev, :logRT_prev]] = zscore(data[[:log_freq, :nr_char, :word_pos, :log_freq_prev, :nr_char_prev, :logRT_prev]])
    elseif experiment=="ET"
        data[:logRTfirstpass] = log.(data[:RTfirstpass])
        data[:log_freq_prev]  = [missing; data[1:end-1,:log_freq]]
        data[:nr_char_prev]   = [missing; data[1:end-1,:nr_char]]

        data[[:log_freq, :nr_char, :word_pos, :log_freq_prev, :nr_char_prev]] = zscore(data[[:log_freq, :nr_char, :word_pos, :log_freq_prev, :nr_char_prev]])
    else
        data[[:log_freq, :nr_char, :word_pos, :baseline]] = zscore(data[[:log_freq, :nr_char, :word_pos, :baseline]])
    end

    data = data[.!data[:reject_data],:]
    data = data[.!data[:reject_word],:]
    for s = subjects[(subjects[:experiment].==experiment) .& (subjects[:reject]),:subj_nr]
        data = data[.!(data[:subj_nr].==s),:]
    end

    data
end

max_stars_repo_name: caurnhammer/gated-rnns

max_issues_repo_name: caurnhammer/gated-rnns

max_forks_repo_name: caurnhammer/gated-rnns


text:
using MultivariateStats
using Base.Test

srand(34568)

## data
n = 10
d = 5
X = randn(d, n)

# step-by-step kernel centralization
for K in [
    reshape(1.:12., 3, 4),
    reshape(1.:9., 3, 3),
    reshape(1.:12., 4, 3),
    rand(n,d),
    rand(d,d),
    rand(d,n) ]

    x, y = size(K)
    I1 = ones(x,x)/x
    I2 = ones(y,y)/y
    Z = K - I1*K - K*I2 + I1*K*I2

    KC = fit(MultivariateStats.KernelCenter, K)
    @test all(isapprox.(Z, MultivariateStats.transform!(KC, copy(K))))
end

# kernel calculations
K = MultivariateStats.pairwise((x,y)->norm(x-y), X, X[:,1:2])
@test size(K) == (n, 2)
@test K[1,1] == 0
@test K[3,2] == norm(X[:,3] - X[:,2])

K = MultivariateStats.pairwise((x,y)->norm(x-y), X[:,1:3], X)
@test size(K) == (3, n)
@test K[1,1] == 0
@test K[3,2] == norm(X[:,2] - X[:,3])

K = similar(X, n, n)
MultivariateStats.pairwise!(K, (x,y)->norm(x-y), X)
@test size(K) == (n, n)
@test K[1,1] == 0
@test K[2,1] == norm(X[:,2] - X[:,1])

KC = fit(MultivariateStats.KernelCenter, K)
Iₙ = ones(n,n)/n
@test MultivariateStats.transform!(KC, copy(K)) ≈ K - Iₙ*K - K*Iₙ + Iₙ*K*Iₙ

K = MultivariateStats.pairwise((x,y)->norm(x-y), X, X[:,1])
@test size(K) == (n, 1)
@test K[1,1] == 0
@test K[2,1] == norm(X[:,2] - X[:,1])

KC = fit(MultivariateStats.KernelCenter, K)
@test all(isapprox.(MultivariateStats.transform!(KC, copy(K)), 0.0, atol=10e-7))

## check different parameters
X = randn(d, n)
M = fit(KernelPCA, X, maxoutdim=d)
M2 = fit(PCA, X, method=:cov, pratio=1.0)
@test indim(M) == d
@test outdim(M) == d
@test abs.(transform(M, X)) ≈ abs.(transform(M2, X))
@test abs.(transform(M, X[:,1])) ≈ abs.(transform(M2, X[:,1]))

M = fit(KernelPCA, X, maxoutdim=3, solver=:eigs)
M2 = fit(PCA, X, method=:cov, maxoutdim=3)
@test indim(M) == d
@test outdim(M) == 3
@test abs.(transform(M, X)) ≈ abs.(transform(M2, X))
@test abs.(transform(M, X[:,1])) ≈ abs.(transform(M2, X[:,1]))

# issue #44
Y = randn(d, 2*n)
@test size(transform(M, Y)) == size(transform(M2, Y))

# reconstruction
@test_throws ArgumentError reconstruct(M, X)
M = fit(KernelPCA, X, inverse=true)
@test all(isapprox.(reconstruct(M, transform(M, X)), X, atol=0.75))

# use rbf kernel
γ = 10.
rbf=(x,y)->exp(-γ*norm(x-y)^2.0)
M = fit(KernelPCA, X, kernel=rbf)
@test indim(M) == d
@test outdim(M) == d

# use precomputed kernel
K = MultivariateStats.pairwise((x,y)->x'*y, X)
@test_throws AssertionError fit(KernelPCA, rand(1,10), kernel=nothing) # symmetric kernel
M = fit(KernelPCA, K, maxoutdim = 5, kernel=nothing, inverse=true) # use precomputed kernel
M2 = fit(PCA, X, method=:cov, pratio=1.0)
@test_throws ArgumentError reconstruct(M, X) # no reconstruction for precomputed kernel
@test abs.(transform(M)) ≈ abs.(transform(M2, X))

@test_throws ArgumentError fit(KernelPCA, rand(1,10), kernel=1)

# fit a Float32 matrix
X = randn(Float32, d, n)
M = fit(KernelPCA, X)
@test indim(M) == d
@test outdim(M) == d
@test eltype(transform(M, X[:,1])) == Float32

## fit a sparse matrix
X = sprandn(100d, n, 0.6)
M = fit(KernelPCA, X, maxoutdim=3, solver=:eigs)
@test indim(M) == 100d
@test outdim(M) == 3


text:
# Transformation of tensors

"""
    rot3(C, a, b, x) -> Cr

Return a rotated tensor `Cr`, the result of rotating in turn the tensor `C`:

1. About the x1 axis by `a` degrees clockwise when looking down the x1 axis
2.   "    "  x2   "   " `b`    "        "       "     "      "   "  x2   "
3.   "    "  x3   "   " `b`    "        "       "     "      "   "  x3   "

Rotations are performed in that order,
"""
rot3(C, a, b, c) = transform(C, RotXYZ(deg2rad(a), deg2rad(b), deg2rad(c)))

"""
    transform(C, M) -> C′

Apply a transformation matrix `M` to a 6x6 matrix `C`, returning the transformed matrix `C′`.
"""
function transform(C, M)
    K = _transformation_matrix(M)
    return K * (C * transpose(K))
end

"""
    _transformation_matrix(M) -> K

Return the matrix `K` which represents the matrix with which to transform a Voigt
6x6 matrix by the 3x3 transformation matrix `M`.
"""
@inline function _transformation_matrix(M)
    K = MMatrix{6, 6, Float64}(undef)
    @inbounds for i = 1:3
        i1 = (i + 1)%3
        i1 = i1 == 0 ? 3 : i1
        i2 = (i + 2)%3
        i2 = i2 == 0 ? 3 : i2
        for j = 1:3
            j1 = (j + 1)%3
            j1 = j1 == 0 ? 3 : j1
            j2 = (j + 2)%3
            j2 = j2 == 0 ? 3 : j2
            K[i,j] = M[i,j]^2
            K[i,j+3] = 2*M[i,j1]*M[i,j2]
            K[i+3,j] = M[i1,j]*M[i2,j]
            K[i+3,j+3] = M[i1,j1]*M[i2,j2] + M[i1,j2]*M[i2,j1]
        end
    end
    K
end


text:
using Interpolations

to11interval(x,a,b) = 2*(x - a)/(b-a)-1

years = to11interval.(2003:2016,-1,1)
year_space = LinRange(-1,1,14)

A = [log(x) for x in years]

# linear interpolation
interp_linear = LinearInterpolation(xs, A)
interp_linear(3) # exactly log(3)
interp_linear(3.1) # approximately log(3.1)

# cubic spline interpolation
interp_cubic = CubicSplineInterpolation(xs, A)
interp_cubic(3) # exactly log(3)
interp_cubic(3.1) # approximately log(3.1)


text:
module StrandbeestRobot

using RigidBodyDynamics
using RigidBodyDynamics.PDControl
using RigidBodyDynamics.Contact
const rbd = RigidBodyDynamics
using LightXML
using StaticArrays
using Optim
using LinearAlgebra
using RigidBodySim

urdfpath() = joinpath(@__DIR__, "..", "data", "Strandbeest.urdf")

function mechanism(::Type{T} = Float64;
                   floating = true,
                   contact_model = default_contact_model(),
                   add_flat_ground = true) where {T}
    mechanism = parse_urdf(T, urdfpath())
    add_loop_joints!(mechanism, urdfpath())

    if floating
        crossbar = findbody(mechanism, "crossbar")
        basejoint = joint_to_parent(crossbar, mechanism)
        floatingjoint = Joint("floating_base", frame_before(basejoint),
            frame_after(basejoint), QuaternionFloating{Float64}())
        replace_joint!(mechanism, basejoint, floatingjoint)
    end

    add_foot_contacts!(mechanism, contact_model)
    if add_flat_ground
        add_environment_primitive!(mechanism,
            HalfSpace3D(
                Point3D(root_frame(mechanism), 0., 0, 0),
                FreeVector3D(root_frame(mechanism), 0., 0, 1)))
    end

    mechanism
end

"""
Parse the Drake-specific `loop_joint` tags from a URDF and use them
to add non-tree joints to the given mechanism.
"""
function add_loop_joints!(mechanism::Mechanism{T}, urdf::AbstractString) where T
    doc = parse_file(urdf)
    xml_root = LightXML.root(doc)
    xml_loops = get_elements_by_tagname(xml_root, "loop_joint")
    for xml_loop in xml_loops
        name = attribute(xml_loop, "name")
        @assert attribute(xml_loop, "type") == "continuous"
        axis = SVector{3}(rbd.parse_vector(T, find_element(xml_loop, "axis"), "xyz", "1 0 0"))
        joint = Joint(name, Revolute(axis))
        xml_link1 = find_element(xml_loop, "link1")
        body1 = findbody(mechanism, attribute(xml_link1, "link"))
        H1 = Transform3D(frame_before(joint), default_frame(body1),
            rbd.parse_pose(T, xml_link1)...)
        xml_link2 = find_element(xml_loop, "link2")
        body2 = findbody(mechanism, attribute(xml_link2, "link"))
        H2 = Transform3D(frame_after(joint), default_frame(body2),
            rbd.parse_pose(T, xml_link2)...)
        attach!(mechanism, body1, body2, joint,
            joint_pose = H1,
            successor_pose = inv(H2))
    end
end

"""
Return a function which maps a configuration vector
for the given mechanism to the total squared error in
the alignment of the robot's non-tree joints.
"""
function loop_joint_error(mechanism::Mechanism)
    let statecache = StateCache(mechanism), resultcache = DynamicsResultCache(mechanism)
        function(q)
            state = statecache[eltype(q)]
            result = resultcache[eltype(q)]
            set_configuration!(state, q)
            zero_velocity!(state)
            gains = rbd.CustomCollections.ConstDict{JointID}(SE3PDGains(PDGains(1, 0), PDGains(1, 0)))
            rbd.constraint_bias!(result, state; stabilization_gains=gains)
            ret = result.constraintbias ⋅ result.constraintbias
        end
    end
end

default_contact_model() = SoftContactModel(hunt_crossley_hertz(k = 500e3), ViscoelasticCoulombModel(1.0, 20e3, 100.))

function add_foot_contacts!(mechanism::Mechanism, contactmodel=default_contact_model())
    for body in bodies(mechanism)
        if occursin("bars_g_h_i", string(body))
            frame = default_frame(body)
            point = Point3D(frame, 0, 0, 0.49)
            add_contact_point!(body, ContactPoint(
                point, contactmodel))

        end
    end
end

function solve_initial_configuration!(state::MechanismState{T}) where T
    # Find an initial state of the robot which correctly aligns all
    # the loop joints.
    mechanism = state.mechanism

    # We do this before adding the contact points to work around
    # https://github.com/JuliaRobotics/RigidBodyDynamics.jl/issues/483

    # First, initialize the state with a configuration within the
    # joint limits:
    lb = similar(configuration(state))
    ub = similar(configuration(state))
    for joint in tree_joints(mechanism)
        lb_joint = lb[joint]
        ub_joint = ub[joint]
        if isfloating(joint)
            # use identity transform for floating joint
            tf = one(Transform3D{T}, frame_after(joint), frame_before(joint))
            set_configuration!(state, joint, tf)
            # Optim can't handle equal bounds, so give it a bit of wiggle room
            lb_joint .= configuration(state, joint) .- 1e-3
            ub_joint .= configuration(state, joint) .+ 1e-3
        else
            lb_joint .= rbd.lower.(position_bounds(joint))
            ub_joint .= rbd.upper.(position_bounds(joint))
        end
    end
    q = configuration(state)
    for i in eachindex(q)
        if isfinite(lb[i]) && isfinite(ub[i])
            q[i] = (lb[i] + ub[i]) / 2
        else
            q[i] = clamp(0, lb[i], ub[i])
        end
    end
    setdirty!(state)

    # Use Optim's Fminbox to minimize the loop joint error within
    # the joint limits:
    cost = loop_joint_error(mechanism)
    result = Optim.optimize(cost, lb, ub,
        Vector(configuration(state)),
        Fminbox(LBFGS()), autodiff=:forward)

    # Verify that we've actually closed all the loops
    @assert Optim.minimum(result) < 1e-9

    set_configuration!(state, Optim.minimizer(result))
    normalize_configuration!(state)
    state
end

end


text:
import ..SamplingConstraints:
    SamplingConstraint,
    SequentialSamplingConstraint,
    OrderedSamplingAlgorithm,
    StartToEnd,
    StrictlyIncreasing,
    strictly_increasing_sequence_exists,
    TruncateMinimum,
    TruncateRange,
    TruncateQuantiles,
    constrain

import IntervalArithmetic:
    interval

"""
    resample(udata::DT, sequential_constraint::StrictlyIncreasing{T};
        quantiles = [0.0001, 0.9999]) where {DT <: AbstractUncertainValueDataset, T <: StartToEnd}

Element-wise resample the uncertain values in the dataset such that each preceding value
is strictly less in magnitude than the next one.

## Arguments:

- **`udata`**: An uncertain dataset.
- **`sequential_constraint`**: An instance of a `StrictlyIncreasing` sequential
    sampling constraint.
- **`ordered_sampling_alg`**: An instance of a `StartToEnd` ordered
    sampling constraint, indicating that the sequence of increasing values should
    be created in one go from start to finish (as opposed to chunking the data set
    first, then gluing partial sequences together afterwards).
- **`quantiles`**: A two-element vector representing a quantile range which is used
    to truncate the supports values in the dataset before drawing the sequence of
    increasing values. This deals with distributions with infinite support.
"""
function resample(udata::DT,
        sequential_constraint::StrictlyIncreasing{T};
        quantiles = [0.0001, 0.9999]) where {DT <: AbstractUncertainValueDataset, T <: StartToEnd}
    n_vals = length(udata)

    # First, constrain all data such that the supports of all furnishing distributions
    # are finite.
    c = TruncateQuantiles(quantiles...,)
    constrained_data = [udata[i] isa CertainValue ? udata[i] : constrain(udata[i], c) for i = 1:length(udata)]

    if !strictly_increasing_sequence_exists(udata, quantiles)
        throw(ArgumentError("No strictly increasing sequence through dataset exists."))
    end

    # Pre-allocate a sample
    sample = Vector{Float64}(undef, n_vals)

    # Sample the first value in a way that ensures a strictly
    # increasing sequence from indices 2:end will exist.
    lo = minimum(constrained_data[1])
    hi = minimum(maximum.(constrained_data[1:end]))

    if lo == hi
        sample[1] = resample(constrained_data[1])
    else
        sample[1] = resample(constrained_data[1], TruncateRange(lo, hi))
    end

    for i = 2:n_vals - 1
        # Find the lower and upper bounds of the support from which we can draw values
        # while still ensuring an increasing sequence of values.
        lo = sample[i - 1]
        hi = minimum([maximum(constrained_data[i]);
                        maximum.(constrained_data[(i + 1):end])])

        # Constrain the support of the furnishing distribution according
        # to the bounds and sample a value from it.
        cd = constrained_data[i]

        if cd isa CertainValue # no need to constrain if value is not uncertain
            sample[i] = resample(cd)
        else
            constrained_val = constrain(cd, TruncateRange(lo, hi))
            sample[i] = resample(constrained_val)
        end
    end

    # The last value doesn't need an upper bound, because there's no
    # value following it we need to consider.
    sample[n_vals] = resample(udata[n_vals], TruncateMinimum(sample[end - 1]))

    return sample
end

function resample(udata::DT,
        constraint::SCT,
        sequential_constraint::StrictlyIncreasing{T};
        quantiles = [0.0001, 0.9999]) where {SCT <: SamplingConstraint,
                DT <: AbstractUncertainValueDataset,
                T <: StartToEnd}

    # Truncate the support of each value in the dataset according to the provided
    # constraint(s) and then resample

    resample(constrain(udata, constraint), sequential_constraint, quantiles = quantiles)
end

function resample(udata::DT,
    constraint::Vector{SCT},
    sequential_constraint::StrictlyIncreasing{T};
    quantiles = [0.0001, 0.9999]) where {SCT <: SamplingConstraint,
        DT <: AbstractUncertainValueDataset,
         T <: StartToEnd}

    # Truncate the support of each value in the dataset according to the provided
    # constraint(s) and then resample
    resample(constrain(udata, constraint), sequential_constraint, quantiles = quantiles)
end

export resample

"""
    resample(udata::AbstractUncertainValueDataset,
        constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
        sequential_constraint::StrictlyIncreasing{OrderedSamplingAlgorithm};
        quantiles = [0.0001, 0.9999])

Draw a sequence of values strictly increasing in magnitude from the dataset, sampling
each of the furnishing distributions once each, after first truncating the supports
of the values in the dataset using the provided `constraint`s.

## Arguments:

- **`udata`**: An uncertain dataset.
- **`constraint`**: Sampling constraint(s) to apply to each of the values in the dataset
    before drawing the sequence of increasing values.
- **`sequential_constraint`**: An instance of a `StrictlyIncreasing` sequential
    sampling constraint. For example, `StrictlyIncreasing(StartToEnd())` indicates
    that a strictly increasing sequence should be created in one go from start to
    finish (as opposed to chunking the data set first, then gluing partial sequences
    together afterwards).
- **`quantiles`**: A two-element vector representing a quantile range which is used
    to truncate the supports values in the dataset before drawing the sequence of
    increasing values. This deals with distributions with infinite support.
"""
resample(udata::AbstractUncertainValueDataset,
    constraint::Union{SamplingConstraint, Vector{SamplingConstraint}},
    sequential_constraint::StrictlyIncreasing{OrderedSamplingAlgorithm};
    quantiles = [0.0001, 0.9999])

"""
function resample(udata::UncertainIndexValueDataset,
    sequential_constraint::StrictlyIncreasing,
    quantiles = [0.0001, 0.9999])

Resample an uncertain index-value dataset by enforcing strictly increasing indices.
"""
function resample(udata::UncertainIndexValueDataset,
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})
    resample(udata.indices, sequential_constraint), resample(udata.values)
end

function resample(udata::UncertainIndexValueDataset,
        constraint::SamplingConstraint,
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        constraint::Vector{SamplingConstraint},
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::SamplingConstraint,
        value_constraint::SamplingConstraint,
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::Vector{<:SamplingConstraint},
        value_constraint::SamplingConstraint,
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::SamplingConstraint,
        value_constraint::Vector{<:SamplingConstraint},
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end

function resample(udata::UncertainIndexValueDataset,
        idx_constraint::Vector{<:SamplingConstraint},
        value_constraint::Vector{<:SamplingConstraint},
        sequential_constraint::StrictlyIncreasing{<:OrderedSamplingAlgorithm})

    inds = resample(constrain(udata.indices, idx_constraint), sequential_constraint)
    vals = resample(constrain(udata.indices, value_constraint))

    inds, vals
end