text:
# ============================================================
# Gas Kinetic Flux Functions
# ============================================================

"""
$(SIGNATURES)

Gas kinetic flux
"""
function flux_gks!(
    face::Interface,
    ctrL::T,
    ctrR::T,
    gas::Gas,
    p,
    dt = 1.0,
) where {T<:ControlVolume}

    dxL, dxR = p[1:2]

    if size(ctrL.w, 1) == 3
        flux_gks!(
            face.fw,
            ctrL.w .+ dxL .* ctrL.sw,
            ctrR.w .- dxR .* ctrR.sw,
            gas.γ,
            gas.K,
            gas.μᵣ,
            gas.ω,
            dt,
            dxL,
            dxR,
            ctrL.sw,
            ctrR.sw,
        )
    else
    end

    return nothing

end

function flux_gks!(
    face::Interface,
    ctrL::T,
    ctrR::T,
    gas::Mixture,
    p,
    dt = 1.0,
) where {T<:ControlVolume}

    dxL, dxR = p[1:2]

    if size(ctrL.w, 1) == 3
        flux_gks!(
            face.fw,
            ctrL.w .+ dxL .* ctrL.sw,
            ctrR.w .- dxR .* ctrR.sw,
            gas.γ,
            gas.K,
            gas.mi,
            gas.ni,
            gas.me,
            gas.ne,
            gas.Kn[1],
            dt,
            dxL,
            dxR,
            ctrL.sw,
            ctrR.sw,
        )
    else
    end

    return nothing

end

"""
$(SIGNATURES)

Gas kinetic scalar flux

* @args: conservative scalars and their slopes
* @args: viscosity
* @args: time step and cell size
* @return: scalar flux
"""
function flux_gks(u, μ, dt, a = 0, su = 0.0)
    prim = ifelse(a == 0, conserve_prim(u), conserve_prim(u, a))

    Mu = gauss_moments(prim)[1]
    tau = 2.0 * μ

    fa = pdf_slope(u, su)
    Δ = -prim[1] * moments_conserve_slope(fa, Mu, 1)
    faT = pdf_slope(u, Δ)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = dt
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]

    # flux related to upwind distribution
    Muv = moments_conserve(Mu, 1)
    Mau = moments_conserve_slope(fa, Mu, 2)
    MauT = moments_conserve_slope(faT, Mu, 1)

    fw =
        Mt[4] * prim[1] * Muv - (Mt[5] + tau * Mt[4]) * prim[1] * Mau -
        tau * Mt[4] * prim[1] * MauT

    return fw / dt
end

"""
$(SIGNATURES)
"""
function flux_gks(
    uL::T,
    uR::T,
    μ,
    dt,
    dxL,
    dxR,
    a = 0,
    suL = 0.0,
    suR = 0.0,
) where {T<:Real}

    primL = ifelse(a == 0, conserve_prim(uL), conserve_prim(uL, a))
    primR = ifelse(a == 0, conserve_prim(uR), conserve_prim(uR, a))

    Mu1, MuL1, MuR1 = gauss_moments(primL)
    Mu2, MuL2, MuR2 = gauss_moments(primR)

    u = primL[1] * moments_conserve(MuL1, 0) + primR[1] * moments_conserve(MuR2, 0)
    prim = ifelse(a == 0, conserve_prim(u), conserve_prim(u, a))
    tau = 2.0 * μ + 1.0 * abs(uL - uR) / (abs(uL) + abs(uR)) * dt

    faL = pdf_slope(uL, suL)
    Δ = -primL[1] * moments_conserve_slope(faL, Mu1, 1)
    faTL = pdf_slope(uL, Δ)

    faR = pdf_slope(uR, suR)
    Δ = -primR[1] * moments_conserve_slope(faR, Mu2, 1)
    faTR = pdf_slope(uR, Δ)

    Mu, MuL, MuR = gauss_moments(prim)
    sw0L = (u - (uL - suL * dxL)) / dxL
    sw0R = ((uR + suR * dxR) - u) / dxR
    gaL = pdf_slope(u, sw0L)
    gaR = pdf_slope(u, sw0R)
    Δ =
        -prim[1] *
        (moments_conserve_slope(gaL, MuL, 1) + moments_conserve_slope(gaR, MuR, 1))
    #sw = (uR - uL) / (dxL + dxR)
    #ga = pdf_slope(u, sw)
    #Δ = -prim[1] .* moments_conserve_slope(ga, Mu, 1)
    gaT = pdf_slope(u, Δ)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, 1)
    MauL = moments_conserve_slope(gaL, MuL, 2)
    MauR = moments_conserve_slope(gaR, MuR, 2)
    #Mau = moments_conserve_slope(ga, MuR, 2)
    MauT = moments_conserve_slope(gaT, Mu, 1)

    fw = Mt[1] * prim[1] * Muv + Mt[2] * prim[1] * (MauL + MauR) + Mt[3] * prim[1] * MauT
    #fw = Mt[1] * prim[1] * Muv + Mt[2] * prim[1] * Mau + Mt[3] * prim[1] * MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, 1)
    MauL = moments_conserve_slope(faL, MuL1, 2)
    MauLT = moments_conserve_slope(faTL, MuL1, 1)

    MuvR = moments_conserve(MuR2, 1)
    MauR = moments_conserve_slope(faR, MuR2, 2)
    MauRT = moments_conserve_slope(faTR, MuR2, 1)

    fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    return fw

end

"""
$(SIGNATURES)

Gas kinetic Navier-Stokes flux

* @args: conservative variables and their left/right slopes
* @args: molecular and thermodynamic parameters
* @args: time step and cell size

Continuous case
"""
function flux_gks!(
    fw::X,
    w::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    sw = zero(w)::Y,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1}}

    prim = conserve_prim(w, γ)
    mus = gauss_moments(prim, inK)
    tau = vhs_collision_time(prim, μᵣ, ω)

    ax = pdf_slope(prim, sw, inK)
    if length(prim) == 3
        ∂ft = -prim[1] .* moments_conserve_slope(ax, mus[1], mus[2], 1)
        at = pdf_slope(prim, ∂ft, inK)

        Muv = moments_conserve(mus[1], mus[2], 1, 0)
        Mau = moments_conserve_slope(ax, mus[1], mus[2], 2)
        Mtu = moments_conserve_slope(at, mus[1], mus[2], 1)
    elseif length(prim) == 4
        ∂ft = -prim[1] .* moments_conserve_slope(ax, mus[1], mus[2], mus[3], 1, 0)
        at = pdf_slope(prim, ∂ft, inK)

        Muv = moments_conserve(mus[1], mus[2], mus[3], 1, 0, 0)
        Mau = moments_conserve_slope(ax, mus[1], mus[2], mus[3], 2, 0)
        Mtu = moments_conserve_slope(at, mus[1], mus[2], mus[3], 1, 0)
    elseif length(prim) == 5
        ∂ft = -prim[1] .* moments_conserve_slope(ax, mus[1], mus[2], mus[3], 1, 0, 0)
        at = pdf_slope(prim, ∂ft, inK)

        Muv = moments_conserve(mus[1], mus[2], mus[3], 1, 0, 0)
        Mau = moments_conserve_slope(ax, mus[1], mus[2], mus[3], 2, 0, 0)
        Mtu = moments_conserve_slope(at, mus[1], mus[2], mus[3], 1, 0, 0)
    end

    @. fw = prim[1] * (Muv - tau * Mau - tau * Mtu)

    return nothing

end

"""
$(SIGNATURES)

1D
"""
function flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mxi1, 1)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mxi2, 1)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mxi, 1) .+
            moments_conserve_slope(gaR, MuR, Mxi, 1)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mxi, 1)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mxi, 2)
    MauR = moments_conserve_slope(gaR, MuR, Mxi, 2)
    # Mau = moments_conserve_slope(ga, MuR, Mxi, 2)
    MauT = moments_conserve_slope(gaT, Mu, Mxi, 1)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mxi1, 1, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mxi1, 2)
    MauLT = moments_conserve_slope(faTL, MuL1, Mxi1, 1)

    MuvR = moments_conserve(MuR2, Mxi2, 1, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mxi2, 2)
    MauRT = moments_conserve_slope(faTR, MuR2, Mxi2, 1)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    return nothing

end

"""
$(SIGNATURES)

Gas mixture
"""
function flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    dxL::Real,
    dxR::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AA{<:FN,2},Y<:AA{<:Real,2}}

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    Mu1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mxi1, 0, 0)
    Muv2 = mixture_moments_conserve(MuR2, Mxi2, 0, 0)

    faL = mixture_pdf_slope(primL, swL, inK)
    mm = mixture_moments_conserve_slope(faL, Mu1, Mxi1, 1)
    sw = similar(swL)
    for j in axes(sw, 2)
        @. sw[:, j] = -primL[1, j] * mm[:, j]
    end
    faTL = mixture_pdf_slope(primL, sw, inK)

    faR = mixture_pdf_slope(primR, swR, inK)
    mm = mixture_moments_conserve_slope(faR, Mu2, Mxi2, 1)
    for j in axes(sw, 2)
        @. sw[:, j] = -primR[1, j] * mm[:, j]
    end
    faTR = mixture_pdf_slope(primR, sw, inK)

    w = similar(wL)
    for j in axes(w, 2)
        @. w[:, j] = primL[1, j] * Muv1[:, j] + primR[1, j] * Muv2[:, j]
    end
    prim = mixture_conserve_prim(w, γ)
    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    Mu, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = mixture_pdf_slope(prim, sw0L, inK)
    gaR = mixture_pdf_slope(prim, sw0R, inK)
    mmL = mixture_moments_conserve_slope(gaL, MuL, Mxi, 1)
    mmR = mixture_moments_conserve_slope(gaR, MuR, Mxi, 1)
    for j in axes(sw, 2)
        @. sw[:, j] = -prim[1, j] * (mmL[:, j] + mmR[:, j])
    end
    gaT = mixture_pdf_slope(prim, sw, inK)

    Mt = zeros(5, 2)
    for j in axes(Mt, 2)
        Mt[4, j] = tau[j] * (1.0 - exp(-dt / tau[j]))
        Mt[5, j] = -tau[j] * dt * exp(-dt / tau[j]) + tau[j] * Mt[4]
        Mt[1, j] = dt - Mt[4, j]
        Mt[2, j] = -tau[j] * Mt[1, j] + Mt[5, j]
        Mt[3, j] = 0.5 * dt^2 - tau[j] * Mt[1, j]
    end

    # flux related to central distribution
    Muv = mixture_moments_conserve(Mu, Mxi, 1, 0)
    MauL = mixture_moments_conserve_slope(gaL, MuL, Mxi, 2)
    MauR = mixture_moments_conserve_slope(gaR, MuR, Mxi, 2)
    MauT = mixture_moments_conserve_slope(gaT, Mu, Mxi, 1)
    for j in axes(fw, 2)
        @. fw[:, j] =
            Mt[1, j] * prim[1, j] * Muv[:, j] +
            Mt[2, j] * prim[1, j] * (MauL[:, j] + MauR[:, j]) +
            Mt[3, j] * prim[1, j] * MauT[:, j]
    end

    # flux related to upwind distribution
    MuvL = mixture_moments_conserve(MuL1, Mxi1, 1, 0)
    MauL = mixture_moments_conserve_slope(faL, MuL1, Mxi1, 2)
    MauLT = mixture_moments_conserve_slope(faTL, MuL1, Mxi1, 1)

    MuvR = mixture_moments_conserve(MuR2, Mxi2, 1, 0)
    MauR = mixture_moments_conserve_slope(faR, MuR2, Mxi2, 2)
    MauRT = mixture_moments_conserve_slope(faTR, MuR2, Mxi2, 1)

    for j in axes(fw, 2)
        @. fw[:, j] +=
            Mt[4, j] * primL[1, j] * MuvL[:, j] -
            (Mt[5, j] + tau[j] * Mt[4, j]) * primL[1, j] * MauL[:, j] -
            tau[j] * Mt[4, j] * primL[1, j] * MauLT[:, j] +
            Mt[4, j] * primR[1, j] * MuvR[:, j] -
            (Mt[5, j] + tau[j] * Mt[4, j]) * primR[1, j] * MauR[:, j] -
            tau[j] * Mt[4, j] * primR[1, j] * MauRT[:, j]
    end

    return nothing

