text:
"""
$(SIGNATURES)

Kinetic central-upwind (KCU) flux
"""
function flux_kcu!(
    face::Interface1F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace1D,
    p,
    dt = 1.0,
) where {T<:ControlVolume1F}

    dxL, dxR = p[1:2]

    flux_kcu!(
        face.fw,
        face.ff,
        ctrL.w + ctrL.sw * dxL,
        ctrL.f + ctrL.sf * dxL,
        ctrR.w - ctrR.sw * dxR,
        ctrR.f - ctrR.sf * dxR,
        vs.u,
        vs.weights,
        gas.K,
        gas.γ,
        gas.μᵣ,
        gas.ω,
        gas.Pr,
        dt,
    )

    return nothing

end

function flux_kcu!(
    face::Interface2F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace1D,
    p,
    dt = 1.0,
) where {T<:ControlVolume2F}

    dxL, dxR = p[1:2]

    if gas isa Mixture
        flux_kcu!(
            face.fw,
            face.fh,
            face.fb,
            ctrL.w + ctrL.sw * dxL,
            ctrL.h + ctrL.sh * dxL,
            ctrL.b + ctrL.sb * dxL,
            ctrR.w - ctrR.sw * dxR,
            ctrR.h - ctrR.sh * dxR,
            ctrR.b - ctrR.sb * dxR,
            vs.u,
            vs.weights,
            gas.K,
            gas.γ,
            gas.mi,
            gas.ni,
            gas.me,
            gas.ne,
            gas.Kn[1],
            dt,
        )
    else
        flux_kcu!(
            face.fw,
            face.fh,
            face.fb,
            ctrL.w + ctrL.sw * dxL,
            ctrL.h + ctrL.sh * dxL,
            ctrL.b + ctrL.sb * dxL,
            ctrR.w - ctrR.sw * dxR,
            ctrR.h - ctrR.sh * dxR,
            ctrR.b - ctrR.sb * dxR,
            vs.u,
            vs.weights,
            gas.K,
            gas.γ,
            gas.μᵣ,
            gas.ω,
            gas.Pr,
            dt,
        )
    end

    return nothing

end

function flux_kcu!(
    face::Interface1F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace2D,
    p,
    dt = 1.0,
) where {T<:ControlVolume1F}

    dxL, dxR, len, n, dirc = p
    swL = @view ctrL.sw[:, dirc]
    swR = @view ctrR.sw[:, dirc]
    sfL = @view ctrL.sf[:, :, dirc]
    sfR = @view ctrR.sf[:, :, dirc]

    flux_kcu!(
        face.fw,
        face.ff,
        local_frame(ctrL.w .+ dxL .* swL, n),
        ctrL.f .+ dxL .* sfL,
        local_frame(ctrR.w .- dxR .* swR, n),
        ctrR.f .- dxR .* sfR,
        vs.u .* n[1] .+ vs.v .* n[2],
        vs.v .* n[1] .- vs.u .* n[2],
        vs.weights,
        gas.K,
        gas.γ,
        gas.μᵣ,
        gas.ω,
        gas.Pr,
        dt,
        len,
    )

    face.fw .= global_frame(face.fw, n)

    return nothing

end

function flux_kcu!(
    face::Interface2F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace2D,
    p,
    dt = 1.0,
) where {T<:ControlVolume2F}

    dxL, dxR, len, n, dirc = p
    swL = @view ctrL.sw[:, dirc]
    swR = @view ctrR.sw[:, dirc]
    shL = @view ctrL.sh[:, :, dirc]
    sbL = @view ctrL.sb[:, :, dirc]
    shR = @view ctrR.sh[:, :, dirc]
    sbR = @view ctrR.sb[:, :, dirc]

    flux_kcu!(
        face.fw,
        face.fh,
        face.fb,
        local_frame(ctrL.w .+ dxL .* swL, n),
        ctrL.h .+ dxL .* shL,
        ctrL.b .+ dxL .* sbL,
        local_frame(ctrR.w .- dxR .* swR, n),
        ctrR.h .- dxR .* shR,
        ctrR.b .- dxR .* sbR,
        vs.u .* n[1] .+ vs.v .* n[2],
        vs.v .* n[1] .- vs.u .* n[2],
        vs.weights,
        gas.K,
        gas.γ,
        gas.μᵣ,
        gas.ω,
        gas.Pr,
        dt,
        len,
    )

    face.fw .= global_frame(face.fw, n)

    return nothing

end

function flux_kcu!(
    face::Interface1F,
    ctrL::T,
    ctrR::T,
    gas::AbstractProperty,
    vs::AbstractVelocitySpace3D,
    p,
    dt = 1.0,
) where {T<:ControlVolume1F}

    dxL, dxR = p[1:2]

    if length(p) == 2 || p[3] == 1
        flux_kcu!(
            face.fw,
            face.ff,
            ctrL.w .+ dxL .* ctrL.sw,
            ctrL.f .+ dxL .* ctrL.sf,
            ctrR.w .- dxR .* ctrR.sw,
            ctrR.f .- dxR .* ctrR.sf,
            vs.u,
            vs.v,
            vs.w,
            vs.weights,
            gas.K,
            gas.γ,
            gas.μᵣ,
            gas.ω,
            gas.Pr,
            dt,
        )
    else

    end

    return nothing

end

"""
$(SIGNATURES)

Kinetic central-upwind (KCU) method

- @args: particle distribution functions and their slopes at left/right sides of interface
- @args: particle velocity quadrature points and weights
- @args: time step and cell size

1D1F1V
"""
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1},Z<:AA{<:Real,1},A<:AA{<:FN,1},B<:AA{<:FN,1}} # 1D1F1V

    # upwind reconstruction
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mxi1, 0, 0)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mxi2, 0, 0)

    w = similar(wL)
    @. w = primL[1] * Muv1 + primR[1] * Muv2

    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g = maxwellian(u, prim)

    fw[1] += Mt[2] * sum(ω .* u .* f)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* f)
    fw[3] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* f))

    @. ff = Mt[1] * u * g + Mt[2] * u * f

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
) where {X<:AA{<:FN,2},Y<:AA{<:FN,2},Z<:AA{<:Real,2},A<:AA{<:FN,2},B<:AA{<:FN,2}}

    # upwind reconstruction
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

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

    # construct interface distribution
    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)

    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)
    #@. tau +=
    #    abs(cellL.prim[1, :] / cellL.prim[end, :] - cellR.prim[1, :] / cellR.prim[end, :]) /
    #    (cellL.prim[1, :] / cellL.prim[end, :] + cellR.prim[1, :] / cellR.prim[end, :]) *
    #    dt *
    #    2.0
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mxi, 1, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    M = mixture_maxwellian(u, prim)

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* f[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* f[:, j])
        fw[3, j] += Mt[2, j] * 0.5 * sum(ω[:, j] .* u[:, j] .^ 3 .* f[:, j])

        @. ff[:, j] = Mt[1, j] * u[:, j] * M[:, j] + Mt[2, j] * u[:, j] * f[:, j]
    end

    return nothing

end

"""
$(SIGNATURES)

1D1F3V
"""
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    uVelo::B,
    vVelo::B,
    wVelo::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len = 1.0,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,3},Z<:AA{<:Real,1},A<:AA{<:FN,3},B<:AA{<:FN,3}} # 1D1F1V

    # upwind reconstruction
    δ = heaviside.(uVelo)
    f = @. fL * δ + fR * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    w = similar(wL)
    @. w = primL[1] * Muv1 + primR[1] * Muv2

    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g = maxwellian(uVelo, vVelo, wVelo, prim)

    fw[1] += Mt[2] * sum(ω .* uVelo .* f)
    fw[2] += Mt[2] * sum(ω .* uVelo .^ 2 .* f)
    fw[3] += Mt[2] * sum(ω .* uVelo .* vVelo .* f)
    fw[4] += Mt[2] * sum(ω .* uVelo .* wVelo .* f)
    fw[5] += Mt[2] * 0.5 * sum(ω .* uVelo .* (uVelo .^ 2 .+ vVelo .^ 2 .+ wVelo .^ 2) .* f)

    @. fw .*= len
    @. ff = (Mt[1] * uVelo * g + Mt[2] * uVelo * f) * len

    return nothing

end

"""
$(SIGNATURES)

1D2F1V
"""
function flux_kcu!(
    fw::X,
    fh::Y,
    fb::Y,
    wL::Z,
    hL::A,
    bL::A,
    wR::Z,
    hR::A,
    bR::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1},Z<:AA{<:Real,1},A<:AA{<:FN,1},B<:AA{<:FN,1}}

    # upwind reconstruction
    δ = heaviside.(u)
    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mxi1, 0, 0)
    Mu2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mxi2, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)

    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mxi, 1, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    Mh = maxwellian(u, prim)
    Mb = Mh .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[2] * sum(ω .* u .* h)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h)
    fw[3] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* h) + sum(ω .* u .* b))

    @. fh = Mt[1] * u * Mh + Mt[2] * u * h
    @. fb = Mt[1] * u * Mb + Mt[2] * u * b

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kcu!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    ω::T5,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
) where {T1<:AA{<:FN,2},T2<:AA{<:FN,2},T3<:AA{<:Real,2},T4<:AA{<:FN,2},T5<:AA{<:FN,2}}

    # upwind reconstruction
    δ = heaviside.(u)

    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

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

    # construct interface distribution
    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)

    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)
    #@. tau +=
    #    abs(cellL.prim[1, :] / cellL.prim[end, :] - cellR.prim[1, :] / cellR.prim[end, :]) /
    #    (cellL.prim[1, :] / cellL.prim[end, :] + cellR.prim[1, :] / cellR.prim[end, :]) *
    #    dt *
    #    2.0
    prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mxi, 1, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    MH = mixture_maxwellian(u, prim)
    MB = similar(MH)
    for j in axes(MB, 2)
        MB[:, j] .= MH[:, j] .* inK ./ (2.0 * prim[end, j])
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* h[:, j])
        fw[3, j] +=
            Mt[2, j] *
            0.5 *
            (sum(ω[:, j] .* u[:, j] .^ 3 .* h[:, j]) + sum(ω[:, j] .* u[:, j] .* b[:, j]))

        @. fh[:, j] = Mt[1, j] * u[:, j] * MH[:, j] + Mt[2, j] * u[:, j] * h[:, j]
        @. fb[:, j] = Mt[1, j] * u[:, j] * MB[:, j] + Mt[2, j] * u[:, j] * b[:, j]
    end

    return nothing

end