end

"""
$(SIGNATURES)

1D1F1V
"""
function flux_gks!(
    fw::T1,
    fh::T2,
    wL::T3,
    wR::T3,
    u::T4,
    inK,
    γ,
    visRef,
    visIdx,
    dt,
    dxL,
    dxR,
    swL::T3,
    swR::T3,
) where {T1<:AA{<:FN,1},T2<:AA{<:FN,1},T3<:AA{<:Real,1},T4<:AA{<:FN,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mxi1, MuL1, _ = gauss_moments(primL, inK)
    Mu2, Mxi2, _, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, visRef, visIdx) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mxi1, 1)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mxi2, 1)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mxi, 1) .+
            moments_conserve_slope(gaR, MuR, Mxi, 1)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mxi, 1)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mxi, 2)
    MauR = moments_conserve_slope(gaR, MuR, Mxi, 2)
    # Mau = moments_conserve_slope(ga, MuR, Mxi, 2)
    MauT = moments_conserve_slope(gaT, Mu, Mxi, 1)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mxi1, 1, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mxi1, 2)
    MauLT = moments_conserve_slope(faTL, MuL1, Mxi1, 1)

    MuvR = moments_conserve(MuR2, Mxi2, 1, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mxi2, 2)
    MauRT = moments_conserve_slope(faTR, MuR2, Mxi2, 1)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    #--- fluxes of distribution functions ---#
    δ = heaviside.(u)
    HL = maxwellian(u, primL)
    HR = maxwellian(u, primR)
    H = maxwellian(u, prim)

    @. fh =
        Mt[1] * u * H +
        Mt[2] * u^2 * (gaL[1] * H + gaL[2] * u * H + 0.5 * gaL[3] * (u^2 * H)) * δ +
        Mt[2] * u^2 * (gaR[1] * H + gaR[2] * u * H + 0.5 * gaR[3] * (u^2 * H)) * (1.0 - δ) +
        Mt[3] * u * (gaT[1] * H + gaT[2] * u * H + 0.5 * gaT[3] * (u^2 * H)) +
        Mt[4] * u * HL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faL[1] * HL + faL[2] * u * HL + 0.5 * faL[3] * (u^2 * HL)) *
        δ + Mt[4] * u * HR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faR[1] * HR + faR[2] * u * HR + 0.5 * faR[3] * (u^2 * HR)) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (faTL[1] * HL + faTL[2] * u * HL + 0.5 * faTL[3] * (u^2 * HL)) *
        δ -
        tau *
        Mt[4] *
        u *
        (faTR[1] * HR + faTR[2] * u * HR + 0.5 * faTR[3] * (u^2 * HR)) *
        (1.0 - δ)

    return nothing

end

"""
$(SIGNATURES)

1D2F1V
"""
function flux_gks!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    wR::T3,
    u::T4,
    inK,
    γ,
    visRef,
    visIdx,
    dt,
    dxL,
    dxR,
    swL::T3,
    swR::T3,
) where {T1<:AA{<:FN,1},T2<:AA{<:FN,1},T3<:AA{<:Real,1},T4<:AA{<:FN,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mxi1, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mxi2, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, visRef, visIdx) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mxi1, 1)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mxi2, 1)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mxi, 1) .+
            moments_conserve_slope(gaR, MuR, Mxi, 1)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mxi, 1)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mxi, 2)
    MauR = moments_conserve_slope(gaR, MuR, Mxi, 2)
    # Mau = moments_conserve_slope(ga, MuR, Mxi, 2)
    MauT = moments_conserve_slope(gaT, Mu, Mxi, 1)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mxi1, 1, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mxi1, 2)
    MauLT = moments_conserve_slope(faTL, MuL1, Mxi1, 1)

    MuvR = moments_conserve(MuR2, Mxi2, 1, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mxi2, 2)
    MauRT = moments_conserve_slope(faTR, MuR2, Mxi2, 1)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    #--- fluxes of distribution functions ---#
    δ = heaviside.(u)

    HL = maxwellian(u, primL)
    BL = HL .* inK ./ (2.0 * primL[end])
    HR = maxwellian(u, primR)
    BR = HR .* inK ./ (2.0 * primR[end])

    H = maxwellian(u, prim)
    B = H .* inK ./ (2.0 * prim[end])

    @. fh =
        Mt[1] * u * H +
        Mt[2] * u^2 * (gaL[1] * H + gaL[2] * u * H + 0.5 * gaL[3] * (u^2 * H + B)) * δ +
        Mt[2] *
        u^2 *
        (gaR[1] * H + gaR[2] * u * H + 0.5 * gaR[3] * (u^2 * H + B)) *
        (1.0 - δ) +
        Mt[3] * u * (gaT[1] * H + gaT[2] * u * H + 0.5 * gaT[3] * (u^2 * H + B)) +
        Mt[4] * u * HL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faL[1] * HL + faL[2] * u * HL + 0.5 * faL[3] * (u^2 * HL + BL)) *
        δ + Mt[4] * u * HR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faR[1] * HR + faR[2] * u * HR + 0.5 * faR[3] * (u^2 * HR + BR)) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (faTL[1] * HL + faTL[2] * u * HL + 0.5 * faTL[3] * (u^2 * HL + BL)) *
        δ -
        tau *
        Mt[4] *
        u *
        (faTR[1] * HR + faTR[2] * u * HR + 0.5 * faTR[3] * (u^2 * HR + BR)) *
        (1.0 - δ)

    @. fb =
        Mt[1] * u * B +
        Mt[2] *
        u^2 *
        (gaL[1] * B + gaL[2] * u * B + 0.5 * gaL[3] * (u^2 * B + Mxi[2] * H)) *
        δ +
        Mt[2] *
        u^2 *
        (gaR[1] * B + gaR[2] * u * B + 0.5 * gaR[3] * (u^2 * B + Mxi[2] * H)) *
        (1.0 - δ) +
        Mt[3] * u * (gaT[1] * H + gaT[2] * u * H + 0.5 * gaT[3] * (u^2 * B + Mxi[2] * H)) +
        Mt[4] * u * BL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faL[1] * BL + faL[2] * u * BL + 0.5 * faL[3] * (u^2 * BL + Mxi[2] * HL)) *
        δ + Mt[4] * u * BR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (faR[1] * BR + faR[2] * u * BR + 0.5 * faR[3] * (u^2 * BR + Mxi[2] * HR)) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (faTL[1] * BL + faTL[2] * u * BL + 0.5 * faTL[3] * (u^2 * BL + Mxi[2] * HL)) *
        δ -
        tau *
        Mt[4] *
        u *
        (faTR[1] * BR + faTR[2] * u * BR + 0.5 * faTR[3] * (u^2 * BR + Mxi[2] * HR)) *
        (1.0 - δ)

    return nothing

end

# ------------------------------------------------------------
# 2D & 3D
# ------------------------------------------------------------

"""
$(SIGNATURES)

2D & 2D
"""
function flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    dy::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mw1, 1, 0, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mw2, 1, 0, 0)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mv, Mw, 1, 0, 0) .+
            moments_conserve_slope(gaR, MuR, Mv, Mw, 1, 0, 0)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw =  -prim[1] .* moments_conserve_slope(ga, Mu, Mv, Mw, 1, 0, 0)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mv, Mw, 2, 0, 0)
    MauR = moments_conserve_slope(gaR, MuR, Mv, Mw, 2, 0, 0)
    # Mau = moments_conserve_slope(ga, Mu, Mv, Mw, 2, 0, 0)
    MauT = moments_conserve_slope(gaT, Mu, Mv, Mw, 1, 0, 0)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mv1, Mw1, 1, 0, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mv1, Mw1, 2, 0, 0)
    MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mw1, 1, 0, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mw2, 1, 0, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mv2, Mw2, 2, 0, 0)
    MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mw2, 1, 0, 0)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    fw .*= dy

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    dxL,
    dxR,
    len,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AA{<:FN,2},Y<:AA{<:Real,2}}

    primL = mixture_conserve_prim(wL, γ)
    primR = mixture_conserve_prim(wR, γ)

    Mu1, Mv1, Mw1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Mu2, Mv2, Mw2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)

    w =
        primL[1] .* mixture_moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0) .+
        primR[1] .* mixture_moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)
    prim = mixture_conserve_prim(w, γ)
    tau = aap_hs_collision_time(prim, mi, ni, me, ne, Kn)
    for i in eachindex(tau)
        tau[i] +=
            2.0 * dt * abs(primL[1, i] / primL[end, i] - primR[1, i] / primR[end, i]) /
            (primL[1, i] / primL[end, i] + primR[1, i] / primR[end, i])
    end
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn) # pseudo primitive variables

    faL = mixture_pdf_slope(primL, swL, inK)
    sw = mixture_moments_conserve_slope(faL, Mu1, Mv1, Mw1, 1, 0)
    for j in axes(sw, 2)
        sw[:, j] .*= -primL[1, j]
    end
    faTL = mixture_pdf_slope(primL, sw, inK)
    faR = mixture_pdf_slope(primR, swR, inK)
    sw = mixture_moments_conserve_slope(faR, Mu2, Mv2, Mw2, 1, 0)
    for j in axes(sw, 2)
        sw[:, j] .*= -primR[1, j]
    end
    faTR = mixture_pdf_slope(primR, sw, inK)

    Mu, Mv, Mw, MuL, MuR = mixture_gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = mixture_pdf_slope(prim, sw0L, inK)
    gaR = mixture_pdf_slope(prim, sw0R, inK)
    sw =
        mixture_moments_conserve_slope(gaL, MuL, Mv, Mw, 1, 0) .+
        mixture_moments_conserve_slope(gaR, MuR, Mv, Mw, 1, 0)
    for j in axes(sw, 2)
        sw[:, j] .*= -prim[1, j]
    end
    gaT = mixture_pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5, axes(fw, 2))
    for j in axes(Mt, 2)
        Mt[4, j] = tau[j] * (1.0 - exp(-dt / tau[j]))
        Mt[5, j] = -tau[j] * dt * exp(-dt / tau[j]) + tau[j] * Mt[4]
        Mt[1, j] = dt - Mt[4, j]
        Mt[2, j] = -tau[j] * Mt[1, j] + Mt[5, j]
        Mt[3, j] = 0.5 * dt^2 - tau[j] * Mt[1, j]
    end

    # flux related to central distribution
    Muv = mixture_moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    MauL = mixture_moments_conserve_slope(gaL, MuL, Mv, Mw, 2, 0)
    MauR = mixture_moments_conserve_slope(gaR, MuR, Mv, Mw, 2, 0)
    MauT = mixture_moments_conserve_slope(gaT, Mu, Mv, Mw, 1, 0)

    for j in axes(fw, 2)
        @. fw[:, j] =
            Mt[1, j] * prim[1, j] * Muv[:, j] +
            Mt[2, j] * prim[1, j] * (MauL[:, j] + MauR[:, j]) +
            Mt[3, j] * prim[1, j] * MauT[:, j]
    end

    # flux related to upwind distribution
    MuvL = mixture_moments_conserve(MuL1, Mv1, Mw1, 1, 0, 0)
    MauL = mixture_moments_conserve_slope(faL, MuL1, Mv1, Mw1, 2, 0)
    MauLT = mixture_moments_conserve_slope(faTL, MuL1, Mv1, Mw1, 1, 0)

    MuvR = mixture_moments_conserve(MuR2, Mv2, Mw2, 1, 0, 0)
    MauR = mixture_moments_conserve_slope(faR, MuR2, Mv2, Mw2, 2, 0)
    MauRT = mixture_moments_conserve_slope(faTR, MuR2, Mv2, Mw2, 1, 0)

    for j in axes(fw, 2)
        @. fw[:, j] +=
            Mt[4, j] * primL[1, j] * MuvL[:, j] -
            (Mt[5, j] + tau[j] * Mt[4, j]) * primL[1, j] * MauL[:, j] -
            tau[j] * Mt[4, j] * primL[1, j] * MauLT[:, j] +
            Mt[4, j] * primR[1, j] * MuvR[:, j] -
            (Mt[5, j] + tau[j] * Mt[4, j]) * primR[1, j] * MauR[:, j] -
            tau[j] * Mt[4, j] * primR[1, j] * MauRT[:, j]
    end

    @. fw .* len

    return nothing

end

"""
$(SIGNATURES)

2D1F2V
"""
function flux_gks!(
    fw::T1,
    ff::T2,
    wL::T3,
    wR::T3,
    u::T4,
    v::T4,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    dy::Real,
    swL::T3,
    swR::T3,
) where {T1<:AV{<:FN},T2<:AM{<:FN},T3<:AV{<:FN},T4<:AM{<:FN}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mxi1, 1, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mxi2, 1, 0)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mv, Mxi, 1, 0) .+
            moments_conserve_slope(gaR, MuR, Mv, Mxi, 1, 0)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mv, Mxi, 1, 0)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mv, Mxi, 2, 0)
    MauR = moments_conserve_slope(gaR, MuR, Mv, Mxi, 2, 0)
    # Mau = moments_conserve_slope(ga, Mu, Mv, Mxi, 2, 0)
    MauT = moments_conserve_slope(gaT, Mu, Mv, Mxi, 1, 0)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mv1, Mxi1, 2, 0)
    MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mxi1, 1, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mv2, Mxi2, 2, 0)
    MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mxi2, 1, 0)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    fw .*= dy

    #--- fluxes of distribution functions ---#
    δ = heaviside.(u)
    HL = maxwellian(u, v, primL)
    HR = maxwellian(u, v, primR)
    H = maxwellian(u, v, prim)

    @. ff =
        Mt[1] * u * H +
        Mt[2] *
        u^2 *
        (gaL[1] * H + gaL[2] * u * H + gaL[3] * v * H + 0.5 * gaL[4] * (u^2 + v^2) * H) *
        δ +
        Mt[2] *
        u^2 *
        (gaR[1] * H + gaR[2] * u * H + gaR[3] * v * H + 0.5 * gaR[4] * (u^2 + v^2) * H) *
        (1.0 - δ) +
        Mt[3] *
        u *
        (gaT[1] * H + gaT[2] * u * H + gaT[3] * v * H + 0.5 * gaT[4] * (u^2 + v^2) * H) +
        Mt[4] * u * HL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faL[1] * HL +
            faL[2] * u * HL +
            faL[3] * v * HL +
            0.5 * faL[4] * (u^2 + v^2) * HL
        ) *
        δ + Mt[4] * u * HR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faR[1] * HR +
            faR[2] * u * HR +
            faR[3] * v * HL +
            0.5 * faR[4] * (u^2 + v^2) * HR
        ) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (
            faTL[1] * HL +
            faTL[2] * u * HL +
            faTL[3] * v * HL +
            0.5 * faTL[4] * (u^2 + v^2) * HL
        ) *
        δ -
        tau *
        Mt[4] *
        u *
        (
            faTR[1] * HR +
            faTR[2] * u * HR +
            faTR[3] * v * HR +
            0.5 * faTR[4] * (u^2 + v^2) * HR
        ) *
        (1.0 - δ)

    ff .*= dy

    return nothing

end

"""
$(SIGNATURES)

2D2F2V
"""
function flux_gks!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    wR::T3,
    u::T4,
    v::T4,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    dxL::Real,
    dxR::Real,
    dy::Real,
    swL::T3,
    swR::T3,
) where {T1<:AV{<:FN},T2<:AM{<:FN},T3<:AV{<:FN},T4<:AM{<:FN}}

    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mxi1, 1, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mxi2, 1, 0)
    faTR = pdf_slope(primR, sw, inK)

    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)
    sw0L = (w .- (wL .- swL .* dxL)) ./ dxL
    sw0R = ((wR .+ swR .* dxR) .- w) ./ dxR
    gaL = pdf_slope(prim, sw0L, inK)
    gaR = pdf_slope(prim, sw0R, inK)
    sw =
        -prim[1] .* (
            moments_conserve_slope(gaL, MuL, Mv, Mxi, 1, 0) .+
            moments_conserve_slope(gaR, MuR, Mv, Mxi, 1, 0)
        )
    # ga = pdf_slope(prim, sw, inK)
    # sw = -prim[1] .* moments_conserve_slope(ga, Mu, Mv, Mxi, 1, 0)
    gaT = pdf_slope(prim, sw, inK)

    # time-integration constants
    Mt = zeros(5)
    Mt[4] = tau * (1.0 - exp(-dt / tau))
    Mt[5] = -tau * dt * exp(-dt / tau) + tau * Mt[4]
    Mt[1] = dt - Mt[4]
    Mt[2] = -tau * Mt[1] + Mt[5]
    Mt[3] = 0.5 * dt^2 - tau * Mt[1]

    # flux related to central distribution
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    MauL = moments_conserve_slope(gaL, MuL, Mv, Mxi, 2, 0)
    MauR = moments_conserve_slope(gaR, MuR, Mv, Mxi, 2, 0)
    # Mau = moments_conserve_slope(ga, Mu, Mv, Mxi, 2, 0)
    MauT = moments_conserve_slope(gaT, Mu, Mv, Mxi, 1, 0)

    fw .=
        Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* (MauL .+ MauR) .+
        Mt[3] .* prim[1] .* MauT
    # fw .= Mt[1] .* prim[1] .* Muv .+ Mt[2] .* prim[1] .* Mau .+ Mt[3] .* prim[1] .* MauT

    # flux related to upwind distribution
    MuvL = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mv1, Mxi1, 2, 0)
    MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mxi1, 1, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mv2, Mxi2, 2, 0)
    MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mxi2, 1, 0)

    @. fw +=
        Mt[4] * primL[1] * MuvL - (Mt[5] + tau * Mt[4]) * primL[1] * MauL -
        tau * Mt[4] * primL[1] * MauLT + Mt[4] * primR[1] * MuvR -
        (Mt[5] + tau * Mt[4]) * primR[1] * MauR - tau * Mt[4] * primR[1] * MauRT
    # @. fw += Mt[4] * primL[1] * MuvL + Mt[4] * primR[1] * MuvR

    fw .*= dy

    #--- fluxes of distribution functions ---#
    δ = heaviside.(u)
    HL = maxwellian(u, v, primL)
    HR = maxwellian(u, v, primR)
    H = maxwellian(u, v, prim)
    BL = HL .* inK ./ (2.0 * primL[end])
    BR = HR .* inK ./ (2.0 * primR[end])
    B = H .* inK ./ (2.0 * prim[end])

    @. fh =
        Mt[1] * u * H +
        Mt[2] *
        u^2 *
        (
            gaL[1] * H +
            gaL[2] * u * H +
            gaL[3] * v * H +
            0.5 * gaL[4] * ((u^2 + v^2) * H + B)
        ) *
        δ +
        Mt[2] *
        u^2 *
        (
            gaR[1] * H +
            gaR[2] * u * H +
            gaR[3] * v * H +
            0.5 * gaR[4] * ((u^2 + v^2) * H + B)
        ) *
        (1.0 - δ) +
        Mt[3] *
        u *
        (
            gaT[1] * H +
            gaT[2] * u * H +
            gaT[3] * v * H +
            0.5 * gaT[4] * ((u^2 + v^2) * H + B)
        ) +
        Mt[4] * u * HL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faL[1] * HL +
            faL[2] * u * HL +
            faL[3] * v * HL +
            0.5 * faL[4] * ((u^2 + v^2) * HL + BL)
        ) *
        δ + Mt[4] * u * HR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faR[1] * HR +
            faR[2] * u * HR +
            faR[3] * v * HR +
            0.5 * faR[4] * ((u^2 + v^2) * HR + BR)
        ) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (
            faTL[1] * HL +
            faTL[2] * u * HL +
            faTL[3] * v * HL +
            0.5 * faTL[4] * ((u^2 + v^2) * HL + BL)
        ) *
        δ -
        tau *
        Mt[4] *
        u *
        (
            faTR[1] * HR +
            faTR[2] * u * HR +
            faTR[3] * v * HR +
            0.5 * faTR[4] * ((u^2 + v^2) * HR + BR)
        ) *
        (1.0 - δ)
    @. fb =
        Mt[1] * u * B +
        Mt[2] *
        u^2 *
        (
            gaL[1] * B +
            gaL[2] * u * B +
            gaL[3] * v * B +
            0.5 * gaL[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) *
        δ +
        Mt[2] *
        u^2 *
        (
            gaR[1] * B +
            gaR[2] * u * B +
            gaR[3] * v * B +
            0.5 * gaR[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) *
        (1.0 - δ) +
        Mt[3] *
        u *
        (
            gaT[1] * B +
            gaT[2] * u * B +
            gaT[3] * v * B +
            0.5 * gaT[4] * ((u^2 + v^2) * B + Mxi[2] * H)
        ) +
        Mt[4] * u * BL * δ -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faL[1] * BL +
            faL[2] * u * BL +
            faL[3] * v * BL +
            0.5 * faL[4] * ((u^2 + v^2) * BL + Mxi[2] * HL)
        ) *
        δ + Mt[4] * u * BR * (1.0 - δ) -
        (Mt[5] + tau * Mt[4]) *
        u^2 *
        (
            faR[1] * BR +
            faR[2] * u * BR +
            faR[3] * v * BR +
            0.5 * faR[4] * ((u^2 + v^2) * BR + Mxi[2] * HR)
        ) *
        (1.0 - δ) -
        tau *
        Mt[4] *
        u *
        (
            faTL[1] * BL +
            faTL[2] * u * BL +
            faTL[3] * v * BL +
            0.5 * faTL[4] * ((u^2 + v^2) * BL + Mxi[2] * HL)
        ) *
        δ -
        tau *
        Mt[4] *
        u *
        (
            faTR[1] * BR +
            faTR[2] * u * BR +
            faTR[3] * v * BR +
            0.5 * faTR[4] * ((u^2 + v^2) * BR + Mxi[2] * HR)
        ) *
        (1.0 - δ)

    fh .*= dy
    fb .*= dy

    return nothing

end

"""
$(SIGNATURES)

Flux reconstruction
"""
function flux_gks!(
    fw::X,
    wL::Y,
    wR::Y,
    inK::Real,
    γ::Real,
    μᵣ::Real,
    ω::Real,
    dt::Real,
    swL = zero(wL)::Y,
    swR = zero(wR)::Y,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1}}
    primL = conserve_prim(wL, γ)
    primR = conserve_prim(wR, γ)

    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)

    w =
        primL[1] .* moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0) .+
        primR[1] .* moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)
    prim = conserve_prim(w, γ)
    tau =
        vhs_collision_time(prim, μᵣ, ω) +
        2.0 * dt * abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end])

    faL = pdf_slope(primL, swL, inK)
    sw = -primL[1] .* moments_conserve_slope(faL, Mu1, Mv1, Mxi1, 1, 0)
    faTL = pdf_slope(primL, sw, inK)

    faR = pdf_slope(primR, swR, inK)
    sw = -primR[1] .* moments_conserve_slope(faR, Mu2, Mv1, Mxi2, 1, 0)
    faTR = pdf_slope(primR, sw, inK)

    MuvL = moments_conserve(MuL1, Mv1, Mxi1, 1, 0, 0)
    MauL = moments_conserve_slope(faL, MuL1, Mv1, Mxi1, 2, 0)
    MauLT = moments_conserve_slope(faTL, MuL1, Mv1, Mxi1, 1, 0)

    MuvR = moments_conserve(MuR2, Mv2, Mxi2, 1, 0, 0)
    MauR = moments_conserve_slope(faR, MuR2, Mv2, Mxi2, 2, 0)
    MauRT = moments_conserve_slope(faTR, MuR2, Mv2, Mxi2, 1, 0)

    @. fw =
        primL[1] * (MuvL - tau * MauL - tau * MauLT) +
        primR[1] * (MuvR - tau * MauR - tau * MauRT)

    return nothing