"""
$(SIGNATURES)

1D4F1V
"""
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    h3L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    h3R::A,
    u::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,1},Z<:AA{<:Real,1},A<:AA{<:FN,1},B<:AA{<:FN,1}}

    # upwind reconstruction
    δ = heaviside.(u)

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)
    h3 = @. h3L * δ + h3R * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)

    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mw, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    g0 = maxwellian(u, prim)
    g1 = Mv[1] .* g0
    g2 = Mw[1] .* g0
    g3 = (Mv[2] + Mw[2]) .* g0

    fw[1] += Mt[2] * sum(ω .* u .* h0)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h0)
    fw[3] += Mt[2] * sum(ω .* u .* h1)
    fw[4] += Mt[2] * sum(ω .* u .* h2)
    fw[5] += Mt[2] * 0.5 * (sum(ω .* u .^ 3 .* h0) + sum(ω .* u .* h3))

    @. fh0 = Mt[1] * u * g0 + Mt[2] * u * h0
    @. fh1 = Mt[1] * u * g1 + Mt[2] * u * h1
    @. fh2 = Mt[1] * u * g2 + Mt[2] * u * h2
    @. fh3 = Mt[1] * u * g3 + Mt[2] * u * h3

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    fh3::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    h3L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    h3R::A,
    u::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    isMHD = false::Bool,
) where {X<:AA{<:FN,2},Y<:AA{<:FN,2},Z<:AA{<:Real,2},A<:AA{<:FN,2},B<:AA{<:FN,2}}

    #--- upwind reconstruction ---#
    δ = heaviside.(u)

    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)
    h3 = @. h3L * δ + h3R * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mw1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mv1, Mw1, 0, 0, 0)
    Mu2, Mv2, Mw2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv2 = mixture_moments_conserve(MuR2, Mv2, Mw2, 0, 0, 0)

    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)
    #@. tau +=
    #    abs(primL[1, :] / primL[end, :] - primR[1, :] / primR[end, :]) /
    #    (primL[1, :] / primL[end, :] + primR[1, :] / primR[end, :]) *
    #    dt *
    #    5.0
    if !isMHD
        prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)
    end

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate fluxes
    Mu, Mv, Mw, MuL, MuR = mixture_gauss_moments(prim, inK)
    Muv = mixture_moments_conserve(Mu, Mv, Mw, 1, 0, 0)

    ## flux from M0
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    g0 = mixture_maxwellian(u, prim)

    g1 = similar(h0)
    g2 = similar(h0)
    g3 = similar(h0)
    for j in axes(g0, 2)
        g1[:, j] .= Mv[1, j] .* g0[:, j]
        g2[:, j] .= Mw[1, j] .* g0[:, j]
        g3[:, j] .= (Mv[2, j] + Mw[2, j]) .* g0[:, j]
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h0[:, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .^ 2 .* h0[:, j])
        fw[3, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h1[:, j])
        fw[4, j] += Mt[2, j] * sum(ω[:, j] .* u[:, j] .* h2[:, j])
        fw[5, j] +=
            Mt[2, j] *
            0.5 *
            (sum(ω[:, j] .* u[:, j] .^ 3 .* h0[:, j]) + sum(ω[:, j] .* u[:, j] .* h3[:, j]))

        @. fh0[:, j] = Mt[1, j] * u[:, j] * g0[:, j] + Mt[2, j] * u[:, j] * h0[:, j]
        @. fh1[:, j] = Mt[1, j] * u[:, j] * g1[:, j] + Mt[2, j] * u[:, j] * h1[:, j]
        @. fh2[:, j] = Mt[1, j] * u[:, j] * g2[:, j] + Mt[2, j] * u[:, j] * h2[:, j]
        @. fh3[:, j] = Mt[1, j] * u[:, j] * g3[:, j] + Mt[2, j] * u[:, j] * h3[:, j]
    end

    return nothing

end

"""
$(SIGNATURES)

2D1F2V
"""
function flux_kcu!(
    fw::X,
    ff::Y,
    wL::Z,
    fL::A,
    wR::Z,
    fR::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,2},Z<:AA{<:Real,1},A<:AA{<:FN,2},B<:AA{<:FN,2}}

    # reconstruct initial distribution
    δ = heaviside.(u)
    f = @. fL * δ + fR * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv * len

    ## flux from f0
    g = maxwellian(u, v, prim)

    fw[1] += Mt[2] * sum(ω .* u .* f) * len
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* f) * len
    fw[3] += Mt[2] * sum(ω .* v .* u .* f) * len
    fw[4] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* f)) * len

    @. ff = (Mt[1] * u * g + Mt[2] * u * f) * len

    return nothing

end

"""
$(SIGNATURES)

2D2F2V
"""
function flux_kcu!(
    fw::T1,
    fh::T2,
    fb::T2,
    wL::T3,
    hL::T4,
    bL::T4,
    wR::T3,
    hR::T4,
    bR::T4,
    u::T5,
    v::T5,
    ω::T5,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {T1<:AA{<:FN,1},T2<:AA{<:FN,2},T3<:AA{<:Real,1},T4<:AA{<:FN,2},T5<:AA{<:FN,2}}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h = @. hL * δ + hR * (1.0 - δ)
    b = @. bL * δ + bR * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv * len

    ## flux from f0
    H = maxwellian(u, v, prim)
    B = H .* inK ./ (2.0 * prim[end])

    fw[1] += Mt[2] * sum(ω .* u .* h) * len
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h) * len
    fw[3] += Mt[2] * sum(ω .* v .* u .* h) * len
    fw[4] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h) + sum(ω .* u .* b)) * len

    @. fh = (Mt[1] * u * H + Mt[2] * u * h) * len
    @. fb = (Mt[1] * u * B + Mt[2] * u * b) * len

    return nothing

end

"""
$(SIGNATURES)

2D3F2V
"""
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    visRef,
    visIdx,
    Pr,
    dt,
    len,
) where {X<:AA{<:FN,1},Y<:AA{<:FN,2},Z<:AA{<:Real,1},A<:AA{<:FN,2},B<:AA{<:FN,2}}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = gauss_moments(primL, inK)
    Muv1 = moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = gauss_moments(primR, inK)
    Muv2 = moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    w = @. primL[1] * Muv1 + primR[1] * Muv2
    prim = conserve_prim(w, γ)
    tau = vhs_collision_time(prim, visRef, visIdx)
    tau +=
        abs(primL[1] / primL[end] - primR[1] / primR[end]) /
        (primL[1] / primL[end] + primR[1] / primR[end]) *
        dt *
        2.0

    Mt = zeros(2)
    Mt[2] = tau * (1.0 - exp(-dt / tau)) # f0
    Mt[1] = dt - Mt[2] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = gauss_moments(prim, inK)

    ## flux from M0
    Muv = moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    @. fw = Mt[1] * prim[1] * Muv

    ## flux from f0
    H0 = maxwellian(u, v, prim)
    H1 = H0 .* prim[4]
    H2 = H0 .* (prim[4]^2 + 1.0 / (2.0 * prim[5]))

    fw[1] += Mt[2] * sum(ω .* u .* h0)
    fw[2] += Mt[2] * sum(ω .* u .^ 2 .* h0)
    fw[3] += Mt[2] * sum(ω .* v .* u .* h0)
    fw[4] += Mt[2] * sum(ω .* u .* h1)
    fw[5] += Mt[2] * 0.5 * (sum(ω .* u .* (u .^ 2 .+ v .^ 2) .* h0) + sum(ω .* u .* h2))

    @. fw *= len
    @. fh0 = (Mt[1] * u * H0 + Mt[2] * u * h0) * len
    @. fh1 = (Mt[1] * u * H1 + Mt[2] * u * h1) * len
    @. fh2 = (Mt[1] * u * H2 + Mt[2] * u * h2) * len

    return nothing

end

"""
$(SIGNATURES)

Mixture
"""
function flux_kcu!(
    fw::X,
    fh0::Y,
    fh1::Y,
    fh2::Y,
    wL::Z,
    h0L::A,
    h1L::A,
    h2L::A,
    wR::Z,
    h0R::A,
    h1R::A,
    h2R::A,
    u::B,
    v::B,
    ω::B,
    inK,
    γ,
    mi,
    ni,
    me,
    ne,
    Kn,
    dt,
    len,
    isMHD = false::Bool,
) where {X<:AA{<:FN,2},Y<:AA{<:FN,3},Z<:AA{<:Real,2},A<:AA{<:FN,3},B<:AA{<:FN,3}}

    # reconstruct initial distribution
    δ = heaviside.(u)
    h0 = @. h0L * δ + h0R * (1.0 - δ)
    h1 = @. h1L * δ + h1R * (1.0 - δ)
    h2 = @. h2L * δ + h2R * (1.0 - δ)

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

    # construct interface distribution
    Mu1, Mv1, Mxi1, MuL1, MuR1 = mixture_gauss_moments(primL, inK)
    Muv1 = mixture_moments_conserve(MuL1, Mv1, Mxi1, 0, 0, 0)
    Mu2, Mv2, Mxi2, MuL2, MuR2 = mixture_gauss_moments(primR, inK)
    Muv2 = mixture_moments_conserve(MuR2, Mv2, Mxi2, 0, 0, 0)

    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)
    #@. tau +=
    #    abs(primL[1, :] / primL[end, :] - primR[1, :] / primR[end, :]) /
    #    (primL[1, :] / primL[end, :] + primR[1, :] / primR[end, :]) *
    #    dt *
    #    2.0
    if !isMHD
        prim = aap_hs_prim(prim, tau, mi, ni, me, ne, Kn)
    end

    Mt = zeros(2, 2)
    @. Mt[2, :] = tau * (1.0 - exp(-dt / tau)) # f0
    @. Mt[1, :] = dt - Mt[2, :] # M0

    # calculate interface flux
    Mu, Mv, Mxi, MuL, MuR = mixture_gauss_moments(prim, inK)

    ## flux from M0
    Muv = mixture_moments_conserve(Mu, Mv, Mxi, 1, 0, 0)
    for j in axes(fw, 2)
        @. fw[:, j] = Mt[1, j] * prim[1, j] * Muv[:, j]
    end

    ## flux from f0
    H0 = mixture_maxwellian(u, v, prim)
    H1 = similar(H0)
    H2 = similar(H0)
    for j in axes(H0, 3)
        H1[:, :, j] = H0[:, :, j] .* prim[4, j]
        H2[:, :, j] .= H0[:, :, j] .* (prim[4, j]^2 + 1.0 / (2.0 * prim[5, j]))
    end

    for j in axes(fw, 2)
        fw[1, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .* h0[:, :, j])
        fw[2, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .^ 2 .* h0[:, :, j])
        fw[3, j] += Mt[2, j] * sum(ω[:, :, j] .* v[:, :, j] .* u[:, :, j] .* h0[:, :, j])
        fw[4, j] += Mt[2, j] * sum(ω[:, :, j] .* u[:, :, j] .* h1[:, :, j])
        fw[5, j] +=
            Mt[2, j] *
            0.5 *
            (
                sum(
                    ω[:, :, j] .* u[:, :, j] .* (u[:, :, j] .^ 2 .+ v[:, :, j] .^ 2) .*
                    h0[:, :, j],
                ) + sum(ω[:, :, j] .* u[:, :, j] .* h2[:, :, j])
            )

        @. fh0[:, :, j] =
            Mt[1, j] * u[:, :, j] * H0[:, :, j] + Mt[2, j] * u[:, :, j] * h0[:, :, j]
        @. fh1[:, :, j] =
            Mt[1, j] * u[:, :, j] * H1[:, :, j] + Mt[2, j] * u[:, :, j] * h1[:, :, j]
        @. fh2[:, :, j] =
            Mt[1, j] * u[:, :, j] * H2[:, :, j] + Mt[2, j] * u[:, :, j] * h2[:, :, j]
    end

    @. fw *= len
    @. fh0 *= len
    @. fh1 *= len
    @. fh2 *= len

    return nothing

end


text:
function HLop!(HL,sbasis,vbasis,quad,linemap::LineMap,normals,components,
    ED,scale)

    @assert length(normals) == length(ED)
    @assert length(components) == length(ED)
    NED = [normals[k]*ED[k] for k = 1:length(normals)]

    detjac = determinant_jacobian(linemap)
    projector = components*components'

    for (idx,(p,w)) in enumerate(quad)
        svals = sbasis(p)
        vvals = vbasis(linemap(p...))
        NP = make_row_matrix(svals,projector)
        for nk in NED
            NK = make_row_matrix(vvals,nk)
            HL .+= NP'*NK*detjac*scale*w
        end
    end
end

function HLop!(HL,sbasis,vbasis,quad,imap::InterpolatingPolynomial,
    normals,components,ED,scale)

    NQ = length(quad)
    @assert length(scale) == NQ
    @assert size(normals)[2] == NQ
    @assert size(components)[2] == NQ
    @assert size(normals)[1] == length(ED)

    for (idx,(p,w)) in enumerate(quad)
        n = normals[:,idx]
        t = components[:,idx]
        projector = t*t'
        svals = sbasis(p)
        vvals = vbasis(imap(p))
        detjac = determinant_jacobian(imap,p)
        NP = make_row_matrix(svals,projector)
        for k in 1:length(ED)
            NED = n[k]*ED[k]
            NK = make_row_matrix(vvals,NED)

            HL .+= NP'*NK*detjac*scale[idx]*w
        end
    end
end