end


text:
using JLD2
using Statistics

mutable struct TTMazeOldError <: LoggerKeyData
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int

    function TTMazeOldError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TTMazeOldEvalSet.jld2") TTMazeOldEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1

        new(zeros(4,num_logged_steps), TTMazeOldEvalSet, logger_init_info[LoggerInitKey.INTERVAL])
    end
end
function lg_start!(self::TTMazeOldError, env, agent)
    num_demons = 4
    true_values = TabularTMazeUtils.get_true_values(env, self.eval_set["ests"])

    for gvf_i in 1:num_demons
        errs = []
        for (sample_i,(state,action)) in enumerate(zip(self.eval_set["states"], self.eval_set["actions"]))
            if self.eval_set["is_action"][gvf_i,sample_i]
                pred = get_demon_prediction(agent, state, Int(action))[gvf_i]
                err = (pred - true_values[gvf_i, sample_i]) .^2
                push!(errs, err)
            end
        end
        self.error[gvf_i,1] = sqrt(mean(errs))
    end
end

function lg_step!(self::TTMazeOldError, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval) + 1
        # Q_est = hcat([predict(agent.demon_learner, agent, agent.demon_weights, state, action) for (state,action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)
        num_demons = 4
        true_values = TabularTMazeUtils.get_true_values(env, self.eval_set["ests"])

        for gvf_i in 1:num_demons
            errs = []
            for (sample_i,(state,action)) in enumerate(zip(self.eval_set["states"], self.eval_set["actions"]))
                if self.eval_set["is_action"][gvf_i,sample_i]
                    pred = get_demon_prediction(agent, state, Int(action))[gvf_i]
                    err = (pred - true_values[gvf_i, sample_i]) .^2
                    push!(errs, err)
                end
            end
            self.error[gvf_i,ind] = sqrt(mean(errs))
        end
    end
end

function lg_episode_end!(self::TTMazeOldError, cur_step_in_episode, cur_step_total)
end

function save_log(self::TTMazeOldError, save_dict::Dict)
    save_dict[:ttmaze_old_error] = self.error
end

mutable struct TTMazeDirectError <: LoggerKeyData
    error::Array{Float64,2}
    eval_set::Dict
    log_interval::Int

    function TTMazeDirectError(logger_init_info)
        eval_set = @load string(pwd(),"/src/data/TTMazeDirectEvalSet.jld2") TTMazeDirectEvalSet
        num_logged_steps = fld(logger_init_info[LoggerInitKey.TOTAL_STEPS], logger_init_info[LoggerInitKey.INTERVAL]) + 1

        new(zeros(4,num_logged_steps), TTMazeDirectEvalSet, logger_init_info[LoggerInitKey.INTERVAL])
    end
end
function lg_start!(self::TTMazeDirectError, env, agent)
    num_demons = 4
    true_values = TabularTMazeUtils.get_true_values(env, self.eval_set["ests"])

    for gvf_i in 1:num_demons
        errs = []
        for (sample_i,(state,action)) in enumerate(zip(self.eval_set["states"], self.eval_set["actions"]))
            if self.eval_set["is_action"][gvf_i,sample_i]
                pred = get_demon_prediction(agent, state, Int(action))[gvf_i]
                err = (pred - true_values[gvf_i, sample_i]) .^2
                push!(errs, err)
            end
        end
        self.error[gvf_i,1] = sqrt(mean(errs))
    end
end

function lg_step!(self::TTMazeDirectError, env, agent, s, a, s_next, r, is_terminal, cur_step_in_episode, cur_step_total)
    if rem(cur_step_total, self.log_interval) == 0
        ind = fld(cur_step_total, self.log_interval) + 1
        # Q_est = hcat([predict(agent.demon_learner, agent, agent.demon_weights, state, action) for (state,action) in zip(self.eval_set["states"], self.eval_set["actions"])]...)
        num_demons = 4
        true_values = TabularTMazeUtils.get_true_values(env, self.eval_set["ests"])

        for gvf_i in 1:num_demons
            errs = []
            for (sample_i,(state,action)) in enumerate(zip(self.eval_set["states"], self.eval_set["actions"]))
                if self.eval_set["is_action"][gvf_i,sample_i]
                    pred = get_demon_prediction(agent, state, Int(action))[gvf_i]
                    err = (pred - true_values[gvf_i, sample_i]) .^2
                    push!(errs, err)
                end
            end
            self.error[gvf_i,ind] = sqrt(mean(errs))
        end
    end
end

function lg_episode_end!(self::TTMazeDirectError, cur_step_in_episode, cur_step_total)
end

function save_log(self::TTMazeDirectError, save_dict::Dict)
    save_dict[:ttmaze_direct_error] = self.error
end

ext: jl

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
#=
This file is auto-generated. Do not edit.
=#

mutable struct HydroFix <: HydroGen
    name::String
    available::Bool
    bus::Bus
    activepower::Float64
    reactivepower::Float64
    tech::TechHydro
    internal::InfrastructureSystemsInternal
end

function HydroFix(name, available, bus, activepower, reactivepower, tech, )
    HydroFix(name, available, bus, activepower, reactivepower, tech, InfrastructureSystemsInternal())
end

function HydroFix(; name, available, bus, activepower, reactivepower, tech, )
    HydroFix(name, available, bus, activepower, reactivepower, tech, )
end

# Constructor for demo purposes; non-functional.

function HydroFix(::Nothing)
    HydroFix(;
        name="init",
        available=false,
        bus=Bus(nothing),
        activepower=0.0,
        reactivepower=0.0,
        tech=TechHydro(nothing),
    )
end

"""Get HydroFix name."""
get_name(value::HydroFix) = value.name
"""Get HydroFix available."""
get_available(value::HydroFix) = value.available
"""Get HydroFix bus."""
get_bus(value::HydroFix) = value.bus
"""Get HydroFix activepower."""
get_activepower(value::HydroFix) = value.activepower
"""Get HydroFix reactivepower."""
get_reactivepower(value::HydroFix) = value.reactivepower
"""Get HydroFix tech."""
get_tech(value::HydroFix) = value.tech
"""Get HydroFix internal."""
get_internal(value::HydroFix) = value.internal


text:
using MNIST

# Extremely ugly quick and dirty
function getE2LSHData(f::IO)
    l1 = "Total time for R-NN query at radius 0.600000 (radius no. 0):\t"
    l2 = "Query point "
    time = 0
    r = Dict{Any,Any}()
    for x in EachLine(f)
        if beginswith(x,l1)
            time = float64(x[endof(l1):end])
        elseif beginswith(x,l2)
            point = int(x[(endof(l2)+1):(search(x,':')-1)])
            if contains(x,"no NNs")
                nn = 0
            else
                nn = int(x[(last(search(x,"found "))+1):(first(search(x," NNs"))-1)])
            end
            push!(r,point,(time,nn))
        end
    end
    r
end

getE2LSHData(filename) = open(getE2LSHData,filename,"r")

function write_e2lsh_points(filename,points::Vector{Vector{Float64}})
    open(filename,"w") do f
        for i = 1:size(points,1)
            for j = 1:size(points[i],1)
                print(f,@sprintf("%0.6lf ",points[i][j]))
            end
            println(f)
        end
    end
end

function write_e2lsh_points(filename,points::Array{Float64,2})
    open(filename,"w") do f
        for i = 1:size(points,1)
            for j = 1:size(points,2)
                print(f,@sprintf("%0.6lf ",points[i,j]))
            end
            println(f)
        end
    end
end

function read_el2sh_parameter_file(filename)
    params = Dict{Symbol,Any}()
    open(filename) do f
        readline(f) # skip the first lne
        while !eof(f)
            x = readline(f)
            if beginswith(x,"R")
                params[:R] = parsefloat(Float64,readline(f))
            elseif beginswith(x,"Success probability")
                params[:P] = parsefloat(Float64,readline(f))
            elseif beginswith(x,"Dimension")
                params[:d] = parseint(readline(f))
            elseif beginswith(x,"Use <u> functions")
                params[:u] = bool(parseint(readline(f)))
            elseif beginswith(x,"k")
                params[:k] = parseint(readline(f))
            elseif beginswith(x,"m")
                params[:m] = parseint(readline(f))
            elseif beginswith(x,"L")
                params[:L] = parseint(readline(f))
            elseif beginswith(x,"R^2") ||
                   beginswith(x,"W") ||
                   beginswith(x,"T") ||
                   beginswith(x,"typeHT")
                readline(f) #skip
            end
        end
        params
    end
end