function HLop(sbasis,vbasis,facequad,facemap,normals,components,
    Dhalf,scale)

    dim = dimension(vbasis)

    NHF = number_of_basis_functions(sbasis)
    NLF = number_of_basis_functions(vbasis)

    sdim = symmetric_tensor_dimension(dim)
    Ek = vec_to_symm_mat_converter(dim)
    ED = [Ek[k]'*Dhalf for k = 1:dim]

    HL = zeros(dim*NHF,sdim*NLF)

    HLop!(HL,sbasis,vbasis,facequad,facemap,normals,components,ED,scale)

    return HL
end

function HUop!(HU,sbasis,vbasis,quad,linemap::LineMap,components,scale,
    nhdofs)

    @assert length(components) == nhdofs
    detjac = determinant_jacobian(linemap)
    projector = components*components'

    for (idx,(p,w)) in enumerate(quad)
        svals = sbasis(p)
        vvals = vbasis(linemap(p...))
        NP = make_row_matrix(svals,projector)
        NI = interpolation_matrix(vvals,nhdofs)

        HU .+= scale*NP'*NI*detjac*w
    end
end

function HUop!(HU,sbasis,vbasis,quad,imap::InterpolatingPolynomial,
    components,scale,nhdofs)

    NQ = length(quad)
    @assert length(scale) == NQ
    @assert size(components) == (nhdofs,NQ)

    for (idx,(p,w)) in enumerate(quad)
        t = components[:,idx]
        projector = t*t'
        svals = sbasis(p)
        vvals = vbasis(imap(p))
        detjac =  determinant_jacobian(imap,p)
        NP = make_row_matrix(svals,projector)
        NI = interpolation_matrix(vvals,nhdofs)

        HU .+= scale[idx]*NP'*NI*detjac*w
    end

end

function HUop(sbasis,vbasis,facequad,facemap,components,scale)

    dim = dimension(vbasis)

    NHF = number_of_basis_functions(sbasis)
    NLF = number_of_basis_functions(vbasis)

    HU = zeros(dim*NHF,dim*NLF)

    HUop!(HU,sbasis,vbasis,facequad,facemap,components,scale,dim)

    return HU

end

function hybrid_local_operator_traction_components(sbasis,vbasis,facequad,
    facemap,normals,components,Dhalf,stabilization,facescale)

    scale = stabilization*facescale

    HL = HLop(sbasis,vbasis,facequad,facemap,normals,components,Dhalf,scale)
    HU = HUop(sbasis,vbasis,facequad,facemap,components,scale)

    return [HL HU]
end

ext: jl


text:
using BenchmarkTools
using LightPropagation

suite = BenchmarkGroup()

suite["Layered_fluence_CW"] = BenchmarkGroup()
suite["Layered_fluence_CW"]["2-layer"] = @benchmarkable fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, N_J0Roots) setup=(ρ = 1.0; μa = (0.1, 0.2); μsp = (10.0, 11.0); n_ext = 1.0; n_med = (1.4, 1.4); l = (1.0, 10.0); a = 10.0; z = 0.0; N_J0Roots = 1000)
suite["Layered_fluence_CW"]["3-layer"] = @benchmarkable fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, N_J0Roots) setup=(ρ = 1.0; μa = (0.1, 0.2, 0.15); μsp = (10.0, 11.0, 12.5); n_ext = 1.0; n_med = (1.4, 1.4, 1.4); l = (1.0, 2.0, 10.0); a = 10.0; z = 0.0; N_J0Roots = 1000)
suite["Layered_fluence_CW"]["4-layer"] = @benchmarkable fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, N_J0Roots) setup=(ρ = 1.0; μa = (0.1, 0.2, 0.15, 0.2); μsp = (10.0, 11.0, 12.5, 12.0); n_ext = 1.0; n_med = (1.4, 1.4, 1.4, 1.4); l = (1.0, 2.0, 2.0, 10.0); a = 10.0; z = 0.0; N_J0Roots = 1000)
suite["Layered_fluence_CW"]["10-layer"] = @benchmarkable fluence_DA_Nlay_cylinder_CW(ρ, μa, μsp, n_ext, n_med, l, a, z, N_J0Roots) setup=(ρ = 1.0; μa = (0.1, 0.2, 0.15, 0.2, 0.1, 0.2, 0.15, 0.2, 0.1, 0.2); μsp = (10.0, 11.0, 12.5, 12.0, 10.0, 11.0, 12.5, 12.0, 10.0, 11.0); n_ext = 1.0; n_med = (1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4, 1.4); l = (1.0, 2.0, 2.0, 1.0, 1.0, 2.0, 2.0, 1.0, 2.0, 2.0); a = 10.0; z = 0.0; N_J0Roots = 1000)

suite["semiinf_fluence_CW"] = BenchmarkGroup()
suite["semiinf_fluence_CW"] = @benchmarkable fluence_DA_semiinf_CW(ρ, μa, μsp) setup=(ρ = 1.0*rand(); μa = 0.1*rand(); μsp = 10.0*rand())

tune!(suite)


text:
module TestTrees

using Phylo
using DataFrames
#using JuliaDB

using Compat.Test
using IterableTables: getiterator

species = ["Dog", "Cat", "Human"]
ntips = 10
df = DataFrame(species = species, count = [10, 20, 3])
observations = ["Dog", "Cat", "Dog", "Dog"]
#jdb = table(@NT(species = observations, count = 1:4))
jdb = DataFrame(species = observations, count = 1:4)

@testset "NamedTree()" begin
    ntn = NamedTree(ntips)
    @test length(nodefilter(isroot, ntn)) == ntips
    @test length(nodefilter(isleaf, ntn)) == ntips
    @test length(nodefilter(isinternal, ntn)) == 0
    @test validate(ntn)
    nt = NamedTree(species)
    @test_throws ErrorException noderoute(nt, "Dog", "Cat")
    @test_throws ErrorException branchroute(nt, "Dog", "Human")
    @test validate(nt)
    n = addnode!(nt)
    b1 = addbranch!(nt, n, "Dog", 2.0)
    b2 = addbranch!(nt, n, "Cat", 2.0)
    @test_throws ErrorException addbranch!(nt, n, "Human", 2.0)
    @test validate(nt)
    r = addnode!(nt)
    @test_throws ErrorException addbranch!(nt, r, "Potato", 2.0)
    b3 = addbranch!(nt, r, "Human", 5.0)
    b4 = addbranch!(nt, r, n, 3.0)
    @test maximum(distances(nt)) ≈ 10.0
    @test validate(nt)
    @test noderoute(nt, "Human", "Dog") == ["Human", r, n, "Dog"]
    @test branchroute(nt, "Human", "Dog") == [b3, b4, b1]
    nb = BinaryTree(DataFrame(name=["Human", "Cat"]))
    @test getleafinfo(nb) === getleafinfo(BinaryTree(nb))
    @test getleafinfo(nb) !== getleafinfo(BinaryTree(nb; copyinfo = true))
    @test_nowarn addnode!(nb, "Dog")
    if VERSION < v"0.7.0-"
        @test_warn "LeafInfo names do not match actual leaves of tree" !validate(nb) || error("validate() should have returned false")
    else
        @test !validate(nb)
    end
    np = PolytomousTree(DataFrame(name=["Human", "Cat"]))
    @test getleafinfo(np) === getleafinfo(PolytomousTree(np))
    @test getleafinfo(np) !== getleafinfo(PolytomousTree(np; copyinfo = true))
    @test_nowarn addnode!(np, "Dog")
    if VERSION < v"0.7.0-"
        @test_warn "LeafInfo names do not match actual leaves of tree" !validate(np) || error("validate() should have returned false")
    else
        @test !validate(np)
    end
    @test_throws ErrorException PolytomousTree(np)
    @test getleafnames(nb) ⊆ getleafnames(np)
    @test getleafinfo(BinaryTree(df)) == df
    @test getleafinfo(PolytomousTree(df)) == df
end

@testset "BinaryTree()" begin
    btn = BinaryTree{DataFrame, Vector{Float64}}(ntips)
    @test length(nodefilter(isroot, btn)) == ntips
    @test length(nodefilter(isleaf, btn)) == ntips
    @test length(nodefilter(isinternal, btn)) == 0

    tmp = BinaryTree(DataFrame(names=species))
    @test validate(tmp)
    @test_nowarn addbranch!(tmp, addnode!(tmp), species[1])
    @test Set(getleafnames(tmp)) == Set(species)
    nt = BinaryTree{DataFrame, Vector{Float64}}(species)
    @test validate(nt)
    n = addnode!(nt)
    @test Set(getleafnames(nt)) == Set(species) ∪ Set([n])
    addbranch!(nt, n, "Dog", 2.0)
    addbranch!(nt, n, "Cat", 2.0)
    @test_throws ErrorException addbranch!(nt, n, "Human", 2.0)
    @test validate(nt)
    r = addnode!(nt)
    @test_throws ErrorException addbranch!(nt, r, "Potato", 2.0)
    addbranch!(nt, r, "Human", 4.0)
    addbranch!(nt, r, n, 2.0)
    @test validate(nt)
    setnoderecord!(nt, "Dog", [1.0])
    @test getnoderecord(nt, "Dog")[1] ≈ 1.0
    @test length(getnoderecord(nt, "Cat")) == 0
    nt2 = BinaryTree(nt)
    @test length(getnoderecord(nt2, "Dog")) == 0
    @test Set(getleafnames(nt2)) == Set(species)
    @test validate(nt)
    nt3 = BinaryTree(nt, empty=false)
    a=IOBuffer()
    b=IOBuffer()
    show(a, nt)
    show(b, nt3)
    @test a.data == b.data

    tdf = BinaryTree(df)
    @test nleaves(tdf) == 3
    tj = BinaryTree(jdb)
    @test nleaves(tj) == 2

    lij = getiterator(getleafinfo(tj, "Dog"))
    @test sum(map(line -> line.count, lij)) == 8

    lij = getiterator(getleafinfo(tj, "Cat"))
    @test sum(map(line -> line.count, lij)) == 2

    lidf = getiterator(getleafinfo(tdf, "Dog"))
    @test sum(map(line -> line.count, lidf)) == 10
end

end


text:
# This code is duplicated from TerminalExtensions
# Maybe they will be reunified in the future

function prepare_display_file(buf;filename="Unnamed file", size=nothing, width=nothing, height=nothing, preserveAspectRatio::Bool=true, inline::Bool=false)
    q = "\e]1337;File="
    options = String[]
    filename != "Unnamed file" && push!(options,"name=" * base64encode(filename))
    size !== nothing && push!(options,"size=" * dec(size))
    height !== nothing && push!(options,"height=" * dec(height))
    width !== nothing && push!(options,"width=" * dec(width))
    preserveAspectRatio !== true && push!(options,"preserveAspectRatio=0")
    inline !== false && push!(options,"inline=1")
    q *= join(options,';')
    q *= ":"
    write(buf,q)
end

function display_file(buf,data::Vector{UInt8}; kwargs...)
    prepare_display_file(buf;kwargs...)
    write(buf,base64encode(data))
    write(buf,'\a')
end

#=
# Incomplete list. Will be extended as necessity comes up
const iterm2_mimes = ["image/png", "image/gif", "image/jpeg", "application/pdf", "application/eps"]

for mime in iterm2_mimes
    @eval begin
        function display(d::InlineDisplay, m::MIME{symbol($mime)}, x)
            prepare_display_file(;filename="image",inline=true)
            buf = IOBuffer()
            pipe = Base64EncodePipe(buf)
            writemime(pipe,m,x)
            close(pipe)
            write(STDOUT, takebuf_array(buf))
            write(STDOUT,'\a')
        end
    end
end
=#


text:
function analyze_icn()
    @info """Start the analyze of files in $datadir("compositions")"""
    df = DataFrame()
    for f in readdir(datadir("compositions"); join = true)
        d = JSON.parsefile(f)
        inds1 = Indices(["icn_time", "nthreads"])
        d1 = view(Dictionary(d),inds1)
        inds2 = Indices(["icn_iterations", "loss_sampler", "domains_size", "memoize", "metric", "sampling", "generations", "search", "population", "concept"])
        d2 = view(Dictionary(d["params"]), inds2)
        d3 = merge(d1, d2)
        d3["loss_sampler"] = string(d3["loss_sampler"])

        if isempty(df)
            df = DataFrame(Dict(pairs(d3)))
        else
            push!(df, Dict(pairs(d3)))
        end
    end
    return df
end