function run_e2lsh_on_data(data, query; e2lshdir="/home/ubuntu/E2LSH-0.1")
    datafile = tempname()
    write_e2lsh_points(datafile,data)
    queryfile = tempname()
    write_e2lsh_points(queryfile,query)
    paramfile = "$datafile.params"
    local resultdata
    cd(e2lshdir) do
        resultdata = getE2LSHData(`bin/lsh 0.6 $datafile $queryfile`)
    end
    params = read_el2sh_parameter_file(paramfile)
    rm(datafile)
    rm(queryfile)
    rm(paramfile)
    (params,resultdata)
end

function normalized(data)
    [ d / norm(d) for d in data]
end

function run_mnist_e2lsh_on_data(data, query)
    run_e2lsh_on_data(
        normalized([ trainfeatures(i) for i in data ]),
        normalized([ testfeatures(i) for i in query ])
    )
end


text:
using PGFPlots

using CSV
using DataFrames
using Missings
using DocOpt
using Multilane

# pushPGFPlotsOptions("scale=1.2")

doc = """
Usage:
    plot_curves.jl [<filename>...]
"""

args = docopt(doc)

function only(a)
    @assert length(a) == 1
    return first(a)
end

names = Dict("qmdp"=>"QMDP",
             "outcome"=>"Naive MDP",
             "pomcpow"=>"POMCPOW",
             "meanmpc"=>"Mean MPC",
             "baseline"=>"Assume normal",
             "omniscient"=>"Omniscient"
            )

# default = Pkg.dir("Multilane", "data/all_gaps_Friday_13_Apr_10_17.csv")
default = Pkg.dir("Multilane", "data/gap_Sunday_29_Apr_11_40.csv")
filenames = args["<filename>"]
if isempty(filenames)
    push!(filenames, default)
end
# @show filenames

data = vcat([CSV.read(f) for f in filenames]...)

data[:time_to_lane] = data[:steps_to_lane].*data[:dt]
data[:brakes_per_step] = data[:nb_brakes]./data[:n_steps]

combined = by(data, [:cor, :solver, :lambda]) do df
    n = size(df, 1)
    success = df[:terminal].=="lane"
    safe = (df[:nb_brakes] .< 1) .& (df[:min_speed] .>= 15.0)
    bpm = df[:nb_brakes]./df[:distance]
    return DataFrame(reached_lane=sum(success)/n,
                     rl_std=std(success),
                     safe=sum(safe)/n,
                     safe_std=std(safe),
                     nb_brakes=mean(df[:nb_brakes]),
                     nb_brakes_std=std(df[:nb_brakes]),
                     bpm=mean(bpm),
                     bpm_std=std(bpm),
                     n=n
                    )
end

for cordf in groupby(combined, [:cor])
    @show sort(cordf, [:solver])
    cor = only(unique(cordf[:cor]))
    plts = Plots.Plot[]
    for df in groupby(cordf, [:solver])
        df = sort(df, [:lambda])
        # if cor == 1.0
        if cor == 0.0 || cor == 0.75
            name = names[only(unique(df[:solver]))]
        else
            name = nothing
        end
        plt = Plots.Linear(df[:safe], df[:reached_lane],
                           legendentry=name,
                           errorBars=ErrorBars(x=df[:safe_std]./sqrt.(df[:n]),
                                               y=df[:rl_std]./sqrt.(df[:n]))
                          )
        push!(plts, plt)
    end
    a = Axis(plts,
             # title="Correlation $cor",
             ylabel="Fraction successful",
             xlabel="Fraction safe",
             ymin=0.0,
             ymax=1.0,
             # style="legend style={at={(0.5,-0.2)}, anchor=north}, legend columns=2"
             style="legend style={at={(0.5,1.05)}}, font=\\footnotesize, anchor=south}, font=\\footnotesize, legend columns=2, grid=both",
             width="3.5in",
             height="2.7in"
            )

    bd = standard_uniform(correlation=cor)
    N = 500
    bs = [rand(Base.GLOBAL_RNG, bd) for i in 1:N]
    v0s_0 = collect(map(b->b.p_idm.v0, bs))
    Ts_0 = collect(map(b->b.p_idm.T, bs))

    s = Axis(Plots.Scatter(v0s_0, Ts_0; markSize=1);
             width="2in",
             height="2in",
             xlabel="\$T\$ (\\si{\\second})",
             ylabel="\$\\dot{x}_0\$ (\\si{\\meter\\per\\second})",
             title="Parameter distribution slice",
             style="xshift=4in, yshift=0.2in, grid=both",
            )

    file = @sprintf("/home/zach/Devel/thesis/media/gaps_%03d.tex", cor*100)
    # save(file, a)
    save(file, [a,s], include_preamble=false)
    # run(`xdg-open $file`)
end


text:
#######################################################################
# Coverage.jl
# Take Julia test coverage results and send them to various renderers
# https://github.com/JuliaCI/Coverage.jl
#######################################################################

using Coverage, Test, LibGit2