"""
icn_iterations, time, composition_number, rsd, mean, std, icn_time, memoize, var, concept, partial_search_limit, sampling, generations
cov, complete_search_limit, median, search, population, selection_rate
{
  "icn_iterations": 16,
  "time": 0.356123463,
  "composition_number": 1,

  "rsd": 2.8284271247461907,
  "maths": "all_different = identity ∘ sum ∘ sum ∘ count_eq_left",

  "mean": 0.037037037037037035,
  "std": 0.10475656017578483,
  "icn_time": 5.040572709,
  "memoize": false,
  "composition": "\"\"\"\n    all_different(x; X = zeros(length(x), 1), param=nothing, dom_size)\n\nComposition `all_different` generated by CompositionalNetworks.jl.\n```\nall_different = identity ∘ sum ∘ sum ∘ count_eq_left\n```\n\"\"\"\nfunction all_different(x; X=zeros(length(x), 1), param=nothing, dom_size)\n    CompositionalNetworks.tr_in(\n        Tuple([CompositionalNetworks.tr_count_eq_left]), X, x, param\n    )\n    for i in 1:length(x)\n        X[i, 1] = CompositionalNetworks.ar_sum(@view X[i, :])\n    end\n    return CompositionalNetworks.ag_sum(@view X[:, 1]) |>\n           (y -> CompositionalNetworks.co_identity(y; param, dom_size, nvars=length(x)))\nend\n",
  "var": 0.010973936899862827,
  "concept": "all_different",
  "partial_search_limit": 256,
  "sampling": 100,
  "generations": 16,
  "cov": 0.010973936899862827,
  "complete_search_limit": 46656,
  "median": 0.0,
  "search": "complete",
  "population": 1024,
  "selection_rate": 0.9375
}

"""

function analyze_composition()
    @info """Start the analyze of files in $datadir("composition_results")"""
    df = DataFrame()
    for f in filter(x -> endswith(x, ".json"),readdir(datadir("composition_results"); join = true))
        d = JSON.parsefile(f)

        inds1 = Indices(["symbols_count","icn_iterations", "time", "composition_number", "mean", "std", "icn_time", "memoize", "var", "concept", "partial_search_limit", "sampling", "generations",
        "cov", "complete_search_limit", "median", "search", "population", "selection_rate"])
        d1 = view(Dictionary(d),inds1)

        if isempty(df)
            df = DataFrame(Dict(pairs(d1)))
        else
            push!(df, Dict(pairs(d1)))
        end
    end
    return df
end

visualize_icn() = Voyager(analyze_icn())

visualize_compositions() = Voyager(analyze_composition())


text:
export linscan_pq, linscan_opq, linscan_lsq, eval_recall, linscan_cq

# Linear scan using PQ codebooks no rotation
function linscan_pq(
  B::Matrix{UInt8},          # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b:: Int,                   # Number of bits per code -- log2(h) * m
  k:: Int = 10000)           # Number of knn results to return

  m, n  = size(B)
  d, nq = size(X)

  @show k, nq
  dists = zeros(Cfloat, k, nq)
  res   = zeros( Cuint, k, nq)

  ccall(("linscan_aqd_query", linscan_aqd), Nothing,
    (Ptr{Cfloat}, Ptr{Cuint}, Ptr{Cuchar}, Ptr{Cfloat},
    Ptr{Cfloat}, Cint, Cuint, Cint, Cint, Cint, Cint, Cint),
    dists, res, B, cat(C..., dims=3), X, Cint(n), Cuint(nq),
    Cint(b), Cint(k), Cint(m), Cint(d), Cint(d/m))

  return dists, (res.+=1)
end

function linscan_pq(
  B::Matrix{T},          # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b:: Int,                   # Number of bits per code -- log2(h) * m
  k:: Int = 10000)  where T <: Integer

  B_uint8 = convert(Matrix{UInt8}, B.-1)
  return linscan_pq(B_uint8, X, C, b, k)
end

"My attempt to implement linscan in julia. Very slow because sortperm is slow and no multithreading"
function linscan_pq_julia(
  B::Matrix{Int16},          # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  knn:: Int = 10000)           # Number of knn results to return

  m, n  = size( B )
  d, nq = size( X )
  subdims = splitarray( 1:d, m )

  @show knn, nq
  dists = zeros( Cfloat, knn, nq )
  idx   = zeros(  Cuint, knn, nq )
  Bt = B'

  # Compute distance tables between queries and codebooks
  tables = Vector{Matrix{Cfloat}}(m)
  for i = 1:m
    # tables[i] = Distances.pairwise(Distances.SqEuclidean(), X[subdims[i],:], C[i])
    tables[i] = Distances.pairwise(Distances.SqEuclidean(), C[i], X[subdims[i],:])
  end

  # Compute approximate distances and sort
  @inbounds for i = 1:nq # Loop over each query

    # println(i)
    xq_dists  = zeros(Cfloat, n)
    p         = zeros(Cuint, n)

     for j = 1:m # Loop over each codebook
      t = tables[j][:,i]

      for k = 1:n # Loop over each code
        xq_dists[k] += t[ Bt[k,j] ]
      end
    end

    # p = sortperm(xq_dists; alg=PartialQuickSort(knn))
    sortperm!(p, xq_dists; alg=PartialQuickSort(knn))
    # @simd for j=1:n; p[j]=j; end
    # sortperm!(p, xq_dists; alg=PartialQuickSort(knn), initialized=true)
    # sort(xq_dists; alg=PartialQuickSort(knn))

    dists[:,i] = xq_dists[ p[1:knn] ]
    idx[:,i]   = p[1:knn]

  end # @inbounds

  return dists, idx

end

# Linear scan using OPQ
function linscan_opq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b::Int,                     # Number of bits per code -- log2(h) * m
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000)             # Number of knn results to return

  # Rotate and call the function as usual
  return linscan_pq( B, R'*X, C, b, k )
end

function linscan_opq(
  B::Matrix{T},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # The cluster centers
  b::Int,                     # Number of bits per code -- log2(h) * m
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000) where T <: Integer # Number of knn results to return

  B_uint8 = convert(Matrix{UInt8}, B.-1)
  return linscan_opq( B_uint8, X, C, b, R, k )
end

# Linear scan using LSQ, with dbnorms not encoded.
function linscan_lsq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # (m-1)-long. The cluster centers
  dbnorms::Vector{Cfloat},   # h-long. Table with database norms
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000)             # Number of knn results to return

  RX = R' * X;

  m, n  = size(B)
  d, nq = size(RX)
  _, h  = size(C[1])

  dists = zeros(Cfloat, k, nq)
  res   = zeros( Cuint, k, nq)

  ccall(("linscan_aqd_query_extra_byte", linscan_aqd_pairwise_byte), Nothing,
    (Ptr{Cfloat}, Ptr{Cint},
    Ptr{Cuchar}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat},
    Cuint, Cint, Cint, Cint, Cint, Cint),
    dists, res,
    B, RX, hcat(C...), dbnorms,
    Cint(nq), Cint(n), Cint(m), Cint(h), Cint(d), Cint(k));

  return dists, res
end

# Linear scan using LSQ, with dbnorms not encoded.
function linscan_lsq(
  B::Matrix{T},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # (m-1)-long. The cluster centers
  dbnorms::Vector{Cfloat},   # h-long. Table with database norms
  R::Matrix{Cfloat},         # Rotation matrix
  k::Int = 10000) where T <: Integer # Number of knn results to return

  B_uint8 = convert(Matrix{UInt8}, B.-1)
  return linscan_lsq(B_uint8, X, C, dbnorms, R, k)
end

# Linear scan using Composite Quantization (Zhang et al, ICML 14)
function linscan_cq(
  B::Matrix{UInt8},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # (m-1)-long. The cluster centers
  k::Int = 10000)             # Number of knn results to return

  m, n  = size( B );
  d, nq = size( X );
  _, h  = size( C[1] );

  dists = zeros( Cfloat, k, nq );
  res   = zeros(  Cuint, k, nq  );

  ccall(("linscan_aqd_cq_query_extra_byte", linscan_aqd_pairwise_byte), Nothing,
    (Ptr{Cfloat}, Ptr{Cint},
    Ptr{Cuchar}, Ptr{Cfloat}, Ptr{Cfloat}, #Ptr{Cfloat},
    Cuint, Cint, Cint, Cint, Cint, Cint),
    dists, res,
    B, X, hcat(C...), #dbnorms,
    Cint(nq), Cint(n), Cint(m), Cint(h), Cint(d), Cint(k) );

  return dists, res
end

# Linear scan using Composite Quantization (Zhang et al, ICML 14)
function linscan_cq(
  B::Matrix{T},           # m-by-n. The database, encoded
  X::Matrix{Cfloat},         # d-by-nq. The queries.
  C::Vector{Matrix{Cfloat}}, # (m-1)-long. The cluster centers
  k::Int = 10000) where T <: Integer # Number of knn results to return

  B_uint8 = convert(Matrix{UInt8},B-1)
  return linscan_cq(B_uint8, X, C, k)
end

# Evaluate ANN search vs ground truth. Produces a recall@N curve.
function eval_recall(ids_gnd::Vector{T}, ids_predicted::Matrix{T}, k::Integer) where T <: Integer

  # Modified from Hervé Jégou's test_compute_stats.m matlab code.

  nquery = size(ids_predicted, 2)
  @assert nquery == length(ids_gnd)

  nn_ranks = zeros(nquery)
  hist_pqc = zeros(k+1)

  for i = 1:nquery
    gnd_ids = ids_gnd[i]
    nn_pos = findall(ids_predicted[:,i] .== gnd_ids)

    if length(nn_pos) == 1
      nn_ranks[i] = nn_pos[1]
    else
      nn_ranks[i] = k+1
    end
  end

  nn_ranks = sort(nn_ranks)

  recall_at_i = zeros(k)

  for i = [1 2 5 10 20 50 100 200 500 1000 2000 5000 10000]
    if i <= k
      recall_at_i[i] = length(findall((nn_ranks .<= i) .& (nn_ranks .<= k) )) ./ nquery * 100
      println("r@$(i) = $(recall_at_i[i])")
    end
  end

  for i = 1:k
    recall_at_i[i] = length(findall((nn_ranks .<= i) .& (nn_ranks .<= k) )) ./ nquery
  end

  return recall_at_i

end


text:
using HttpServer
using WebSockets
using JSON

struct User
    name::String
    client::WebSocket
end
#global Dict to store open connections in
global connections = Dict{String,User}()

function decodeMessage( msg )
    JSON.parse(String(copy(msg)))
end

wsh = WebSocketHandler() do req, client
    global connections
    while true
        msg = read(client)
        msg = decodeMessage(msg)
        id = msg["id"]
        if haskey(msg,"userName") && !haskey(connections,id)
            uname = msg["userName"]
            println("SETTING USERNAME: $(uname)")
            connections[id] = User(uname,client)
        end
        if haskey(msg,"say")
            content = msg["say"]
            println("EMITTING MESSAGE: $(content)")
            for (k,v) in connections
                if k != id
                    write(v.client, (v.name * ": " * content))
                end
            end
        end
    end
end

httph = HttpHandler() do req::Request, res::Response
    onepage = readstring(Pkg.dir("WebSockets","examples","chat-client.html"))
    Response(onepage)
end

server = Server(httph, wsh)
println("Chat server listening on 8000...")
run(server,8000)


text:
using ResumableFunctions, SimJulia, BenchmarkTools

@resumable function fibonnaci(sim::Simulation)
  a = 0.0
  b = 1.0
  while true
    @yield timeout(sim, 1)
    a, b = b, a+b
  end
end

function run_test()
  sim = Simulation()
  @process fibonnaci(sim)
  run(sim, 10)
end

run_test()
@btime run_test()


text:
module adapters_postgresql_copy_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo

Repo.debug_sql()

include("options.jl")

Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    Options.for_postgresql...
)

# https://discourse.julialang.org/t/postgresql-in-julia-libpq-jl/9379/4
Repo.execute(Raw("CREATE TABLE IF NOT EXISTS a_uts (ut text, id serial)"))

using LibPQ.libpq_c: PQputCopyData, PQputCopyEnd
loader = Repo.current_loader()
jl_conn = loader.current_conn()

Repo.execute(Raw("COPY a_uts (ut, id) FROM STDIN"))
d = [["W:000060362500001",1], ["W:000060362500002",2], ["W:000070603200027",3]]
for (ut, id) in d
    buf = string(join((ut, id), '\t'), '\n')
    PQputCopyData(jl_conn.conn, pointer(buf), Cint(length(buf)))
end
PQputCopyEnd(jl_conn.conn, C_NULL)

struct UTS
end
Schema.model(UTS, table_name="a_uts")

df = Repo.query(UTS)
@test size(df) == (3,)

Repo.execute([DROP TABLE :a_uts])

Repo.disconnect()

end # module adapters_postgresql_copy_test


text:
using ConfParser
using Random
using GLM
using Distributions
using StatsBase
using Plots

include("dictionaries.jl")
include("structs.jl")

function generate_config(conpath::String="none")
    """
    Takes in the path of a config file and produces sim_dict, a dictionary of
    experimental parameters.
    """
    if conpath == "none" # DEFAULT config is sim.cfg, currently in src
        print("using default config...")
        conpath = string(@__DIR__) * "\\sim.cfg"
    end
    conf = ConfParse(conpath)
    parse_conf!(conf)

    sim_dict = Dict{String,Any}()
    sim_dict["random_seed"] = parse(Int64,retrieve(conf, "RANDOM_SEED"))
    sim_dict["genes"] = parse(Int64,retrieve(conf, "GENES"))
    sim_dict["ancestor_connectivity"] = parse(Float64,retrieve(conf, "ANCESTOR_CONNECTIVITY"))
    sim_dict["interaction_type"] = retrieve(conf,"INTERACTION_TYPE")
    sim_dict["mu_gain"]= parse(Float64,retrieve(conf, "GAIN_MUTATION_RATE"))
    sim_dict["mu_loss"] = parse(Float64,retrieve(conf, "LOSS_MUTATION_RATE"))
    sim_dict["mu_change"] = parse(Float64,retrieve(conf, "CHANGE_MUTATION_RATE"))
    sim_dict["treatment"] = retrieve(conf,"TREATMENT")
    sim_dict["max_convergence_time"] = parse(Int64,retrieve(conf, "MAX_CONVERGENCE_TIME"))
    sim_dict["activation"] = parse(Float64,retrieve(conf, "INIT_COND_PROB_ACTIVATED"))

    sim_dict["experiment_type"] = retrieve(conf,"EXPERIMENT_TYPE")
    sim_dict["selection_method"] = retrieve(conf,"SELECTION_METHOD")
    sim_dict["sigma"] = parse(Float64,retrieve(conf, "SIGMA"))

    sim_dict["population_size"] = parse(Int64,retrieve(conf, "POPULATION_SIZE"))

    sim_dict["walk_length"] = parse(Int64,retrieve(conf,"WALK_LENGTH"))
    sim_dict["reverse_signal"] = parse(Bool,retrieve(conf,"REVERSE_SIGNAL"))

    sim_dict["parent_folder"] = retrieve(conf,"PARENT_FOLDER")

    # for the purpose of testing
    sim_dict["test_param"] = parse(Bool,retrieve(conf,"TEST_PARAM"))

    return sim_dict
end

function assess_stability(test_grn::Array{AbstractFloat,2},test_S_0::Array{AbstractFloat,1},max_convergence_time::Int)
    """
    Goal: This function measures whether a given GRN results in a stable final
    gene state vector for a given initial gene state vector.

    Inputs:
        1. test_grn: gene regulatory network matrix to be tested
        2. test_S_0: vector of initial gene states

    Outputs:
        1. boolean variable on whether the grn is stable
        2. the vector of final gene states
        3. the number of steps it took to reach a stable state
    """

    function filter_function(x::AbstractFloat)
        """
        Goal: This function normalizes the value of a gene's state to either
        0.0 or 1.0.
        """

        if x > 0.5
            return 1.0
        else
            return 0.0
        end
    end

    S_prev = deepcopy(test_S_0)
    for t=1:max_convergence_time
        S_next_temp = test_grn*S_prev
        S_next = map(filter_function,S_next_temp)

        #If stable, return stability, final gene state, and number of iterations
        #it took to reach state
        if S_next == S_prev
            return true, S_next, t
        else
            S_prev = deepcopy(S_next)
        end
    end

    #If not stable, return 0 vector for stable state
    return false,zeros(length(test_S_0)),-1
end

function generate_genotype(genes::Int,interaction_type::String,
    ancestor_connectivity::Real)
    """
    returns a geneotype, based on the specified size, connectivity, and
    interaction type.

    interaction_types:
    binary - all interacting genes have interaction of 1 or -1
    normal - normal distribution mean 0, sd 1
    uniform - uniform distribution over [-1,1]
    """

    grn = rand(interaction_dict[interaction_type],genes,genes)

    num_unconnected = convert(Int64, floor((1-ancestor_connectivity)*genes^2))
    idx = sample(1:genes^2,num_unconnected,replace=false,ordered=true)
    for id in idx
        grn[id] = 0.0
    end
    return Genotype(grn, Array{Mutation,1}())
end

function generate_condition(genes::Int,connectivity::AbstractFloat,deterministic::Bool=false)
    """
    NEED to rename
    generates an initial or final condition, based on the genome size and
    connectivity given. If deterministic, the number of 1s and 0s uses the
    floor of genes*connectivity
    """
    if deterministic
        n = floor(genes*connectivity)
        condition = [0.0 for _ in 1:(genes-n)]
        push!(condition,[1.0 for _ in 1:n]...)
        return sample(condition,10,replace=false)
    else
        return sample([0.0,1.0],Weights([connectivity,1-connectivity]),10)
    end
end

# function generate_genotype(genes::Int,sample_distribution::Sampleable,
#     ancestor_connectivity::Real)
#     """
#     returns a genotype based on the specified size and connectivity, using the
#     speicified sampling method.
#     """
#     grn = rand(sample_distribution,genes,genes)
#
#     num_unconnected = convert(Int64, (1-ancestor_connectivity)*genes^2)
#     idx = sample(1:genes^2,num_unconnected,replace=false,ordered=true)
#     for id in idx
#         grn[id] = 0.0
#     end
#     return Genotype(grn, Array{Mutation,1}())
# end

function hamming_distance(A::Array{AbstractFloat,1},B::Array{AbstractFloat,1})
    """
    returns the hamming distance of two vectors
    """
    distance = 0.0
    for i in 1:length(A)
        if A[i] != B[i]
            distance += 1.0
        end
    end
    return distance
end

function gen_grn_from_states(S_0::Array{AbstractFloat,1},target_eq::Array{AbstractFloat,1},
    interaction_type::String,ancestor_connectivity::Real,max_convergence_time::Int)
    """
    generates a genotype given an initial state, that converges to the target state
    """

    genes = length(S_0)
    @assert length(target_eq) == genes

    grn = nothing
    stable = false
    current_eq = target_eq
    while stable==false
        grn = generate_genotype(genes, interaction_type, ancestor_connectivity).grn
        stable, current_eq, t = assess_stability(grn,S_0,max_convergence_time)
    end

    p_change = 0.5
    p_gain = ancestor_connectivity*0.5
    p_loss = (1.0 - ancestor_connectivity)*0.5

    grn_found = false
    distance = hamming_distance(current_eq,target_eq)
    total_steps = 0
    while grn_found == false
        total_steps += 1
        new_grn = simple_random_mutate(grn,[p_change,p_gain,p_loss],interaction_type,
        sample(1:10))
        stable, next_eq, t = assess_stability(new_grn,S_0,max_convergence_time)
        if stable
            next_distance = hamming_distance(next_eq,target_eq)
            accept = false
            if next_distance < distance
                accept = true
            elseif next_distance == distance
                accept = sample([true,false])
            else
                p_accept = (distance/next_distance)^2
                accept = sample([true,false],Weights([p_accept,1-p_accept]))
            end
            if accept
                grn = new_grn
                current_eq = next_eq
                distance = next_distance
            end
        end
        if distance == 0.0
            grn_found = true
        end
    end
    println(total_steps)
    return grn
end

function simple_random_mutate(oldgrn::Array{AbstractFloat,2},mutation_probs::Array{AbstractFloat,1},
        interaction_type::String,mutations::Int)
    """
    finds a new mutation without altering the old grn, and without updating the
    mutation record of any genotype. iterating this function will produce a
    steady-state connectivity of p_gain/p_change
    """
    grn = deepcopy(oldgrn)
    genes = size(grn,1)

    for _ in 1:mutations
        mutation_found = false
        loc = -1
        newval = 0
        while mutation_found == false
            loc = sample(1:genes^2)
            mutation_type = sample(["change","gain","loss"],Weights(mutation_probs))

            if mutation_type == "change" && grn[loc] != 0.0

                if interaction_type == "binary"
                    newval = grn[loc]*-1
                else
                    newval =  rand(interaction_dict[interaction_type])
                end
            elseif mutation_type == "gain" && grn[loc] == 0.0
                newval = rand(interaction_dict[interaction_type])
            elseif mutation_type == "loss" && grn[loc] != 0.0
                newval = 0
            else
                continue
            end
            mutation_found = true
        end
        grn[loc] = newval
    end
    return grn
end

function mutate_genotype!(genotype::Genotype,mutation_type::String,interaction_type::String)
    """
    makes a single mutation on a grn of the specified type, at a random
    valid location
    """
    mutation_found = false
    grn = genotype.grn
    genes = size(grn,1)
    loc = -1
    oldval = 0
    newval = 0
    while mutation_found == false
        loc = sample(1:genes^2)
        if mutation_type == "loss" && grn[loc] != 0.0
            newval = 0
        elseif mutation_type == "change" && grn[loc] != 0.0

            if interaction_type == "binary"
                newval = grn[loc]*-1
            else
                newval =  rand(interaction_dict[interaction_type])
            end
        elseif mutation_type == "gain" && grn[loc] == 0.0
            newval = rand(interaction_dict[interaction_type])
        else
            continue
        end

        oldval = grn[loc]
        mutation_found = true
    end
    grn[loc] = newval
    push!(genotype.fixed_mutations, Mutation(loc, oldval, newval))
    return genotype
end

function pi(X, duplicates)
    """
    Produces a vector with the given indices duplicated
    """
    output = Array{Float64,1}()
    for i in 1:length(X)
        push!(output, X[i])
        if i in duplicates
            push!(output, X[i])
        end
    end
    return output
end

function pi_mat(A, duplicates)
    new_A = []
    for i in 1:length(A[1,:])
        push!(new_A,transpose(pi(A[i,:],duplicates)))
        if i in duplicates
            push!(new_A,transpose(pi(A[i,:],duplicates)))
        end
    end
    return vcat(new_A...)
end

# function mutate_genotype!(genotype::Genotype,mutation_type::String,sample_distribution::Sampleable)
#     """
#     Mutates a grn with using the specified mutation type, with the specified
#     distribution
#     """
#     mutation_found = false
#     grn = genotype.grn
#     genes = size(grn,1)
#     loc = -1
#     oldval = 0
#     newval = 0
#     while mutation_found == false
#         loc = sample(1:genes^2)
#         if mutation_type == "loss" && grn[loc] != 0.0
#             newval = 0
#         elseif mutation_type == "change" && grn[loc] != 0.0
#             # currently just inverts interaction
#             newval = grn[loc]*-1
#         elseif mutation_type == "gain" && grn[loc] == 0.0
#             newval = rand(sample_distribution)
#         else
#             continue
#         end
#         oldval = grn[loc]
#         mutation_found = true
#     end
#     grn[loc] = newval
#     push!(genotype.fixed_mutations, Mutation(loc, oldval, newval))
#     return genotype
# end

function binom_mutate!(genotype::Genotype,config::Dict)
    """
    determines the number of mutations for an individual using the given
    probabilities of interaction loss, gain, and change
    """
    genes = size(genotype.grn,1)
    empty = count(i->(i==0.0),genotype.grn)

    gains = rand(Binomial(empty,config["mu_gain"]))
    losses = rand(Binomial(genes^2 - empty,config["mu_loss"]))
    changes = rand(Binomial(genes^2 - empty,config["mu_change"]))

    for i in 1:gains
        mutate_genotype!(genotype,"gain",config["interaction_type"])
    end
    for i in 1:losses
        mutate_genotype!(genotype,"loss",config["interaction_type"])
    end
    for i in 1:changes
        mutate_genotype!(genotype,"change",config["interaction_type"])
    end
    return gains + losses + changes
end