@testset "Coverage" begin
# set up base system ENV vars for testing
withenv(
    "DISABLE_AMEND_COVERAGE_FROM_SRC" => nothing,
    "COVERALLS_TOKEN" => "token_name_1",
    "COVERALLS_URL" => nothing,
    "CODECOV_URL" => nothing,
    "CODECOV_TOKEN" => nothing,
    "CODECOV_NAME" => nothing,
    "CODECOV_FLAGS" => nothing,
    "TRAVIS" => nothing,
    "TRAVIS_BRANCH" => nothing,
    "TRAVIS_COMMIT" => nothing,
    "TRAVIS_PULL_REQUEST" => nothing,
    "TRAVIS_BUILD_NUMBER" => nothing,
    "TRAVIS_JOB_ID" => nothing,
    "TRAVIS_JOB_NUMBER" => nothing,
    "TRAVIS_REPO_SLUG" => nothing,
    "APPVEYOR" => nothing,
    "APPVEYOR_PULL_REQUEST_NUMBER" => nothing,
    "APPVEYOR_ACCOUNT_NAME" => nothing,
    "APPVEYOR_PROJECT_SLUG" => nothing,
    "APPVEYOR_BUILD_VERSION" => nothing,
    "APPVEYOR_REPO_BRANCH" => nothing,
    "APPVEYOR_REPO_COMMIT" => nothing,
    "APPVEYOR_REPO_NAME" => nothing,
    "APPVEYOR_BUILD_NUMBER" => nothing,
    "APPVEYOR_BUILD_ID" => nothing,
    "APPVEYOR_JOB_ID" => nothing,
    "JENKINS" => nothing,
    "BUILD_ID" => nothing,
    "CI_PULL_REQUEST" => nothing,
    "GIT_BRANCH" => nothing
    ) do

    @testset "codecovio.jl" begin
        # these methods are only used for testing the token generation for local repos
        # and CI whilst not breaking the current API
        construct_uri_string_local(dir=pwd(); kwargs...) = Coverage.Codecov.construct_uri_string(
            Coverage.Codecov.add_local_to_kwargs(dir; kwargs...))

        construct_uri_string_ci(; kwargs...) = Coverage.Codecov.construct_uri_string(
            Coverage.Codecov.add_ci_to_kwargs(; kwargs...))

        # empty file coverage for testing
        fcs = FileCoverage[]

        # test local submission process (but only if we are in a git repo)
        _dotgit = joinpath(dirname(@__DIR__), ".git")
        if isdir(_dotgit) || isfile(_dotgit)
            LibGit2.with(LibGit2.GitRepoExt(pwd())) do repo
                # default values
                codecov_url = construct_uri_string_local()
                @test occursin("codecov.io", codecov_url)
                @test occursin("commit=", codecov_url)
                @test occursin("branch=", codecov_url)
                @test !occursin("service", codecov_url)

                # env var url override
                withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                    codecov_url = construct_uri_string_local()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("commit=", codecov_url)
                    @test occursin("branch=", codecov_url)
                    @test !occursin("service", codecov_url)

                    # function argument url override
                    codecov_url = construct_uri_string_local(codecov_url = "https://enterprise-codecov-2.com")
                    @test occursin("enterprise-codecov-2.com", codecov_url)
                    @test occursin("commit=", codecov_url)
                    @test occursin("branch=", codecov_url)
                    @test !occursin("service", codecov_url)
                    @test !occursin("name", codecov_url)
                    @test !occursin("flags", codecov_url)

                    # env var token
                    withenv( "CODECOV_TOKEN" => "token_name_1",
                             "CODECOV_NAME" => "cv_name",
                             "CODECOV_FLAGS" => "cv_flags" ) do

                        codecov_url = construct_uri_string_local()
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_1", codecov_url)
                        @test occursin("name=cv_name", codecov_url)
                        @test occursin("flags=cv_flags", codecov_url)
                        @test !occursin("service", codecov_url)

                        # function argument token url override
                        codecov_url = construct_uri_string_local(token="token_name_2")
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_2", codecov_url)
                        @test !occursin("service", codecov_url)
                    end
                end
            end
        else
            @warn "skipping local repo tests for Codecov, since not a git repo"
        end

        # test faulty non-CI submission

        @test_throws(ErrorException("No compatible CI platform detected"),
                     Coverage.Codecov.submit(fcs; dry_run = true))

        # test travis-ci submission process

        # set up travis env
        withenv(
            "TRAVIS" => "true",
            "TRAVIS_BRANCH" => "t_branch",
            "TRAVIS_COMMIT" => "t_commit",
            "TRAVIS_PULL_REQUEST" => "t_pr",
            "TRAVIS_BUILD_NUMBER" => "t_build_num",
            "TRAVIS_JOB_ID" => "t_job_id",
            "TRAVIS_REPO_SLUG" => "t_slug",
            "TRAVIS_JOB_NUMBER" => "t_job_num",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=travis-org", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("job=t_job_id", codecov_url)
            @test occursin("slug=t_slug", codecov_url)
            @test occursin("build=t_job_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(;codecov_url = "https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=travis-org", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_job_id", codecov_url)
                @test occursin("slug=t_slug", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_2", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_job_id", codecov_url)
                    @test occursin("slug=t_slug", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)
                end
            end
        end

        # test appveyor submission process

        # set up appveyor env
        withenv(
            "APPVEYOR" => "true",
            "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
            "APPVEYOR_ACCOUNT_NAME" => "t_account",
            "APPVEYOR_PROJECT_SLUG" => "t_slug",
            "APPVEYOR_BUILD_VERSION" => "t_version",
            "APPVEYOR_REPO_BRANCH" => "t_branch",
            "APPVEYOR_REPO_COMMIT" => "t_commit",
            "APPVEYOR_REPO_NAME" => "t_repo",
            "APPVEYOR_BUILD_ID" => "t_build_id",
            "APPVEYOR_JOB_ID" => "t_job_num",
            ) do

                # default values
                codecov_url = construct_uri_string_ci()
                @test occursin("codecov.io", codecov_url)
                @test occursin("service=appveyor", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                @test occursin("build=t_job_num", codecov_url)

                # env var url override
                withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # function argument url override
                    codecov_url = construct_uri_string_ci(codecov_url = "https://enterprise-codecov-2.com")
                    @test occursin("enterprise-codecov-2.com", codecov_url)
                    @test occursin("service=appveyor", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                    @test occursin("build=t_job_num", codecov_url)

                    # env var token
                    withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                        codecov_url = construct_uri_string_ci()
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_1", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)

                        # function argument token url override
                        codecov_url = construct_uri_string_ci(token="token_name_2")
                        @test occursin("enterprise-codecov-1.com", codecov_url)
                        @test occursin("token=token_name_2", codecov_url)
                        @test occursin("branch=t_branch", codecov_url)
                        @test occursin("commit=t_commit", codecov_url)
                        @test occursin("pull_request=t_pr", codecov_url)
                        @test occursin("job=t_account%2Ft_slug%2Ft_version", codecov_url)
                        @test occursin("build=t_job_num", codecov_url)
                    end
                end
            end

        # test circle ci submission process

        # set up circle ci env
        withenv(
            "CIRCLECI" => "true",
            "CIRCLE_PR_NUMBER" => "t_pr",
            "CIRCLE_PROJECT_USERNAME" => "t_proj",
            "CIRCLE_BRANCH" => "t_branch",
            "CIRCLE_SHA1" => "t_commit",
            "CIRCLE_PROJECT_REPONAME" => "t_repo",
            "CIRCLE_BUILD_URL" => "t_url",
            "CIRCLE_BUILD_NUM" => "t_num",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=circleci", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("pull_request=t_pr", codecov_url)
            @test occursin("build_url=t_url", codecov_url)
            @test occursin("build=t_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(codecov_url="https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=circleci", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("pull_request=t_pr", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=circleci", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("pull_request=t_pr", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)
                end
            end
        end

        # test Jenkins ci submission process

        # set up Jenkins ci env
        withenv(
            "JENKINS" => "true",
            "GIT_BRANCH" => "t_branch",
            "GIT_COMMIT" => "t_commit",
            "JOB_NAME" => "t_job",
            "BUILD_ID" => "t_num",
            "BUILD_URL" => "t_url",
            "JENKINS_URL" => "t_jenkins_url",
            ) do

            # default values
            codecov_url = construct_uri_string_ci()
            @test occursin("codecov.io", codecov_url)
            @test occursin("service=jenkins", codecov_url)
            @test occursin("branch=t_branch", codecov_url)
            @test occursin("commit=t_commit", codecov_url)
            @test occursin("build_url=t_url", codecov_url)
            @test occursin("build=t_num", codecov_url)

            # env var url override
            withenv( "CODECOV_URL" => "https://enterprise-codecov-1.com" ) do

                codecov_url = construct_uri_string_ci()
                @test occursin("enterprise-codecov-1.com", codecov_url)
                @test occursin("service=jenkins", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # function argument url override
                codecov_url = construct_uri_string_ci(codecov_url="https://enterprise-codecov-2.com")
                @test occursin("enterprise-codecov-2.com", codecov_url)
                @test occursin("service=jenkins", codecov_url)
                @test occursin("branch=t_branch", codecov_url)
                @test occursin("commit=t_commit", codecov_url)
                @test occursin("build_url=t_url", codecov_url)
                @test occursin("build=t_num", codecov_url)

                # env var token
                withenv( "CODECOV_TOKEN" => "token_name_1" ) do

                    codecov_url = construct_uri_string_ci()
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("token=token_name_1", codecov_url)
                    @test occursin("service=jenkins", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)

                    # function argument token url override
                    codecov_url = construct_uri_string_ci(token="token_name_2")
                    @test occursin("enterprise-codecov-1.com", codecov_url)
                    @test occursin("service=jenkins", codecov_url)
                    @test occursin("branch=t_branch", codecov_url)
                    @test occursin("commit=t_commit", codecov_url)
                    @test occursin("build_url=t_url", codecov_url)
                    @test occursin("build=t_num", codecov_url)
                end
            end
        end

        # test codecov token masking
        withenv(
            "APPVEYOR" => "true",
            "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
            "APPVEYOR_ACCOUNT_NAME" => "t_account",
            "APPVEYOR_PROJECT_SLUG" => "t_slug",
            "APPVEYOR_BUILD_VERSION" => "t_version",
            "APPVEYOR_REPO_BRANCH" => "t_branch",
            "APPVEYOR_REPO_COMMIT" => "t_commit",
            "APPVEYOR_REPO_NAME" => "t_repo",
            "APPVEYOR_BUILD_ID" => "t_build_id",
            "APPVEYOR_JOB_ID" => "t_job_num",
            "CODECOV_URL" => "https://enterprise-codecov-1.com",
            "CODECOV_TOKEN" => "token_name_1"
            ) do
                codecov_url = construct_uri_string_ci()
                masked = Coverage.Codecov.mask_token(codecov_url)
                @test !occursin("token_name_1", masked)
                @test occursin("token=<HIDDEN>", masked)
        end

        # in the case above, the token is at the end. Let's test explicitly,
        # that this also works if the token occurs earlier in the url
        url = "https://enterprise-codecov-1.com/upload/v2?token=token_name_1&build=t_job_num"
        masked = Coverage.Codecov.mask_token(url)
        @test masked == "https://enterprise-codecov-1.com/upload/v2?token=<HIDDEN>&build=t_job_num"
    end

    @testset "coveralls" begin
        # NOTE: this only returns actual content if this package is devved.
        # Hence the test is basically on this function returning something
        # (rather than erroring)
        git = Coverage.Coveralls.query_git_info()
        @test git["remotes"][1]["name"] == "origin"
        @test haskey(git["remotes"][1], "url")

        # for testing submit_***()
        fcs = FileCoverage[]

        # an error should be raised if there is no coveralls token set
        withenv("COVERALLS_TOKEN" => nothing,
                "REPO_TOKEN" => nothing,  # this is deprecrated, use COVERALLS_TOKEN
                "APPVEYOR" => "true",  # use APPVEYOR as an example to make the test reach the repo token check
                "APPVEYOR_BUILD_NUMBER" => "my_job_num",
                "APPVEYOR_BUILD_ID" => "my_build_id",
                "APPVEYOR_JOB_ID" => "my_job_id") do
            @test_throws(ErrorException("Coveralls submission requires a COVERALLS_TOKEN environment variable"),
                         Coverage.Coveralls.prepare_request(fcs, false))
        end

        # test error if not local and no CI platform can be detected from ENV
        @test_throws(ErrorException("No compatible CI platform detected"),
                     Coverage.Coveralls.prepare_request(fcs, false))

        # test local submission, when we are local
        _dotgit = joinpath(dirname(@__DIR__), ".git")
        if isdir(_dotgit) || isfile(_dotgit)
            request = Coverage.Coveralls.prepare_request(fcs, true)
            @test request["repo_token"] == "token_name_1"
            @test isempty(request["source_files"])
            @test haskey(request, "git")
            @test request["git"]["remotes"][1]["name"] == "origin"
            @test !haskey(request, "service_job_id")
            @test request["service_name"] == "local"
            @test !haskey(request, "service_pull_request")
        else
            @warn "skipping local repo tests for Coveralls, since not a git repo"
        end

        # test custom (environment variables)
        withenv("COVERALLS_SERVICE_NAME" => "My CI",
                "COVERALLS_PULL_REQUEST" => "c_pr",
                "COVERALLS_SERVICE_NUMBER" => "ci_num",
                "COVERALLS_SERVICE_JOB_NUMBER" => "ci_job_num",
                "COVERALLS_SERVICE_JOB_ID" => "ci_job_id",
                "COVERALLS_FLAG_NAME" => "t_pr") do
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_number"] == "ci_num"
                @test request["service_job_number"] == "ci_job_num"
                @test request["service_job_id"] == "ci_job_id"
                @test request["service_name"] == "My CI"
                @test request["service_pull_request"] == "c_pr"
        end

        # test APPVEYOR
        withenv("APPVEYOR" => "true",
                "APPVEYOR_PULL_REQUEST_NUMBER" => "t_pr",
                "APPVEYOR_BUILD_NUMBER" => "my_build_num",
                "APPVEYOR_BUILD_ID" => "my_build_id") do
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_job_number"] == "my_build_num"
                @test request["service_job_id"] == "my_build_id"
                @test request["service_name"] == "appveyor"
                @test request["service_pull_request"] == "t_pr"
                @test !haskey(request, "parallel")
        end

        # test Travis
        withenv("TRAVIS" => "true",
                "TRAVIS_BUILD_NUMBER" => "my_job_num",
                "TRAVIS_JOB_ID" => "my_job_id",
                "TRAVIS_PULL_REQUEST" => "t_pr",
                "COVERALLS_PARALLEL" => "true") do
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_number"] == "my_job_num"
                @test request["service_job_id"] == "my_job_id"
                @test request["service_name"] == "travis-ci"
                @test request["service_pull_request"] == "t_pr"
                @test request["parallel"] == "true"
        end

        # test Jenkins
        withenv("JENKINS" => "true",
                "BUILD_ID" => "my_job_id",
                "CI_PULL_REQUEST" => true,
                "COVERALLS_PARALLEL" => "not") do
                my_git_info = Dict("remote_name" => "my_origin")
                request = Coverage.Coveralls.prepare_request(fcs, false)
                @test request["repo_token"] == "token_name_1"
                @test request["service_job_id"] == "my_job_id"
                @test request["service_name"] == "jenkins-ci"
                @test !haskey(request, "service_pull_request")
                @test !haskey(request, "parallel")

                withenv("CI_PULL_REQUEST" => "false",
                        "GIT_BRANCH" => "my_remote/my_branch",
                        "COVERALLS_SERVICE_NUMBER" => "ci_num",
                        "COVERALLS_SERVICE_JOB_NUMBER" => "ci_job_num",
                        "COVERALLS_SERVICE_JOB_ID" => "ci_job_id"
                       ) do
                        request = Coverage.Coveralls.prepare_request(fcs, false)
                        @test haskey(request, "git")
                        @test request["git"]["branch"] == "my_branch"
                        @test request["service_number"] == "ci_num"
                        @test request["service_job_number"] == "ci_job_num"
                        @test request["service_job_id"] == "ci_job_id"
                end
        end

        # test git_info (only works with Jenkins & local at the moment)
        withenv("JENKINS" => "true",
                "BUILD_ID" => "my_job_id",
                "CI_PULL_REQUEST" => true) do
                # we can pass in our own function, that returns a dict
                my_git_info() = Dict("test" => "test")
                request = Coverage.Coveralls.prepare_request(fcs, false, my_git_info)
                @test haskey(request, "git")

                # or directly a dict
                request = Coverage.Coveralls.prepare_request(fcs, false, Dict("test" => "test"))
                @test haskey(request, "git")
        end
    end

end # of withenv( => nothing)

end # of @testset "Coverage"


text:
using ApproxFun, IntervalSets, SpecialFunctions, LinearAlgebra, Random, Test
    import ApproxFun: HeavisideSpace, PointSpace, DiracSpace, PiecewiseSegment

@testset "Singularities" begin
    @testset "sqrt" begin
        x=Fun(identity);
        @test sqrt(cos(π/2*x))(.1) ≈ sqrt(cos(.1π/2))

        x=Fun(identity,-2..2)
        u=sqrt(4-x^2)/(2π)

        @test u(.1) ≈ sqrt(4-0.1^2)/(2π)
        @test sum(u) ≈ 1

        #this call threw an error, which we check
        @test length(values(u)) == 1

        f = Fun(x->x*cot(π*x/2))
        x = Fun(identity)
        u = Fun(JacobiWeight(1.,1.,ChebyshevInterval()), (f/(1-x^2)).coefficients)
        @test 1/(0.1*cot(π*.1/2)) ≈ (1/u)(.1)

        @test (x/u)(.1) ≈ tan(π*.1/2)

        f=Fun(x->exp(-x^2),Line(0.,0.,-.5,-.5),400)
        @test sum(f) ≈ sqrt(π)

        f=Fun(x->exp(x)/sqrt(1-x.^2),JacobiWeight(-.5,-.5))
        @test f(.1) ≈ (x->exp(x)/sqrt(1-x.^2))(.1)

        @test norm(Fun(exp,Legendre(0..1))+sqrt(Fun(0..1))) ≈ 2.491141949903508

        # sampling Chebyshev
        x=Fun(identity)
        f = exp(x)/sqrt(1-x^2)
        g = cumsum(f)
        @test abs(g(-1)) ≤ 1E-15
        @test g'(0.1) ≈ f(0.1)
    end

    @testset "JacobiWeight Derivative" begin
        S = JacobiWeight(-1.,-1.,Chebyshev(0..1))

        # Checks bug in Derivative(S)
        @test typeof(ConstantSpace(0..1)) <: Space{ClosedInterval{Int},Float64}

        D=Derivative(S)
        f=Fun(S,Fun(exp,0..1).coefficients)
        x=0.1
        @test f(x) ≈ exp(x)*x^(-1)*(1-x)^(-1)/4
        @test (D*f)(x) ≈ -exp(x)*(1+(x-3)*x)/(4*(x-1)^2*x^2)

        S=JacobiWeight(-1.,0.,Chebyshev(0..1))
        D=Derivative(S)

        f=Fun(S,Fun(exp,0..1).coefficients)
        x=.1
        @test f(x) ≈ exp(x)*x^(-1)/2
        @test (D*f)(x) ≈ exp(x)*(x-1)/(2x^2)
    end

    @testset "Jacobi singularity" begin
        x = Fun(identity)
        f = exp(x)/(1-x.^2)

        @test f(.1) ≈ exp(.1)/(1-.1^2)
        f = exp(x)/(1-x.^2).^1
        @test f(.1) ≈ exp(.1)/(1-.1^2)
        f = exp(x)/(1-x.^2).^1.0
        @test f(.1) ≈ exp(.1)/(1-.1^2)

        ## 1/f with poles

        x=Fun(identity)
        f=sin(10x)
        g=1/f

        @test g(.123) ≈ csc(10*.123)
    end

    @testset "Ray and Line" begin
        @testset "Ray" begin
            @test Inf in Ray()   # this was a bug
            @test Space(0..Inf) == Chebyshev(Ray())
            f=Fun(x->exp(-x),0..Inf)
            @test f(0.1) ≈ exp(-0.1)
            @test f'(.1) ≈ -f(.1)

            x=Fun(identity,Ray())
            f=exp(-x)
            u=integrate(f)
            @test (u(1.)-u(0)-1) ≈ -f(1)

            x=Fun(identity,Ray())
            f=x^(-0.123)*exp(-x)
            @test integrate(f)'(1.) ≈ f(1.)

            @test ≈(sum(Fun(sech,0..Inf)),sum(Fun(sech,0..40));atol=1000000eps())
            @test Line() ∪ Ray() == Line()
            @test Line() ∩ Ray() == Ray()

            f=Fun(sech,Line())
            Fun(f,Ray())(2.0) ≈ sech(2.0)
            Fun(f,Ray(0.,π))(-2.0) ≈ sech(-2.0)
            Fun(sech,Ray(0.,π))(-2.0) ≈ sech(-2.0)
        end

        @testset "Ei (Exp Integral)" begin
            y=Fun(Ray())
            q=integrate(exp(-y)/y)
            @test (q-last(q))(2.) ≈ (-0.04890051070806113)

            ## Line
            f=Fun(x->exp(-x^2),Line())

            @test f'(0.1) ≈ -2*0.1exp(-0.1^2)
            @test (Derivative()*f)(0.1) ≈ -2*0.1exp(-0.1^2)
        end

        @testset "PeriodicLine" begin
            d=PeriodicLine()
            D=Derivative(d)

            f = Fun(x->sech(x-0.1),d,200)
            @test f(1.) ≈ sech(1-0.1)

            f=Fun(x->sech(x-0.1),d)
            @test f(1.) ≈ sech(1-0.1)

            @test ≈((D*f)(.2),-0.0991717226583897;atol=100000eps())
            @test ≈((D^2*f)(.2),-0.9752522555114987;atol=1000000eps())
        end
    end

    @testset "LogWeight" begin
        x=Fun(identity,-1..1)
        f=exp(x+1)-1
        @test log(f)(0.1) ≈ log(f(0.1))

        x=Fun(identity,0..1)
        f=exp(x)-1
        @test log(f)(0.1) ≈ log(f(0.1))

        x=Fun(identity,0..1)
        @test Fun(exp(x)/x-1/x,Chebyshev)(0.1) ≈ (exp(0.1)-1)/0.1

        x=Fun(identity,0..1)
        f=1/x
        p=integrate(f)
        @test (p-p(1.))(0.5) ≈ log(0.5)

        f=1/(1-x)
        p=integrate(f)
        @test (p-p(0.))(0.5) ≈ -log(1-0.5)
    end

    @testset "Complex domains sqrt" begin
        a=1+10*im; b=2-6*im
        d=Curve(Fun(x->1+a*x+b*x^2))

        x=Fun(d)
        w=sqrt(abs(leftendpoint(d)-x))*sqrt(abs(rightendpoint(d)-x))

        @test sum(w/(x-2.))/(2π*im) ≈ (-4.722196879007759+2.347910413861846im)
        @test linesum(w*log(abs(x-2.)))/π ≈ (88.5579588360686)

        a=Arc(0.,1.,0.,π/2)
        ζ=Fun(identity,a)
        f=Fun(exp,a)*sqrt(abs((ζ-1)*(ζ-im)))
    end

    @testset "DiracDelta and PointSpace" begin
        a,b=DiracDelta(0.),DiracDelta(1.)
        f=Fun(exp)
        g=a+0.2b+f
        @test components(g)[2](0.) ≈ 1.
        @test g(.1) ≈ exp(.1)
        @test sum(g) ≈ (sum(f)+1.2)

        #Checks prevoius bug
        δ=DiracDelta()
        x=Fun()
        w=sqrt(1-x^2)
        w+δ

        ## PointSpace

        @test eltype(domain(PointSpace([0,0.1,1])) ) == Float64

        f=Fun(x->(x-0.1),PointSpace([0,0.1,1]))
        @test roots(f) == [0.1]

        a=Fun(exp,space(f))
        @test f/a == Fun(x->(x-0.1)*exp(-x),space(f))

        f = Fun(space(f),[1.,2.,3.])

        g = f + Fun(2..3)
        @test f(0.0) ≈ g(0.0) ≈ 1.0
        @test f(0.1) ≈ g(0.1) ≈ 2.0
        @test f(1.0) ≈ g(1.0) ≈ 3.0

        @test g(2.3) ≈ 2.3

        h = a + Fun(2..3)

        # for some reason this test is broken only on Travis
        @test_skip g/h ≈ f/a + Fun(1,2..3)
    end

    @testset "DiracDelta integration and differentiation" begin
        δ = DiracDelta()
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0,Inf])
        @test h(-2) == 0
        @test h(2) == 1

        δ = 0.3DiracDelta(0.1) + 3DiracDelta(2.3)
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0.1,2.3,Inf])
        @test h(-2) == 0
        @test h(2) == 0.3
        @test h(3) == 3.3

        δ = (0.3+1im)DiracDelta(0.1) + 3DiracDelta(2.3)
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0.1,2.3,Inf])
        @test h(-2) == 0
        @test h(2) == 0.3+1im
        @test h(3) == 3.3+1im
    end

    @testset "DiracDelta sampling" begin
        δ = 0.3DiracDelta(0.1) + 3DiracDelta(2.3)
        Random.seed!(0)
        for _=1:10
            @test sample(δ) ∈ [0.1, 2.3]
        end
        Random.seed!(0)
        r = sample(δ, 10_000)
        @test count(i -> i == 0.1, r)/length(r) ≈ 0.3/(3.3) atol=0.01
    end

    @testset "Multiple roots" begin
        x=Fun(identity,-1..1)

        @test (1/x^2)(0.1) ≈ 100.
        @test (1/x^2)(-0.1) ≈ 100.

        fc=x*(1+x)^2
        @test (1/fc)(0.1) ≈ 1/fc(0.1)

        fc=x*(1-x)^2
        @test (1/fc)(0.1) ≈ 1/fc(0.1)
    end

    @testset "special function singularities" begin
        x=Fun(0..1)
        @test erf(sqrt(x))(0.1) ≈ erf(sqrt(0.1))
        @test erfc(sqrt(x))(0.1) ≈ erfc(sqrt(0.1))

        ## roots of log(abs(x-y))
        x=Fun(-2..(-1))
        @test space(abs(x)) == Chebyshev(-2 .. (-1))

        @test roots(abs(x+1.2)) ≈ [-1.2]

        f = sign(x+1.2)
        @test space(f) isa PiecewiseSpace
        @test f(-1.4) == -1
        @test f(-1.1) == 1

        f=abs(x+1.2)
        @test abs(f)(-1.3) ≈ f(-1.3)
        @test abs(f)(-1.1) ≈ f(-1.1)
        @test norm(abs(f)-f)<10eps()

        @test norm(sign(f)-Fun(1,space(f)))<10eps()

        @test log(f)(-1.3) ≈ log(abs(-1.3+1.2))
        @test log(f)(-1.1) ≈ log(abs(-1.1+1.2))
    end

    @testset "#393" begin
        x = Fun(0..1)
        f = exp(x)*sqrt(x)*log(1-x)
        @test f(0.1) ≈ exp(0.1)*sqrt(0.1)*log(1-0.1)
    end

    @testset "Jacobi conversions" begin
        S1,S2=JacobiWeight(3.,1.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(0.,1.))
        f=Fun(S1,[1,2,3.])
        C=Conversion(S1,S2)
        Cf=C*f
        @test Cf(0.1) ≈ f(0.1)

        S1,S2=JacobiWeight(3.,2.,Jacobi(1.,1.)),JacobiWeight(1.,1.,Jacobi(0.,0.))
        f=Fun(S1,[1,2,3.])
        C=Conversion(S1,S2)
        Cf=C*f
        @test Cf(0.1) ≈ f(0.1)
    end

    @testset "Derivative operator for HeavisideSpace" begin
        H = HeavisideSpace([-1.0,0.0,1.0])
        @test Fun(H, [1.0])(1.0) == 0.0
        @test Fun(H, [0.0,1.0])(1.0) == 1.0

        H=HeavisideSpace([1,2,3])
        D=Derivative(H)
        @test domain(D)==PiecewiseSegment([1,2,3])
        @test D[1,1]==-1
        @test D[1,2]==1

        H=HeavisideSpace([1,2,3,Inf])
        D=Derivative(H)
        @test domain(D)==PiecewiseSegment([1,2,3,Inf])
        @test D[1,1]==-1
        @test D[2,2]==-1
        @test D[1,2]==1

        S = HeavisideSpace([-1.0,0.0,1.0])
        @test Derivative(S) === Derivative(S,1)
    end