function poisson_mutate!(genotype::Genotype,config::Dict)
    """
    Mutates using a poisson distribution, rather than binomial
    """

    gains = rand(Poisson(config["mu_gain"]))
    losses = rand(Poisson(config["mu_loss"]))
    changes = rand(Poisson(config["mu_change"]))

    for i in 1:gains
        mutate_genotype!(genotype,"gain",config["interaction_type"])
    end
    for i in 1:losses
        mutate_genotype!(genotype,"loss",config["interaction_type"])
    end
    for i in 1:changes
        mutate_genotype!(genotype,"change",config["interaction_type"])
    end
    return gains + losses + changes
end

function generate_population(population_size,genes,interaction_type,
        ancestor_connectivity,init_conditions::Array{Float64,1},
        multiplicities::Array{Int,1})
    """
    Generates a population given the specifications of the config file.
    Each member of the population will be viable under the init_conditions specified

    Improvement: take in list of clonal sizes
    TOTAL REWORK NECESSARY
    Assess stability needs fix (max convergence time)
    """
    individuals = []
    popsize = population_size
    for individual in 1:length(multiplicities)
        viable = false
        genotype = nothing
        while viable == false
            genotype = generate_genotype(genes,interaction_type,ancestor_connectivity)
            viable = assess_stability(genotype.grn,init_conditions,100)[1]
        end
        for _ in 1:multiplicities[individual]
            push!(individuals,deepcopy(genotype))
        end
    end
    return population
end

function generate_population(config,init_conditions::Array{AbstractFloat,1},
        multiplicities::Array{Int,1})
    """
    Generates a population given the specifications of the config file.
    Each member of the population will be viable under the init_conditions specified

    Improvement: take in list of clonal sizes
    TOTAL REWORK NECESSARY
    """
    individuals = []
    popsize = config["population_size"]
    for individual in 1:length(multiplicities)
        viable = false
        genotype = nothing
        while viable == false
            genotype = generate_genotype(config["genes"],config["interaction_type"],
                config["ancestor_connectivity"])
            viable = assess_stability(genotype.grn,init_conditions,
                config["max_convergence_time"])[1]
        end
        for _ in 1:multiplicities[individual]
            push!(individuals,deepcopy(genotype))
        end
    end
    return individuals
end

function wagner_experiment(genes,k_values,trials,subtrials)
    """
    framework for performing the experiment in Wagner's 1994 paper, which
    introduced the model
    """
    activation = 0.5
    connectivity = 1
    S_0_list = []
    S_eq_list = []
    grn_list = []
    for i in 1:trials
        push!(S_0_list,generate_condition(genes,activation,false))
        push!(S_eq_list,generate_condition(genes,activation,false))
        push!(grn_list,gen_grn_from_states(S_0_list[i],S_eq_list[i],"binary",connectivity,100))
    end

    different = zeros(length(k_values))
    stables = zeros(length(k_values))
    for (i,k) in enumerate(k_values)
        for trial in 1:trials
            for subtrial in 1:subtrials
                duplicates = sample(1:genes,k,replace=false,ordered=true)
                # println(duplicates)
                new_grn = pi_mat(grn_list[trial],duplicates)
                stable, new_S_eq, t = assess_stability(new_grn,pi(S_0_list[trial],duplicates),100)
                # print(stable)
                if stable
                    stables[i] += 1.0
                    # println(new_S_eq,pi(S_eq_list[trial],duplicates))
                    if new_S_eq != pi(S_eq_list[trial],duplicates)
                        different[i] += 1.0
                    end
                else
                    print(k)
                end
            end
        end
        different[i] = different[i]/stables[i]
    end
    print(stables)
    return different,S_0_list,S_eq_list,grn_list
end

function epistasis_experiment(config, generations)
    population_size = 500
    S_0 = generate_condition(config["genes"], config["activation"], true)
    population = generate_population(config,S_0,[population_size])
    for generation in 1:generations
        new_population = []
        for individual in 1:population_size
            individual_found = false
            new_individual = nothing
            while individual_found == false
                parent = sample(1:population_size)
                new_individual = deepcopy(population[parent])
                mutations = poisson_mutate!(new_individual,config)
                if mutations == 0
                    individual_found = true
                elseif assess_stability(new_individual.grn,S_0,config["max_convergence_time"])[1]
                    println("here!")
                    individual_found = true
                else
                    println("here2")
                end
            end
            push!(new_population,new_individual)
        end
        population = new_population
    end
    return population
end

# code for testing
config = generate_config()

S_0 = generate_condition(10,0.5,true)
S_eq = generate_condition(10,0.5,true)
grn = gen_grn_from_states(S_0,S_eq,config["interaction_type"], config["ancestor_connectivity"],20)
grn2 = generate_genotype(config["genes"],config["interaction_type"],
    config["ancestor_connectivity"])
# init_cond = rand([-1.0,1.0],config["genes"])
# oldgrn = deepcopy(grn)
poisson_mutate!(grn2, config)
# population = generate_population(config,init_cond,false)
# m =  pi_mat(grn,[1,2,3])
#
#
results,S_0_list,S_eq_list,grn_list = wagner_experiment(10,[0,1,2,3,4,5,6,7,8,9,10],100,100)
plot([0,1,2,3,4,5,6,7,8,9,10],results)

# population = epistasis_experiment(config, 100)


text:
foo() = 1
const bar = 2
# const_b has a function scope where the const on the previous line is
# We test here that we don't confuse the two files
include("const_b.jl")


text:
struct ConservationLawsProblem{islinear,isstochastic,MeshType,F,F3,F4,F5} <: AbstractConservationLawProblem{islinear,isstochastic,MeshType}
 f0::F5
 f::F
 CFL::F3
 tspan::Tuple{F4,F4}
 numvars::Int
 mesh::MeshType
end

isinplace{islinear,isstochastic,MeshType}(prob::AbstractConservationLawProblem{islinear,isstochastic,MeshType}) = false

function ConservationLawsProblem(f0,f,CFL,tend,mesh)
 numvars = size(f0(cell_faces(mesh)[1]),1)
 islinear = false
 isstochastic = false
 ConservationLawsProblem{islinear,isstochastic,typeof(mesh),typeof(f),typeof(CFL),typeof(tend),typeof(f0)}(f0,f,CFL,(0.0,tend),numvars,mesh)
end

### Displays
Base.summary(prob::AbstractConservationLawProblem{islinear,isstochastic,mType}) where {islinear, isstochastic, mType} = string("ConservationLawsProblem"," with mType ",mType)

function Base.show(io::IO, A::AbstractConservationLawProblem)
  println(io,summary(A))
  print(io,"timespan: ")
  show(io,A.tspan)
  println(io)
  print(io,"num vars: ")
  show(io, A.numvars)
end

max_stars_repo_name: Paulms/DiscontinuousGalerkin

max_issues_repo_name: Paulms/DiscontinuousGalerkin

max_forks_repo_name: Paulms/DiscontinuousGalerkin


text:
using OneVision
using OneVision.Examples
using OneVision.Examples.Car1DExample
using OneVision.Examples.Car2DExamples

using StaticArrays
using Plots
using Test
using Statistics: mean
using SimpleTypePrint: config_type_display
config_type_display()  # for simpler type error

@testset "Discretize a double integrator" begin
    A, B = [0.0 1.0; 0.0 0.0], colvec([0.0; 1.0])
    let dt = 1.23
        @show A′, B′ = discretize(A, B, dt)
        x0 = 1.1
        v0 = 1.4
        a = 1.45
        x1, v1 = A′ * [x0 v0]' + B′ * [a]'
        @test x1 == x0 + v0 * dt + a * dt * dt / 2
        @test v1 == v0 + a * dt
    end
end

let
    dt = 0.1
    H = 20  # horizon
    ΔT = 5  # control interval
    HΔT = H * ΔT
    sys = car_system(dt)

    @testset "Path tracking test" begin
        x_weights = CarX(5.0, 2.0, 0.0)
        u_weights = CarU(0.1)

        prob = TrajPlanningProblem(
            H,
            ΔT,
            sys,
            x_weights,
            u_weights,
        )
        x0 = CarX(-1.0, -1.0, 0.0)

        x_path = SVector{HΔT}(fill(CarX(1.0, 0.0, 0.0), HΔT))
        u_path = SVector{HΔT}(fill(CarU(0.0), HΔT))
        t0 = 0
        u, obj = plan_trajectory(prob, x0, x_path, u_path, t0)
        x = x_path_from_u(x0, t0, u, prob.dy, Val(ΔT))
        u_times = (1:H)*(dt * ΔT)
        x_times = (1:HΔT) * dt
        x_data = hcat(to_matrix(x), to_matrix(x_path))
        plot(
            plot(u_times, first.(u); label = "u"),
            plot(x_times, x_data; label = ["x" "v" "u'" "x*" "v*" "u'*"]),
            layout = (2, 1),
            size = (500, 300 * 3),
        ) |> display
        @test true
    end

    struct NoObs <: ObsDynamics end
    OneVision.obs_forward(::NoObs, x, z, t::𝕋) = z

    struct RendezvousControl <: CentralControlStateless{CarU{ℝ}} end

    OneVision.control_one(::RendezvousControl, xs, zs, t::𝕋, id::ℕ) = begin
        target = mean(xs)
        CarU(sum(target - xs[id]))
    end

    @testset "Forward prediciton" begin
        N = 2
        world = WorldDynamics(fill((sys, NoObs()), 2))
        s1 = (CarX(0.0, 0.0, 0.0), CarZ(false, 0.0))
        s2 = (CarX(2.0, 0.0, 0.0), CarZ(false, 0.0))
        prob = ForwardPredictProblem(world, RendezvousControl();
            X = CarX{ℝ}, Z = CarZ{ℝ}, Hf = HΔT)
        u_traj, x_traj = forward_predict!(prob, [s1,s2], nothing, 0, ΔT)
        plot(
            plot(1:HΔT, to_matrix(x_traj[:,1]); label = ["x1" "v1" "u1'"]),
            plot(1:HΔT, to_matrix(u_traj[:,1]); label = ["u1"]),
            layout = (2, 1),
            size = (500, 300*2),
        ) |> display
        @test true
    end

    @testset "Self estimation" begin
        s0 = CarX(0.0, 0.0, 0.0), 0
        times = 1:1000
        u_history = [CarU(sin(t/100)) for t in times]
        xs = self_estimate(sys, s0, u_history)
        plot(times * dt, to_matrix(xs); label = ["x" "v" "u'"]) |> display
        @test true
    end
end # Double integrater let

int_tests = let
    setting = ExampleSetting(
        time_end = 2, freq = 100, sensor_noise = 0.01, noise = 0.01,
        delays = default_delays, H = 20, model_error = 0.05)
    [
    () -> Car1DExample.run_example(;setting, plot_result = false)
    () -> Car2DExamples.tracking_example(;setting, plot_result = false)
    () -> Car2DExamples.formation_example(;setting, plot_result = false)
    ]
end
@testset "Integration test $i" for (i, f) in enumerate(int_tests)
    f()
    @test true
end


text:
# Stuff (mostly adjoints) to make things play nice with Zygote

nograd(f) = f()
Flux.Zygote.@nograd nograd

Flux.Zygote.@nograd mutate
Flux.Zygote.@adjoint! function dispatch!(lm::LazyMutable, m::ResetLazyMutable, x...)
     dispatch!(lm, m, x...), Δ -> nothing
end

# Basically a workaround as the pullback for the recursive implementation in NaiveNASlib which
#       1) uses get! which does not have a pullback function
#        and
#       2) was in some cases extremely slow or even stalled completely
Flux.Zygote.@adjoint function output!(memo::AbstractDict, v::AbstractVertex)
    return Flux.Zygote.pullback(__context__, output_loop!, memo, v)
end

function output_loop!(memo, v)
    vs = nograd() do
        # ancestors returns all input ancestors to v in topological order
        # We also provide all vertices for which we have the output already in memo
        # so we don't do unnecessary calculations.
        ancestors(v, collect(AbstractVertex, keys(memo)))[length(memo)+1:end]
    end

    for vn in vs
        vnins = inputs(vn) # Types don't seem to be inferred if put in map
        inpt = map(iv -> memo[iv], vnins)
        memo[vn] = vn(inpt...)
    end
    return memo[v]
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/BoundingSphere.jl/blob/master/LICENSE

__precompile__()
module BoundingSphere

if VERSION < v"0.7-"
else
    using Random
    using LinearAlgebra
    using Statistics: middle
end

include("api.jl")
include("boundary.jl")
include("geometry.jl")
include("welzl.jl")
include("ritter.jl")
include("util.jl")

end # module


text:
# Added from Oscar Dowson code for JuMP at https://github.com/jump-dev/JuMP.jl/pull/2484
generate_precompile = true
package_path = "C:\\Users\\wilhe\\Desktop\\Package Development"
example_path = ""

module Foo

using JuMP
using EAGO

function stress_precompile()
    for file in readdir("precompiles")
        if !endswith(file, ".jl")
            continue
        end
        include(file)
    end
    return
end

if generate_precompile
    using SnoopCompile
    tinf = @snoopi_deep Foo.stress_precompile()
    ttot, pcs = SnoopCompile.parcel(tinf)
    SnoopCompile.write("precompiles", pcs)
    for file in readdir("precompiles")
        if !endswith(file, ".jl")
            continue
        end
        src = joinpath("precompiles", file)
        m = match(r"precompile\_(.+)\.jl", file)
        modules = split(m[1], ".")
        modules = vcat(modules[1], "src", modules[2:end])
        if !(modules[1] in ["EAGO"])
            continue
        end
        dest = joinpath(package_path, modules..., "precompile.jl")
        @show dest
        cp(src, dest; force = true)
    end
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
struct Sorted{T}
    idx::UInt32
    value::T
end

Base.print(io::IO, s::Sorted) = print(io, s.value)
Base.isless(s1::Sorted, s2::Sorted) = isless((s1.idx, s1.value), (s2.idx, s2.value))

struct Renamer{U, L}
    uniquevalues::U
    labels::L
end

function renamer(p::Pair, ps::Pair...)
    pairs = (p, ps...)
    k, v = map(first, pairs), map(last, pairs)
    return Renamer(k, v)
end

function (r::Renamer)(x)
    i::UInt32 = findfirst(isequal(x), r.uniquevalues)
    label = r.labels[i]
    return Sorted(i, label)
end

struct NonNumeric{T}
    x::T
end

nonnumeric(x) = NonNumeric(x)

Base.print(io::IO, n::NonNumeric) = print(io, n.x)
Base.isless(n1::NonNumeric, n2::NonNumeric) = isless(n1.x, n2.x)


text:
"""
Depricated file
"""

function layerBackProp(
    cLayer::FCLayer,
    model::Model,
    actFun::SoS,
    labels::AbstractArray,
) where {SoS<:Union{Type{softmax},Type{σ}}}

    lossFun = model.lossFun
    dlossFun = Symbol("d", lossFun)
    A = cLayer.A
    Y = labels
    dZ = eval(:($dlossFun($A, $Y)))

    cLayer.dA = cLayer.W'dZ
    return dZ
end #softmax or σ layerBackProp

function layerBackProp(
    cLayer::Activation,
    model::Model,
    actFun::SoS,
    labels::AbstractArray,
) where {SoS<:Union{Type{softmax},Type{σ}}}

    lossFun = model.lossFun
    dlossFun = Symbol("d", lossFun)
    A = cLayer.A
    Y = labels
    dZ = eval(:($dlossFun($A, $Y)))

    return dZ
end #softmax or σ layerBackProp

function layerBackProp!(
    cLayer::FCLayer,
    model::Model,
    dA::AoN = nothing;
    labels::AoN = nothing,
    kwargs...,
) where {AoN<:Union{AbstractArray,Nothing}}
    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    m = size(cLayer.A)[end]

    A = cLayer.A

    regulization, λ = model.regulization, model.λ

    actFun = cLayer.actFun

    dZ = []
    if cLayer.actFun == model.outLayer.actFun

        dZ = layerBackProp(cLayer, model, eval(:($actFun)), labels)

    elseif dA != nothing
        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dA .*= D
            dA ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dA .* eval(:($dActFun($Z)))

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dA = []
        for nextLayer in cLayer.nextLayers
            try
                dA .+= nextLayer.dA
            catch e #in case first time DimensionMismatch
                dA = nextLayer.dA
            end
        end #for

        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dA .*= D
            dA ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dA .* eval(:($dActFun($Z)))

    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.dW = dZ * cLayer.prevLayer.A' ./ m

    if regulization > 0
        if regulization == 1
            cLayer.dW .+= (λ / 2m)
        else
            cLayer.dW .+= (λ / m) .* cLayer.W
        end
    end

    cLayer.dB = 1 / m .* sum(dZ, dims = 2)

    cLayer.dA = cLayer.W'dZ

    cLayer.backCount += 1

    return nothing

end #function layerBackProp(cLayer::FCLayer)

function layerBackProp!(
    cLayer::Activation,
    model::Model,
    dA::AoN = nothing;
    labels::AoN = nothing,
    kwargs...,
) where {AoN<:Union{AbstractArray,Nothing}}
    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    m = size(cLayer.A)[end]

    regulization, λ = model.regulization, model.λ

    actFun = cLayer.actFun

    dZ = []
    if cLayer.actFun == model.outLayer.actFun

        dZ = layerBackProp(cLayer, model, eval(:($actFun)), labels)

    elseif dA != nothing
        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dA .*= D
            dA ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dA .* eval(:($dActFun($Z)))

    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dA = []
        for nextLayer in cLayer.nextLayers
            try
                dA .+= nextLayer.dA
            catch e #if first time DimensionMismatch
                dA = nextLayer.dA
            end #try/catch
        end #for

        try
            keepProb = cLayer.keepProb
            if keepProb < 1.0 #to save time of multiplication in case keepProb was one
                D = rand(cLayer.channels, 1) .< keepProb
                dA .*= D
                dA ./= keepProb
            end
        catch e

        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.prevLayer.A #this is the Activation layer

        dZ = dA .* eval(:($dActFun($Z)))

    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.dA = dZ

    cLayer.backCount += 1

    return nothing
end #function layerBackProp!(cLayer::Activation

function layerBackProp!(
    cLayer::AddLayer,
    model::Model,
    dA::AoN = nothing;
    labels::AoN = nothing,
    kwargs...,
) where {AoN<:Union{AbstractArray,Nothing}}

    # m = size(cLayer.A)[end]
    #
    # A = cLayer.A
    #
    # regulization, λ = model.regulization, model.λ
    if dA != nothing
        keepProb = cLayer.keepProb
        if keepProb < 1.0 #to save time of multiplication in case keepProb was one
            D = rand(cLayer.channels, 1) .< keepProb
            dA .*= D
            dA ./= keepProb
        end

        dActFun = Symbol("d", cLayer.actFun)

        Z = cLayer.Z

        dZ = dA .* eval(:($dActFun($Z)))
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dA = []
        for nextLayer in cLayer.nextLayers
            try
                dA .+= nextLayer.dA
            catch e
                dA = nextLayer.dA #to initialize dA
            end #try/catch
        end #for

        cLayer.dA = dA

    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.backCount += 1

    return nothing
end #function layerBackProp!(cLayer::AddLayer

function layerBackProp!(
    cLayer::Input,
    model::Model,
    dA::AoN = nothing;
    labels::AoN = nothing,
    kwargs...,
) where {AoN<:Union{AbstractArray,Nothing}}

    # m = size(cLayer.A)[end]
    #
    # A = cLayer.A
    #
    # regulization, λ = model.regulization, model.λ

    if dA != nothing
        cLayer.dA = dA
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dA = []
        for nextLayer in cLayer.nextLayers
            try
                dA .+= nextLayer.dA
            catch e
                dA = nextLayer.dA #need to initialize dA
            end #try/catch
        end #for

        cLayer.dA = dA

    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.backCount += 1

    return nothing
end #function layerBackProp!(cLayer::Input

function layerBackProp!(
    cLayer::BatchNorm,
    model::Model,
    dA::AoN = nothing;
    labels::AoN = nothing,
    kwargs...,
) where {AoN<:Union{AbstractArray,Nothing}}
    prevLayer = cLayer.prevLayer
    lossFun = model.lossFun

    m = size(cLayer.A)[end]

    A = cLayer.A

    normDim = cLayer.dim

    regulization, λ = model.regulization, model.λ
    dZ = []
    if dA != nothing
        dZ = cLayer.W .* dA
    elseif all(
        i -> (i.backCount == cLayer.nextLayers[1].backCount),
        cLayer.nextLayers,
    )
        dA = []
        for nextLayer in cLayer.nextLayers
            try
                dA .+= nextLayer.dA
            catch e
                dA = nextLayer.dA #need to initialize dA
            end #try/catch
        end #for

        Z = cLayer.Z

        dZ = cLayer.W .* dA

    else #in case not every next layer done backprop
        return nothing

    end #if all(i->(i.backCount==cLayer.nextLayers[1].backCount), cLayer.nextLayers)

    cLayer.dW = sum(dZ .* cLayer.Z, dims = 1:normDim)

    if regulization > 0
        if regulization == 1
            cLayer.dW .+= (λ / 2m)
        else
            cLayer.dW .+= (λ / m) .* cLayer.W
        end
    end

    cLayer.dB = sum(dZ, dims = 1:normDim)

    N = prod(size(dZ)[1:normDim])

    varep = cLayer.var .+ cLayer.ϵ

    dẐ =
        dZ ./ sqrt.(varep) .*
        (.-(cLayer.Ai_μ_s) .* ((N - 1) / N^2) ./ (varep) .^ 2 .+ 1) .*
        (1 - 1 / N)

    cLayer.dA = dẐ

    cLayer.backCount += 1

    return nothing

end #function layerBackProp(cLayer::BatchNorm)

export layerBackProp!


text:
function WrappedGaussian(x, μ, σ)
    p = 0
    K = 10
    for k = -K:K
        p += exp(-0.5 * ((x - μ - 2pi * k) / σ)^2)
    end
    return p / (sqrt(2π) * σ)
end

function ω2k(ω, h)
    g = 9.81
    α = ω^2*(h/g)
    out = α*(tanh(α))^(-0.5)/h
    return out .* !isnan(out)
end


text:
# License for this file: MIT
# Copyright 2018-2019, DLR Institute of System Dynamics and Control
#
"""
   module Fluid1_Components

Date: Jan. 22, 2019

Fluid component models with ModiaMedia based on the new approach described in

- Zimmer, Bender, Pollok (2018): Robust Modeling of Directed Thermofluid Flows in Complex Networks.
  2nd Japanese Modelica Conference, Tokyo, Proceedings, pp 39-48. Download:
  https://www.modelica.org/events/modelica2018japan/conference-proceedings/modelica-final-proceedings-2018-Japan.pdf

# Main Authors of this module
```

      /|      Dirk Zimmer, Martin Otter
  ___/_|___   German Aerospace Center (DLR)
 /  /  /  /   Institute of System Dynamics and Control
/__/__/__/    https://www.dlr.de/sr/en/
   | /
   |/

```
"""
module Fluid1_Components

using  Modia
using  ModiaMedia

using  ModiaMedia.StaticArrays # included via ModiaMedia, to avoid requirement to add it in the standard environment
using  ModiaMedia.Unitful      # included via ModiaMedia, to avoid requirement to add it in the standard environment

export add!
export BaseProperties
export StaticPressure, InertialPressure, SpecificEnthalpy, Temperature, MassFlowRate
export InPlug, OutPlug, FixedSource_pT, FixedSink, ShortPipe, Splitter, Junction, Compressor, ClosedVolume

# Variable types
StaticPressure(  ; args...) = Variable(; start=1.0e5, nominal=1.0e5, size = (), T = u"Pa"  , info = "Pressure for static mass-flow", args...)
InertialPressure(; args...) = Variable(; start=0.0  , nominal=1.0e5, size = (), T = u"Pa"  , info = "Pressure for accelerating the fluid", args...)
SpecificEnthalpy(; args...) = Variable(; start=1.0e4, nominal=1.0e4, size = (), T = u"J/kg", info = "Specific enthalpy", args...)
Temperature(     ; args...) = Variable(; start=300.0, nominal=300.0, size = (), T = u"K"   , info = "Temperature", args...)
MassFlowRate(    ; args...) = Variable(; start=0.0  ,                size = (), T = u"kg/s", info = "Mass flow rate", args...)

const MediumVariable() = Var(size=())
const MediumState()    = Var(size=())

const mediumFunctionMap = Dict{String, Symbol}("T" => :temperature,
                                               "d" => :density)