end


text:
using MarkovTransitionMatrices

# state grids
grid_x1 = 0.0:1.0:10.0
grid_x2 = -1.5:0.5:15.0

# probabilities must all be greater than minp
minp = 1e-8

# Correlated random walk
μ(s) = [s...]
Σ(s) = [1.0 0.5; 0.5 1.0]
P = markov_transition(μ, Σ, minp, grid_x1, grid_x2)

# Markov-switching process

# regime transition matrix (NOT tranpsosed! - rows sum to 1.0)
πswitch = [.9 .1; .4 .6]

μswitch(r::Int, s) = r==1 ? [s...] : [s...] .+ ones(2)
Σswitch(r::Int, s) = r==1 ? Matrix(1.0I,2,2) : [1.0 0.5; 0.5 1.0]
Pswitch = markovswitching_transition(μswitch, Σswitch, πswitch, 1e-8, grid_x1, grid_x2)


text:
# This file contains functions for the Enhanced Russell Graph Slack Based Measure DEA model
"""
    EnhanceRussellGraphDEAModel
An data structure representing an Enhanced Russell Graph Slack Based Measure DEA model.
"""
struct EnhancedRussellGraphDEAModel <: AbstractTechnicalDEAModel
    n::Int64
    m::Int64
    s::Int64
    rts::Symbol
    dmunames::Union{Vector{String},Nothing}
    eff::Vector
    beta::Vector
    tX::Matrix
    tY::Matrix
    slackX::Matrix
    slackY::Matrix
    lambda::SparseMatrixCSC{Float64, Int64}
    Xtarget::Matrix
    Ytarget::Matrix