"""
    add!(result, portName, var::Vector{String})

Add the variables of `var` to port `portName` and store the vectors in `result`.
"""
function add!(result, portName::AbstractString, var::Vector{String})
    Medium = result[portName * ".Medium"][1]
    p = result[portName * ".p"]
    h = result[portName * ".h"]

    for j in eachindex(var)
        fc = getfield(ModiaMedia, mediumFunctionMap[var[j]])
        v  = fill(0.0, length(p))
        for i in eachindex(p)
            v[i] = fc(Medium, setState_ph(Medium,p[i],h[i]))
        end
        result[portName * "." * var[j]] = v
    end

    return nothing
end

###### Define BaseProperties models supported here (user can provide further definitions)

@model SimpleMedium_BaseProperties begin
    p_start = 1.0e5                # parameter
    T_start = 300.0                # parameter

    Medium = MediumVariable()      # input
    p      = StaticPressure(  start=p_start)  # input
    h      = SpecificEnthalpy(start=1.0)      # input
    T      = Temperature(     start=T_start)  # output
    d      = Float(start=1.0)      # output
    u      = Float(start=1.0)      # output
    der_d  = Float(start=1.0)      # output
    der_u  = Float(start=1.0)      # output
    state  = MediumState()         # output

@equations begin
    d     = density(Medium)
    u     = specificInternalEnergy_T(Medium,T)
    state = setState_pT(Medium,p,T)
    h     = specificEnthalpy(Medium,state)
    der_d = 0.0
    der_u = specificInternalEnergy_T_der_2(Medium,T)*der(T)
    end
end

@model SimpleIdealGasMedium_BaseProperties begin
    p_start = 1.0e5                # parameter
    T_start = 300.0                # parameter

    Medium = MediumVariable()      # input
    p      = StaticPressure(  start=p_start)  # input
    h      = SpecificEnthalpy(start=1.0)      # input
    T      = Temperature(     start=T_start)  # output
    d      = Float(start=1.0)      # output
    u      = Float(start=1.0)      # output
    der_d  = Float(start=1.0)      # output
    der_u  = Float(start=1.0)      # output
    state  = MediumState()         # output

@equations begin
    d     = density_pT(Medium,p,T)
    u     = specificInternalEnergy_T(Medium,T)
    h     = specificEnthalpy_T(Medium,T)
    der_d = density_pT_der_2(Medium,p,T)*der(p) + density_pT_der_3(Medium,p,T)*der(T)
    der_u = specificInternalEnergy_T_der_2(Medium,T)*der(T)
    state = setState_pT(Medium,p,T)
    end
end

BaseProperties(Medium::ModiaMedia.SimpleMedium        ; p_start=1e5, T_start=300.0, args...) = SimpleMedium_BaseProperties(        Medium=Medium, p_start=p_start, T_start=T_start, args...)
BaseProperties(Medium::ModiaMedia.SimpleIdealGasMedium; p_start=1e5, T_start=300.0, args...) = SimpleIdealGasMedium_BaseProperties(Medium=Medium, p_start=p_start, T_start=T_start, args...)

####### Define Fluid components

"""
    connector InPlug - Uni-directional flow of fluid (mass flow is entering the component)
"""
@model InPlug(:connector) begin
    Medium = MediumVariable()
    m_flow = MassFlowRate(flow = true)
    r      = InertialPressure()
    p      = StaticPressure()
    h      = SpecificEnthalpy()
end

"""
    connector OutPlug - Uni-directional flow of fluid (mass flow is leaving the component)
"""
@model OutPlug(:connector) begin
    Medium = MediumVariable()
    m_flow = MassFlowRate(flow = true)
    r      = InertialPressure()
    p      = StaticPressure()
    h      = SpecificEnthalpy()
end

"""
    model FixedSource_pT - Fixed pressure and fixed specific enthalpy at source
"""
@model FixedSource_pT begin
    Medium  = MediumVariable()    # Medium MUST be redefined when instanting Source, otherwise error
    outPlug = OutPlug()
    p0 = StaticPressure(variability=parameter, info = "Fixed static pressure at source")
    T0 = Temperature(   variability=parameter, info = "Fixed temperature at source")
    @equations begin
        outPlug.Medium = Medium
        outPlug.r = 0.0
        outPlug.p = p0
        outPlug.h = specificEnthalpy(Medium, setState_pT(Medium,p0,T0))
    end
end

"""
    model FixedSink - Fixed pressure at sink
"""
@model FixedSink begin
    inPlug = InPlug()
    p0 = StaticPressure(variability=parameter, info = "Fixed static pressure at sink")
    @equations begin
        inPlug.p + inPlug.r = p0
    end
end

"""
    model ShortPipe - Simple pipe (without volume)
"""
@model ShortPipe begin
    inPlug  = InPlug()
    outPlug = OutPlug()

    k = Parameter(1.0e5, info="Quadratic pressure drop coefficient")
    l = Parameter(1.0  , info="Length of pipe")
    A = Parameter(1.0  , info="Area of pipe")    # should be taken from a global world model or propagated through the connection structure

    dp     = StaticPressure(  info = "delta pressure in flow direction")
    dr     = InertialPressure(info = "delta inertial pressure in flow direction")
    dh     = SpecificEnthalpy(info = "delta in spec. enthalpy in flow direction")
    m_flow = MassFlowRate()

    state = Var(size=())
    eta   = Float(size=())
    @equations begin
        # Medium propagation
        outPlug.Medium = inPlug.Medium

        # mass flow balance
        m_flow = inPlug.m_flow
        inPlug.m_flow + outPlug.m_flow = 0

        # Propagation of specific quantities
        outPlug.p  = inPlug.p + dp;
        outPlug.r  = inPlug.r + dr;
        outPlug.h  = inPlug.h + dh;

        dp = -m_flow*abs(m_flow) * k    # pressure drop
        dh = 0.0                        # isenthalpic process

        # inertial component of pressure
        der(m_flow)/A*l = -dr

        # compute dynamic viscosity to test dispatch of function
        state = setState_ph(inPlug.Medium, inPlug.p, inPlug.h)
        eta   = dynamicViscosity(inPlug.Medium, state)
    end
end

"""
    model Splitter - one inflow, two outflows
"""
@model Splitter begin
    inPlugA  = InPlug()
    outPlugB = OutPlug()
    outPlugC = OutPlug()

    eps = 0.000001

    @equations begin
        # Medium propagation
        outPlugB.Medium = inPlugA.Medium
        outPlugC.Medium = inPlugA.Medium

        # mass flow balance
        inPlugA.m_flow + outPlugB.m_flow + outPlugC.m_flow = 0.0

        # propagation of flow and pressure balance
        outPlugB.r = inPlugA.r
        outPlugB.p = inPlugA.p
        outPlugB.h = inPlugA.h

        outPlugC.r = inPlugA.r
        outPlugC.p = inPlugA.p
        outPlugC.h = inPlugA.h
    end
end

"""
    model Junction - two inflows, one outflow
"""
@model Junction begin
    inPlugA  = InPlug()
    inPlugB  = InPlug()
    outPlugC = OutPlug()

    eps = 0.000001

    @equations begin
        # Medium propagation
        inPlugA.Medium = outPlugC.Medium
        inPlugB.Medium = outPlugC.Medium

        # mass flow balance
        inPlugA.m_flow + inPlugB.m_flow + outPlugC.m_flow = 0.0

        # total pressure balance
        inPlugA.p + inPlugA.r = outPlugC.p + outPlugC.r
        inPlugB.p + inPlugB.r = outPlugC.p + outPlugC.r

        # mixing of fluid
            # formulation safeguarded against zero-massflow and backward flow
            (abs(outPlugC.m_flow) + 2*eps)*outPlugC.h = (abs(inPlugA.m_flow) + eps)*inPlugA.h + (abs(inPlugB.m_flow) + eps)*inPlugB.h

            # for simplicity m_flow is taken for the weighted mean of p but it should be V_flow
            (abs(outPlugC.m_flow) + 2*eps)*outPlugC.p = (abs(inPlugA.m_flow) + eps)*inPlugA.p + (abs(inPlugB.m_flow) + eps)*inPlugB.p
    end
end

"""
    model Compressor - Simple compressor
"""
@model Compressor begin
    inPlug  = InPlug()
    outPlug = OutPlug()

    pressureRatio = Parameter(2.0, info="Fixed pressure ratio")
    efficiency    = Parameter(0.9)
    kappa         = Parameter(1.4)
    R             = Parameter(8.314, info="Gas constant")

    k = Parameter(1.0e5, info="Quadratic pressure drop coefficient")
    l = Parameter(1.0  , info="Length of compressor")
    A = Parameter(1.0  , info="Area of compressor")    # should be taken from a global world model or propagated through the connection structure

    dr     = InertialPressure(info = "delta inertial pressure in flow direction")
    dh     = SpecificEnthalpy(info = "delta in spec. enthalpy in flow direction")
    m_flow = MassFlowRate()
    y      = Float(info="Specific energy inserted into the mass-flow")
    k_frac = Float()

    @equations begin
        # Medium propagation
        outPlug.Medium = inPlug.Medium

        k_frac = kappa/(kappa-1.0)

        # mass flow balance
        m_flow = inPlug.m_flow
        inPlug.m_flow + outPlug.m_flow = 0

        # Propagation of specific quantities
        outPlug.p = pressureRatio*inPlug.p
        outPlug.r = inPlug.r + dr
        outPlug.h = inPlug.h + dh
        y  = k_frac*R*(pressureRatio^(1/k_frac)-1)
        dh = y/efficiency;

        # inertial component of pressure
        der(m_flow)/A*l = -dr
    end
end

"""
    model ClosedVolume - Simple volume
"""
@model ClosedVolume begin
    Medium  = AbstractMedium    # Medium MUST be redefined when instanting ClosedVolume, otherwise error

    inPlug  = InPlug()
    outPlug = OutPlug()

    V  = 1.0       # Fixed volume
    p0 = Medium.infos.p_default     # initial pressure
    T0 = Medium.infos.T_default     # initial temperature

    M     = Float(info="Mass of volume"      , start=1.0)
    der_M = Float(info="Total mass flow rate", start=0.0)

    #medium = getBaseProperties(Medium)(p_start=p0, T_start=T0)
    medium = BaseProperties(Medium; p_start=p0, T_start=T0)

    h = Float(size=())
@equations begin
    # Medium propagation
    # outPlug.Medium = inPlug.Medium
    outPlug.Medium = Medium

    # Mass balance
    der_M          = inPlug.m_flow + outPlug.m_flow
    medium.der_d*V = der_M

    # Energy balance
    M = medium.d*V
    der_M*medium.u + M*medium.der_u = inPlug.m_flow*(inPlug.h + inPlug.r/medium.d) + outPlug.m_flow*outPlug.h

    # Properties at outPlug
    outPlug.p = medium.p
    outPlug.h = medium.h
    outPlug.r = 0.0

    # Properties of inPlug
    inPlug.p + inPlug.r = medium.p

    # For testing
    h = specificEnthalpy(Medium, medium.state)
    end
end

end


text:
function F_mmek(x, parms)
  S, E, SE, P = x
  k1, k2, k3 = parms

  binding = k1 * S * E
  dissociation = k2 * SE
  conversion = k3 * SE

  return [binding, dissociation, conversion]
end

function gillespie_mmek(;
  S = 301, E = 120, SE = 0, P = 0,
  k1 = 0.00166, k2 = 0.0001, k3 = 0.1)
  # initial conditions
  x0 = [S, E, SE, P]
  # stoichiometries
  nu = [
    [-1 -1 1 0]; # binding, S + E --> SE
    [1 1 -1 0];  # dissociation, SE --> S + E
    [0 1 -1 1]   # conversion, SE --> P + E
  ]
  # parameter vector
  parms = [k1, k2, k3]

  return x0, nu, parms
end

x0, nu, parms = gillespie_mmek()
F = F_mmek


text:
__precompile__()

module DBAPI

using Reexport

include("DBAPIBase.jl")
include("arrays.jl")

@reexport using .DBAPIBase
import .DBAPIBase: NotImplementedError, NotSupportedError
import .ArrayInterfaces: ColumnarArrayInterface
export ColumnarArrayInterface

end # module