end

"""
    deaerg(X, Y)
Compute data envelopment analysis Enhanced Russell Graph Slack Based Measure for inputs `X` and outputs `Y`.

# Optional Arguments
- `rts=:CRS`: chooses constant returns to scale. For variable returns to scale choose `:VRS`.
- `Xref=X`: Identifies the reference set of inputs against which the units are evaluated.
- `Yref=Y`: Identifies the reference set of outputs against which the units are evaluated.
- `names`: a vector of strings with the names of the decision making units.

# Examples
```jldoctest
julia> X = [2; 4; 8; 12; 6; 14; 14; 9.412];

julia> Y = [1; 5; 8; 9; 3; 7; 9; 2.353] ;

julia> deaerg(X, Y, rts = :VRS)
Enhanced Russell Graph Slack Based Measure DEA Model
DMUs = 8; Inputs = 1; Outputs = 1
Orientation = Graph; Returns to Scale = VRS
───────────────────────────────────────
   efficiency    beta  slackX1  slackY1
───────────────────────────────────────
1    1.0       1.0     0.0        0.0
2    1.0       1.0     0.0        0.0
3    1.0       1.0     0.0        0.0
4    1.0       1.0     0.0        0.0
5    0.4       0.6     2.0        2.0
6    0.47619   1.0     7.33333    0.0
7    0.857143  1.0     2.0        0.0
8    0.2       0.4706  5.412      2.647
───────────────────────────────────────
```
"""
function deaerg(X::Union{Matrix,Vector}, Y::Union{Matrix,Vector};
    rts::Symbol = :CRS, slack::Bool = true,
    Xref::Union{Matrix,Vector,Nothing} = nothing, Yref::Union{Matrix,Vector,Nothing} = nothing,
    names::Union{Vector{String},Nothing} = nothing,
    optimizer::Union{DEAOptimizer,Nothing} = nothing)::EnhancedRussellGraphDEAModel

    # Check parameters
    nx, m = size(X, 1), size(X, 2)
    ny, s = size(Y, 1), size(Y, 2)

    if Xref === nothing Xref = X end
    if Yref === nothing Yref = Y end

    nrefx, mref = size(Xref, 1), size(Xref, 2)
    nrefy, sref = size(Yref, 1), size(Yref, 2)

    if nx != ny
        throw(DimensionMismatch("number of rows in X and Y ($nx, $ny) are not equal"));
    end
    if nrefx != nrefy
        throw(DimensionMismatch("number of rows in Xref and Yref ($nrefx, $nrefy) are not equal"));
    end
    if m != mref
        throw(DimensionMismatch("number of columns in X and Xref ($m, $mref) are not equal"));
    end
    if s != sref
        throw(DimensionMismatch("number of columns in Y and Yref ($s, $sref) are not equal"));
    end

    # Default optimizer
    if optimizer === nothing
        optimizer = DEAOptimizer(:LP)
    end

    # Compute efficiency for each DMU
    n = nx
    nref = nrefx

    effi = zeros(n)
    betai = zeros(n)
    tXi = zeros(n, m)
    tYi = zeros(n, s)
    lambdaeff = spzeros(n, nref)

    for i=1:n
        # Value of inputs and outputs to evaluate
        x0 = X[i,:]
        y0 = Y[i,:]

        # Create the optimization model
        deamodel = newdeamodel(optimizer)

        @variable(deamodel, beta >= 0)
        @variable(deamodel, tX[1:m] >= 0)
        @variable(deamodel, tY[1:s] >= 0)
        @variable(deamodel, mu[1:nref] >= 0)

        @objective(deamodel, Min, beta - 1 / m * sum(tX[t] / x0[t] for t in 1:m))

        @constraint(deamodel, beta + 1 / s * sum(tY[t] / y0[t] for t in 1:s) == 1)
        @constraint(deamodel, [j in 1:m], sum(Xref[t,j] * mu[t] for t in 1:nref) == beta * x0[j] - tX[j])
        @constraint(deamodel, [j in 1:s], sum(Yref[t,j] * mu[t] for t in 1:nref) == beta * y0[j] + tY[j])

        # Add return to scale constraints
        if rts == :CRS
            # No contraint to add for constant returns to scale
        elseif rts == :VRS
            @constraint(deamodel, sum(mu) == beta)
        else
            throw(ArgumentError("`rts` must be :CRS or :VRS"));
        end

        # Optimize and return results
        JuMP.optimize!(deamodel)

        effi[i]  = JuMP.objective_value(deamodel)
        betai[i] = JuMP.value(beta)
        tXi[i,:] = JuMP.value.(tX)
        tYi[i,:] = JuMP.value.(tY)
        mui = JuMP.value.(mu)
        lambdaeff[i,:] = mui ./ betai[i]

        # Check termination status
        if (termination_status(deamodel) != MOI.OPTIMAL) && (termination_status(deamodel) != MOI.LOCALLY_SOLVED)
            @warn ("DMU $i termination status: $(termination_status(deamodel)). Primal status: $(primal_status(deamodel)). Dual status: $(dual_status(deamodel))")
        end

    end

    slackX = tXi ./ betai
    slackY = tYi ./ betai

    # Get X and Y targets
    Xtarget = X - slackX
    Ytarget = Y + slackY

    return EnhancedRussellGraphDEAModel(n, m, s, rts, names, effi, betai, tXi, tYi, slackX, slackY, lambdaeff, Xtarget, Ytarget)

end

function Base.show(io::IO, x::EnhancedRussellGraphDEAModel)
    compact = get(io, :compact, false)

    n = nobs(x)
    m = ninputs(x)
    s = noutputs(x)
    dmunames = names(x)

    eff = efficiency(x)
    beta = x.beta
    slackX = slacks(x, :X)
    slackY = slacks(x, :Y)
    hasslacks = ! isempty(slackX)

    if !compact
        print(io, "Enhanced Russell Graph Slack Based Measure DEA Model \n")
        print(io, "DMUs = ", n)
        print(io, "; Inputs = ", m)
        print(io, "; Outputs = ", s)
        print(io, "\n")
        print(io, "Orientation = Graph")
        print(io, "; Returns to Scale = ", string(x.rts))
        print(io, "\n")
        show(io, CoefTable(hcat(eff, beta, slackX, slackY), ["efficiency"; "beta"; ["slackX$i" for i in 1:m ]; ["slackY$i" for i in 1:s ]], dmunames))
    end

end

function efficiency(model::EnhancedRussellGraphDEAModel, type::Symbol)::Vector

    if type == :beta return model.beta end

    throw(ArgumentError("$(typeof(model)) has no efficiency $(type)"));

end


text:
function test_laplace(T = Float64, N = 128)
    M = 2N
    splinedegree = 1
    obstacle = Kite{T}(2)
    param = parameterization(obstacle)
    paramdomain = domain(param)
    SingleLayerPotential = Laplace_SLP_2D{T}()
    BIO = BoundaryIntegralOperator(SingleLayerPotential, obstacle, param)
    if paramdomain == Interval{:closed,:open,T}(0,1)
        basis = BSplineTranslatesBasis{T}(N, splinedegree)
    else
        a, b = extrema(paramdomain)
        m = mapto(a..b, 0..1)
        splines = BSplineTranslatesBasis{T}(N, splinedegree)
        basis = MappedDict(splines, m)
    end
    basis_obstacle = BasisFunctions.ParamDict(basis, param, obstacle)
    coll_points = PeriodicEquispacedGrid(M, paramdomain)
    coll_points_obstacle = map_grid(param, coll_points)
    sampling_col = GridSampling(coll_points, T)
    sampling_col_obstacle = GridSampling(coll_points_obstacle, T)
    sampling_gal = ProjectionSampling(basis, measure(BIO))
    bcond = laplace_param_bcond_harmonic(param)
    bcond_field = laplace_bcond_harmonic()
    BEM_col = (sampling_col * BIO) * basis
    BEM_gal = (sampling_gal * BIO) * basis
    quad_qbf = QuadQBF(splinedegree; oversamplingfactor=2)
    quad_adaptive = QuadAdaptive()
    bemquad_qbf_graded = BEMQuadQBF_Graded(quad_qbf)
    bemquad_qbf_adaptive = BEMQuadQBF_Adaptive(quad_qbf, quad_adaptive)

    z1 = compute_BEM_entry(BEM_col, 2, 2, bemquad_qbf_graded)
    @test abs(z1 - (-1.8800039020044172)) < 1e-3

    assemble!(BEM_col, bemquad_qbf_graded)
    assemble!(BEM_gal, bemquad_qbf_adaptive)
    A_col = copy(matrix(BEM_col))
    b_col = sampling_col * bcond
    coef_col = A_col \ b_col
    density_col = Expansion(basis, coef_col)
    A_gal = copy(matrix(BEM_gal))
    b_gal = sampling_gal * bcond
    coef_gal = A_gal \ b_gal
    density_gal = Expansion(basis, coef_gal)

    point = SVector(0.05, -0.2)
    z_exact = bcond_field(point...)

    z_col = eval_field(BEM_col, density_col, point)
    @test abs(z_col-z_exact) / abs(z_exact) < 1e-3

    z_gal = eval_field(BEM_gal, density_gal, point)
    @test abs(z_gal - z_exact) / abs(z_exact) < 1e-5
end